CAPÍTULO 23 Diferenciación numérica En el capítulo 4 ya se introdujo la noción de diferenciación numérica. Recuerde que se emplearon las expansiones en serie de Taylor Taylor para obtener las la s aproximaciones de las derivadas en diferencias divididas finitas. En el mismo capítulo se desarrollaron las aproximaciones en diferencias divididas hacia adelante, hacia atrás y centradas para la primer pri mer derivada y las derivadas der ivadas de orden superior. Recuerde que, en el mejor de los casos, dichas estimaciones tenían errores que fueron O (h2); es decir, sus errores eran era n proporcionales al cuadrado del tamaño tama ño de paso. Este nivel de exactitud exactitud se debe al número de términos térmi nos de la serie de Taylor Taylor que se utilizaron utiliza ron durante la deducción de esas fórmulas. Ahora se mostrará cómo desarrollar fórmulas de mayor mayor exactitud utilizando más términos. tér minos.
23.1 FÓRMULAS DE DIFERENC DIFERENCIACIÓN IACIÓN CON ALT ALTA A EXACTIT EXACTITUD UD Como se indica antes, se pueden generar fórmulas por diferencias divididas de alta exactitud tomando términos adicionales en la expansión de la serie de Taylor. Por ejemplo, la expansión de la serie de Taylor hacia adelante se escribe como [ecuación (4.21)]: f ( x x i+1) = f ( x x i) + f ′( x x i)h +
f ′′( x i )
2
h2 +···
(23.1)
de la que se despeja f ′ x i) = ′( x
f ( xi +1 ) − f ( x i )
−
f ′′( x i )
h + O(h2) (23.2) 2 En el capítulo 4 truncamos este resultado al excluir los términos de la segunda derivada en adelante y nos quedamos con un resultado final,
′( x f ′ x i) =
h
f ( xi +1 ) − f ( x i ) h
+ O(h)
(23.3)
Ahora retendremos, en cambio, el término de la segunda derivada sustituyendo la siguiente aproximación de la segunda derivada [recuerde la ecuación (4.24)] f ′′ x i) = ′′( x
f ( x i +2 ) − 2 f ( xi +1 ) + f ( x i ) h
2
+ O(h)
(23.4)
en la ecuación (23.2) para dar f ′( x i ) =
f ( x i +1 ) − f ( x i ) h
−
f ( x i +2 ) − 2 f ( x i +1 ) + f ( x i )
2h
2
h + O(h 2 )
669
23.1 FÓRMULAS DE DIFERENCIACIÓN CON ALTA EXACTITUD
o, al agrupar términos, f ′( x i ) =
− f ( xi+2 ) + 4 f ( xi +1 ) − 3 f ( x i ) 2h
+ O(h 2 )
(23.5)
Observe que al incluir el término de la segunda derivada mejora la exactitud a O (h2). Es posible desarrollar versiones similares mejoradas para las fórmulas hacia adelante y centradas, así como para las aproximaciones de derivadas de orden superior. Las fórmulas se resumen en las figuras 23.1 a 23.3, junto con todos los resultados del capítulo 4. El siguiente ejemplo ilustra la utilidad de esas fórmulas para la estimación de las derivadas.
FIGURA 23.1 Fórmulas de diferencias divididas finitas hacia delante: se presentan dos versiones para cada derivada. La última versión emplea más términos de la expansión de la serie de Taylor y, en consecuencia, es más exacta. Primera derivada
Error
f (x i +1) − f (x i ) h −f (x i +2 ) + 4f (x i +1) − 3f (x i ) f ′(x i ) = 2h f ′(x i ) =
O (h) O (h2)
Segunda derivada
f (x i +2 ) − 2f (x i +1) + f (x i ) h2 −f (x i +3 ) + 4f (x i +2 ) − 5f (xi +1) + 2f (x i ) f ″ (x i ) = h2 f ″ (x i ) =
O (h)
O (h2)
Tercera derivada
f (x i +3 ) − 3f (x i +2 ) + 3f (x i +1) − f (x i ) h3 −3f (x i + 4 ) + 14f (x i+3 ) − 24f (x i +2 ) + 18f (x i +1) − 5f (x i ) f ″ (x i ) = 2h3 f ″ (x i ) =
O (h) O (h2)
Cuarta derivada
f (x i + 4 ) − 4f (x i+3 ) + 6f (x i +2 ) − 4f (x i +1) + f (x i ) h4 −2f (x i +5 ) + 11f (x i+ 4 ) − 24f (x i +3 ) + 26f (x i+ 2 ) − 14f (x i+ 1) + 3f (x i ) f ″″ (x i ) = h4 f ″″ (x i ) =
O (h)
O (h2)
670
DIFERENCIACIÓN NUMÉRICA Primera derivada
Error
f (x i ) − f (x i −1) h 3f (x i ) − 4f (x i −1) + f (x i − 2) f ′(x i ) = 2h f ′(x i ) =
O (h) O (h2)
Segunda derivada
f (x i ) − 2f (x i −1) + f (x i − 2 ) h2 2f (x i ) − 5f (x i −1) + 4f (x i −2 ) − f (x i − 3 ) f ″ (x i ) = h2 f ″ (x i ) =
O (h) O (h2)
Tercera derivada
FIGURA 23.2 Fórmulas de diferencias divididas fi nitas hacia atrás: se presentan dos versiones para cada derivada. La última versión emplea más términos de la expansión de la serie de Taylor y, en consecuencia, es más exacta.
FIGURA 23.3 Fórmulas de diferencias divididas fi nitas centradas: se presentan dos versiones para cada derivada. La última versión emplea más términos de la expansión de la serie de Taylor y, en consecuencia, es más exacta.
f (x i ) − 3f (x i −1) + 3f (x i −2 ) − f (x i − 3 ) h3 5f (x i ) − 18f (x i −1) + 24f (x i −2 ) − 14f (x i −3 ) + 3f (x i − 4 ) f ″ (x i ) = 2h3 f ″ (x i ) =
O (h) O (h2)
Cuarta derivada
f (x i ) − 4f (x i −1) + 6f (x i −2 ) − 4f (x i −3 ) + f (x i − 4 ) h4 3f (x i ) − 14f (x i −1) + 26f (x i −2 ) − 24f (x i −3 ) + 11f (x i −4) − 2f (x i −5 ) f ″″ (x i ) = h4 f ″″ (x i ) =
O (h)
O (h2)
Primera derivada
f (x i +1) − f (x i − 1) 2h −f (x i +2 ) + 8f (x i +1) − 8f (x i −1) + f (x i − 2 ) f ′(x i ) = 12h
Error
f ′(x i ) =
O (h2)
O (h4)
Segunda derivada
f (x i +1) − 2f (x i ) + f (x i − 1) h2 −f (x i +2 ) + 16f (x i +1) − 30f (x i ) + 16f (x i −1) − f (x i − 2 ) f ″ (x i ) = 12h2 f ″ (x i ) =
O (h2)
O (h4)
Tercera derivada f (x i +2 ) − 2f (x i +1) + 2f (x i −1) − f (x i − 2 ) 2h 3 −f (x i +3 ) + 8f (x i +2) − 13f (x i +1) + 13f (x i −1) − 8f (x i −2 ) + f (x i −3 ) f ′′′(x i ) = 8h 3 f ′′′(x i ) =
O (h2)
O (h4)
Cuarta derivada
f (x i +2 ) − 4f (x i +1) + 6f (x i ) − 4f (x i −1) + f (x i − 2 ) h4 −f (x i +3 ) + 12f (x i +2 ) + 39f (x i +1) + 56f (x i ) − 39f (x i −1) +12f (x i −2 ) + f (x i −3 ) f ″″ (x i ) = 6h 4 f ″″ (x i ) =
O (h2)
O (h4)
671
23.1 FÓRMULAS DE DIFERENCIACIÓN CON ALTA EXACTITUD
EJEMPLO 23.1
Fórmulas de diferenciación con alta exactitud Planteamiento del problema. Recuerde que en el ejemplo 4.4 estimamos la derivada de f ( x ) = –0.1 x 4 – 0.15 x 3 – 0.5 x 2 – 0.25 x + 1.2
en x = 0.5 usando diferencias divididas finitas y un tamaño de paso: h = 0.25,
Estimación et (%)
Hacia adelante O (h )
Hacia atrás O (h )
Centrada O (h 2)
–1.155 –26.5
–0.714 21.7
–0.934 –2.4
donde los errores fueron calculados basándose en el valor verdadero: –0.9125. Repita este cálculo, pero ahora emplee las fórmulas con alta exactitud a partir de las figuras 23.1 a 23.3.
Solución. Los datos necesarios para este ejemplo son x i–2 = 0
f ( x i–2) = 1.2
x i–1 = 0.25
f ( x i–1) = 1.103516
x i = 0.5
f ( x i) = 0.925
x i+1 = 0.75
f ( x i+1) = 0.6363281
x i+2 = 1
f ( x i+2) = 0.2
La diferencia hacia adelante de exactitud O (h2) se calcula como sigue (figura 23.1): f ′(0.5) =
− 0.2 + 4(0.6363281) − 3(0.925 ) 2(0.25)
= −0.859375
ε t = 5.82%
La diferencia hacia atrás de exactitud O (h2) se calcula como (figura 23.2): f ′( 0.5) =
3(0 .925) − 4(1.035156) + 1 .2 = −0.878125 2(0.25)
ε t = 3.77%
La diferencia centrada de exactitud O (h 4) se calcula como (figura 23.3): f ′( 0.5) =
−0.2 + 8(0.6363281) − 8(1.035156) + 1.2 12(0.25)
= −0.9125
ε t = 0%
Como se esperaba, los errores para las diferencias hacia adelante y hacia atrás son considerablemente menores y los resultados más exactos que los del ejemplo 4.4. Sin embargo, de manera sorprendente, la diferencia centrada da un resultado perfecto. Esto es porque las fórmulas se basan en la serie de Taylor, y son equivalentes a polinomios que pasan a través de los puntos.
672
DIFERENCIACIÓN NUMÉRICA
23.2 EXTRAPOLACIÓN DE RICHARDSON Hasta aquí hemos visto que hay dos formas para mejorar la estimación obtenida al emplear diferencias divididas finitas: 1. disminuir el tamaño de paso o 2. usar una fórmula de grado superior que emplee más puntos. Un tercer procedimiento, basado en la extrapolación de Richardson, utiliza dos estimaciones de la derivada para calcular una tercera aproximación más exacta. Recuerde de la sección 22.1.1 que la extrapolación de Richardson constituye un medio para obtener una mejor estimación de la integral I por medio de la fórmula [ecuación (22.4)] 1 (23.6) I ≅ I (h2 ) + [ I (h2 ) − I (h1 )] (h1 / h2 ) 2 − 1 donde I (h1) e I (h2) son estimaciones de la integral obtenidas usando dos tamaños de paso, h1 y h2. Debido a su conveniencia cuando se expresa como un algoritmo computacional, esta fórmula usualmente se escribe para el caso en que h2 = h1 /2, como 4 1 (23.7) I ≅ I (h2) – I (h1) 3 3 De manera similar, la ecuación (23.7) se escribirá para las derivadas como 4 1 (23.8) D ≅ D(h2) – D(h1) 3 3 Para aproximaciones por diferencia centrada con O (h2), la aplicación de esta fórmula dará una nueva estimación de la derivada de O (h4).
EJEMPLO 23.2
Extrapolación de Richardson Planteamiento del problema. Utilizando la misma función que en el ejemplo 23.1, estime la primera derivada en x = 0.5 empleando tamaños de paso h1 = 0.5 y h2 = 0.25. Después, con la ecuación (23.8) calcule una mejor estimación con la extrapolación de Richardson. Recuerde que el valor verdadero es –0.9125.
Solución. Las estimaciones de la primera derivada se calculan con diferencias centradas como sigue: D(0.5) =
0.2 − 1.2 = –1.0 et = –9.6% 1
y 0.6363281 − 1.103516 = –0.934375 et = –2.4% 0.5 Se determina una mejor estimación aplicando la ecuación (23.8) al obtener D(0.25) =
4 1 ( −0.934375) – ( −1) = −0.9125 3 3 que, en este caso, es un resultado perfecto. D =
23.3 DERIVADAS DE DATOS IRREGULARMENTE ESPACIADOS
673
El ejemplo anterior dio un resultado perfecto debido a que la función analizada era un polinomio de cuarto grado. El resultado perfecto se debió al hecho de que la extrapolación de Richardson, en realidad, es equivalente a ajustar un polinomio de grado superior a los datos y después evaluar las derivadas con diferencias divididas centradas. Así, este caso concuerda, precisamente, con la derivada del polinomio de cuarto grado. Para las otras funciones que no son polinomios, por supuesto, esto no ocurrirá y nuestra estimación de la derivada será mejor, pero no perfecta. En consecuencia, como en el caso de la aplicación de la extrapolación de Richardson, el procedimiento puede aplicarse de manera iterativa usando un algoritmo de Romberg, hasta que el resultado se halle por debajo de un criterio de error aceptable.
23.3 DERIVADAS DE DATOS IRREGULARMENTE ESPACIADOS Los procedimientos analizados hasta ahora se diseñaron principalmente para determinar la derivada de una función dada. Para las aproximaciones por diferencias divididas finitas de la sección 23.1, los datos deben estar igualmente espaciados. Para la técnica de extrapolación de Richardson de la sección 23.2, los datos deben estar igualmente espaciados y generados por sucesivas divisiones a la mitad de los intervalos. Para tener un buen control del espaciamiento de datos, con frecuencia, sólo es posible cuando se utiliza una función para generar la tabla de valores. Sin embargo, la información empírica (es decir, datos a partir de experimentos o de estudios de campo) con frecuencia se obtiene a intervalos desiguales. Tal información no puede analizarse con las técnicas estudiadas hasta aquí. Una manera de emplear datos irregularmente espaciados consiste en ajustar un polinomio de interpolación de Lagrange de segundo grado [recuerde la ecuación (18.23)] a cada conjunto de tres puntos adyacentes. Recuerde que estos polinomios no requieren que los puntos estén igualmente espaciados. Si se deriva analíticamente el polinomio de segundo grado se obtiene 2 x − xi − x i +1 2 x − xi −1 − x i +1 + f ( x i ) ( xi−1 − xi )( xi −1 − x i +1 ) ( xi − xi −1 )( xi − x i +1 ) 2 x − xi −1 − x i + f ( x i +1 ) ( xi +1 − xi −1 )( xi +1 − x i )
f ′( x ) = f ( x i −1 )
(23.9)
donde x es el valor en el cual se quiere estimar la derivada. Aunque esta ecuación es más complicada que las aproximaciones de la primera derivada de las figuras 23.1 a 23.3, tiene importantes ventajas. Primero, sirve para estimar la derivada en cualquier punto dentro de un intervalo determinado por los tres puntos. Segundo, los puntos no tienen que estar igualmente espaciados y tercero, la estimación de la derivada tiene la misma exactitud que la diferencia centrada [ecuación (4.22)]. De hecho, con puntos igualmente espaciados, la ecuación (23.9) evaluada en x = x i se reduce a la ecuación (4.22).
EJEMPLO 23.3
Diferenciación de datos irregularmente espaciados Planteamiento del problema. Como se muestra en la figura 23.4, un gradiente de temperatura puede medirse abajo del suelo. El flujo de calor en la interfaz suelo-aire puede calcularse mediante la ley de Fourier,
674
DIFERENCIACIÓN NUMÉRICA
Aire
10
12
13.5 T ( C)
Suelo 1.25
3.75
z, cm
FIGURA 23.4 Temperatura contra la profundidad bajo el suelo.
dT q( z = 0) = − kρ C dz z =0
donde q = flujo de calor (W/m2), k = coeficiente de difusividad térmica en el suelo ( ≅ 3.5 × 10–7 m2 /s), r = densidad del suelo ( ≅ 1 800 kg/m3) y C = calor específico del suelo (≅ 840 J/(kg · ºC)). Observe que un valor positivo del flujo indica que el calor se transfiere del aire al suelo. Utilice diferenciación numérica para evaluar el gradiente en la interfaz suelo-aire y emplee dicha estimación para determinar el flujo de calor bajo el suelo.
Solución. La ecuación (23.9) se utiliza para calcular la derivada como sigue 2( 0) − 1.25 − 3.75 2(0) − 0 − 3.75 + 12 (0 − 1.25)(0 − 3.75) (1.25 − 0)(1.25 − 3.75) 2(0) − 0 − 1.25 +10 (3.75 − 0)(3.75 − 1.25) = −14.4 + 14.4 − 1.333333 = −1.333333 °C/cm
f ′( x ) = 13.5
que al ser sustituida se obtiene (advierta que 1 W = 1 J/s), q( z = 0) = –3.5 × 10
–7
m2 s
1 800 kg 840 J −133.3333 °C m 3 kg ⋅°C m
= 70.56 W/m2
23.4 DERIVADAS E INTEGRALES PARA DATOS CON ERRORES Además de tener espaciados irregulares, otro problema en la diferenciación de datos empíricos es que generalmente se presentan errores de medición. Una desventaja de la
675
23.4 DERIVADAS E INTEGRALES PARA DATOS CON ERRORES
diferenciación numérica es que tiende a amplificar los errores de los datos. La figura 23.5a muestra datos uniformes sin errores, que al diferenciarse numéricamente producen un resultado adecuado (figura 23.5c). En cambio, la figura 23.5b usa los mismos datos, pero con algunos puntos ligeramente por arriba y otros por abajo. Esta pequeña modificación es apenas notoria en la figura 23.5b. Sin embargo, el efecto resultante en la figura 23.5d es significativo, ya que el proceso de diferenciación amplifica los errores. Como era de esperarse, el principal procedimiento para determinar derivadas de datos imprecisos consiste en usar regresión por mínimos cuadrados, para ajustar una función suave y diferenciable a los datos. Si no se tiene alguna otra información, una regresión polinomial de grado inferior podría ser una buena elección. Obviamente, si se conoce la verdadera relación funcional entre las variables dependiente e independiente, esta relación deberá ser la base para el ajuste por mínimos cuadrados.
23.4.1 Diferenciación versus integración de datos inciertos Así como las técnicas para el ajuste de curvas, como la regresión, se utilizan para diferenciar datos inciertos, se emplea un proceso similar para la integración. No obstante, debido a la diferencia en estabilidad entre diferenciación e integración, esto rara vez se hace. Como se ilustró en la figura 23.5, la diferenciación tiende a ser inestable; es decir, amplifica los errores. En cambio, el hecho de que la integración sea un proceso de suma
FIGURA 23.5 Ilustración de cómo los pequeños errores en los datos se amplifi can mediante la diferencición numérica: a) datos sin error, b ) datos modifi cados ligeramente, c ) resultado de la diferenciación numérica que se obtiene de la curva a), y d ) la diferenciación resultante de la curva b ) que manifi esta un aumento en la variabilidad. En cambio, la operación inversa de integración [moviéndose de d ) a b ) y tomando el área bajo d )] tienden a suavizar o atenuar los errores en los datos.
y
y
t
a)
b)
D
dy dt
t
o d n a i c n e r e f i
D
dy dt
t
c)
o d n a i c n e r e f i
t
d )
676
DIFERENCIACIÓN NUMÉRICA
tiende a hacerlo muy estable respecto a datos inciertos. En esencia, conforme los puntos se suman para formar una integral, los errores aleatorios positivos y negativos tienden a compensarse. En cambio, debido a que la diferenciación es sustractiva, los errores aleatorios positivos y negativos tienden a sumarse.
23.5 INTEGRACIÓN/DIFERENCIACIÓN NUMÉRICAS CON BIBLIOTECAS Y PAQUETES DE SOFTWARE Las bibliotecas y los paquetes de software tienen muchas capacidades para la integración y la diferenciación numérica. En esta sección le daremos una muestra de algunas de las más útiles.
23.5.1 MATLAB MATLAB tiene varias funciones prediseñadas que permiten integrar y diferenciar funciones y datos. El siguiente ejemplo ilustra cómo se utilizan algunas de ellas.
EJEMPLO 23.4
Uso de MATLAB para integración y diferenciación Planteamiento del problema. Explore cómo se utiliza MATLAB para integrar y diferenciar la función f ( x ) = 0.2 + 25 x – 200 x 2 + 675 x 3 – 900 x 4 + 400 x 5
desde a = 0 hasta b = 0.8. De los capítulos 21 y 22 recuerde que el valor verdadero de la integral analíticamente se determina igual a 1.640533.
Solución. Primero, usaremos la función quad del MATLAB para integrar la función. Para usar quad, primero desarrollamos un archivo M que contendrá la función. Con un editor de textos creamos el siguiente archivo: function y=fx(x) y=0.2+25*x–200*x.^2+675*x.^3–900*x.^4+400*x.^5;
Éste se guarda en el directorio de MATLAB como fx.m. Después de entrar a MATLAB, llamamos a quad tecleando >> Q=quad(‘fx’,0,.8)
donde la segunda y tercera entradas son los límites de integración. El resultado es Q= 1.6405
Así, MATLAB proporciona una estimación exacta de la integral. Ahora investiguemos cómo se manipula en MATLAB las integrales de datos tabulados. Para ello, repetiremos el ejemplo 21.7, donde muestreamos la función en diferentes intervalos (recuerde la tabla 21.3). Es posible generar la misma información en MATLAB definiendo primero los valores de la variable independiente, >> x=[0 .12 .22 .32 .36 .4 .44 .54 .64 .7 .8];
23.5 INTEGRACIÓN/DIFERENCIACIÓN NUMÉRICAS CON BIBLIOTECAS
677
Después, se genera un vector y que contiene los valores correspondientes de la variable dependiente llamando a fx, >> y=fx(x) y = Columns 1 through 7 0.2000 1.3097 1.3052 2.8430 Columns 8 through 11 3.5073 3.1819 2.3630
1.7434
2.0749
2.4560
0.2320
Se integran estos valores llamando a la función trapz, >> integral=trapz(x,y) integral = 1.5948
Como su nombre lo indica, trapz aplica la regla del trapecio a cada intervalo y suma los resultados para obtener la integral total. Por último, se diferencian los datos irregularmente espaciados en x y y. Para ello se utiliza la función diff , que sólo determina las diferencias entre los elementos adyacentes de un vector, por ejemplo, >> diff(x) ans = Columns 1 0.1200 0.1000 Columns 8 0.1000
through 7 0.1000 0.1000
0.0400
0.0400
0.0400
through 10 0.0600 0.1000
El resultado representa las diferencias entre cada par de elementos de x . Para calcular aproximaciones por diferencias divididas de la derivada, sólo realizamos una división vectorial de las diferencias de y entre las diferencias de x tecleando >> d=diff(y)./diff(x)
que da d = Columns 1 9.2477 6.6431 Columns 8 –3.2537
through 7 –0.0449 4.3815 through 10 –13.6488
8.2877
9.5274
9.6746
–21.3100
Éstas representan estimaciones burdas de las derivadas en cada intervalo. Tal procedimiento se detallará utilizando espaciamientos más finos.
678
DIFERENCIACIÓN NUMÉRICA
23.5.3 IMSL IMSL tiene varias rutinas para la integración y la diferenciación (tabla 23.1). En el presente análisis, nos concentraremos en la rutina QDAG. Dicha rutina integra una función por medio de un esquema globalmente adaptable basado en las reglas de Gauss-Kronrod. QDAG se implementa con la siguiente declaración CALL: CALL QDAG (F, A, B, ERRABS, ERRREL, IRULE, RESULT, ERREST)
donde F = Función que introduce el usuario para que sea integrada. La forma es F(X), donde X es la variable independiente. Observe que F se debe declarar como EXTERNAL en el programa principal. A = Límite inferior de integración. (Entrada) B = Límite superior de integración. (Entrada) ERRABS = Exactitud absoluta deseada. (Entrada) ERRREL = Exactitud relativa deseada. (Entrada) IRULE = Selección de la regla de cuadratura. (Entrada). IRULE = 2 se recomienda para la mayoría de las funciones. Si la función tiene una singularidad, use IRULE = 1; si la función es oscilatoria, IRULE = 6.
TABLA 23.1 Rutinas IMSL para integrar y diferenciar. Categoría
Rutinas
Capacidad
QDAGS QDAG QDAGP QDAGI QDAWO QDAWF QDAWS QDAWC QDNG
Adaptativa de propósito general con singularidad en puntos extremos Adaptativa de propósito general Adaptativa de propósito general con puntos de singularidad Adaptativa de propósito general con intervalos infinitos Adaptativa con oscilación ponderada (trigonométrica) Adaptativa de Fourier ponderada (trigonométrica) Adaptativa algebraica ponderada con singularidad en puntos extremos Adaptativa de Cauchy ponderada con valor principal No adaptiva de propósito general
TWODQ QAND
Cuadratura bidimensional (integral iterada) Adaptativa cuadratura N-dimensional sobre un hiperrectángulo
GQRUL GQRCF RECCF RECQR FQRUL
Regla de cuadratura de Gauss para pesos clásicos Regla de cuadratura de Gauss a partir de coeficientes de recurrencia Coeficientes de recurrencia para pesos clásicos Coeficientes de recurrencia a partir de la regla de cuadratura Regla de cuadratura de Fejer
DERIV
Aproximación a la primera, segunda o tercera derivadas
Cuadratura univariada
Cuadratura multidimensional
Reglas de Gauss y recurrencias de tres términos
Diferenciación
679
PROBLEMAS
RESULT = Estimación de la integral de A a B de F (Salida) ERREST = Estimación del valor absoluto del error. (Salida)
EJEMPLO 23.5
Uso de IMSL para integrar una función Planteamiento del problema. Utilice QDAG para determinar la integral de f ( x ) = 0.2 + 25 x – 200 x 2 + 675 x 3 – 900 x 4 + 400 x 5
desde a = 0 hasta b = 0.8. De los capítulos 21 y 22 recuerde que el valor exacto de la integral analíticamente se determina igual a 1.640533.
Solución. Un ejemplo del programa principal en Fortran 90 y de una función usando QDAG para resolver este problema se describirá como sigue PROGRAM Integrate USE mimsl IMPLICIT NONE INTEGER::irule=1 REAL::a=0.,b=0.8,errabs=0.0.errrel=0.001 REAL::errest,res,f EXTERNAL f CALL QDAG (f,a,errabs,errrel,irule,res,errest) PRINT ‘(‘‘ Computed =’’,F8.4)’, res PRINT ‘(“ Error estimate =“,1PE10.3)’, errest END PROGRAM FUNCTION f(x) IMPLICIT NONE REAL::x,f f=0.2+25.*X–200.*X**2+675.*X**3.–900.*X**4+400.*X**5 END FUNCTION Output: Computed = 1.6405 Error estimate = 5.000E–05
PROBLEMAS
Calcule las aproximaciones por diferencias hacia delante y hacia atrás, de O(h) y O(h2), y aproximaciones por diferencia central de O(h2) y O(h4) para la primera derivada de y = cos x , en x = p /4, con el uso de un valor de h = p /12. Estime el error relativo porcentual verdadero et para cada aproximación. 23.1
Repita el problema 23.1, pero para y = log x evaluada en x = 25 con h = 2. 23.3 Use aproximaciones por diferencias centradas para estimar las derivadas primera y segunda de y = e x en x = 2 para h = 0.1. Emplee las dos fórmulas de O(h2) y O(h4) para hacer sus estimaciones. 23.2
680
DIFERENCIACIÓN NUMÉRICA
23.4 Emplee
la extrapolación de Richardson para estimar la primera derivada de y = cos x en x = p /4, con el uso de tamaños de paso de h1 = p /3 y h2 = p /6. Utilice diferencias centradas de O(h2) para las estimaciones iniciales. 23.5 Repita el problema 23.4, pero para la primera derivada de ln x en x = 5, con h1 = 2 y h2 = 1. 23.6 Emplee la ecuación (23.9) para determinar la primera derivada de y = 2 x 4 – 6 x 3 – 12 x – 8 en x = 0, con base en los valores de x 0 = -0.5, x 1 = 1 y x 2 = 2. Compare este resultado con el valor verdadero y con una estimación obtenida con el uso de una aproximación por diferencias centradas con base en h = 1. 23.7 Demuestre que para puntos de datos equidistantes, la ecuación (23.9) se reduce a la ecuación (4.22) en x = x i. 23.8 Calcule las aproximaciones por diferencia central de primer orden de O(h4) para cada una de las funciones siguientes en la ubicación y con el tamaño de paso que se especifica: a) b) c) d ) e)
y = x 3 + 4 x – 15 y = x 2 + cos x /3) y = tan ( x y = sen (0.5√ x )/ x x y = e + x
en x = 0, en x = 0.4, en x = 3, en x = 1, en x = 2,
h = 0.25 h = 0.1 h = 0.5 h = 0.2 h = 0.2
x f (x )
–2
–1.5
0
25
50
75
100
125
y , km
0
32
58
78
92
100
Use diferenciación numérica para estimar la velocidad y aceleración del cohete en cada momento. 23.10 Desarrolle un programa amigable para el usuario a fin de aplicar el algoritmo de Romberg para estimar la derivada de una función dada. 23.11 Desarrolle un programa amigable para el usuario, que obtenga estimaciones de la primera derivada para datos irregularmente espaciados. Pruébelo con los datos siguientes: x
1
1.5
1.6
2.5
3.5
f (x )
0.6767
0.3734
0.3261
0.08422
0.01596
donde f ( x ) = 5e x . Compare sus resultados con las derivadas verdaderas. 23.12 Recuerde que para el problema del paracaidista que desciende, la velocidad está dada por gm c
(1 − e − ( / ) ) c m t
y la distancia recorrida se obtiene con
t
∫ (1 − e
− ( c / m ) t
) dt
0
(P23.12b)
0
0.5
0.35207 0.39894
0.35207
1
1.5
2
0.24197 0.12952
0.5399
a)
Utilice MATLAB para integrar estos datos de x = –1 a 1 y –2 a 2, con la función trap. b) Emplee MATLAB para estimar los puntos de inflexión de estos datos. 23.15 Emplee IMSL para integrar la distribución normal (véase el problema 23.13) de x = –1 a 1, de –2 a 2, y de –3 a 3. 23.16 Escriba un programa en MATLAB para integrar π / 2
∫ 0
23.17
cos(cos x ) dx
Escriba un programa en MATLAB que integre 2/ π sen t
∫ 0
t
dt
con el uso de las funciones quad y quadL. Para aprender más acerca de quadL, escriba help quadL
–2 x
(t ) =
c
–0.5
0.05399 0.12952 0.24197
t , s
gm
Dadas g = 9.81, m = 70 y c = 12, a) Use MATLAB para integrar la ecuación (P23.12 a), de t = 0 a 10. b) Integre la ecuación (P23.12 b) en forma analítica, con la condición inicial de que d = 0 y t = 0. Evalúe el resultado en t = 10 para confirmar el inciso a). c) Emplee MATLAB para diferenciar la ecuación (P23.12a) en t = 10. d ) Diferencie en forma analítica la ecuación (P23.12 a) en t = 10 para confirmar el inciso c). 23.13 La distribución normal se define como 1 − x 2 / 2 f ( x ) = e 2π a) Utilice MATLAB para integrar esta función de x = –1 a 1, y de –2 a 2. b) Use MATLAB para determinar los puntos de inflexión de esta función. 23.14 Los datos siguientes se generaron a partir de la distribución normal:
–1
Compare sus resultados con las soluciones analíticas. 23.9 Para un cohete, se recabaron los datos siguientes de la distancia recorrida versus el tiempo:
v
d (t ) =
(P23.12a)
en la barra de MATLAB. 23.18 Use el comando diff(y) en MATLAB y calcule la aproximación por diferencia finita de la primera y segunda derivadas en cada valor de x de los que se muestran en la siguiente tabla, excepto los dos puntos extremos. Use aproximaciones por diferencias finitas que sean correctas en el segundo orden, O(∆ x 2). x
0
1
2
3
4
5
6
7
y
1.4
2.1
3.3
4.5
6.8
6.6
8.6
7.5
8
9
10
8.9 10.9 10
681
PROBLEMAS El objetivo de este problema es comparar las aproximaciones por diferencias finitas de segundo orden exactas hacia delante, atrás y centradas, de la primera derivada de una función con el valor real de la derivada. Esto se hará para 23.19
f ( x ) = e–2 x – x a)
Use el cálculo para determinar el valor correcto de la derivada en x = 2. b) Para evaluar las aproximaciones por diferencias finitas centradas, comience con x = 0.5. Así, para la primera evaluación, los valores de x para la aproximación por diferencias centradas será x = 2 ± 0.5 o x = 1.5 y 2.5. Entonces, disminuya en pasos de 0.01 hacia abajo hasta un valor mínimo de ∆ x = 0.01. c) Repita el inciso b) para las diferencias de segundo orden hacia delante y hacia atrás. (Observe que esto se puede hacer al mismo tiempo que la diferencia centrada se calcula en el lazo.) d ) Grafique los resultados de b) y c) versus x . Para efectos de comparación, incluya el resultado exacto en la gráfica. 23.20 Use una expansión en series de Taylor para obtener una aproximación a la tercera derivada que tenga una exactitud de segundo orden, por diferencias finitas centradas. Para hacer esto, tendrá que usar cuatro expansiones diferentes para los puntos x i-2, x i-1, x i+1 y x i+2. En cada caso, la expansión será alrededor del punto x i. El intervalo ∆ x se usará en cada caso de i –1 e i + 1, y 2∆ x se empleará en cada caso de i – 2 e i + 2. Las cuatro ecuaciones deben combinarse en forma que se elimine las derivadas primera y segunda. Utilice suficientes términos en cada expansión para evaluar el primer término que se truncará a fin de determinar el orden de la aproximación. 23.21 Use los datos siguientes para encontrar la velocidad y aceleración en t = 10 segundos: Tiempo, t , s
0
2
4
6
8
10
12
14
16
Posición, x , m
0
0.7
1.8
3.4
5.1
6.3
7.3
8.0
8.4
Emplee los métodos de diferencias finitas correctas de segundo orden a) centradas, b) hacia delante, c) hacia atrás. 23.22 Un avión es seguido por radar, y se toman datos cada segundo en coordenadas polares q y r . t , s
200
202
204
206
208
210
q ,
0.75
0.72
0.70
0.68
0.67
0.66
5 120
5 370
5 560
5 800
6 030
6 240
(rad)
r , m
A los 206 segundos, utilice diferencias finitas centradas (correctas de segundo orden) para encontrar las expresiones vectoriales para la velocidad u→, y aceleración → a. La velocidad y aceleración en coordenadas polares son:
v
⋅ r + rθ ⋅ eθ = re
⋅ 2 ) e + (rθ¨ + 2r⋅ θ ⋅) e r θ
a = (r¨ − rθ
y
23.23 Desarrolle
un programa de macros en Excel VBA para leer en columnas adyacentes de una hoja de cálculo los valores de x y y. Evalúe las derivadas en cada punto con el uso de la ecuación 23.9, y muestre los resultados en una tercera columna que se construya en la hoja, adyacente a las de los valores x y y. Pruebe su programa aplicándolo para evaluar las velocidades para los valores tiempo-posición del problema 23.21. 23.24 Use regresión para estimar la aceleración en cada momento para los datos siguientes con polinomios de segundo, tercero y cuarto orden. Grafique los resultados. t
1
2
3.25
4.5
6
7
8
8.5
9.3
10
v
10
12
11
14
17
16
12
14
14
10
Usted tiene que medir la tasa de flujo de agua a través de un tubo pequeño. Para ello, coloque una boquilla en la salida del tubo y mida el volumen a través de ella como función del tiempo, según se ha tabulado a continuación. Estime la tasa de flujo en t = 7 s. 23.25
Tiempo, s
0
1
5
8
Volumen, cm3
0
1
8
16.4
Se mide la velocidad v (m/s) del aire que fluye por una superficie plana a distintas distancias, y (m) de la superficie. Determine el esfuerzo cortante t (N/m2) en la superficie ( y = 0). 23.26
τ
= µ
d v dy
Suponga un valor de viscosidad dinámica µ = 1.8 × 10–5 N · s/m2. y , m
0
0.002
0.006 0.0012 0.018
0.024
v , m/s
0
0.287
0.899
4.299
23.27
1.915
3.048
Es frecuente que las reacciones químicas sigan este modelo: dc dt
= − kcn
donde c = concentración, t = tiempo, k = tasa de reacción, y n = orden de reacción. Es posible evaluar valores dados de c y dc/dt , k y n, por medio de regresión lineal del logaritmo de esta ecuación: log −
dc
= log k + n log c
dt
Use este enfoque y los datos que siguen para estimar los valores de k y n: t
10
20
30
40
50
60
c
3.52
2.48
1.75
1.23
0.87
0.61
