Derivación e integración numéricas

Derivación numérica

Mediante f'[x], D[f[x],{x,1}] o D[f[x],x] podemos conseguir la derivada de la función f, de la variable independiente x. Análogamente, con f'[x], f''[x] se logran las de órdenes dos y tres y con D[f[x],{x,n}] obtenemos la derivada n-ésima. Sin embargo, a veces es necesario emplear aproximaciones numéricas de la derivada en un punto (por ejemplo, cuando sólo se conoce de f los valores que toma en un conjunto discreto de puntos). A ello es a lo que dedicamos esta práctica.

Ejemplo de construcción de una fórmula de derivación en un punto a, usando los valores de a, a+h y a+2h

Podemos usar las funciones incorporadas por Mathematica para encontrar la fórmula exacta en P2, espacio de los polinomios de grado menor o igual que dos,  para aproximar f'(a) utilizando el valor de la  función en los puntos a, a+h, a+3h, por ejemplo.

Clear[a, h, x]\(pol[x_] := 
      Simplify[InterpolatingPolynomial[{{a, f[a]}, {a + h, f[a + h]}, 

            			\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ 	{a + 2\ h, f[a + 2\ h]}}, x]]; 

    \(pol'\)[a] // Simplify\)

Fórmulas de derivación numérica y errores asociados

De la misma manera podemos establecer las siguientes fórmulas para aproximar la derivada de f en a,
              der1=(f[a]-f[a+h])/h
              der2=(f[a+h]-f[a-h])/(2h),
que serán utilizadas a continuación.

Consideramos la función f(z)=(z-0.5)/(1+z^2). Vamos a estimar su derivada en el punto a=1 para mvalores de h=10^(-j) , j=1,2,3,4, y de las fórmulas indicadas anteriormente.

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

n = 4; derivadas1 = 
      Table[(f[1 + (10^(-j))] - f[1])/ (10^(-j)), {j, n}]errores1 = derivadas1 - f'[1]

Operamos ahora de la misma manera pero utilizando la fórmula der2.

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

Integración numérica

Otro problema que hay que abordar con frecuencia es el que consiste en estimar el valor de la integral definida de una cierta función. También podemos utilizar el procedimiento consistente en calcular fórmulas exactas en espacios de polinomios deseados, como en el caso de la derivación, además de emplear herramientas específicas proporcionadas por Mathematica.

Ejemplo de construcción de una fórmula de integración numérica
de tipo interpolatorio exacta en P2 que utiliza los puntos a, (a+b)/2 y b

Vamos a establecer la fórmula de Simpson utilizando la facilidad de cálculo del polinomio de interpolación. Comenzamos calculando el polinomio de interpolación en los tres nodos indicados.

Clear[a, b, f]\(\(pol = 
      Simplify[InterpolatingPolynomial[{{a, f[a]}, {\(a + b\)\/2, 
              f[\(a + b\)\/2]}, 

            				\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ {b, f[b]}}, x]]; \)\)

Lo Integramos en [a,b].

Simplify[(∫_a^b) pol \[DifferentialD]x]

Ejemplos de uso de fórmulas compuestas asociadas

Fórmula del rectángulo y su compuesta.

La fórmula del rectángulo a la derecha (que aproxima la integral en [a,b] por (b-a) f(b) ) aplicada a la partición del intervalo [a,b] en n subintervalos de igual longitud es la fórmula del rectángulo a la derecha compuesta.  Por ejemplo, podemos estimar la integral de la función (x^3-x+2) E^x en [0,2] subdividiendo este intervalo en 200 partes iguales y aplicando la fórmula del rectángulo a la derecha..

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]]]

Podemos comparar los resultados obtenidos con el valor exacto,  que es 3+3 E^2.

error1 = Integrate[f[x], {x, a, b}] - aprox1

Fórmulas compuestas basadas en las del trapecio y Simpson

Vamos a operar de la misma forma empleando como fórmulas base las del trapecio y Simpson. Comenzamos con la del trapecio

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

El error cometido es

error2 = Integrate[f[x], {x, a, b}] - aprox2

Con la fórmula de Simpson se trabaja de manera análoga.

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

El error que se comete es

error3 = Integrate[f[x], {x, a, b}] - aprox3

Mathematica y la integración numérica.

Mathematica incluye una función para la integración numérica cuya sintaxis es 
      NIntegrate[función,{variable, extremo inferior, extremo superior}]
Se trata de una orden muy precisa, construida a partir de un algoritmo de tipo adaptativo.

NIntegrate[f[x], {x, 0, 2}]

error = Integrate[f[x], {x, a, b}] - NIntegrate[f[x], {x, 0, 2}]

Ejercicios

1.- Obtenga mediante interpolación en el espacio P2 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 ejercicio1, halle una tabla de aproximaciones y errores de f''(2.5), siendo f(x)=x^x, para h=10^(-i), i=1,...,5

3.- Sea f(x)=(x^2+dm +ds)/(x+(dm x)^(1/2)+d1). Calcule una tabla que recoja las derivadas de f en a=3, utilizando la fórmula der1 con h=10^(-3) y mostrando al mismo tiempo el error cometido en cada punto. Repita el ejercicio con la fórmula der2 y 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 de rectángulo, Simpson y trapecio compuestas para aproximar la integral en dicho intervalo de  f(x)=x^x.  Compare los resultados obtenidos con los que se calculan mediante el comando NIntegrate.

5.- Repita el ejercico 4 param  f(x)=(x^2+dm +ds)/(x+(dm x)^(1/2)+d1).

6.- Deseamos calcular el área de una región plana acotada por un contorno curvilíneo del que se conocen 10 puntos: (0,0), (2,5), (4,12), (6,8), (8,0), (6,-10), (4,-8), (2,-4) y (0,0). Estos puntos están ordenados en el sentido de las agujas del reloj. 
a) Calcule  dos splines cuadráticos de clase 2, uno para los puntos de ordenada positiva y otros para los de ordenada negativa, sabiendo que, en ambos casos, la tangente en el punto (0,0) tiene pendiente nula.
b) Represente el dominio aproximado obtenido.
c) Halle  el área aproximada del recinto considerado utilizando la fórmula de Simpson compuesta para lo cual, en cada punto medio entre dos nodos consecutivos, se aproxima la ordenada en el punto por la imagen,  mediante el spline calculado, de la abscisa correspondiente.