1. No category

PRÁCTICAS DE LABORATORIO­MÉTODOS NUMÉRICOS PRACTICAS DE LABORATORIO MÉTODOS NUMÉRICOS PRÁCTICA 19­20: Resumen estadística. Nociones de errores. En la práctica 19 tocaremos algunos conceptos de estadística que no hemos visto en la prácticas anteriores. En la práctica 20 haremos una introducción a lo métodos numéricos y trabajaremos con los diferente errores que se generan con el uso de aproximaciones para representar cantidades y/u operaciones, conocido con el nombre de Teoría de errores. NOTA: estás prácticas no se realizarán en el presente curso académico.
ISIDORO PONTE­E.S.M.C.127 PRÁCTICAS DE LABORATORIO­MÉTODOS NUMÉRICOS PRACTICAS DE LABORATORIO MÉTODOS NUMÉRICOS PRÁCTICA 21­22: Resolución numérica de ecuaciones. En estas prácticas aprenderemos a buscar raíces(soluciones) de ecuaciones no lineales, para ello usaremos diferentes métodos de aproximación a dichas raíces: bisección, regla falsa, Newton­ Raphson, método de la secante y punto fijo. Trataremos también en algunas de las prácticas siguientes de crear “pequeños programas”(explicándose convenientemente) que efectúen automáticamente los cálculos aproximativos, con el simple cambio de los datos y, por otra parte, cuando el problema no tenga solución nos lo advierta. PROGRAMA DE BISECCIÓN. Este programa utiliza el método de la bisección para aproximar la raíz de f ( x ) en un intervalo [ a , b ] ( recordemos que este método se basa en el teorema de valor medio Sea f ( x ) continua en un intervalo [a , b ] y supongamos que f ( a ) < f ( b ) . Entonces para cada z tal que f ( a ) < z < f ( b ) , existe un x 0 Î (a , b ) tal que f ( x 0 ) = z . La misma conclusión se obtiene para el caso que f ( a ) > f ( b ) . Básicamente el Teorema del Valor medio nos dice que toda función continua en un intervalo cerrado, una vez que alcanzó ciertos valores en los extremos del intervalo, entonces debe alcanzar todos los valores intermedios. En particular, si f (a ) y f (b ) tienen signos opuestos, entonces un valor intermedio es precisamente z = 0 , y por lo tanto, el Teorema del Valor Intermedio nos asegura que debe existir x 0 Î (a , b ) tal que f ( x 0 ) = 0 , es decir, debe haber por lo menos una raíz de f ( x ) en el intervalo ( a , b ) .
ISIDORO PONTE­E.S.M.C.128 PRÁCTICAS DE LABORATORIO­MÉTODOS NUMÉRICOS El método de bisección sigue los siguientes pasos: Sea f ( x ) continua, i) Encontrar valores iniciales x a , x b tales que f ( x a ) y f ( x b ) tienen signos opuestos, es decir, ii) La primera aproximación a la raíz se toma igual al punto medio entre x a y x b : iii) Evaluar f ( x r ) . Forzosamente debemos caer en uno de los siguientes casos: En este caso, tenemos que f ( x a ) y f ( x r ) tienen signos opuestos, y por lo tanto la raíz se encuentra en el intervalo [x a , x r ] . Y asi sucesivamente hasta donde queramos aproximar) Veamos un ejemplo con MATHEMATICA; introducimos la función In[1]:= f@x_D =
Exp@-xD - Log@xD; introducimos los extremos del intervalo In[2]:= x0 = 1; y 0 = 2; introducimos el numero de iteraciones In[4]:= k =
8; introducimos el error inicial en porcentaje In[5]:= e0 = 100;
ISIDORO PONTE­E.S.M.C.129 PRÁCTICAS DE LABORATORIO­MÉTODOS NUMÉRICOS calculamos la primera aproximación a la raíz In[6]:= z0 = Hx0 + y0L • 2; le decimos a matemática que use la fórmula de la bisección en el caso que proceda y en este caso, que nos despliegue una tabla con las aproximaciones calculadas. En caso negativo que nos mande un mensaje de error. En la tabla aparece la numeración de las raíces, el valor y el error In[7]:= IfASign@f@x0DD
!= Sign@f@y0DD, ForAi = 0, i < k, i++, IfASign@f@xiDD != Sign@f@ziDD, xi+1 = x i ; yi+1 = zi; zi+1 = Hxi+1 + yi+1L • 2, xi+1 = z i; yi+1 = yi; zi+1 = Hxi+1 + yi+1L • 2EE; Do@ei+1 = Abs@Hzi+1 - ziL • zi+1D * 100, 8i, 0, k<D; Table@8xri+1 , N@zi, 20D, N@ei, 20D<, 8i, 0, k<D ••
TableForm, Print@"El método de la bisección no se puede aplicar ya que fHx0L y fHy0L tienen el mismo signo"DE
Out[7]//TableForm= xr1 xr2 xr3 xr4 xr5 xr6 xr7 xr8 xr9 1.5000000000000000000 1.2500000000000000000 1.3750000000000000000 1.3125000000000000000 1.2812500000000000000 1.2968750000000000000 1.3046875000000000000 1.3085937500000000000 1.3105468750000000000 100.00000000000000000 20.000000000000000000 9.0909090909090909091 4.7619047619047619048 2.4390243902439024390 1.2048192771084337349 0.59880239520958083832 0.29850746268656716418 0.14903129657228017884
ISIDORO PONTE­E.S.M.C.130 PRÁCTICAS DE LABORATORIO­MÉTODOS NUMÉRICOS podemos comprobar que si la función no tiene raíces , MATHEMATICA nos lo advierte veamos un ejemplo como el anterior In[9]:= f@x_D =
Exp@-xD - Log@xD; In[10]:= x0 = 2; y0 = 3; In[12]:= k = 8; In[13]:= e 0 = 100; In[14]:= z0 = Hx0 + y0L • 2; In[15]:= IfASign@f@x 0DD
!= Sign@f@y 0DD, ForAi = 0, i < k, i++, IfASign@f@xiDD != Sign@f@ziDD, xi+1 = x i ; yi+1 = zi; zi+1 = Hxi+1 + yi+1L • 2, xi+1 = zi; yi+1 = yi; zi+1 = Hx i+1 + y i+1L • 2EE; Do@ei+1 = Abs@Hzi+1 - ziL • zi+1D * 100, 8i, 0, k<D; Table@8xri+1 , N@zi, 20D, N@ei, 20D<, 8i, 0, k<D ••
TableForm, Print@"El método de la bisección no se puede aplicar ya que fHx0L y fHy0L tienen el mismo signo"DE
El método de la bisección no se puede aplicar ya que fHx0L y fHy0L tienen el mismo signo estos datos los podríamos comprobar gráficamente, usando el comando que dibuja las funciones, y ver que hay una raíz en el intervalo [1 , 2 ] y ninguna en el intervalo [ 2 , 3 ] ISIDORO PONTE­E.S.M.C.131 PRÁCTICAS DE LABORATORIO­MÉTODOS NUMÉRICOS In[16]:= Plot@8 Exp@-xD -
Log@xD<, 8x, 1, 3<D
0.4 0.2 1.5 2 2.5 3 ­0.2 ­0.4 ­0.6 ­0.8 ­1 Out[16]= … Graphics …
METODO DE LA REGLA FALSA También este programa se usa para aproximar raíces de f ( x ) en un intervalo [ a , b ] ( es un método bueno para considerar si la raíz se encuentra cerca de los extremos. Consideremos una gráfica del tipo: Donde hemos agregado la línea recta que une los puntos extremos de la gráfica en el intervalo [a , b ] . Es claro que si en lugar de considerar el punto medio del intervalo, tomamos el punto donde cruza al eje x esta recta, nos aproximaremos mucho más rápido a la raíz; ésta es en sí, la idea central del método de la regla falsa y ésta es realmente la única diferencia con el método de bisección, puesto que en todo lo demás los dos métodos son prácticamente idénticos.
ISIDORO PONTE­E.S.M.C.132 PRÁCTICAS DE LABORATORIO­MÉTODOS NUMÉRICOS Supongamos que tenemos una función f ( x ) que es contínua en el intervalo [x a , x b ] y además, f ( x a ) y f ( x b ) tienen signos opuestos. Calculemos la ecuación de la línea recta que une los puntos ( x a , f ( x a )) , ( x b , f ( x b )) . Sabemos que la pendiente de esta recta esta dada por: Por lo tanto la ecuación de la recta es: Para obtener el cruce con el eje x , hacemos y = 0 : Multiplicando por xb - x a nos da: Finalmente, de aquí despejamos x : Este punto es el que toma el papel de x r en lugar del punto medio del método de bisección. Así pues, el método de la regla falsa sigue los siguientes pasos: Sea f ( x ) contínua, Encontrar valores iniciales x a , x b tales que f ( x a ) y f ( x b ) tienen signos opuestos, es decir,
ISIDORO PONTE­E.S.M.C.133 PRÁCTICAS DE LABORATORIO­MÉTODOS NUMÉRICOS ii) La primera aproximación a la raíz se toma igual a: iii) Evaluar f ( x r ) . Forzosamente debemos caer en uno de los siguientes casos: En este caso, tenemos que f ( x a ) y f ( x r ) tienen signos opuestos, y por lo tanto la raíz se encuentra en el intervalo [x a , x r ] . En este caso, tenemos que f ( x a ) y f ( x r ) tienen el mismo signo, y de aquí que f ( x r ) y f ( x b ) tienen signos opuestos. Por lo tanto, la raíz se encuentra en el intervalo [xr , x b ] . En este caso se tiene que f ( x r ) = 0 y por lo tanto ya localizamos la raíz. El proceso se vuelve a repetir con el nuevo intervalo, hasta que: Resolvamos con MATHEMATICA un ejercicio Introducimos la función In[1]:= f@x_D =
Sin@xD - 0.5 x; introducimos los extremos del intervalo
ISIDORO PONTE­E.S.M.C.134 PRÁCTICAS DE LABORATORIO­MÉTODOS NUMÉRICOS In[2]:= x0 = 1; y 0 = 2; introducimos el numero de iteraciones In[4]:= k =
3; introducimos el error inicial en porcentaje In[5]:= e0 = 100; calculamos la primera aproximación a la raíz In[6]:= z0 = y0 - f@y0D * Hx0 - y0L • Hf@x0D - f@y0DL; le decimos a matemática que verifique el método. En caso afirmativo, le pedimos que aplique reiteradamente hasta el número de iteraciones introducidas, y que nos despliegue una tabla con los datos obtenidos. En caso negativo que nos mande un mensaje de error. En la tabla aparece la numeración de las raíces, el valor y el error In[7]:= If@Sign@f@x0DD
¹ Sign@f@y0DD, For@i = 0, i < k, i++, If@Sign@f@xiDD ¹ Sign@f@ziDD, xi+1 = xi; yi+1 = zi; zi+1 = yi+1 - f@yi+1D * Hxi+1 - yi+1L • Hf@xi+1D - f@yi+1DL, xi+1 = zi; yi+1 = yi; zi+1 = yi+1 - f@yi+1D * Hxi+1 - yi+1L • Hf@xi+1D - f@yi+1DLDD; Do@ei+1 = Abs@Hzi+1 - ziL • zi+1D * 100, 8i, 0, k<D; Table@8xri+1 , N@zi, 7D, N@ei, 7D<, 8i, 0, k<D •• TableForm, Print@"El metodo de la regla falsa no se puede aplicar ya que fHx0L y fHy0L tienen el mismo signo"DD
Out[7]//TableForm= xr 1 xr 2 xr 3 xr 4 1.79012 1.88912 1.89513 1.89547 100. 5.24032 0.317292 0.0179524
ISIDORO PONTE­E.S.M.C.135 PRÁCTICAS DE LABORATORIO­MÉTODOS NUMÉRICOS METODO DE NEWTON­RAPHSON Este método, el cual es un método iterativo, es uno de los más usados y efectivos. A diferencia de los métodos anteriores, el método de Newton­Raphson no trabaja sobre un intervalo sino que basa su fórmula en un proceso iterativo. Supongamos que tenemos la aproximación f ( x ) , x i a la raíz x r de Trazamos la recta tangente a la curva en el punto ( x i , f ( x i ) ) ; ésta cruza al eje x en un punto aproximación a la raíz x r . x i + 1 que será nuestra siguiente Para calcular el punto x i + 1 , calculamos primero la ecuación de la recta tangente. Sabemos que tiene pendiente Y por lo tanto la ecuación de la recta tangente es: Hacemos y = 0 : Y despejamos x :
ISIDORO PONTE­E.S.M.C.136 PRÁCTICAS DE LABORATORIO­MÉTODOS NUMÉRICOS Que es la fómula iterativa de Newton­Raphson para calcular la siguiente aproximación: , si Note que el método de Newton­Raphson no trabaja con intervalos donde nos asegure que encontraremos la raíz, y de hecho no tenemos ninguna garantía de que nos aproximaremos a dicha raíz. Desde luego, existen ejemplos donde este método no converge a la raíz, en cuyo caso se dice que el método diverge. Sin embargo, en los casos donde si converge a la raíz lo hace con una rapidez impresionante, por lo cual es uno de los métodos preferidos por excelencia. También observe que en el caso de que f ¢( x i ) = 0 , el método no se puede aplicar. De hecho, vemos geométricamente que esto significa que la recta tangente es horizontal y por lo tanto no intersecta al eje x en ningún punto, a menos que coincida con éste, en cuyo caso x i mismo es una raíz de f ( x ) . Veamos un ejemplo con MATHEMATICA; introducimos la función In[1]:= f@x_D =
Exp@-xD - Log@xD; calculamos la derivada de nuestra función f ( x ) y la llamamos g ( x ) In[2]:= g@x_D =
D@f@xD, xD; introducimos el valor inicial para el proceso iterativo In[3]:= x0 = 1; introducimos el número de iteraciones In[4]:= k =
4; introducimos un error inicial, digamos del 100%
ISIDORO PONTE­E.S.M.C.137 PRÁCTICAS DE LABORATORIO­MÉTODOS NUMÉRICOS In[5]:= e0 = 100; le pedimos a MATHEMATICA que aplique la fórmula de Newton­ Raphson In[6]:= DoAxi+1 = x i - f@xiD “
gAx i E, 8i, 0, k<E; también le pedimos que calcule los errores aproximados y que nos despliegue una tabla en la cual aparezcan el número de aproximación, el valor de la aproximación y el error aproximado correspondiente In[7]:= Do@ei+1 = Abs@Hxi+1 - xiL • xi+1D * 100, 8i, 0, k<D; Table@8x , N@x ri i, 20D, N@ei, 20D<, 8i, 0, k + 1<D ••
TableForm Out[8]//TableForm= xr0 xr1 xr2 xr3 xr4 1.0000000000000000000 1.2689414213699951207 1.3091084032740159137 1.3097993886689735953 1.3097995858041344422 100.00000000000000000 21.194155761708544507 3.0682701144966407035 0.052755055540212554120 0.000015050788149843876043 xr 5 1.3097995858041504777 1.2242709553333811460 ´ 10-12 METODO DE LA SECANTE Este método se basa en la fórmula de Newton­Raphson, pero evita el cálculo de la derivada usando la siguiente aproximación: Sustituyendo en la fórmula de Newton­Raphson, obtenemos:
ISIDORO PONTE­E.S.M.C.138 PRÁCTICAS DE LABORATORIO­MÉTODOS NUMÉRICOS Que es la fórmula del método de la secante. Nótese que para poder calcular el valor de , necesitamos conocer los dos valores anteriores y . Obsérvese también, el gran parecido con la fórmula del método de la regla falsa. La diferencia entre una y otra es que mientras el método de la regla falsa trabaja sobre intervalos cerrados, el método de la secante es un proceso iterativo y por lo mismo, encuentra la aproximación casi con la misma rapidez que el método de Newton­Raphson. Claro, corre el mismo riesgo de éste último de no converger a la raíz, mientras que el método de la regla falsa va a la segura. Introducimos la función In[1]:= f@x_D =
ArcSin@xD - Exp@-xD; Introducimos los dos valores iniciales In[2]:= x0 = 0; In[3]:= x1 = 0.5; el número de iteraciones In[4]:= k =
4; el error inicial In[5]:= e0 = 100; usamos un iterador para que MATHEMATICA aplique la fórmula de la secante In[6]:= Do@xn+1 = xn - f@xnD * Hxn - xn-1L • Hf@xnD - f@xn-1DL, 8n, 1, k<D; ahora finalmente nos presentará la tabla en la que aparecerán las aproximaciones y los errores aproximados
ISIDORO PONTE­E.S.M.C.139 PRÁCTICAS DE LABORATORIO­MÉTODOS NUMÉRICOS In[7]:= Do@ei+1 = Abs@Hxi+1 - xiL • xi+1D * 100, 8i, 0, k<D; Table@8x , N@x ri i, 10D, N@f@x iD, 10D, N@ei, 10D<, 8i, 0, k<D ••
TableForm Out[8]//TableForm= xr 0 xr 1 xr 2 xr 3 0.
0.5
0.545216
0.546954 -1. 0.0000116128 100. 100. 8.29319 0.31781 xr 4 0.546947
-1.99027 ´ 10-9 0.00119758 -0.0829319 -0.0030702 METODO DEL PUNTO FIJO Este programa usa el método de iteración del punto fijo para aproximar la raíz de una ecuación. Debemos de introducir una función g ( x ) tal que se desee resolver la ecuación g ( x ) = x . Introducimos la función In[1]:= g@x_D = H4 -
Log@xDL • 2; Introducimos el valor inicial In[2]:= x0 = 1.5; el número de iteraciones In[3]:= k =
4; el error inicial In[4]:= e0 = 100; usamos un iterador para que MATHEMATICA aplique la fórmula del punto fijo In[5]:= Do@xi+1 = g@xiD, 8i, 0, k<D; ahora finalmente nos presentará la tabla en la que aparecerán las aproximaciones y los errores aproximados
ISIDORO PONTE­E.S.M.C.140 PRÁCTICAS DE LABORATORIO­MÉTODOS NUMÉRICOS In[6]:= Do@e i+1 = Abs@Hx i+1 - x iL • xi+1D * 100, 8i, 0, k<D; Table@8x , N@x ri i, 8D, N@e i, 8D<, 8i, 0, k<D •• TableForm Out[7]//TableForm= xr0 xr1 xr2 xr3 xr4 1.5 1.79727 1.70687 1.73267 1.72517 100. 16.54 5.29632 1.48927 0.434877
ISIDORO PONTE­E.S.M.C.141 PRÁCTICAS DE LABORATORIO­MÉTODOS NUMÉRICOS Ejercicios complementarios personalizados 1A Usa el método de la bisección para aproximar la raíz de f ( x ) = x 2 +
1 - tg x comenzando en el intervalo [0 . 75 , 1 ] y hasta d 5 + 1 que el error sea menor que el 1%. 1B Usa el método de la regla falsa para aproximar la raíz de f ( x ) = 4 - x 2 - x 3 comenzando en el intervalo [1 + 0 . d 4 , 2 + 0 . d 7 ] y hasta que el error sea menor que el 1%. 2A Usa el método de la Newton­Raphson para aproximar la raíz de f ( x ) = 1 . d 2 - x 2 - arc tg x comenzando con x 0 = 0 . 5 y hasta que el error sea menor que el 1%. 2B Usa el método de la secante para aproximar la raíz de f ( x ) = e - x - x comenzando con x0 = 0 . d 1 y x1 = 0 . d 3 hasta que el error sea menor que el 1%. 2C Usa el método de iteración del punto fijo para aproximar la raíz de f ( x ) = sen x + x - 1 comenzando con x0 = 0 . 5 d 8 y hasta que el error sea menor que el 1%.
ISIDORO PONTE­E.S.M.C.142 PRÁCTICAS DE LABORATORIO­MÉTODOS NUMÉRICOS PRACTICAS DE LABORATORIO MÉTODOS NUMÉRICOS PRÁCTICA 23­24: Interpolación. En estas prácticas estudiaremos el importantísimo tema de la interpolación de datos. Veremos dos tipos de interpolación: la interpolación polinomial (a la que dedicaremos casi todas las prácticas) y la interpolación segmentaria (splines). Comencemos dando la definición general. Definición. Dados n + 1 puntos que corresponden a los datos: y los cuales se representan gráficamente como puntos en el plano cartesiano, Si existe una función f ( x ) definida en el intervalo [x0 , x n ] (donde x0 < x 1 < L < x n ), tal que f ( x i ) = y i para suponemos que i = 0 , 1 , 2 , L, n , entonces a f ( x ) se le llama una función de interpolación de los datos, cuando es usada para aproximar valores dentro del intervalo [x 0 , x n ] , y se le llama función de extrapolación de los datos, cuando está definida y es usada para aproximar valores fuera del intervalo.
ISIDORO PONTE­E.S.M.C.143 PRÁCTICAS DE LABORATORIO­MÉTODOS NUMÉRICOS Evidentemente pueden existir varios tipos de funciones que interpolen los mismos datos; por ejemplo, funciones trigonométricas, funciones exponenciales, funciones polinomiales, combinaciones de éstas, etc. El tipo de interpolación que uno elige, depende generalmente de la naturaleza de los datos que se están manejando, así como de los valores intermedios que se están esperando. Un tipo muy importante es la interpolación por funciones polinomiales. Puesto que evidentemente pueden existir una infinidad de funciones polinomiales de interpolación para una misma tabla de datos, se hace una petición extra para que el polinomio de interpolación, sea único. Definición. Un polinomio de interpolación es una función polinomial que además de interpolar los datos, es el de menor grado posible. DIFERENCIAS DIVIDIDAS FINITAS DE NEWTON Las diferencias divididas finitas de Newton, se define de la siguiente manera: f [ x i , x j ] = f ( x i ) - f ( x j ) x i - x j f [ x i , x j , x k ] = f [ x i , x j ] - f [ x j , x k ] x i - x k · ·
·
f [ x n , x n -1 , L , x 1 , x 0 ] =
f [ x n , L , x 1 ] - f [ x n -1 , L, x 0 ] x n - x 0 ISIDORO PONTE­E.S.M.C.144 PRÁCTICAS DE LABORATORIO­MÉTODOS NUMÉRICOS A manera de ejemplo citemos el siguiente caso específico : f [ x 3 , x 2 , x 1 , x 0 ] =
f [ x 3 , x 2 , x 1 ] - f [ x 2 , x 1 , x 0 ] x 3 - x 0 donde a su vez: f [ x 3 , x 2 , x 1 ] =
f [ x 3 , x 2 ] - f [ x 2 , x 1 ] x 3 - x 1 y f [ x 2 , x 1 , x 0 ] =
f [ x 2 , x 1 ] - f [ x 1 , x 0 ] x 2 - x 01 Y donde a su vez: f [ x 3 , x 2 ] =
f ( x 3 ) - f ( x 2 ) x 3 - x 2 etc. Podemos ahora definir nuestro primer tipo de polinomio de interpolación. POLINOMIO DE INTERPOLACIÓN DIFERENCIAS DIVIDIDAS DE NEWTON CON Dados n + 1 datos: ­ El polinomio de interpolación de Newton se define de la siguiente manera:
f ( x ) = b 0 + b 1 ( x - x 0 ) + b 2 ( x - x 0 )( x - x 1 ) + L + b n (x - x 0 )( x - x 1 )L ( x - x n -1 ) donde :
b0 = f ( x 0 ) b1 = f [ x 1 , x 0 ] b2 = f [x 2 , x 1 , x 0 ] M
b n = f [x n ,L
, x 0 ] ISIDORO PONTE­E.S.M.C.145 PRÁCTICAS DE LABORATORIO­MÉTODOS NUMÉRICOS Para calcular los coeficientes b0 , b 1 , L , b n , es conveniente construir una tabla de diferencias divididas como la siguiente : Obsérvese que los coeficientes del polinomio de interpolación de Newton, se encuentran en la parte superior de la tabla de diferencias divididas. Veamos un ejemplo con MATHEMATICA, para una tabla de datos con cinco puntos. (2, 8 ), (- 2 , 16 ), (- 6 , 10 ), (4 , -2 ), (3 , 0 . 4055 ) Introducimos primero los valores de las abscisas y luego los de las ordenadas In[1]:= x 0 = 2; x1 = -2; x 2 = -6; x3 = 4; x 4 = 3; In[6]:= y0 = 8; y1 = 16; y2 = 10; y3 = -2; y4 = 0.4055; Ahora introducimos el grado del polinomio de interpolación que buscamos In[11]:= n = 4; calculamos la tabla de diferencias divididas finitas de Newton
ISIDORO PONTE­E.S.M.C.146 PRÁCTICAS DE LABORATORIO­MÉTODOS NUMÉRICOS In[12]:= Do@B i = Hy i - y i-1L•Hx i - x i-1L, 8i, 1, n<D; Do@C i = HB i+1 - BiL•Hx i+1 - xi-1L, 8i, 1, n - 1<D; Do@Di = HCi+1 - CiL•Hxi+2 - xi-1L, 8i, 1, n - 2<D; Do@Ei = HDi+1 - DiL•Hxi+3 - xi-1L, 8i, 1, n - 3<D; calculamos el polinomio de interpolación de newton correspondiente a la tabla de datos In[16]:= f@x_D = y 0 + B 1 * Hx - x 0L + C1 * Hx - x 0L * Hx - x1L +
D1* Hx - x0L * Hx- x1L * Hx - x2L +
E 0L * Hx- x 1L * Hx - x 2L * Hx - x 3L; 1 * Hx - x por fin , pedimos a MATHEMATICA que nos escriba en pantalla el polinomio y que nos haga una gráfica con los puntos dados y el polinomio de interpolación de Newton obtenido
ISIDORO PONTE­E.S.M.C.147 PRÁCTICAS DE LABORATORIO­MÉTODOS NUMÉRICOS In[17]:= Print@"El polinomio de Newton para los datos dados es:"D
f@xD
Print@
"la grafica de los puntos y el polinomio de newton se ven como sigue:"D
datos:= 88x 0, y 0<, 8x1, y1<, 8x 2, y 2<, 8x 3, y 3<, 8x4, y4<<; puntos:=
ListPlot@datos, PlotStyle -> [email protected], DisplayFunction -> IdentityD; curva:= Plot@f@xD, 8x, Min@x 0, x 1, x2, x3, x4D, Max@x 0, x 1, x 2, x 3, x 4D<, DisplayFunction -> IdentityD; Show@puntos, curva, DisplayFunction ® $DisplayFunctionD; El polinomio de Newton para los datos dados es: 7 1
H- 2 + xL H2 + xL H- 2 + xL H2 + xL H6 + xL +
16
160
0.0694611 H-4 + xL H- 2 + xL H2 + xL H6 + xL
Out[18]= 8 - 2 H- 2 + xL -
la grafica de los puntos y el polinomio de newton se ven como sigue: 20 15 10 5 ­6 ­4 ­2 2 4 POLINOMIO DE INTERPOLACION DE LAGRANGE. Nuevamente tenemos los datos : El polinomio de interpolación de Lagrange se plantea como sigue:
ISIDORO PONTE­E.S.M.C.148 PRÁCTICAS DE LABORATORIO­MÉTODOS NUMÉRICOS P ( x ) = y 0l 0 ( x ) + y 1 l 1 ( x ) + L + y nl n ( x ) Donde los polinomios l i ( x ) se llaman los polinomios de Lagrange, correspondientes a la tabla de datos. Como se debe satisfacer que P ( x 0 ) = y 0 , esto se cumple si l 0 ( x 0 ) = 1 y l i ( x 0 ) = 0 para toda i ¹ 0 . Como se debe satisfacer que P ( x 1 ) = y 1 , esto se cumple si l 1 ( x 1 ) = 1 y l i ( x 1 ) = 0 para toda i ¹ 1 . Y así sucesivamente, veremos finalmente que la condición Pn ( x n ) = y n se cumple si l n ( x n ) = 1 y l i ( x n ) = 0 para toda i ¹ n . Esto nos sugiere como plantear los polinomios de Lagrange. Para ser más claros, analicemos detenidamente el polinomio l 0 ( x ) . De acuerdo al análisis anterior vemos que deben cumplirse las siguientes condiciones para l 0 ( x ) : l 0 ( x 0 ) = 1 y l 0 ( x j ) = 0 , para toda j ¹ 0
Por lo tanto, planteamos l 0 ( x ) como sigue:
lo ( x ) = c ( x - x 1 )( x - x 2 )L ( x - x n ) Con esto se cumple la segunda condición sobre l 0 ( x ) . La constante c se determinará para hacer que se cumpla la primera condición:
l0 ( x 0 ) = 1 Þ 1 = c ( x 0 - x 1 )( x 0 - x 2 )L( x 0 - x n ) Þ c =
1 ( x 0 - x 1 )( x 0 - x 2 )L ( x 0 - x n ) Por lo tanto el polinomio l 0 ( x ) queda definido como:
l0 ( x ) = ( x - x 1 )( x - x 2 )L ( x - x n )
( x 0 - x 1 )( x 0 - x 2 )L ( x 0 - x n ) Análogamente se puede deducir que:
l j ( x ) = ( x - x i ) Õ
i j ¹
( x j - x i ) Õ
i j ¹
, para j = 1, K, n ISIDORO PONTE­E.S.M.C.149 PRÁCTICAS DE LABORATORIO­MÉTODOS NUMÉRICOS Veamos un ejemplo con MATHEMATICA, para una tabla de datos con cuatro puntos. (- 3, 112 ), (- 1 , -120 ), (4 , -70 ), (5 , 0 ) Introducimos primero los valores de las abscisas y luego los de las ordenadas, In[1]:= x 0 = -3; x 1 = -1; x 2 = 4; x 3 = 5; In[5]:= y0 = 112; y1 = -120; y2 = -70; y3 = 0; introducimos el grado del polinomio In[9]:= n = 3; calculamos los polinomios de Lagrange In[10]:= L 0@x_D =
L 1@x_D =
L2@x_D =
L3@x_D =
Hx - x 1L * Hx - x 2L * Hx- x 3L
Hx 0 - x 1L * Hx0 - x 2L * Hx 0 - x3L
Hx - x 0L * Hx - x 2L * Hx- x 3L
Hx 1 - x 0L * Hx1 - x 2L * Hx 1 - x3L
Hx - x 0L * Hx - x 1L * Hx- x 3L
Hx 2 - x 0L * Hx 2 - x 1L * Hx 2 - x 3L
Hx - x 0L * Hx - x 1L * Hx- x 2L
Hx 3 - x 0L * Hx 3 - x 1L * Hx 3 - x 2L
; ; ; ; calculamos el polinomio de interpolación correspondiente a la tabla de datos de Lagrange n In[14]:= p@x_D = âyi* L i@xD; i=0 finalmente pedimos a MATHEMATICA que nos escriba en pantalla el polinomio de interpolación de Lagrange, así como una gráfica
ISIDORO PONTE­E.S.M.C.150 PRÁCTICAS DE LABORATORIO­MÉTODOS NUMÉRICOS donde veamos los puntos de la tabla de datos junto con el polinomio de Lagrange In[15]:= Print@"El polinomio de Lagrange para los datos dados es:"D
p@xD
Print@
"La grafica de los puntos y el polinomio de interpolacion de lagrange se ven como sigue:"D
datos:= 88x0, y 0<, 8x 1, y 1<, 8x 2, y 2<, 8x3, y3<<; puntos:=
ListPlot@datos, PlotStyle -> [email protected], DisplayFunction -> IdentityD; curva:= Plot@p@xD, 8x, Min@x 0, x 1, x 2, x 3D, Max@x 0, x 1, x 2, x 3D<, DisplayFunction -> IdentityD; Show@puntos, curva, DisplayFunction ® $DisplayFunctionD; El polinomio de Lagrange para los datos dados es: Out[16]= -H- 5 + xL H-4 + xL H1 + xL -
2 H-5 + xL H- 4 + xL H3 + xL + 2 H- 5 + xL H1 + xL H3 + xL
La grafica de los puntos y el polinomio de interpolacion de lagrange se ven como sigue: 100 50 ­2 2 4 ­50 ­100 ­150 INTERPOLACION POR SPLINES CUBICAS. Terminamos este capítulo, estudiando un tipo de interpolación que ha demostrado poseer una gran finura, y que inclusive es usado para el diseño por computadora, por ejemplo, de tipos de letra.
ISIDORO PONTE­E.S.M.C.151 PRÁCTICAS DE LABORATORIO­MÉTODOS NUMÉRICOS Esta interpolación se llama interpolación segmentaria o interpolación por splines. La idea central es que en vez de usar un solo polinomio para interpolar los datos, podemos usar segmentos de polinomios y unirlos adecuadamente para formar nuestra interpolación. Cabe mencionar que entre todas, las splines cúbicas han resultado ser las más adecuadas para aplicaciones como la mencionada anteriormente. Así pues, podemos decir de manera informal, que una funcion spline está formada por varios polinomios, cada uno definido en un intervalo y que se unen entre si bajo ciertas condiciones de continuidad. Definición. (Splines de grado k) Dada nuestra tabla de datos, donde suponemos que x0 < x 1 < L < x n , y dado k un número entero positivo, una función de interpolación spline de grado k, para la tabla de datos, es una función s ( x ) tal que : i) s( x i ) = y i , para toda i = 0 , 1 , K, n . ii) s ( x ) es un polinomio de grado £ k en cada subintervalo [xi -1 , x i ] . iii ) s ( x ) tiene derivada continua hasta de orden k - 1 en [x 0 , x n ] . FUNCIONES SPLINES CUBICAS Para hacer más firme el entendimiento, escribimos la definición correspondiente a este caso (k=3). Dados los n + 1 datos: Una spline cúbica que interpola estos datos, es una función s ( x ) definida como sigue :
ISIDORO PONTE­E.S.M.C.152 PRÁCTICAS DE LABORATORIO­MÉTODOS NUMÉRICOS ì s 0 ( x ) si x Î [x 0 , x 1 ]
ï s ( x ) si x Î [x , x ]
ï
1 2 s ( x ) = í 1 M ï
ïîs n -1 ( x ) si x Î [x n -1 , x n ]
donde cada s i (x ) es un polinomio cúbico; si ( x i ) = y i , para toda i = 0 , 1 , K, n y tal que s ( x ) tiene primera y segunda derivadas contínuas en [x0 , x n ] . Ejemplo. Interpolar los siguientes datos mediante una spline cúbica : Solución. Definimos un polinomio cúbico en cada uno de los intervalos que se forman:
ì a 1 x 3 + b 1 x 2 + c 1 x + d 1 s ( x ) = í
3 2 îa 2 x + b 2 x + c 2 x + d 2 si x Î [2 , 3 ]
si x Î [3 , 5 ]
A continuación, hacemos que se cumpla la condición de que la spline debe pasar por los puntos dados en la tabla. Así, tenemos que:
s(2 ) = -1 Þ 8 a 1 + 4 b 1 + 2 c 1 + d 1 = -1 s(3 ) = 2 Þ 27 a 1 + 9 b 1 + 3 c 1 + d 1 = 2 s (5 ) = -7 Þ 125 a 2 + 25 b 2 + 5 c 2 + d 2 = -7 Ahora calculamos la primera derivada de s ( x ) :
ì 3 a 1 x 2 + 2 b 1 x + c 1 s ¢ ( x ) = í
2 î3 a 2 x + 2 b 2 x + c 2 si x Î [2 , 3 ]
si x Î [3 , 5 ]
Al igual que en el caso de las splines cuadráticas, se presentan ecuaciones que pueden presentar discontinuidad en los cambios de intervalo; las posibles discontinuidades son los puntos donde se
ISIDORO PONTE­E.S.M.C.153
PRÁCTICAS DE LABORATORIO­MÉTODOS NUMÉRICOS x = 3 . Para evitar esta cambia de intervalo, en este caso discontinuidad, evaluamos x = 3 en los dos polinomios e igualamos:
3a 1 (3 ) + 2 b 1 (3 ) + c 1 = 3 a 2 (3 ) + 2 b 2 (3 ) + c 2 2 2 o lo que es lo mismo: 27a 1 + 6 b 1 + c 1 = 27 a 2 + 6 b 2 + c 2 Análogamente procedemos con la segunda derivada :
ì 6 a 1 x + 2 b 1 si x Î [2 , 3 ]
î6 a 2 x + 2 b 2 si x Î [3 , 5 ]
s ¢ ¢( x ) = í
Para lograr que s ¢¢( x ) sea continua :
6 a1 (3 ) + 2 b 1 = 6 a 2 (3 ) + 2 b 2 \18a 1 + 2 b 1 = 18 a 2 + 2 b 2 En este punto contamos con 6 ecuaciones y 8 incognitas, por lo tanto tenemos 2 grados de libertad; en general, se agregan las siguientes 2 condiciones:
s ¢¢ ( x 0 ) = 0 s ¢¢( x n ) = 0 De lo cual vamos a obtener :
s ¢¢ (2 ) = 0 Þ 6 a 1 (2 ) + 2 b 1 = 0 \ 12 a1 + 2 b 1 = 0 s ¢¢ (5 ) = 0 Þ 6 a 2 (5 ) + 2 b 2 = 0 \ 30 a 2 + 2 b 2 = 0 Con lo cual, hemos completado un juego de 8 ecuaciones vs. 8 incógnitas, el cual es el siguiente:
ISIDORO PONTE­E.S.M.C.154 PRÁCTICAS DE LABORATORIO­MÉTODOS NUMÉRICOS 8 a 1 + 4 b 1 + 2 c 1 + d 1 = -1 27 a 1 + 9 b 1 + 3 c 1 + d 1 = 2 27 a 2 + 9 b 2 + 3 c 2 + d 2 = 2 125 a 2 + 25 b 2 + 5 c 2 + d 2 = -7 27 a 1 + 6 b 1 + c 1 = 27 a 2 + 6 b 2 + c 2 18 a 1 + 2 b 1 = 18 a 2 + 2 b 2 12 a 1 + 2 b 1 = 0 30 a 2 + 2 b 2 = 0 Cuya forma matricial es la siguiente :
é 8 ê 27 ê
ê 0 ê
ê 0 ê 27 ê
ê 18 ê12 ê
ëê 0 0 ù é a 1 ù é - 1 ù
3 1 0 0 0 0 úú êê b 1 úú êê 2 úú
0 0 27 9 3 1 ú ê c 1 ú ê 2 ú
úê ú ê
ú
0 0 125 25 5 1 ú ê d 1 ú ê - 7 ú
=
1 0 - 27 - 6 - 1 0 ú ê a 2 ú ê 0 ú
úê ú ê
ú
0 0 - 18 - 2 0 0 ú ê b 2 ú ê 0 ú
0 0 0 0 0 0 ú ê c 2 ú ê 0 ú
úê ú ê
ú
0 0 30 2 0 0 ûú ëê d 2 ûú êë 0 úû
4 2 1 9 0 0 6 2 2 0 0 0 0 Usando MATHEMATICA, obtenemos la siguiente solución: a 1
b 1 c 1 d 1 a 2 b 2 c 2 d 2 = - 1 . 25 =
7 . 5 = - 10 . 75 =
0 . 5 =
0 . 625 = - 9 . 375 = 39 . 875 = - 50 . 125 Sustituyendo estos valores en nuestra función inicial, vemos que la spline cúbica para la tabla de datos dada, queda definida como sigue:
si x Î [2 , 3 ]
î0 . 625 x - 9 . 375 x + 39 . 875 x - 50 . 125 si x Î [3 , 5 ]
ì
s( x ) = í
- 1 . 25 x 3 + 7 . 5 x 2 - 10 . 75 x + 0 . 5 3 2 ISIDORO PONTE­E.S.M.C.155 PRÁCTICAS DE LABORATORIO­MÉTODOS NUMÉRICOS Mostramos la gráfica correspondiente a este ejercicio, creada tambien en MATHEMATICA. Obsérvese la finura con la que se unen los polinomios cúbicos que conforman a la spline. Prácticamente ni se nota que se trata de dos polinomios diferentes. Esto es debido a las condiciones que se impusieron sobre las derivadas de la función. Esta finura casi artística, es la que permite aplicar las splines cúbicas, para cuestiones como el diseño de letras por computadoras, o bien a problemas de aplicación donde la interpolación que se necesita es de un carácter bastante delicado, como podría tratarse de datos médicos sobre algún tipo de enfermedad. Veamos un ejemplo con MATHEMATICA: Este programa nos calcula y grafica la spline cúbica correspondiente a una tabla de datos con cuatro puntos. Los datos son los siguientes: In[1]:= datos = 88-1, - 1<, 81, 1<, 82, 5<, 84, - 2<<
Out[1]= 88-1, - 1<, 81, 1<, 82, 5<, 84, - 2<<
Estos se guardan como un objeto gráfico como sigue: In[2]:= mi =
Min@Table@datos@@i, 1DD, 8i, 4<DD
Out[2]= -1
ISIDORO PONTE­E.S.M.C.156 PRÁCTICAS DE LABORATORIO­MÉTODOS NUMÉRICOS In[3]:= mf =
Max@Table@datos@@i, 1DD, 8i, 4<DD
Out[3]= 4 In[4]:= Mi =
Min@Table@datos@@i, 2DD, 8i, 4<DD
Out[4]= -2 In[5]:= Mf =
Max@Table@datos@@i, 2DD, 8i, 4<DD
Out[5]= 5 In[6]:= puntos =
ListPlot@datos, PlotStyle ® [email protected], PlotRange ® 88 mi - 5, mf + 5<, 8 Mi - 5, Mf + 5<<, AxesOrigin ® 80, 0<, DisplayFunction ® IdentityD
Out[6]= … Graphics …
Para calcular la función spline cúbica, primero definimos los polinomios cúbicos: 3
2
In[7]:= DoAsi@x_D = ai * x + bi * x + ci * x + di, 8i, 3<E
Enseguida, definimos las ecuaciones que se forman por las condiciones de la spline cúbica: In[8]:= eq1 = s1@datos@@1, 1DDD ==
datos@@1, 2DD
Out[8]= -a1 + b1 - c1 + d1 == - 1 In[9]:= eq2 = s1@datos@@2, 1DDD ==
datos@@2, 2DD
Out[9]= a1 + b1 + c1 + d1 == 1 In[10]:= eq3 = s2@datos@@2, 1DDD ==
datos@@2, 2DD
Out[10]= a2 + b2 + c2 + d2 == 1 In[11]:= eq4 = s2@datos@@3, 1DDD ==
datos@@3, 2DD
Out[11]= 8 a2 + 4 b2 + 2 c2 + d2 == 5 In[12]:= eq5 = s3@datos@@3, 1DDD ==
datos@@3, 2DD
Out[12]= 8 a3 + 4 b3 + 2 c3 + d3 == 5
ISIDORO PONTE­E.S.M.C.157 PRÁCTICAS DE LABORATORIO­MÉTODOS NUMÉRICOS In[13]:= eq6 = s3@datos@@4, 1DDD ==
datos@@4, 2DD
Out[13]= 64 a3 + 16 b3 + 4 c3 + d3 == - 2 In[14]:= eq7 = s1 '@datos@@2, 1DDD == s2 '@datos@@2, 1DDD
Out[14]= 3 a1 + 2 b1 + c1 == 3 a2 + 2 b2 + c2 In[15]:= eq8 = s2 '@datos@@3, 1DDD == s3 '@datos@@3, 1DDD
Out[15]= 12 a2 + 4 b2 + c2 == 12 a3 + 4 b3 + c3 In[16]:= eq9 = s1 ''@datos@@2, 1DDD == s2 ''@datos@@2, 1DDD
Out[16]= 6 a1 + 2 b1 == 6 a2 + 2 b2 In[17]:= eq 10 = s2 ''@datos@@3, 1DDD == s 3 ''@datos@@3, 1DDD
Out[17]= 12 a2 + 2 b2 == 12 a3 + 2 b3 In[18]:= eq 11 = s1 ''@datos@@1, 1DDD Š 0 Out[18]= -6 a1 + 2 b1 ==
0 In[19]:= eq 12 = s3 ''@datos@@4, 1DDD Š 0 Out[19]= 24 a3 + 2 b3 == 0 Resolvemos el sistema de ecuaciones: In[20]:= Solve@Table@eq i, 8i, 12<DD
51 153 21 297 , b1 ®
, a2 ® , b2 ®
, 140 140 10 35 24 288 89 473 a3 ®
, b3 ® , c1 ®
, c2 ® , 35 35 140 70 1867 153 48 732 c3 ®
, d1 ® , d2 ®
, d3 ® >>
70 140 35 35
Out[20]= ::a1 ®
A continuación vaciamos la información obtenida en el paso anterior e introducimos la regla de correspondencia de la spline:
ISIDORO PONTE­E.S.M.C.158 PRÁCTICAS DE LABORATORIO­MÉTODOS NUMÉRICOS In[21]:= splinecubica =
WhichAdatos@@1, 1DD £ x £ datos@@2, 1DD,
2 51 x 3 153 89 x 153 x +
+
+
, 140 140
140
140 datos@@2, 1DD £ x £ datos@@3, 1DD, 2 21 x 3 48 473 x 297 x +
, 35
70
35
10 datos@@3, 1DD £ x £ datos@@4, 1DD,
2 24 x 3 732 1867 x 288 x +
+
E
35
70
35
35
La cual se guarda como un objeto gráfico como sigue: In[22]:= graficaspline:=
Plot@splinecubica, 8x, mi, mf<, DisplayFunction ® IdentityD
Finalmente, le pedimos a MATHEMATICA que nos escriba la regla de correspondencia de la spline cúbica, así como la grafica correspondiente: In[23]:= PrintA"La función spline cúbica que se obtuvo es: ",
2 51 x 3 153 89 x 153 x +
+
+
, " si ", 140 140
140
140 datos@@1, 1DD £ x £ datos@@2, 1DDE; 2 21 x 3 48 473 x 297 x +
, PrintA" ", 35
70
35
10 " si ", datos@@2, 1DD £ x £ datos@@3, 1DDE; 2 24 x 3 732 1867 x 288 x +
+
, " si ", PrintA35
70
35
35 datos@@3, 1DD £ x £ datos@@4, 1DDE; Print@" La gráfica correspondiente es:"D
Show@ puntos, graficaspline, Ticks -> 88-1, 1, 2, 4<, 8-1, 1, 5, -2<<, DisplayFunction -> $DisplayFunctionD; La función spline cúbica que se obtuvo es:
2 51 x 3 153 89 x 153 x +
+
+
si
-1 £ x £ 1
140
140
140
140 ISIDORO PONTE­E.S.M.C.159 PRÁCTICAS DE LABORATORIO­MÉTODOS NUMÉRICOS 2 21 x 3 48 473 x 297 x +
35
70
35
10 -
2 24 x 3 732 1867 x 288 x +
+
35
70
35
35 si si 1£x£2
2 £ x £ 4 La gráfica correspondiente es: 5 1 ­1 ­1 ­2 1 2 4 Podríamos hacer el programa en dos pasos, si quisiesemos: In[1]:= datos = 88-1, -1<, 81, 1<, 82, 5<, 84, -2<<; mi = Min@Table@datos@@i, 1DD, 8i, 4<DD; mf = Max@Table@datos@@i, 1DD, 8i, 4<DD; Mi = Min@Table@datos@@i, 2DD, 8i, 4<DD; Mf = Max@Table@datos@@i, 2DD, 8i, 4<DD; puntos = ListPlot@datos, PlotStyle ® [email protected], PlotRange ® 88 mi - 5, mf + 5<, 8 Mi - 5, Mf + 5<<, AxesOrigin ® 80, 0<, DisplayFunction ® IdentityD; 3
2
DoAs i@x_D = a i * x + b i * x + c i * x + d i, 8i, 3<E; eq 1 = s 1@datos@@1, 1DDD == datos@@1, 2DD; eq 2 = s 1@datos@@2, 1DDD == datos@@2, 2DD; eq 3 = s 2@datos@@2, 1DDD == datos@@2, 2DD; eq4 = s2@datos@@3, 1DDD == datos@@3, 2DD; eq 5 = s 3@datos@@3, 1DDD == datos@@3, 2DD; eq6 = s3@datos@@4, 1DDD == datos@@4, 2DD; eq 7 = s 1 '@datos@@2, 1DDD == s 2 '@datos@@2, 1DDD; eq8 = s2 '@datos@@3, 1DDD == s3 '@datos@@3, 1DDD; eq 9 = s 1 ''@datos@@2, 1DDD == s 2 ''@datos@@2, 1DDD; eq10 = s2 ''@datos@@3, 1DDD == s3 ''@datos@@3, 1DDD; eq 11 = s 1 ''@datos@@1, 1DDD Š 0; eq12 = s3 ''@datos@@4, 1DDD Š 0; Solve@Table@eq i, 8i, 12<DD
ISIDORO PONTE­E.S.M.C.160 PRÁCTICAS DE LABORATORIO­MÉTODOS NUMÉRICOS 51 153 21 297 , b1 ®
, a2 ® , b2 ®
, 140 140 10 35 24 288 89 473 a3 ®
, b3 ® , c1 ®
, c2 ® , 35 35 140 70 1867 153 48 732 c3 ®
, d1 ® , d2 ®
, d3 ® >>
70 140 35 35
Out[10]= ::a1 ®
incorporar los datos obtenidos, para finalizar el programa In[11]:= splinecubica =
WhichAdatos@@1, 1DD £ x £ datos@@2, 1DD,
2 51 x 3 153 89 x 153 x +
+
+
, 140 140
140
140 datos@@2, 1DD £ x £ datos@@3, 1DD, 2 21 x 3 48 473 x 297 x +
, 35
70
35
10 datos@@3, 1DD £ x £ datos@@4, 1DD,
2 24 x 3 732 1867 x 288 x +
+
E; 35
70
35
35
graficaspline := Plot@splinecubica, 8x, mi, mf<, DisplayFunction ® IdentityD; PrintA"La función spline cúbica que se obtuvo es: ",
2 51 x 3 153 89 x 153 x +
+
+
, " si ", 140 140
140
140 datos@@1, 1DD £ x £ datos@@2, 1DDE; 2 21 x 3 48 473 x 297 x PrintA" ", +
, 35
70
35
10 " si ", datos@@2, 1DD £ x £ datos@@3, 1DDE; 2 24 x 3 732 1867 x 288 x PrintA+
+
, " si ", 35
70
35
35 datos@@3, 1DD £ x £ datos@@4, 1DDE; Print@" La gráfica correspondiente es:"D
Show@ puntos, graficaspline, Ticks -> 88-1, 1, 2, 4<, 8-1, 1, 5, -2<<, DisplayFunction -> $DisplayFunctionD;
ISIDORO PONTE­E.S.M.C.161 PRÁCTICAS DE LABORATORIO­MÉTODOS NUMÉRICOS La función spline cúbica que se obtuvo es:
2 51 x 3 153 89 x 153 x +
+
+
si
- 1 £ x £ 1 140
140
140
140 2 21 x 3 48 473 x 297 x +
35
70
35
10 -
2 24 x 3 732 1867 x 288 x +
+
35
70
35
35 si si 1£x£2
2 £ x £ 4 La gráfica correspondiente es: 5 1 ­1 ­1 ­2
1 2 4 ISIDORO PONTE­E.S.M.C.162 PRÁCTICAS DE LABORATORIO­MÉTODOS NUMÉRICOS Ejercicios complementarios personalizados 1A Calcula el polinomio de interpolación de Newton para los siguientes datos (2, 0 . 5 d 7 ), (- 2 , -3 ), (1 , 2 . 4 d 3 ), (4 , 7 . 8 ) . 1B Calcula el polinomio de interpolación de Newton para los siguientes datos (0. 3 d 6 , -3 ), (0 . 6 , 0 ), (0 . 9 d 4 , -6 ), (1 . 2 , 9 ), (1 . 5 , -12 ) 2A Calcula el polinomio de interpolación de Lagrange para los siguientes datos (1, 1 . 56 ), (- 2 , 3 . 54 ), (3 , -2 . 57 ), (- 5 , -8 . 9 d 2 ) 2B Calcula el polinomio de interpolación de Lagrange para los siguientes datos (1. 5 d 2 , 9 ), (- 0 . 5 d 4 , -2 ), (1 , 5 ), (- 2 , 33 ), (- 4 , 0 ) 2C Calcula las splines cúbicas para los siguientes datos
(- 2, 40 ), (1 , -5 ), (3 , -20 ) ISIDORO PONTE­E.S.M.C.163 PRÁCTICAS DE LABORATORIO­MÉTODOS NUMÉRICOS PRACTICAS DE LABORATORIO MÉTODOS NUMÉRICOS PRÁCTICA 25­26: Derivación e integración Integración aproximada. numerica. En los cursos de Cálculo Integral, nos enseñan como calcular una integral definida de una función continua mediante una aplicación del Teorema Fundamental del Cálculo: Teorema Fundamental del Cálculo : Sea una función continua en el intervalo primitiva de y sea una . Entonces: El problema en la práctica, se presenta cuando nos vemos imposibilitados de encontrar la primitiva requerida, aún para integrales aparentemente sencillas como: la cual simplemente es imposible de resolver con el Teorema Fundamental del Cálculo. Incluso casos en los que no se conoce la primitiva, o en que casos en los que solo se conoce una serie de valores , En estas prácticas estudiaremos diversos métodos numéricos que nos permitirán obtener aproximaciones bastante exactas a integrales como la mencionada anteriormente. Esencialmente, veremos dos tipos de integración numérica: las fórmulas de Newton­Cotes y el algoritmo de Romberg. Las fórmulas de Newton­Cotes están conformadas por las bien conocidas reglas del trapecio y de Simpson (regla de un tercio y de tres octavos). El algoritmo de Romberg forma parte de un método conocido como método de extrapolación de Richardson.
ISIDORO PONTE­E.S.M.C.164 PRÁCTICAS DE LABORATORIO­MÉTODOS NUMÉRICOS Haciendo uso de MATHEMATICA es posible discernir sobre las cualidades y defectos de cada uno de los métodos mencionados arriba. REGLA DEL TRAPECIO Dada una función , dividimos el intervalo misma longitud una función continua en el intervalo en subintervalos, todos de la . Sea la partición que se forma al hacer dicha subdivisión. ya que todos los subintervalos tienen la misma longitud h, tenemos la fórmula de los trapecios que es: Sustituyendo el valor de h y usando la notación sigma, tenemos finalmente: Esta es la regla del trapecio para n subintervalos. Obviamente, cuantos más subintervalos usemos, mejor será la aproximación a la integral. Este programa calcula usando la Regla del Trapecio con n intervalos todos de la misma longitud con el programa MATHEMATICA Introducimos el intervalo de integración: In[1]:= a = 0; b = 6; Introducimos la función a integrar:
ISIDORO PONTE­E.S.M.C.165 PRÁCTICAS DE LABORATORIO­MÉTODOS NUMÉRICOS In[3]:= f@x_D =
Cos@xD
; x + 1 Introducimos el número de intervalos a usar: In[4]:= n = 100000; Calculamos la longitud de cada intervalo: In[5]:= h = Hb- aL • n; Calculamos los puntos de la partición que se genera: In[6]:= Do@xi = a + h* i, 8i, 0, n<D; Aplicamos la regla del trapecio para n intervalos: In[7]:= Trapecio = Hb- aL•H2 * nL *
i
in-1 y
y
k
ki=1 {
{
f@x0D + 2* âf@xiD + f@xnD ; Calculamos el valor exacto de la integral, para comparar: b In[8]:= Verdadero = à f@xD âx; a Le pedimos a Mathematica que nos despliegue la información en pantalla: In[9]:= Print@"El valor aproximado de la integral usando la regla del trapecio es:"D
Print@N@Trapecio, 24DD
Print@"El valor verdadero de la integral es:"D
Print@N@Verdadero, 24DD
El valor aproximado de la integral usando la regla del trapecio es: 0.287037911124493970829599 El valor verdadero de la integral es: 0.287037910818397614552909 0.287038
ISIDORO PONTE­E.S.M.C.166 PRÁCTICAS DE LABORATORIO­MÉTODOS NUMÉRICOS REGLA DE SIMPSON DE UN TERCIO Dada una función , dividimos el intervalo misma longitud una función continua en el intervalo en subintervalos, todos de la , sea x mi el punto medio del subintervalo
[x i -1 , x i ] . Sea la partición que se forma al hacer dicha subdivisión. ya que todos los subintervalos tienen la misma longitud h, tenemos la fórmula de los trapecios que es: Sustituyendo el valor de h y usando la notación sigma, tenemos finalmente: Esta es la regla de Simpson de un tercio para n subintervalos. Obviamente, cuantos más subintervalos usemos, mejor será la aproximación a la integral. Calcula la integral usando la Regla de Simpson de 1/3 para n intervalos, todos de la misma longitud con el programa MATHEMATICA Introducimos el intervalo de integración: In[1]:= a = -2; b = 1; Introducimos la función a integrar: In[3]:= f@x_D = ArcTan@x^3 + 2D; Introducimos el número de intervalos a usar:
ISIDORO PONTE­E.S.M.C.167 PRÁCTICAS DE LABORATORIO­MÉTODOS NUMÉRICOS In[4]:= n = 100000; Calculamos la longitud de los intervalos: In[5]:= h = Hb- aL • n; Calculamos los elementos de la partición que se genera: In[6]:= Do@xi = a + h* i, 8i, 0, n<D; Calculamos los puntos medios de cada intervalo: In[7]:= Do@yi = Hxi + xi-1L • 2, 8i, 1, n<D; Aplicamos la Regla de Simpson de 1/3: In[8]:= Simpson13 = Hb- aL•H6 * nL *
i
i n y
in-1 y
y
k
ki=1 {
ki=1 {
{
f@x 0D + 4* âf@y iD + 2* â f@xiD + f@xnD ; Calculamos el valor exacto de la integral: b In[9]:= Verdadero = à f@xD âx; a Le pedimos a MATHEMATICA que nos despliegue en pantalla la información obtenida: In[10]:= Print@
"El valor aproximado de la integral usando la regla del simpson es:"D
Print@ N@Simpson13, 12DD
Print@"El valor verdadero de la integral es:"D
Print@ N@ Verdadero, 12DD
El valor aproximado de la integral usando la regla del simpson es: 1.56093 El valor verdadero de la integral es: 1.56093 + 3.33067 ´ 10-16 ä
ISIDORO PONTE­E.S.M.C.168 PRÁCTICAS DE LABORATORIO­MÉTODOS NUMÉRICOS REGLA DE SIMPSON DE TRES OCTAVOS Dada una función , dividimos el intervalo misma longitud una función continua en el intervalo en subintervalos, todos de la , sea el subintervalo [x i -1 , x i ] , lo dividimos en tres partes iguales y lo puntos intermedios los llamamos y i y z i .l Sea la partición que se forma al hacer dicha subdivisión. ya que todos los subintervalos tienen la misma longitud h, tenemos la fórmula de los trapecios que es: Sustituyendo el valor de h, tenemos usando el polinomio de interpolación de Lagrange y el método de integración por partes: Esta es la regla de Simpson de tres octavos para n subintervalos iguales. Obviamente, cuantos más subintervalos usemos, mejor será la aproximación a la integral. Introducimos el intervalo de integración: In[1]:= a =
0; b = 2; Introducimos la función a integrar: In[3]:= f@x_D =
x * Exp@xD; Introducimos el número de intervalos a usar: In[4]:= n =
2; Calculamos la longitud de cada intervalo: In[5]:= h = Hb - aL • n;
ISIDORO PONTE­E.S.M.C.169 PRÁCTICAS DE LABORATORIO­MÉTODOS NUMÉRICOS Calculamos los elementos de la partición que se genera: In[6]:= Do@xi = a + h * i, 8i, 0, n<D; Calculamos los puntos que dividen en tres partes iguales a cada intervalo: In[7]:= Do@y i = H2 * xi-1 + x iL • 3, 8i,
1, n<D; In[8]:= Do@z i = Hx i-1 + 2 * x iL • 3, 8i,
1, n<D; Aplicamos la Regla de Simpson de 3/8 para n intervalos: In[9]:= Simpson38 =
H b - aL • H8 * nL *
in
y
i n-1 y
y
f@x0D + 3 * â Hf@yiD + f@ziDL + 2 * â f@xiD + f@xnD ; k
ki=1 {
ki=1 {
{
i
Calculamos el valor exacto de la integral: b In[10]:= Verdadero = à f@xD â x; a Le pedimos a MATHEMATICA que nos despliegue en pantalla la información obtenida: In[11]:= Print@"El valor aproximado de la integral usando la regla del simpson es:"D
Print@ N@Simpson38, 14DD
Print@"El valor verdadero de la integral es:"D
Print@ N@ Verdadero, 14DD
El valor aproximado de la integral usando la regla del simpson es: 8.39411 El valor verdadero de la integral es: 8.38906 De lo visto resulta evidente que la regla de Simpson de tres octavos, es más exacta, que la de un tercio, que a su vez, aproxima
ISIDORO PONTE­E.S.M.C.170 PRÁCTICAS DE LABORATORIO­MÉTODOS NUMÉRICOS mejor que la del trapecio. Se pueden establecer cotas para los errores que se comenten en cada uno de los métodos. los siguientes resultados se mencionan para completar la información, pero omitimos las demostraciones correspondientes. REGLA F O R M U L A E R R O R CON. . . Trapecio Simpson Simpson
ISIDORO PONTE­E.S.M.C.171 PRÁCTICAS DE LABORATORIO­MÉTODOS NUMÉRICOS Ejercicios complementarios personalizados 1A Usar la regla del trapecio para aproximar, i) Dividiendo en un solo intervalo. ii) Dividiendo en 6+d4 intervalos. 1B Usar la regla de Simpson 1/3 para aproximar, i) Dividiendo en un solo intervalo. ii) Dividiendo en 4+d2 intervalos. 2A. Usar la regla de Simpson 3/8 para aproximar, i) Dividiendo en un solo intervalo. ii) Dividiendo en 4+d5 intervalos. 2B. Integrar la siguiente tabla de datos: i) 2C. Integrar la siguiente tabla de datos: ii)
ISIDORO PONTE­E.S.M.C.172 PRÁCTICAS DE LABORATORIO­MÉTODOS NUMÉRICOS PRACTICAS DE LABORATORIO MÉTODOS NUMÉRICOS PRÁCTICA 27­28: Métodos numéricos de resolución de ecuaciones diferenciales. ECUACI ON ES DI FERENCI ALES En esta práctica, haremos un breve estudio de los métodos numéricos básicos que se usan para aproximar soluciones de algunas ecuaciones diferenciales. Recordamos rápidamente, que una ecuación diferencial (ordinaria) es aquella que involucra una variable independiente, una variable dependiente y la derivada (ó derivadas ) de esta última. En una ecuación diferencial, la incógnita es la variable dependiente y se espera encontrarla como función de la variable independiente, de tal forma que si se sustituye dicha variable dependiente, así como las derivadas que aparecen en la ecuación diferencial, la igualdad que resulta es verdadera. De cursos anteriores de ecuaciones diferenciales, sabemos que en general, existen una infinidad de funciones (curvas) que resuelven una misma ecuación diferencial. Por ejemplo, la ecuación: tiene como solución general: donde c es una constante arbitraria que puede ser cualquier número real (y de aquí la infinidad de curvas solución que mencionamos arriba). En este curso, estudiaremos solamente ecuaciones diferenciales de primer orden del tipo:
ISIDORO PONTE­E.S.M.C.173 PRÁCTICAS DE LABORATORIO­MÉTODOS NUMÉRICOS donde es una función de dos variables. Cuando se desea que la curva solución pase por algún punto específico, digamos , entonces se dice que se trata de una ecuación diferencial con una condición inicial dada. Así, estudiaremos ecuaciones diferenciales de la forma con la condición inicial . Obviamente, la importancia de los métodos numéricos radica en la aparición de ecuaciones diferenciales que no pueden resolverse por métodos tradicionales, y de ahí la necesidad de implementar algún método de aproximación. Veremos tres métodos numéricos:
·
·
·
El método de Euler.
El método de Euler mejorado.
El método de Runge­Kutta de orden 4. En todos estos métodos se busca aproximar el valor donde es un valor cercano a (el de la condición inicial dada). Comencemos con el primer método que como siempre, no es el más exacto, pero si el más sencillo y simple de explicar, así como el que marca la pauta para desarrollar los otros métodos. MÉTODO DE EULER La idea del método de Euler es muy sencilla y está basada en el significado geométrico de la derivada de una función en un punto dado. Supongamos que tuviéramos la curva solución de la ecuación diferencial y trazamos la recta tangente a la curva en el punto dado por la condición inicial.
ISIDORO PONTE­E.S.M.C.174 PRÁCTICAS DE LABORATORIO­MÉTODOS NUMÉRICOS Debido a que la recta tangente aproxima a la curva en valores cercanos al punto de tangencia, podemos tomar el valor de la recta tangente en el punto como una aproximación al valor deseado . Así, calculemos la ecuación de la recta tangente a la curva solución de la ecuación diferencial dada en el punto . De los cursos de Geometría Analítica, sabemos que la ecuación de la recta es: donde m es la pendiente. En este caso, sabemos que la pendiente de la recta tangente se calcula con la derivada: Por lo tanto, la ecuación de la recta tangente es : Ahora bien, suponemos que tanto estará dado como siguiente aproximación:
es un punto cercano a , y por lo . De esta forma, tenemos la ISIDORO PONTE­E.S.M.C.175 PRÁCTICAS DE LABORATORIO­MÉTODOS NUMÉRICOS De aquí, tenemos nuestra fórmula de aproximación: Esta aproximación puede ser suficientemente buena, si el valor de h es realmente pequeño, digamos de una décima ó menos. Pero si el valor de h es más grande, entonces podemos cometer mucho error al aplicar dicha fórmula. Una forma de reducir el error y obtener de hecho un método iterativo, es dividir la distancia en n partes iguales (procurando que estas partes sean de longitud suficientemente pequeña) y obtener entonces la aproximación en n pasos, aplicando la fórmula anterior n veces de un paso a otro, con la nueva h igual a . En una gráfica, tenemos lo siguiente: Ahora bien, sabemos que: Para obtener únicamente hay que pensar que ahora el papel de lo toma el punto , y por lo tanto, si sustituimos los datos adecuadamente, obtendremos que: De aquí se ve claramente que la fórmula recursiva general, está dada por: Esta es la conocida fórmula de Euler que se usa para aproximar el valor de aplicándola sucesivamente desde hasta en pasos de longitud h.
ISIDORO PONTE­E.S.M.C.176 PRÁCTICAS DE LABORATORIO­MÉTODOS NUMÉRICOS Este programa usa el método de Euler para aproximar dada la ecuación diferencial de la forma: y'=f(x,y) con y donde se supone que es cercano a . Introducimos la función f(x,y): In[1]:= f@x_, y_D =
Log@x + yD; Introducimos los valores dados por la condición inicial: In[2]:= x0 = 1; y 0 = 1.5; Introducimos el valor de h: In[4]:= h =
0.1; Introducimos el número de iteraciones: In[5]:= k =
5; Aplicamos la fórmula de Euler: In[6]:= Do@xn+1 = xn + h, 8n,
0, k<D; Do@yn+1 = yn + h* f@xn, ynD, 8n, 0, k<D; Le pedimos a MATHEMATICA que nos despliegue una tabla con los datos obtenidos: In[8]:= Print@
"La tabla de valores obtenidos mendiante el método de Euler es:"D
Table@8x i, yi<, 8i, 0, k<D •• TableForm La tabla de valores obtenidos mendiante el método de Euler es: Out[9]//TableForm= 1 1.1 1.2 1.3 1.4 1.5 1.5 1.59163 1.69064 1.79679 1.90983 2.02952
ISIDORO PONTE­E.S.M.C.177 PRÁCTICAS DE LABORATORIO­MÉTODOS NUMÉRICOS MÉTODO DE EULER MEJORADO Este método se basa en la misma idea del método anterior, pero hace un refinamiento en la aproximación, tomando un promedio entre ciertas pendientes. La fórmula es la siguiente: donde Para entender esta fórmula, analicemos el primer paso de la aproximación, con base en la siguiente gráfica: En la gráfica, vemos que la pendiente promedio corresponde a la pendiente de la recta bisectriz de la recta tangente a la curva en el punto de la condición inicial y la “recta tangente” a la curva en el punto , donde es la aproximación obtenida con la primera fórmula de Euler. Finalmente, esta recta bisectriz se traslada paralelamente hasta el punto de la condición inicial, y se considera el valor de esta recta en el punto como la aproximación de Euler mejorada.
ISIDORO PONTE­E.S.M.C.178 PRÁCTICAS DE LABORATORIO­MÉTODOS NUMÉRICOS Este programa usa el método de Euler mejorado igual que en el caso anterior Introducimos la función f(x,y): In[1]:= f@x_, y_D =
Log@x + yD; Introducimos los valores dados por la condición inicial: In[2]:= x0 = 1; y 0 = 1.5; Introducimos el valor de h: In[4]:= h =
0.1; Introducimos el número de iteraciones: In[5]:= k =
5; Aplicamos la fórmula de Euler: In[6]:= Do@xn+1 = xn + h, 8n,
0, k<D; DoA9zn+1 = yn + h * f@xn, ynD, f@xn, ynD + f@xn+1, zn+1D
yn+1 = yn + h * J
N=, 8n, 0, k<E; 2
pedimos a MATHEMATICA una tabla con los datos obtenidos: In[8]:= Print@
"La tabla de valores obtenidos mendiante el método de Euler Mejorado es:"D
Table@8x i, y i<, 8i, 0, k<D •• TableForm La tabla de valores obtenidos mendiante el método de Euler Mejorado es: Out[9]//TableForm= 1 1.1 1.2 1.3 1.4 1.5 1.5 1.59532 1.69804 1.80788 1.9246 2.04794
ISIDORO PONTE­E.S.M.C.179 PRÁCTICAS DE LABORATORIO­MÉTODOS NUMÉRICOS MÉTODO DE RUNGE – KUTTA Sin entrar en mucho detalle, mencionamos solamente que el método de Runge­Kutta cambia la dirección en el sentido de que no sigue la misma línea de los métodos de Euler. De hecho está basado en una aplicación de los polinomios de Taylor. Comentamos sin embargo, que el método de Runge­Kutta si contiene como casos especiales los de Euler. Las fórmulas donde Se conocen como las reglas o fórmulas de Runge­Kutta de orden cuatro para la ecuación diferencial: Este programa usa el método de Runge­Kutta para aproximar y(x1), dada la ecuación diferencial de la forma: y'=f(x,y) con la condición inicial: y(x0)=y0 y donde se supone que x1 es cercano a x0 .
ISIDORO PONTE­E.S.M.C.180 PRÁCTICAS DE LABORATORIO­MÉTODOS NUMÉRICOS Introducimos la función f(x,y): In[1]:= f@x_, y_D =
Log@xD +
1 ; y Introducimos los valores dados por la condición inicial: In[2]:= x0 = 4; y0 = 5; Introducimos el valor de h: In[4]:= h =
0.1; Introducimos el número de iteraciones: In[5]:= k =
5; Aplicamos la fórmula de Euler: In[6]:= Do@xn+1 = xn + h, 8n, 0, k<D; Do@8k1 = h * f@xn, ynD, k2 = h * f@xn + 0.5 * h, yn + 0.5 * k1D, k3 = h * f@xn + 0.5 * h, yn + 0.5 * k2D, k4 = h * f@xn + h, yn + k3D, yn+1 = yn + H1 • 6L * Hk1 + 2 * k2 + 2 * k3 + k4L<, 8n, 0, k - 1<D; Le pedimos a MATHEMATICA que nos despliegue una tabla con los datos obtenidos: In[8]:= Print@
"La tabla de valores obtenidos mendiante el método de Runge-Kutta es:"D
Table@8x i, y i<, 8i, 0, k<D •• TableForm La tabla de valores obtenidos mendiante el método de Runge-Kutta es: Out[9]//TableForm= 4 4.1 4.2 4.3 4.4 4.5 5 5.15956 5.32095 5.48415 5.64913 5.81587
ISIDORO PONTE­E.S.M.C.181 PRÁCTICAS DE LABORATORIO­MÉTODOS NUMÉRICOS Ejercicios complementarios personalizados: 1A. Dada la ecuación diferencial: Usa el método de Euler para aproximar y ( 2 . 3 d 6 ) tomando en cada paso del proceso iterativo. 1B. Dada la ecuación diferencial: Usa el método de Euler para aproximar y ( 1 . 3 d 1 ) tomando en cada paso del proceso iterativo. 2A. Dada la ecuación diferencial: Usa el método de Euler mejorado para aproximar y ( 2 . 3 d 7 ) tomando en cada paso del proceso iterativo. 2B. Dada la ecuación diferencial: Usa el método de Euler mejorado para aproximar y ( 3 . 3 d 6 ) tomando en cada paso del proceso iterativo. 2C. Dada la ecuación diferencial: Usa el método de Runge­Kutta para aproximar tomando en cada paso del proceso iterativo.
y ( 4 . 3 d 1 ) ISIDORO PONTE­E.S.M.C.182