Derivación e integración numéricas

Obtención de fórmulas de tipo interpolatorio

Para obtener fórmulas de derivación o integración numéricas a partir de la interpolación polinómica necesitamos calcular, en primer lugar, el polinomio de interpolación y después la derivada correspondiente en un punto a o la integral en un intervalo [a, b]. En lo que respecta a la derivación numérica, consideraremos distintos órdenes de derivación. Los datos de interpolación para las fórmulas de derivación numérica serán lagrangianos y para la integración numérica serán lagrangianos y de tipo Hermite. El polinomio interpolador se obtiene, en todos los casos, con la orden InterpolatingPolynomial.

Obtención de fórmulas de derivación numérica

En la derivación numérica, los nodos de la fórmula, es decir los puntos de interpolación de la función f, serán el punto a y puntos a su derecha e izquierda, de la forma a - h, a + h, a + 2  h, etc. La primera fórmula que vamos a obtener es la que utiliza los nodos a y a + h para la derivada primera.

Clear[a, h, f, x, datos]\(datos = {{a, f[a]}, {a + h, f[a + h]}}; 

    pol = Simplify[InterpolatingPolynomial[datos, x]]; 

    Simplify[∂\_x pol /. x -> a]\)

Obsérve cómo se ha expresado la derivada del polinomio respecto de su variable.

A continuación vamos a usar los nodos a - h, a y a + h para la derivada primera.

\(datos = {{a - h, f[a - h]}, {a, f[a]}, {a + \ h, f[a + \ h]}}; 

    pol = Simplify[InterpolatingPolynomial[datos, x]]; 

    Simplify[∂\_x pol /. x -> a]\)

Es una fórmula con sólo dos nodos (el coeficiente de f(a) vale cero). Por último, para la derivada primera, vamos a obtener una fórmula (que no se suele usar) que utiliza los nodos a, a + h y a + 3  h.

\(datos = {{a, f[a]}, {a + h, f[a + h]}, {a + 3  h, f[a + 3  h]}}; 

    pol = Simplify[InterpolatingPolynomial[datos, x]]; 

    Simplify[∂\_x pol /. x -> a]\)

Ahora vamos a obtener una fórmula para la derivada segunda, concretamente la que usa cinco nodos simétricos, a - 2  h, a - h, a, a + h y a + 2  h

\(datos = Table[{a + i\ h, f[a + i\ h]}, {i, \(-2\), 2}]; 

    pol = Simplify[InterpolatingPolynomial[datos, x]]; 

    Simplify[D[pol, {x, 2}] /. x -> a]\)

De igual forma se obtendría cualquier otra fórmula. A continuación vamos a analizar el orden de precisión de la última fórmula obtenida. Tenemos que ver hasta qué potencia del tipo (x^i) es exacta. De antemano, para i = 0,  1,  2,  3 y 4 está garantizada la exactitud por ser la fórmula de tipo interpolatorio. Al probar con i = 5 observamos que la fórmula da como resultado 20 (a^3), es decir el valor exacto de la derivada. Para i = 6 el valor proporcionado es 30 (a^4) - 8  (h^4), diferente de 30 (a^4). Luego el orden de precisión es 5.

f[x_] := (x^5)Simplify[
      -((30 f[a] + f[a - 2 h] - 16 f[a - h] - 16 f[a + h] + 
              f[a + 2 h])/(12 (h^2)))]f[x_] := (x^6)Simplify[
      -((30 f[a] + f[a - 2 h] - 16 f[a - h] - 16 f[a + h] + 
              f[a + 2 h])/(12 (h^2)))]

Comportamiento del error en las fórmulas de derivación numérica

Como es conocido, aunque sea muy simple, toda fórmula de derivación numérica que se obtuvo en el desarrollo teórico es convergente. Los errores tienden a cero cuando h tiende a cero. Sin embargo, en la computación aparecen otros errores que pueden desaconsejar el uso de valores muy pequeños de h si se quiere una buena aproximación. Es lo que ocurre con el eror de cancelación en las fórmulas de derivación numérica. Veámoslo. Vamos a dar el valor de las derivadas numéricas de una función f en el punto a = 2, con la fórmula ((f (a + h) - f (a))/h) para valores de h de la forma h = (10^(-j)) con j ≥ 1.Veremos que llega un momento en el que la derivada numérica deja de converger, con unas oscilaciones no controlables, para después pasar a valer cero, porque el numerador vale cero. Esto último ocurre cuando h es tan pequeño que para el ordenador a + h = h.

f[z_] := ((z - 0.5)/(1 + (z^2)))

\(n = 20; a = 2.; 

    derivadasnumericas = 
      Table[\(f[a + 10\^\(-j\)] - f[a]\)\/10\^\(-j\), {j, n}]\)

El valor exacto de la derivada de f en 2 es

f'[2.]

errores = f'[a] - derivadasnumericas

Los redondeos (a seis cifras) de la derivada numérica para j = 6, 7, 8  coinciden con el valor exacto. Al principio, a medida que h disminuye el error disminuye, después empieza un comportamiento extraño donde los errores pueden ser muy grandes y finalmente la aproximación que da es cero. Puede cambiar el valor de a y comprobarlo para otros valores y otras funciones. Por ejemplo el comportamiento del error en a = -1 se muestra a continuación:

\(a = \(-1.\); 

    derivadasnumericas = 
      Table[\(f[a + 10\^\(-j\)] - f[a]\)\/10\^\(-j\), {j, n}]\)errores = f'[a] - derivadasnumericas

El error mínimo en valor absoluto se alcanza cuando j = 8; después empieza a aumentar hasta j = 16. A partir de j = 17, el error se hace constante porque la derivada numérica es cero. Cambiamos ahora a una fórmula de mayor precisión y volvemos a calcular las derivadas numéricas y los errores.

derivadas2 = 
      Table[(1/2) ((f[a + (10^(-j))] - f[a - (10^(-j))])/(10^(-j))), {j, 
          n}]f'[a] - derivadas2

El comportamiento es similar al de la primera fórmula, sólo que el error es menor. Se consiguen más cifras exactas, once, sin que h tenga que disminuir tanto como antes. Ahora con j = 5 se han conseguido once cifras. Antes el error mínimo se alcanzó con j = 8 y fueron ocho cifras.

Integración numérica

Vamos a obtener las fórmulas del Trapecio y de Simpson.

Clear[a, b, f]\(pol = Simplify[InterpolatingPolynomial[{{a, f[a]}, {b, f[b]}}, x]]; 

    Simplify[∫\_a\%b pol \[DifferentialD]x]\)

\(datos = {{a, f[a]}, {\(a + b\)\/2, f[\(a + b\)\/2]}, {b, f[b]}}; 

    pol = Simplify[InterpolatingPolynomial[datos, x]]; 

    Simplify[∫\_a\%b pol \[DifferentialD]x]\)

Un método alternativo, pero equivalente, consiste en imponer exactitud. Veámoslo para la fórmula de Simpson.

\(B = Table[∫\_a\%b\( x\^i\) \[DifferentialD]x, {i, 0, 2}]; 

    A = Table[{a\^i, \((\(a + b\)\/2)\)\^i, b\^i}, {i, 0, 2}]; 

    sol = Simplify[LinearSolve[A, B]]\)

Simplify[sol.{f[a], f[((a + b)/2)], f[b]}]

Por último, vamos a obtener la fórmula del trapecio corregida, es decir, la que usa como datos el valor de la función y el de la derivada en cada extremo.

\(datos = {{a, {f[a], \(f'\)[a]}}, {b, {f[b], \(f'\)[b]}}}; 

    pol = Simplify[InterpolatingPolynomial[datos, x]]; 

    Simplify[∫\_a\%b pol \[DifferentialD]x]\)

Ejemplos de uso de fórmulas compuestas

Fórmula del rectángulo compuesta.

Vamos a aplicar varias fórmulas a la integración en [0, 2] de la función ((x^3) - x + 2) (e^x). Su valor exacto, en este caso, se puede obtener con la orden Integrate[función, {variable, límite inferior, límite superior}]. En aprox1 guardamos la aproximación obtenida con la fórmula del rectángulo derecha, en aprox2 guardamos el resultado del trapecio y en aprox3 el de Simpson. En todos los casos el número de subintervalos es el mismo, n = 200.

Integrate[((x^3) - x + 2) (E^x), {x, 0, 2}]

N[%]

f[x_] := ((x^3) - x + 2) (E^x)\(a = 0.; 
b = 2.; 
n = 200; 
h = \(b - a\)\/n; 

    nod = Table[a + i\ h, {i, 0, n}]; 

    aprox1 = h\ 
        \(∑\+\(i = 1\)\%n f[
            nod[[i + 1]]]\)\)

error1 = 3. + 3. (E^2.) - aprox1

Fórmulas compuestas basadas en las del trapecio y Simpson

aprox2 = 
      h ((1/2) f[nod[[1]]] + 
            (∑_(i = 2)^n) f[
                nod[[i]]] + 
            (1/2) f[nod[[n + 1]]])

error2 = 3. + 3. (E^2.) - aprox2

aprox3 = 
      (1/6) h (
          f[nod[[1]]] + 
            2 (∑_(i = 2)^n) f[
                  nod[[i]]] + 
            4 (∑_(i = 1)^n) f[
                  nod[[i]] + (h/2)] + 
            f[nod[[n + 1]]])

error3 = 3. + 3. (E^2.) - aprox3

Observamos la gran precisión de la fórmula de Simpson, que es una de las más usadas en la práctica.

Fórmulas gaussianas

Las fórmulas gaussianas usan como nodos aquellos que conducen a exactitud máxima. Si vamos a construir una fórmula de integración numérica con n nodos  sabemos que podemos elegirlos, de manera única, de forma que la fórmula correspondiente alcanza exactitud 2  n - 1. Por tanto, para construir fórmulas gaussianas hemos de calcular los nodos y los coeficientes. Con Mathematica tenemos varias alternativas. Una de ellas es plantear el sistema no lineal de las 2  n incógnitas que resulta de imponer exactitud para 1,  x,  ....,   (x^(-1 + 2 n)). Lo haremos a continuación, puesto que Mathematica lo resuelve sin dificultad para valores de n bajos.

La fórmula gaussiana con dos nodos

\(a = \(-1\); 
b = 1; 
n = 2; 
exactitud = 2\ n - 1; 

    nodo = Table[p[i], {i, n}]\)coefic = Table[c[i], {i, n}]incognitas = nodo ⋃ coefic

Observe cómo la lista incognitas se forma mediante la unión de las listas nodo y coefic, mediante el símbolo ⋃.

\(sistema = {}; 

    For[i = 0, i ≤ exactitud, \(i++\), 
      aux = coefic.Table[
            x\^i /. x -> 
                nodo[[j]], {j, n}]; 
      
AppendTo[sistema, 
        aux == \(b\^\(i + 1\) - a\^\(i + 1\)\)\/\(i + 1\)]]\)MatrixForm[sistema]

\(sol = Solve[sistema, incognitas]; 
TableForm[sol]\)

nodosGauss1 = 
      nodo /. sol[[1
          ]]coeficientesGauss1 = 
      coefic /. 
        sol[[1]]

Clear[f]coeficientesGauss1.Table[
        f[nodosGauss1[[i]]], {i, n}]

Uno de los inconvenientes de las fórmulas gaussianas es que los nodos suelen ser números trascencentes, y su aproximación introduce un primer error.

Mathematica y la integración numérica.

Mathematica tiene una orden de integración numérica, es NIntegrate, que nos da la respuesta con 6 dígitos, pero podemos recuperar hasta la precisión interna estándar (unas 15 cifras) si lo pedimos a través de WorkingPrecision.

NIntegrate[((x^3) - x + 2) (E^x), {x, 0, 2}, 
      WorkingPrecision -> 20]

No debe confundir NIntegrate con Integrate. La primera lleva programada una integración numérica y puede aplicarla a cualquier función definida y acotada en un intervalo [a, b]; sin embargo, la segunda requiere que la función tenga una primitiva, a la que aplica la regla de Barrow. Intente integrar en el intervalo [1, 2] la función ((e^x)/x) usando ambas órdenes.

Ejercicios

1.- Obtenga mediante interpolación en el espacio (Subscript[ℙ, 2]) una fórmula para aproximar f'' (a) del tipo combinación de f(a - h), f(a) y f(a + h).

2.- Con la fórmula obtenida en el ejercicio 1, halle una tabla de aproximaciones y errores de \(TraditionalForm\`\(f\_1\%'\)(2.5)\), siendo (Subscript[f, 1])(x) = (x^x), para h = (10^(-i)), i = 1,  ... , 5.

3.- Sea (Subscript[f, 2])(
          x) = (((x^2) + 40)/(x + (Sqrt[5  x]) + 7)). Calcule una tabla que recoja las derivadas de (Subscript[f, 2]) en (Subscript[x, i]) = 1, 2,  ... , 10, utilizando la primera fórmula de derivación numérica obtenida al inicio de la práctica, con h = (10^(-3)), y muestre al mismo tiempo el error cometido en cada punto. Repita el ejercicio con la segunda  fórmula obtenida para la derivada primera y, finalmente, para la obtenida en el ejercicio 1 (con respecto a la segunda derivada).

4.- Divida el intervalo [1, 2] en 100 partes iguales y aplique las fórmulas del rectángulo, Simpson y trapecio compuestas para aproximar la integral en dicho intervalo de (Subscript[f, 1]). Compare los resultados obtenidos con los que se calculan mediante el comando NIntegrate.

5.- Repita el ejercico 4 para (Subscript[f, 2]).

6.- Sea (Subscript[f, 3])(x) = (x^15) (e^x) en [0, 2]. Vamos a dividir el intervalo en 10  (2^n) subintervalos, es decir, 10, 20, 40, 80, ....y a aplicar la fórmula de Simpson compuesta hasta que la diferencia  entre dos aproximaciones consecutivas (por ejemplo, podrían ser con 20 y 40 subintervalos) sea menor que (10^(-2)), dando en tal caso por buena la última aproximación obtenida. Programe y calcule dicha aproximación. Compare el resultado con el obtenido mediante el comando NIntegrate y compare ambas aproximaciones con el valor exacto.

7.- Calcule las fórmulas gaussianas con 2 y 3 nodos,en el intervalo [-1, 1]. Aplíquelas para aproximar x (e^x) en [-1, 1] y compare los resultados con los obtenidos con NIntegrate y con el valor exacto obtenido con Integrate (organizando los cálculos de forma adecuada).