Derivadas e Integrales Numéricas


Los temas cubiertos en este capítulo están basados en el libro
Numerical Analysis, de Richard L. Burden y J. Douglas Faires, PWS-Kent, 3rd Edition (1985)

Derivadas
Los métodos de aproximación vistos anteriormente, son especialmente útiles para derivar e integrar.
Utilizando la notación vista en las interpolaciones de funciones, podemos escribir en orden 1: 


f(x) = P1(x) + (x-x0)(x-x1)/2! f''(ξ)

donde P1(x) es el polinomio de primer grado construido con polinomios de Lagrange 


P1(x) = f(x0)L1,0(x) + f(x1)L1,1(x)
      = f(x0)(x - x1)/(x0 - x1) + f(x1)(x - x0)/(x1 - x0) 
      = f(x0)(x - x1)/(-h) + f(x1)(x - x0)/(h) 
      = f(x0)(x - x1)/(-h) + f(x1)(x - x0)/(h) 
      = f(x0)(x - x1)/(-h) + f(x0+h)(x - x0)/(h)


Diferenciando tenemos 


f'(x) = ( f(x0 + h)- f(x0) )/h + [(x-x0)(x-x1)/2! f''(ξ)]'

Para realizar la última derivada, necesitamos el conocimiento de la derivada segunda de la función, por lo tanto, en general, no podemos evaluar el error de truncamiento. Sin embargo, para el caso en que x-x0: 


f'(x0) = ( f(x0 + h)- f(x0) )/h - f''(ξ)h/2

Ejercicios:

  • Calcular la derivada de f(x) = ln x , y determinar el error de truncado, para x0=1.8 y h=0.1, 0.01, 0.001.

    Interpolación de Richardson
    Existen muchas variantes para hacer la derivación numérica de una función, utilizando ordenes mas altos de los polinomios interpolantes, y combinándolos en forma diferente. Dado que la bibliografía es muy abundante para estos casos, sólo veremos una forma de producir mejores resultados utilizando expansiones de orden bajo.
    Si bien este método fue publicado en 1927, se supone que la idea básica es muy anterior, remontándose a Arquímedes (200 a.c.).
    Supongamos que tenemos una fórmula de aproximación N(h) que nos da como resultado un valor 

    
M = N(h) + K1 h2

    si en lugar de calcular el valor buscado M con un paso h, lo hiciésemos con un paso h/2 tendríamos: 

    
M = N(h/2) + K1 h2/4

    La segunda aproximación se supone que es mas precisa, pero aún podemos hacer un mejor trabajo, aprovechando estos dos calculos. Multiplicamos la segunda ecuación por 4 y restamos a la primera, obteniendo: 

    
3M = 4N(h/2) -  N(h) 

de donde tenemos un nuevo valor de M

               4N(h/2) -  N(h)
          M = ----------------
                      3

    Del mismo modo, podemos continuar dividiendo h/2 una vez mas por 2, y nos quedaría: 

    
                    4j-1Nj-1(h/2) -  Nj-1(h)
          M = Nj = --------------------------
                           4j-1 - 1

    Ejercicios:

  • Calcular la derivada de f(x) = x ex , para x0=2.0 y h=0.2, 0.1 y 0.05. Utilizar la extrapolación de Richardson para mejorar los resultados, y comparar con el valor exacto.

    Integrales
    Del mismo modo en que sugerimos leer la abundante bibliografía para las derivadas numéricas, en la mayoría de los libros de cálculo numérico existe un capítulo de integración.
    Entre los numerosos métodos de integración utilizando polinomios interpoladores, revisaremos el Método de Romberg.
    La primer aproximación de la integral por el método de los trapecios se puede escribir tomando todo el intervalo h1=(b-a) de la siguiente manera: 

    
 Intabf(x) = R1 = [ f(a) + f(b)] h1/2 - h1/12 h12 f''(ξ)

    Si partimos el intervalo en 2, con h2=(b-a)/2 tenemos aproximadamente: 

    
 Intabf(x) ~ R2 = [ f(a) + f(b) + 2f(a + h2) ] h2/2 
                = (1/2)[ R1 + h1 f(a + h1/2) ]

    y en forma general 

    
 Intabf(x) ~ Rk =  (1/2)[ Rk-1 + hk-1 (f(a + hk-1/2) + f(a + hk-13/2) + 
                + f(a + hk-13/2) + ... + f(a + hk-1(2(k-1)-1/2) ) ]


    La integración de Romberg, utiliza todos los cálculos parciales de la integración para encontrar un valor final. Si ahora denominamos a los Ri en la forma Ri,1, podemos ir mejorando los resultados iterativamente: 

    
                    4j-1Ri,j-1 -  Ri-1,j-1
         Ri,j = --------------------------
                           4j-1 - 1

    Ejercicios:

  • Mostrar que la integral entre 0 y π de sinx con las diferentes subdivisiones de los intervalos es:
    1. R1 = 0
    2. R2 = 1.57079633
    3. R3 = 1.89611890
    4. R4 = 1.97423160
    5. R5 = 1.99357034
    6. R6 = 1.99839336
  • Utilizar el método de Romberg para obtener un valor mejorado de esta integral. Comparar con el valor exacto.
  • Escribir un programa en Fortran que calcule una integral utilizando el método de Romberg. Utilizar este programa para calcular la integral de xex entre 0 y 1.

    Cuadratura Gaussiana
    El método de Gauss utiliza la posibilidad de elegir los puntos en los que será calculada la función, de manera de minimizar el error en la integración.
    Este método se basa en la utilización de polinomios ortogonales, como por ejemplo los polinomios de Legendre: 

    
P0(x)= 1
P1(x)= x
Pk(x)=  Pk-1(x)(x - Bk) - Pk-2(x)Ck
siendo
Bk = (Intg x(Pk-1(x))2 dx) /(Int (Pk-1(x))2 dx)

y

Ck = (Int x Pk-1(x) Pk-2(x)  dx) /(Int (Pk-2(x))2 dx)

y las integrales se hacen entre [-1,1].

    Los polinomios Pn son ortogonales y generan a cualquier polinomio de orden n, en el rango [-1,1].
    Supongamos que calcular (en P[-1,1]) la integral de una función f(x), que se puede escribir como 

    
f(x) = Q(x)Pn(x) + R(x)

    donde Pn(x) es el polinomio de Legendre de orden n.
    Como Q se puede escribir como combinación lineal de los Pn, y estos son ortogonales, tenemos que 

    
Int f(x) dx = Int R(x) dx 
 
y expandiendo R en polinomios de Lagrange obtenemos

Int f(x) dx =  Sum ( R(xi) Int Li dx ) = Sum ( R(xi) ci)

    Si nos fijamos en los xi que son raíces de los polinomios de Legendre, 

    
f(xi) = Q(xi)Pn(xi) + R(xi) = 0 + R(xi) = R(xi)

    por lo que nos queda 

    
Int f(x) dx = Sum ( f(xi) ci)


    Cualquier integral en un rango definido [a,b] , se puede convertir en una integral entre [-1,1] mediante la transformación 

    
t = 1/(b-a)  (2x - a - b) 

por lo que 

Intab f(x) dx = Int-11 f([(b-a)t + b + a ]/2) (b-a)/2 dt


    Ejercicios:

  • Calcular la integral e-(x2), entre [1,1.5] por el método de integración Gaussiana con polinomios de Legendre de grado 1,2 y 3.
  • Utilizar el programa gaussian.for para calcular la integral anterior, con mayor precisión.
  • Comparar los métodos de Romberg y Gauss integrando algunos casos particulares.
  • Calcular la integral multidimensional ln(x + 2y), con x entre 1.4 y 2.0 e y entre 1.0 y 1.5. Utilizar un método de Simpson y comparar con la integración Gaussiana de orden 3. Cuántas veces se evalúa la función en cada método?


    Darío Mitnik