Práctica 10. Derivación e integración numéricas.
Aunque Matlab cuenta con comandos para realizar derivadas e integrales simbólicamente, en ocasiones no se conoce la expresión de la función a derivar o integrar y sólo se conocen algunos valores de dicha función. En otras ocasiones, aunque se dispone de la expresión de la función, resulta costoso hallar su función derivada o una primitiva de la función. Incluso, a veces, o no existe una primitiva de la función o bien la primitiva de la función no puede expresarse usando funciones elementales.
Contents
En esta práctica veremos algunas de las fórmulas que permiten aproximar el valor de la derivada de una función en un punto (fórmulas de derivación numérica) y el valor de la integral definida de una función (fórmulas de integración numérica).
Las fórmulas de derivación y las de integración numérica suelen requerir los valores de la función a derivar o a integrar en varios puntos. Estos valores pueden estar recogidos en un vector o hallarse evaluando la función.
Derivación numérica.
La fórmula de diferencia centrada en dos nodos permite aproximar la derivada de una función en un punto a partir de dos evaluaciones de dicha función.
difc=@(x,h) (f(x+h)-f(x-h))/(2*h)
difc =
@(x,h)(f(x+h)-f(x-h))/(2*h)
Ej 1. Emplee la fórmula de diferencia centrada en dos nodos para aproximar la derivada de la función f(x)=sqrt(1-x^2) en el punto x=sqrt(2)/2 para h=0.01 (Sol. -1.0002)
Ej 2. El script malp.m usa la fórmula de diferencia centrada en dos nodos para aproximar la derivada de f(x)=sqrt(1-x^2) en el punto x=sqrt(2)/2 para distintos valores h=10^(-i). Modifique el programa para que pueda representarse en el eje horizontal los valores de i y en el eje vertical el valor absoluto del error absoluto cometido. Vuelva a modificar el programa para que en el eje vertical se muestre log10(abs(y)) donde y=(-1)-derivadaap es el error absoluto cometido. Observe que un valor de h más pequeño no da necesariamente menor error. De los valores h empleados, ¿cuáles dan el error más pequeño? (Sol. h=10^(-6), h=10^(-7).)
Ej 3. Modifique el script malp.m para aproximar el valor de la derivada de f(x)=exp(x) en el punto x=0. Sabiendo que el resultado exacto es f '(0)=exp(0)=1 ¿qué valor de h da el error más pequeño? (Cap 6. Diez lecciones de cálculo numérico. Sanz-Serna) (Sol. h=10^(-5) )
Ej 4. Para aproximar f'(c) se considera la fórmula de derivación numérica en cinco nodos siguiente:
f(x-2h) -2*f(x-h) 2*f(x+h) -f(x+2h) --------- + --------- + --------- + --------- 12h 3h 3h 12h
Implemente dicha fórmula. Use la fórmula para hallar la derivada de la función f(x)=exp(x) en x=1 para h=10^(-3). Sabiendo que el valor exacto de la derivada es exp(1) ¿qué error se comete al considerar el valor que da la fórmula de cinco nodos anterior? (Sol. 5.0182e-013 )
Integración numérica.
A continuación vamos a trabajar con las fórmulas del trapecio y de Simpson y sus respectivas fórmulas compuestas.
Fórmula del trapecio
La fórmula del trapecio aproxima el valor de la integral definida entre x=a y x=b por el área del trapecio determinado por el eje OX, la recta x=a, la recta x=b y la recta que une los puntos (a,f(a)) y (b,f(b)).
f=@(x) x.^2+x+2; a=2; b=3; trapecio=(b-a)*(f(a)+f(b))/2
trapecio =
11
Fórmula del trapecio compuesta
Si dividimos un intervalo [a,b] en n intervalos de longitud h=(b-a)/n con nodos xi=a+i*h, podemos aplicar a cada uno de esos intervalos la fórmula del trapecio y obtener la fórmula del trapecio compuesta (h/2)*(f(a) + 2*f(a+h) + 2*f(a+2h) + 2*f(a+3h) + ... + 2*f(a+(n-1)h) + f(b))
En la función trapecioc.m se encuentra una implementación del método.
Ej 5. Incluya en la función trapecioc.m un comentario en la segunda linea indicando el objetivo de la función y otros comentarios explicando el significado de las variables a, b y n.
Ej 6. Usando la fórmula del trapecio compuesta, aproxime el valor de la integral definida entre x=-1 y x=1 de la función f(x)=sqrt(1-x^2) usando n=2, n=100 y n=1000 nodos. ¿El área de qué figura se está calculando cuando n=2? Sabiendo que el valor exacto de la integral es pi/2, halle el error cometido en las aproximaciones obtenidas. (Sol: 1, 1.5691, 1.5707. Errores absolutos 0.5708, 0.0017, 5.2588e-005 )
Fórmula de Simpson
En la fórmula de Simpson, el valor de la integral definida se aproxima mediante la integral definida de una parábola que interpola la función integrando en tres puntos.
f=@(x) x.^2+x+2; a=-1; b=1; simpson=(b-a)*(f(a)+ 4*f( (a+b)/2 )+f(b))/6
simpson =
4.6667
Ej 7. Defina una función anónima con argumentos a y b que aproxime la integral de una función f (definida previamente) entre x=a y x=b usando la fórmula de Simpson. Use dicha fórmula para aproximar la integral definida de cos(x) entre 0 y pi/2. (Sol. 1.0023)
Fórmula de Simpson compuesta.
Si dividimos un intervalo [a,b] en n intervalos de longitud h=(b-a)/n con nodos xi=a+i*h, podemos aplicar a cada uno de esos intervalos la fórmula de Simpson y obtener la fórmula de Simpson compuesta.
Ej 8. Para h=(b-a)/n y xi=a+i*h, implemente la fórmula de Simpson compuesta. (Puede inicializar un acumulador como s=h*f(a)/6, y vaya sumando 2*h*f( a+(i-1/2)*h )/3 + h*f( a+i*h )/3 para i=1,2,..., n-1; para completar el último intervalo, sume 2*h*f( a+(n-1/2)*h )/3 + h*f(b)/6. Compruebe la fórmula aproximando el valor de la integral de sen(x) entre 0 y pi considerando n=10. (Sol. 2.0000)
Ej 9. Usando la fórmula de integración del trapecio compuesta aproxime la integral de la función x(t)=t^3-t^2+t+2 entre t=0 y t=2.2 considerando n=10. Aproxime la misma integral mediante la fórmula de Simpson compuesta para n=5. Sabiendo que el valor exacto de la integral es 68453/7500, ¿cuál de los dos resultados se ha aproximado más al valor exacto? (Sol. Trapecio compuesta = 9.1679. Simpson compuesta = 9.1271.)
Ej 10. Para hallar el trabajo de expansión irreversible isotérmica de un mol de dióxido de carbono se necesita hallar el valor de la integral de la función P(V)=R*T/(V-b)-a/V^2 desde V1=0.0729762 litros hasta V2=26.4115 litros. Considere a=3.5924 L^2*atm/mol^2, b=0.04267 L/mol, R=0.08205746 (atm*L)/(mol*K) y T=323.15 K . (Sol. aprox. 130.3907 L*atm (usando Simpson compuesta con 5000 subintervalos))
Derivación e integración numéricas mediante aproximación.
Otra manera de hallar el valor de la derivada de una función en un punto (o el valor de la integral definida de una función) a partir de los valores de dicha función consiste en hallar un polinomio que aproxime tales valores y evaluar en dicho punto la derivada de ese polinomio (o hallar la integral definida de dicho polinomio).
En el siguiente ejercicio se aproxima una derivada a partir de un ajuste parabólico de unos datos.
Ej 11. Si f(t) es la concentración de peróxido de hidrógeno (H2O2) en función del tiempo, la derivada -df(t)/dt es la velocidad de la reacción de descomposición del peróxido de hidrógeno
H2O2 ---> H2O + (1/2) O2
Se quiere estimar dicha velocidad, sin conocer el orden de la reacción, a partir de un ajuste parabólico de los siguientes datos (Tabla 15.1, Química General, Petrucci)
tiempo en segundos 0 200 400 600 1200 1800 3000
[H2O2] en moles/litro 2.32 2.01 1.72 1.49 0.98 0.62 0.25
x=[ 0 200 400 600 1200 1800 3000]; y=[2.32 2.01 1.72 1.49 0.98 0.62 0.25];
i) Represente los datos de la tabla considerando los tiempos en el eje de abscisas y las concentraciones en el eje de ordenadas.
ii) Halle un ajuste parabólico de los valores de y respecto a los de x, y represente la parábola obtenida junto con los datos mencionados en i). (Sol: 2.2828e-007*x^2 - 0.0014*x + 2.2683 )
iii) Derive (en papel) la parábola obtenida. Evalúe la expresión de dicha derivada en x=0 y en x=1400 para hallar la velocidad de reacción para t=0 y para t=1400. (Observe que aunque en la gráfica de la parábola obtenida la pendiente de la tangente a cualquier punto es negativa, la velocidad de reacción, por definición, es positiva) (Sol: Para t=0, velocidad=13.5*10^(-4) Moles/segundo. Para t=1400, velocidad=7.13*10^(-4) Moles/segundo.)
iv) Represente el vector x en en el eje horizontal y el vector log(y) en el eje vertical. (Observe que la gráfica es similar a una recta) (Sabiendo que la reacción es de primer orden, se cumple que ln f(t) = -k*t + f(0), con lo que los datos obtenidos en los apartados ii) y iii) son una aproximación, pero no representan un modelo del proceso.)
% Matemáticas II. Grado en Ingeniería Química. (A. Palomares)
disp(date)
19-Apr-2017