Práctica 6: Derivación numérica y extrapolación de Richardson.

1 Introducción

Recordemos que posiblemente la más sencilla de las fórmulas de derivación numérica es la diferencia dividida

f'(x) =

f(x+h)-f(x)

+ e(x) , e(x)=-

f''( x )

h , (1)

que es de primer orden (el error cometido es proporcional a h).

A su vez, la diferencia simétrica

f'(x) =

f(x+h)-f(x-h)

+ e(x) , e(x)=-

f'''( x )

h² , (2)

es una fórmula de segundo orden (el error cometido es proporcional a h²).

Para aproximar la segunda derivada podemos usar

f'' (x) =

f(x+h)-2f(x)+f(x+h)

h²

+ e(x) , e(x)=-

f^(IV)( x )

, (3)

que es también de segundo orden.

Sin embargo, recordemos que todas estas fórmulas sufren del problema de cancelación cuando tomamos un paso h muy pequeño. Se ha visto en teoría que el valor de h cercano al óptimo es el que satisface

|e(x)|=

2 e

, (4)

donde e es la cota del error absoluto en redondeo, que se puede coger igual a un ulp.

2 Fórmulas de derivación numérica

Vamos a intentar reproducir algunos de los resultados numéricos que aparecen en el libro, y de paso estudiar experimentalmente las fórmulas anteriores. Para ello tomemos como función de prueba alguna para la cual las derivadas sean fáciles de calcular analíticamente. Por ejemplo, sea

f(t) = sen (t)

y el punto donde calcularemos la derivada, el x=p/3.2. Puesto que |f(x)| » 1, es razonable tomar e en la fórmula (4) igual al épsilon de la máquina.

Usando la fórmula (4) estime el orden del valor óptimo de h para cada una de las fórmulas (1)-(3).
Cree un script con el nombre derivada.m, que contenga las instrucciones para los siguientes pasos:
- Cree el vector h= 2^-1, 2^-2, ..., 2^-30 de los pasos que se utilizarán para las fórmulas (1)-(3).
- Cree los vectores f1 y f2, que contengan la evaluación de f en los elementos de x+h y x-h, respectivamente.
- Cree finalmente los vectores der1 y der2 que contengan las aproximación dada por las fórmulas (1) y (2), respectivamente.
- Usando el valor exacto de la derivada, calcular los vectores err1 y err2 del error absoluto cometido.
- Dibuje en escala logarítmica las gráficas de los errores como función del paso h. Puede emplear para ello la función loglog de MATLAB.
- Haga lo mismo para la segunda derivada (3) y añádela a la gráfica anterior.
- Agregue al gráfico la leyenda (función legend) en la esquina superior izquierda, que especifique a qué fórmula corresponde cada gráfica (consulte la ayuda del MATLAB si lo precisa).
Interprete los resultados, ¿corroboran o no el orden del error teórico? ¿Cuál es el paso óptimo en cada caso? ¿Confirman el óptimo teórico que hemos calculado en el paso 1?

3 Extrapolación de Richardson

Un paso del algoritmo de extrapolación de Richardson, aplicado a una función G(h), puede venir dado por la fórmula

G₁(h)=

2^-k G(h) - G(h/2)

2^-k-1

donde k es el orden de convergencia esperado de la función G(h).

Considere los elementos del vector der1, calculado anteriormente, como valores de la función G(h) y aplíquele el algoritmo de Richardson.
En una ventana diferente a la anterior (use el comando figure(2) del MATLAB) dibuje en escala logarítmica las gráficas de los errores de G y G₁ como funciones del paso h.
Diga si considera que se ha logrado elevar el orden de convergencia a 2, tal como lo afirma el teorema de Richardson.

This document was translated from L^AT_EX by H^EV^EA.
Profesor: Andrei Martínez Finkelshtein