UNIDAD UNID AD 5 Y 6 MÉTODOS NUMERICOS M.C. BRENDA EDITH MORALES FERNANDEZ IBARRA AYALA LUIS GUSTAVO HERNANDEZ ALVARADO MIGUEL ANGEL ING. ELECTRICA
INSTITUTO TECNOLOGICO DE VERACRUZ
Derivación numérica. Consideremos Consideremos una función f(x) de la cual se conoce un conjunto discreto de valores (x0, f0), (x 1, f 1),...,(xn, f n). El problema que vamos a abordar es el de calcular la derivada de la función en un punto x que en principio no tiene por qué coincidir con alguno de los que figuran en los dat os de que disponemos. disponemos. La forma más sencilla de resolver el problema de la diferenciación diferenciación numérica consiste en estimar la derivada utilizando fórmulas obtenidas mediante la aproximación de Taylor, que se denominan fórmulas de diferencias finitas. Es importan te tener en cuenta que el proceso de diferenciación numérica es inestable. Los errores que tengan los datos, por ejemplo los cometidos en la adquisición adquisición de los mismos o los debidos al redondeo aumentan en el proceso de diferenciación diferenciación como veremos a lo largo de éste capítulo. Fórmulas de diferencias de dos puntos
Este proceso de paso al límite presenta distintos problemas para ser realizado en situaciones prácticas donde no se conozca la forma explícita de f′(x). En primer lugar un límite no puede calcula rse de modo aproximado en un computador donde los números que se manejan son
finitos. A pesar de todo es de esperar que si la función f(x) no se comporta mal y h0 es un número finito pero pequeño se cumpla:
Es más, la misma definición de la derivada impli ca que si f′(x) existe, entonces hay algún h0 a partir del cual nuestra aproximación dista menos
de una cantidad δ del valor real para la derivada. El problema es que esto 101
sólo es cierto con precisión infinita ya que h0 puede ser tan t an pequeño que no pueda representarse presentarse en el ordenador o que la diferencia d iferencia f(x + h0) − f(x) esté seriamente afectada por el error de redondeo. La ecuación (2.1) es la forma más sencilla de aproximar una derivada conocidas f(x) f(x) y f(x+h0). El siguiente siguiente teorema nos proporciona proporciona información sobre la precisión de esta aproximación. Teorema. Sea f(x) ∈ C1 (a, b) y existe f′′(x) en (a, b), entonces se cumple que:
Demostración. Demostración. Escribamos la aproximación de Taylor para la función en un punto x+h:
Reordenando Reordenando la expresión anterior queda demostrado el teorema. El teorema anterior nos indica que el error cometido al aproximar la derivada primera por su fórmula de diferencia adelantada adelantada es una función lineal de h. Cuanto menor sea h (o sea al tomar valores de f(x) más cercanos) la derivada numérica será más precisa. Este error se denomina error de truncación o discretización discretización y puede acotarse fácilmente,
obteniéndose que: E ≤ máx(x,x+h) |f′′(z)|. En realidad, para datos obtenidos a partir de una tabla esta acotación no es de gran utilidad directa ya que si no se conoce la derivada primera menos aún se conocerá la segunda pero al menos nos permite conocer el orden de aproximación de la fórmula.
102
Geométricamente el error O(h) procede del hecho de aproximar la derivada por la pendiente de la cuerda que une los puntos f(x) y f(x + h), Por otro lado, si existe la derivada deben existir las derivadas laterales y entonces
Un problema que presenta esta fórmula es que la precisión de la misma es baja y por lo tanto en situaciones donde sólo dispongamos de un muestreo de baja precisión de f(x), como ocurre en ensayos, datos experimentales, etc., será conveniente utilizar otras fórmulas de derivación más precisas. Fórmulas de orden superior El error de truncación de la fórmula de diferencia adelantada de dos puntos varía linealmente con h, de manera que es necesario usar valores de h muy pequeños para reducir suficientemente los errores de truncación. Es posible deducir fórmulas para las derivadas con errores de truncación más pequeños. Por ejemplo, tomemos
donde z1 ∈ (x, x + h) y z2 ∈ (x − h, x). Restando (2.1a) y (2.1a) obtenemos
que nos proporciona una siguiente aproximación para la derivada con un término de error de truncación que depende cuadráticamente de h. 103
Usando el teorema del valor intermedio f′′′(z) = (f′′′(z 1 ) + f′′′(z 2 ))/2, y entonces, si f es suficientemente derivable.
Usando los desarrollos de Taylor de f(x+h) y f(x+2h) se encuentra la llamada fórmula de diferencia adelantada de tres puntos que es:
Reemplazando h por −h en ( 2.3) obtenemos una fórmula de diferencias retrasadas de tres puntos
De todas estas, la fórmula de diferencia centrada es la que tiene, en principio, menor error de truncación y la que requiere menos evaluaciones de la función, siendo por lo tanto más efic iente desde el punto de vista computacional. Utilizando el valor de la función en más puntos se construyen fórmulas más precisas para las derivadas. Alguna de ellas se muestra en la tabla siguiente junto con las que hemos deducido ya.
104
Derivadas de orden superior El mismo procedimiento que se ha seguido al deducir fórmulas para calcular numérica- mente las derivadas primeras puede usarse para construir derivadas de orden superior partiendo del desarrollo de Taylor y eliminando las derivadas primeras. Consideremos por ejemplo las expresiones:
Sumando las ecuaciones anteriores y despejando se encuentra que:
Procediendo de la misma forma es posible encontrar aproximaciones que usen diferentes puntos y aproximaciones para derivadas de orden superior. La tabla siguiente presenta algunas de las fórmulas más comunes para calcular derivadas de orden superior.
105
Integración numérica: Método del trapecio, Métodos de Simpson 1/3 y 3/8. 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 f(x) una función continua en el intervalo [a,b] y sea F(x) una anti derivada def(x). Entonces:
106
El problema en la práctica, se presenta cuando nos vemos imposibilitados de encontrar la antiderivada requerida, aún para integrales aparentemente sencillas como:
la cual simplemente es imposible de resolver con el Teorema Fundamental del Cálculo.
REGLA DEL TRAPECIO Corresponde al caso donde n=1 , es decir :
donde f 1(x) es un polinomio de interpolación (obviamente de grado 1) para los datos: x a b y f(a) f(b) sabemos que este polinomio de interpolación es:
Integrando este polinomio, tenemos que:
107
Que es la conocida Regla del Trapecio. Este nombre se debe a la interpretación geométrica que le podemos dar a la fórmula. El polinomio de interpolación para una tabla que contiene dos datos, es una línea recta. La integral, corresponde al área bajo la línea recta en el intervalo [a,b], que es precisamente el área del trapecio que se forma.
REGLA DE SIMPSON DE UN TERCIO Suponemos que tenemos los datos: a f(a)
Xm F(Xm)
B F(b)
donde xm es el punto medio entre a y b. En este caso se tiene que:
Donde f 2(x) es el polinomio de interpolación para los datos en la tabla anterior. Usaremos el polinomio de LaGrange. Así, tenemos que:
Si denotamos h= (b-a)/2 = x m-a = b-x m, entonces:
108
Simplificando términos:
Vemos que cada uno de los términos anteriores, es esencialmente de la misma forma, es decir, una constante por (x- α)(x-β). Así, calculamos la siguiente integral por partes:
Obteniendo, por lo tanto,
Usamos esta fórmula para calcular la integral de cada uno de los tres términos de f 2(x) y obteniendo como resultado final
Debido al factor h/3 se le conoce como la regla de Simpson de un tercio . En la práctica, sustituimos el valor de h = (b-a)/ 2 para obtener nuestra fórmula final
REGLA DE SIMPSON DE TRES OCTAVOS 109
Este caso corresponde a n=3 , es decir,
donde f 3(x) es un polinomio de interpolación para los siguientes datos: X0 X1 X2 X3 f(X0) f(X1) fX2) F(X3) Y donde a= x0, b= x3 y x1, x2 son los puntos que dividen en tres partes iguales al intervalo [a,b]. Igual que en el caso anterior, se usa el polinomio de interpolación de Lagrange, y usando el método de integración por partes se llega a la siguiente fórmula:
donde h = (b-a) / 3 . Debido al factor 3h / 8 es que se le dio el nombre deRegla de Simpson de 3/8 . En la práctica, se sustituye el valor de h para obtener:
Integración con intervalos desiguales. Cuando la longitud de los subintervalos no es igual, se usa una combinación de la regla Trapezoidal y las reglas de Simpson, procurando seguir el siguiente orden jerárquico:
1 .- Simpson 3/8 Esta se aplica, si contamos con 4 puntos igualmente espaciados.
110
2 .- Simpson 1/3 Esta se aplica si falla (1) y contamos con 3 puntos igualmente espaciados.
3 .- Regla Trapezoidal Solo se aplica si no se cumple (1) y (2)
Ejemplo Evaluar
usando la siguiente tabla: x f(x)
0 0
0.1 6.84
0.3 4
0.5 4.2
0.7 5.51
0.95 5.77
1.2 1
Solución. Vemos que en el intervalo [0,0.01] podemos aplicar la regla del trapecio, en el intervalo [0.1,0.7] la regla de Simpson de 3/8 y en el intervalo [0.7,1.2] la regla de Simpson de 1/3. Así, tenemos las siguientes integrales:
Finalmente, la integral buscada es la suma de las tres integrales anteriores:
111
5.4 Aplicaciones. ·
Método del trapecio
Ejemplo 1: Utilizar la regla del trapecio para aproximar la integral:
Solución. Usamos la fórmula directamente con los siguientes datos:
Por lo tanto tenemos que:
·
Método de Simpson 1/3
Ejemplo 1. Usar la regla de Simpson de 1/3 para aproximar la siguiente integral:
Solución. Aplicamos la fórmula directamente, con los siguientes datos:
112
Por lo tanto, tenemos que:
·
Método de Simpson 3/8
Ejemplo 1. Aproximar la siguiente integral:
aplicando la regla de Simpson de 3/8, y subdiviendo en 3 intervalos.
Solución
Identificamos n=3 y la partición correspondiente:
Al considerar los puntos que dividen en tres partes iguales a cada subintervalo, tenemos los siguientes datos:
Sustituyendo todos los datos en la fórmula, obtenemos:
113
De acuerdo a los ejemplos vistos, resulta evidente que la regla de Simpson de 3/8, es más exacta que la de 1/3 y a su vez, ésta es más exacta que la regla del trapecio. En realidad, pueden establecerse cotas para los errores que se cometen en cada uno de estos métodos
UNIDAD VI. Ecuaciones Diferenciales Ordinarias 6.1 Fundamentos Matemáticos Las ecuaciones diferenciales tienen importancia fundamental en las ingenierías, debido a que muchas leyes y relaciones físicas se expresan matemáticamente mediante estas relaciones. Las siguientes ecuaciones son ejemplos de ecuaciones diferenciales: d 2 y dx
2
4y
0
dy dx
2 xy
x
2
z x
2
z y
3 xy 4 z sec y 0
Las dos primeras ecuaciones contienen derivadas ordinarias y por la forma en que están escritas vemos que y = f ( x ); la tercera contiene derivadas parciales y podemos ver que z = f ( x , y ). El orden de una ecuación diferencial es el máximo orden de las derivadas que contiene En esta unidad desarrollaremos métodos numéricos para encontrar la solución de ecuaciones diferenciales ordinarias a partir de valores iniciales. Un problema de valor inicial consiste en una ecuación diferencial, y en una condición que debe satisfacer la solución (o varias condiciones que se refieren al mismo valor de x , si la ecuación es de orden superior.)
114
dy dx
y 0 = y ( x 0)
f x y , ,
6.2 Métodos de un Solo Paso Su aplicación parte de y 0 = y ( x 0) y se avanza por pasos. En el primer paso se calcula un valor aproximado de y 1 de la solución y en x = x 0 + h, en el segundo paso se calcula un valor aproximado de y 2 en x = x 0 + 2h, y así sucesivamente. En cada paso, los cálculos e llevan a cabo mediante la misma fórmula, y en ellas h es un valor fijo.
6.2.1 Forma General para Métodos de un Solo Paso Deducción a partir de la serie de Taylor f xi 1 f xi f xi h
1 2!
f xi h 2
1 3!
f xi h3 ...
donde h = x i+1 - x i.
f ( x i+1) = f ( x i) + f ´ ( x i) h + (0) h 2 Si truncamos la serie de Taylor a partir del término
f ( x i+1) =
f ( x i)
+ f ´ ( x i)
1 2!
h
f xi h 2
+ (o) h 2
Valor Actual = Valor Anterior + Pendiente Tamaño del Paso + Error
Si hacemos = f ´ ( x i)
y i+1 = y i+ h
f ( x ) y i+1 = y i + h
115
y i+1
pendiente =
tamaño de
paso h= x i+1 - x i
Figura 6.1 Método de un solo paso
6.2.1.1 Método de Euler La primera derivada proporciona una aproximación directa a la pendiente en xi = f ( x i, y i)
donde f ( x i, y i) es la ecuación diferencial evaluada en ( x i, y i)
y i+1 = y i+ f ( x i, y i) h A esta fórmula se le conoce como método de Euler, o método de EulerCauchy o de pendiente puntual.
Ejemplo 6.1 Hallar el valor de f ( x ) en x =2, sí y (0) = 1 dy dx
yx 2
y
Analíticamente dy dx
yx 2
dy dx
y
y x 2
y x 2
1
1
dy y
x
d ln y
2
x
1 dx 2
1
116
ln y
x 3 3
x
x C
x 3
3
x
y
e
3
x
y
Ae 3
Entonces y (2) = 1.947734
Dado que y (0) = 1, A = 1
Numéricamente Por el método de Euler, usando h = 0.5, y (0) = 1. Ecuación del método
y i+1 = y i+ f ( x i , y i) h
x i = 0 y i = 1 f ( x i, y i) = y i ( x i 2 - 1) f (0, 1) = 1 (02 - 1) = -1 y i+1 = 1+ (-1) (0.5) = 0.5 x i+1 = x i+ h = 0 + 0.5 = 0.5 x i = 0.5 y i = 0.5 f ( x i, y i) = f (0.5, 0.5) = 0.5 [(0.5)2 – 1]
x i = 1 y i = 0.3125 f ( x i, y i) = f (1, 0.3125) = 0.3125 [(1) 2 – 1] = 0
y i+1 = 0.3125 + (0) (0.5) = 0.3125 x i+1 = x i+ h = 1 + 0.5 = 1.5 x i = 1.5 y i = 0.3125 f ( x i, y i) = f (1.5, 0.3125) = 0.3125 [(1.5) 2 – 1] = 0.390625
= -0.375
y i+1 = 0.5 + (-0.375) (0.5) = 0.3125 x i+1 = x i+ h = 0.5 + 0.5 = 1
y i+1 = 0.3125 + 0.390625 (0.5) = 0.5078125 x i+1 = x i+ h = 1.5 + 0.5 = 2
Aplicando el mismo procedimiento para h = 0.25 y h = 0.125 se obtiene 117
Tabla 6.1 Valores de y para distintos valores de h con Euler
h = 0.5 h = 0.25 h = 0.125 x y x y x y 0.0 1.0000 0.00 1.0000 0.000 1.000 0 0.125 0.875 0 0.25 0.7500 0.250 0.767 3 0.375 0.677 4 0.5 0.5000 0.50 0.5742 0.500 0.604 6 0.625 0.548 0 0.75 0.4666 0.750 0.506 2 0.875 0.478 5 1.0 0.3125 1.00 0.4155 1.000 0.464 5 1.125 0.464 5 1.25 0.4155 1.250 0.479 9 1.375 0.513 7 1.5 0.3125 1.50 0.4740 1.500 0.570 9 1.625 0.660 1 1.75 0.6221 1.750 0.795 4
118
1.875
2.0
0.507 8
2.00
1.000 5 0.942 2.000 1.315 8 1
6.2.1.2 Método de Heun (Euler-Gauss) Un método para mejorar la aproximación a la pendiente implica el cálculo de dos derivadas del intervalo, una en el punto inicial y la otra en el punto final. Enseguida se promedian las dos derivadas y se obtiene una aproximación mejorada de la pendiente en el intervalo completo. En el método el Euler, la pendiente al principio del intervalo es
y i´ = f ( x i, y i) se usa para extrapolar linealmente a y i+1 en x i+1 y i+1 = y i+ f ( x i, y i)
(Ecuación predictora)h Pero al final del intervalo se puede calcular una pendiente aproximada
y i+1´ = f ( x i+1, y i+1) Por la tanto se pueden combinar las dos pendientes y obtener una pendiente promedio en el intervalo: y
yi ´ yi 1´ 2
f xi , yi f xi 1 , yi 1 2
por lo que yi
1
yi
f xi , yi f xi 1 , yi
2
1
h
(Ecuación correctora)
Por ello, el método de Heun es un esquema predictor-corrector. Nótese que la ecuación correctora tiene el término y i+1 a ambos lados de la igualdad, y puede aplicarse para corregir en un esquema iterativo hasta que se obtenga una y i+1 mejorada para una tolerancia preestablecida.
119
Ejemplo 6.2 Resolver el ejemplo anterior utilizando el método de Heun, y valores de h de 0.5, 0.25 y 0.125. f ( x i, y i ) = y i ( x i2 - 1)
h = 0.5 x i = 0, y i = 1
1
f ( x i, y i) = 1 [02 - 1] = -1 y i+1 = y i+ f ( x i, y i) h = 1+ (-1) (0.5) = 0.5
2
0.5 0.65625
yi
1
1
1
0.4923 2
0.5
0.62695
Corrector segunda iteración yi
Corrector x i+1 = 0.5
1
1
yi
1
1
= 0.5 [(0.5)2 – 1] = -0.375 1
y i
- 1 - 47021 2
0.5 0.63245
Corrector tercera iteración
y i+1´ = f ( x i+1, y i+1) = f (0.5, 0.5)
0.375
Corrector primera iteración
Predictor
y i
1
f xi , y i f xi 1 , y i
2
1
- 1 - 0.4743 2
0.5
0.63142
Corrector cuarta iteración h yi
1
1
- 1 - 0.4736 2
0.5
0.63161
Aplicando el mismo procedimiento para h = 0.25 y h = 0.125 se obtiene Tabla 6.2 Valores de y para distintos valores de h con Heun h = 0.5 h = 0.25 h = 0.125 x y x y x y 0.0 1.0000 0.00 1.0000 0.00 1.0000 0 0.12 0.8832 5
120
0.5
0.6316
1.0
0.5132
1.5
0.7464
2.0
3.9174
0.25 0.7832 0.25 0.7830 0 0.37 0.6995 5 0.50 0.6322 0.50 0.6323 0 0.62 0.5806 5 0.75 0.5432 0.75 0.5436 0 0.87 0.5211 5 1.00 0.5135 1.00 0.5134 0 1.12 0.5221 5 1.25 0.5523 1.25 0.5501 0 1.37 0.6030 5 1.50 0.7006 1.50 0.6905 0 1.62 0.8296 5 1.75 1.0915 1.75 1.0499 0 1.87 1.4064 5 2.00 2.1965 2.00 2.0031 0
6.2.1.3 Métodos de Runge-Kutta En los métodos de Euler y Heun se aplica la fórmula de recurrencia: y i+i = y i + ( x i, y i) h Donde
121
( x i, y i) = f ( x i, y i) ( x i, y i) = [ f ( x i, y i) + f ( x i+1, y i+1)]
método de Euler método de Heun
1
2
Estos dos métodos tienen los siguientes puntos comunes: 1. Son métodos de un paso, para determinar y i+1 se necesita conocer únicamente los valores de x i y y i del punto anterior. 2. No requiere evaluar ninguna derivada, sino únicamente los valores de la función que representa a la ecuación diferencial. Estas características dan origen a una gran variedad de métodos conocidos como de Runge Kutta. La diferencia entre ellos consiste en la forma como se define la función ( x i, y i):
Segundo Orden. (Método de Ralston) y i+1 = y i + ( k 1 + donde k 1 = f ( x i , y i) ; h k 1) 1
2
3
3
k 2) h k 2 = f ( x i+
3 4
h, y i+
3 4
Tercer Orden. y i+1 = y i + [ (k 1 + 4k 2 + k 3)] h 1
6
donde k 1 = f ( x i, y i) ; y i - hk 1 +2hk 2)
k 2 = f ( x i+
h, y i + h k 1);
1
1
2
2
k 3 = f ( x i+ h,
Cuarto Orden. y i+1 = y i + [ (k 1 + 2k 2 + 2k 3 + k 4)] h 1
6
donde k 1 = f ( x i , y i) ; k 3 = f ( x i + h, y i + hk 2); k 2 = f ( x i+ h, y i + h k 1); k 4 = f ( x i + h, y i + h k 3); Ejemplo 6.3 Resuelva la ecuación diferencial de los ejemplos anteriores por los métodos de Runge-Kutta de segundo, tercer y cuarto orden; sí y (0) = 1, utilizando h = 0.5 1
1
2
2
1
2
1
2
dy dx
y x 2
1
Segundo Orden
y i+1 = y i + ( k 1 + k 2) h 1
2
3
3
122
y i = y (0) = 1 k 1 = f ( x i, y i) = f (0, 1) = 1(0 2-1) =-1 k 2 = f ( x i+ h, y i+ h k 1) x i+ h = 0+ (0.5) = 0.375 y i+ h k 1 = 1 + (0.5)(-1) = 0.625 k 2 = f (0.375, 0.625) = 0.625(0.375 2 - 1) = -0.5371 y i+1 = y (0.5) = 1 + ( (-1) + (-0.5371)) (0.5) = 0.6543
y i = y (1) = 0.5358 k 1 = f (1, 0.5358) = 0 x i+ h = 1.375 y i+ h k 1 = 0.5358 k 2 = f (1.375, 0.5358) = 0.4772 y i+1 = y (1.5) = 0.6948 y i = y (1.5) = 0.6948 k 1 = f (1.5, 0.6948) = 0.8685 x i+ h = 1.875 y i+ h k 1 = 1.0205 k 2 = f (1.875, 1.0205) = 2.5673 y i+1 = y (2) = 1.6953
y i = y (0.5) = 0.6543 k 1 = f (0.5, 0.6543) = -0.4907 x i+ h = 0.875 y i+ h k 1 = 0.4703 k 2 = f (0.875, 0.4703) = -0.1102 y i+1 = y (1) = 0.5358 Tercer Orden
y i+1 = y i + [
(k 1 + 4k 2 + k 3)] h
y i = y (0) = 1 k 1 = f ( x i, y i) = f (0, 1) = 1(0 2-1) = -
k 3 = f (0.5, 0.7969) = -0.5977 y i+1 = y (0.5)
1
= 1 + (-1 + 4(-0.7031) + 0.7969) (0.5) = 0.6325
k 2 = f ( x i+ h, y i+ h k 1) x i+ h = 0+ (0.5) = 0.25 y i+ h k 1 = 1 + (0.5)(-1) = 0.75 k 2 = f (0.25, 0.75) = 0.625(0.375 2 1) = -0.7031 k 3 = f ( x i+h, y i - hk 1 +2hk 2) x i+h = 0 + 0.5 = 0.5 y i - hk 1 +2hk 2 = 1–0.5(-1)+2(0.5)(-0.7031) = 0.7969
y i = y (0.5) = 0.6325 k 1 = f (0.5, 0.6325) = -0.4744 x i+ h = 0.75 y i+ h k 1 = 0.5139 k 2 = f (0.875, 0.4703) = -0.2248 x i+h = 1 y i - hk 1 +2hk 2 = 0.6448 k 3 = f (1, 0.6448) = 0 y i+1 = y (1) = 0.5180 123
y i = y (1) = 0.6995 k 1 = f (1, 0.6995) = 0.8743 x i+ h = 1.75 y i+ h k 1 = 0.9180 k 2 = f (1.75, 0.9180) = 1.8934 x i+h = 2 y i - hk 1 +2hk 2 = 2.1557 k 3 = f (2, 2.1557) = 6.4672 y i+1 = y (1) = 1.9424
y i = y (1) = 0.5180 k 1 = f (1, 0.5180) = 0 x i+ h = 1.25 y i+ h k 1 = 0.5180 k 2 = f (0.875, 0.4703) = -0.2914 x i+h = 1.5 y i - hk 1 +2hk 2 = 0.8094 k 3 = f (1, 0.6448) = 1.0117 y i+1 = y (1) = 0.6995 Cuarto Orden
y i+1 = y i + [
(k 1 + 2k 2 + 2k 3 + k 4)] h
y i = y (0) = 1 k 1 = f ( x i, y i) = f (0, 1) = 1(0 2-1) = 1
k 2 = f ( x i+ h, y i+ h k 1) x i+ h = 0+ (0.5) = 0.25 y i+ h k 1 = 1 + (0.5)(-1) = 0.75 k 2 = f (0.25, 0.75) = 0.625(0.375 2 1) = -0.7031 k 3 = f ( x i + h, y i + hk 2) x i+ h = 0+ (0.5) = 0.25 y i+ h k 2 = 1 + (0.5)(-0.7031) = 0.8242 k 3 = f (0.25, -0.7031) = -0.7727 k 4 = f ( x i + h, y i + h k 3) x i+ h = 0+ (0.5) = 0.5 y i+ h k 3 = 1 + (0.5)(-0.7727) = 0.6136 k 4 = f (0.5, 0.6136) = -0.4602 y i+1 = y (0.5) = 1 + [-1 + 2(-0.7031) + 2(-0.7727) + 4(-0.4602)] (0.5) = 0.6323
y i = y (0.5) = 0.6323 k 1 = f (0.5, 0.6323) = -0.4743 x i+ h = 0.75 y i+ h k 1 = 0.5138 k 2 = f (0.75, 0.5138) = -0.2248 x i+ h = 0.75 y i+ h k 2 = 0.5761 k 3 = f (0.75, -0.5761) = -0.2521 x i+ h = 1 y i+ h k 3 = 0.5063 k 4 = f (1, 0.5063) = 0 y i+1 = y (1) = 0.5133
124
y i = y (1) = 0.5133 k 1 = f (1, 0.5133) = 0 x i+ h = 1.25 y i+ h k 1 = 0.5133 k 2 = f (0.75, 0.5138) = -0.2889 x i+ h = 1.25 y i+ h k 2 = 0.5855 k 3 = f (0.75, -0.5761) = 0.3294 x i+ h = 1.5 y i+ h k 3 = 0.6780 k 4 = f (1, 0.5063) = 0.8475 y i+1 = y (1.5) = 0.6870
y i = y (1.5) = 0.6870 k 1 = f (1.5, 0.6870) = 0.8587 x i+ h = 1.75 y i+ h k 1 = 0.9017 k 2 = f (1.75, 0.9017) = 1.8597 x i+ h = 1.75 y i+ h k 2 = 1.1519 k 3 = f (1.75, 1.1519) = 2.3758 x i+ h = 2 y i+ h k 3 = 1.8749 k 4 = f (2, 1.8749) = 5.6248 y i+1 = y (2) = 1.933
6.3 Métodos de Pasos Múltiples Una técnica alterna para resolver Ecuaciones Diferenciales Ordinarias se puede desarrollar conociendo información de la función en varios puntos, tomando estos como base para predecir el valor de la función en los puntos subsiguientes. dy dx
f x y ,
separando variables e integrando entre los límites i e i+1 y i 1 y i
xi 1
xi
f x, y dx
resolviendo la integral se pueden encontrar los valores de y i+1 conociendo y i. Método de Milne Este es un método predictor-corrector que utiliza información en los primeros cuatro Esta información se puede obtener aplicando alguno de los métodos vistos anteriormente. Para resolver la integral se usa las formulas de integración numérica vistas en la unidad anterior
125
Predictor y i
1
y i
3
4h 3
2 f ( xi , y i ) f ( xi 1 , y i 1 ) 2 f ( xi
2
, yi 2 )
Corrector y i
1
y i
1
h 3
2 f ( xi 1 , y i 1 ) 4 f ( x i , y i ) 2 f ( xi 1 , y i 1 )
Un tipo de fórmulas que tienen la forma general descrita anteriormente son las fórmulas de Adams Fórmula abierta de n-ésimo orden
(Adams-Bashforth) n 1
yi 1 yi h k f i k o h n 1
k 0
Fórmula abierta de n-ésimo orden
(Adams-Moulton) n 1
y i 1 y i h k f i 1 k o h n 1
k 0
donde k son coeficientes reportados en la bibliografía. Combinando estas dos fórmulas en un esquema de predictor corrector se puede desarrollar un método para encontrar la solución de las ecuaciones diferenciales ordinarias. La información de los puntos necesarios para iniciar el procedimiento se obtiene generalmente a partir de un método de un solo paso, con un orden suficiente para que esta información sea confiable.
Ejemplo 6.4 Resolver por el método de Adams de cuarto orden la ecuación. dy dx
y x 2
1
y i+1 = y i + h ( 0 f i-0 + 1 f i-1 + 2 f i-2 + 3 f i-3)
Predictor
59 37 9 55 f i f i 1 f i 2 f i 3 24 24 24 24
y i 1 y i h
Corrector
y i+1 = y i + h ( 0 f i+1-0 + 1 f i+1-1 + 2 f i+1-2 + 3 f i+1-3)
126
19 5 9 f i 1 f i f i 1 24 24 24
y i 1 y i h
1 24
f i 2
Cálculo de los puntos iniciales por el método de Runge-Kutta de cuarto orden
Primer Paso x i = 0, y i = 1 x i - 1 = -0.5, y i - 1 = 1.581052 x i - 2 = -1, y i - 2 = 1.947028 x i - 3 = -1.5, y i - 3 = 1.453834
Segundo Paso x i = 0.5, y i = 0.6379742 x i - 1 = 0, y i - 1 = 1 x i - 2 = -0.5, y i - 2 = 1.581052 x i - 3 = -1, y i - 3 = 1.947028
Predictor y (0.5) = 0.9709569 Corrector y (0.5) = 0.5911456 y (0.5) = 0.6445565 y (0.5) = 0.6370456 y (0.5) = 0.6381018 y (0.5) = 0.6379533 y (0.5) = 0.6379742
Predictor y (1) = 0.735548 Corrector y (1) = 0.5319068
6.4 Sistemas de Ecuaciones Diferenciales Ordinarias. Todo sistema de generalmente como
ecuaciones
dy1 dx dy2 dx dyn dx
diferenciales
puede
representarse
f 1 x, y1 , y2 ,... yn f 2 x, y1 , y2 ,... yn
f n x, y1 , y2 ,... yn
La solución de este sistema requiere de n condiciones iniciales conocidas para un valor inicial de x .
127
Una ecuación diferencial de orden superior puede escribirse como un sistema de ecuaciones diferenciales de primer orden. Escsriba la ecuación diferencial ordinaria y (n) = f ( x , y , y’ , y´´ , ..., y (n - 1)) como un sistema de ecuaciones de primer orden haciendo las sustituciones
y 1 = y , y 2 = y’ , ..., y n = y (n - 1) Entonces:
y´ 1 = y 2 y´ 2 = y 3 y’ n = f ( x , y 1, y 2, y 3, ..., y n ) es un sistema de n ecuaciones diferenciales ordinarias. Por ejemplo, considere el problema de valor inicial.
y´´´ -3 y’’ – y’y = 0 y (0) = 0
y´ (0) = 1 y ´´ (0) = -1
Despeje en la ecuación diferencial, para su derivada de mayor orden escribiendo y´´ ´ en términos de x y de sus derivadas de orden menor y´´´ = 3 y´´ + y´y . Si hacemos las sustituciones
y 1 = y y 2 = y’ y 3 = y´´ entonces
y´ 1 = y 2 y´ 2 = y 3 y 3’ = 3 y 3 + y 2 y 1 con las condiciones iniciales
y 1 (0) = 0 y 2 (0) = 1 y 3 (0) = -1
128
Ejemplo 6.5 Resolver el problema de valores en la frontera definido por la ecuación: d 2 y dx 2
y
0
si y (0) = 1, y (0) = 2; y calcular el valor de y (1).
Analíticamente Teorema.
“Si la ecuación auxiliar m 2 + bm +c = 0 tiene las raíces complejas s
ti, entonces la solución general de y + by + cy = 0 es y = esx (c1 cos tx + c 2 sen tx )” En el ejemplo, para la ecuación auxiliar b = 0 y c = 1 m =i Por ello, s = 0 y t = 1, y la solución general queda:
m 2 + 1 = 0
y = e(0)x (c1 cos (1) x + c 2 sen (1) x ) y = c1 cos x + c 2 sen x y = c 2 cos x – c1 sen x Sustituyendo las condiciones en la frontera
y (0) = c1 cos (0) + c 2 sen (0) = 1 y (0) = c 2 cos (0) – c1 sen (0) = 2
c1 = 1 c 2 = 2
y = cos x + 2sen x y (1) = cos (1) + 2sen (1) = 2.223244
Utilizando el paquete Polymath, para x =1, y = 2.2232
129
Numéricamente Usando el método de Runge-Kutta de segundo orden (método de Ralston) con h = 0.5, y (0) = 1, y (0) = 2; d 2 y
y 0 2 dx d dy1 y1 0 dx dx dy1 y2 dx dy2 dy2 y1 0 y1 dx dx
Ecuaciones del método:-
y j, i+1 = y j, i + (
k 1, j +
k 2, j ) h
k 1, j = f j ( x i, y 1, i , y 2, i ,..., y n, i ); k 2, j = f j ( x i+ h, y 1, i + h k 1, 1, y 2, i + h k 1, 2,..., y n, i + h k 1, n,)
x i = 0; y 1, i = 1; y 2, i = 2 k 1, 1 = f 1 (0, 1, 2) = 2 k 1, 2 = f 2 (0, 1, 2) = -1 x i+ h = 0 + (0.5) = 0.375 y 1, i + h k 1, 1 = 1 + (0.5)(2) = 1.75 y 2, i + h k 1, 2 = 2 + (0.5)(-1) = 1.625
k 2,
y 1 (0.5) = 1 + ( (2) + (1.625) (0.5) = 1.875 y 2 (0.5) = 2 + ( (-1) + (-1.75) (0.5) = 1.25
= f 1 (0.375, 1.75, 1.625) = 1.625 k 2, 2 = f 2 (0.375, 1.75, 1.625) = 1.75 1
130
x i = 0.5; y 1, i = 1.875; y 2, i = 1.25 k 1, 1 = f 1 (0.5, 1.875, 1.25) = 1.25 k 1, 2 = f 2 (0.5, 1.875, 1.25) = -1.875 x i+ h = 0.5 + (0.5) = 0.875 y 1, i + h k 1, 1 = 1.875 +
= 0.546875 k 2, 2 = f 2 (0.875, 0.546875) = -2.34375
(0.5)(1.25) = 2.34375 y 2 , i + h k 1, 2 = 1.25 + (0.5)(1.875) = 0.546875
k 2,
= f 1 (0.875, 0.546875) 1
2.34375,
y 1 (1) = 1.875 + [( (1.25) + (0.546875)](0.5) = 2.265625 y 2 (1) = 1.25 + [ (-1.875) + 2.34375)](0.5) = 0.15625
2.34375,
(-
131
6.5 Aplicaciones
Ejemplo 6.6 Usando el método de Runge-Kutta de cuarto orden con h = 0.5, y las mismas condiciones iniciales Ecuaciones del método: y j, i+1 = y j, i + (k 1, j + 2 k 2, j + 2 k 3, j , + k 4, j ) h k 1, j = f j ( x i, y 1, i , y 2, i ,..., y n, i ); k 2, j = f j ( x i+
h, y 1, i +
h k 1, 1, y 2, i +
h k 1,2,... , y n, i + h k 1, n)
k 3, j = f j ( x i+
h, y 1, i +
h k 2, 1, y 2, i +
h k 2, 2j ,... , y n, i + h k 2, n)
k 4, j = f j ( x i+ h, y 1, i + h k 3, 1, y 2, i +h k 3, 2,... , y n, i + h k 3, n)
x i = 0; y 1, i = 1; y 2, i = 2 k 1, 1 = f 1 (0, 1, 2) = 2 k 1, 2 = f 2 (0, 1, 2) = -1
k 3, 2 = f 2 (0.25, 1.4375, 1.625) = 1.4375
x i + h = 0 + (0.5) = 0.25 y 1, i + h k 1, 1 = 1 + (0.5)(2) = 1.5
y 2, i + h k 1, 2 = 2 + (0.5)(-1) = 1.75 k 2, 1 = f 1 (0.25, 1.5, 1.75) = 1.75 k 2, 2 = f 2 (0.25, 1.5, 1.75) = -1.5
x i + h = 0.25 y 1, i + h k 2, 1 = 1 +
(0.5)(1.75) = 1.4375 y 2, i + h k 2, 2 = 2 + (0.5)(-1.5) = 1.625 k 3, 1 = f 1 (0.25, 1.4375, 1.625) = 1.625
x i + h = 0.5 y 1, i + h k 3, 1 = 1 + (0.5)(1.625) = 1.8125 y 2, i +h k 3, 2 = 2 + (0.5)(-1.4375) = 1.28125 k 4, 1 = f 1 (0.5, 1.8125, 1.28125 = 1.28125 k 4, 2 = f 2 (0.5, 1.8125, 1.28125) = 1.8125
y 1 (0.5) =1+ [2+2(1.75)+2(1.625)+1.28125](0 .5) = 1.835938 y 2 (0.5)
132
=2+ [-1+2(-1.5)+2(-1.4375)1.8125](0.5) = 1.276042
x i = 0.5; y 1, i = 1.835938; y 2, i = 1.276042 k 1,1 = f 1 (0.5, 1.835938, 1.276042) = 1.276042 k 1,2 = f 2 (0.5, 1.835938, 1.276042) = -1.835938
x i + h = 0.5 + (0.5) = 0.75 y 1,i + hk 1,1 =1.835938 + (0.5)(1.276042) = 2.154948 y 2,i hk 1,2 + =1.276042+ (0.5)(-1.835938) = 0.817057 k 2, 1 = f 1 (0.75, 2.154948, 0.817057) = 0.817057 k 2, 2 = f 2 (0.75, 2.154948, 0.817057) = -2.154948
x i + h = 0.75 y 1,i + hk 2,1
x i + h = 1 y 1, i + h k 3, 1 = 1.835938 + (0.5)(0.737305) = 2.204590 y 2, i +h k 3, 2 = 1.276042 + (0.5)(2.040202) = 0.255941 k 4, 1 = f 1 (1, 2.204590, 0.255941) = 0.255941 k 4, 2 = f 2 (1, 2.204590, 0.255941) = -2.204590
y 1 (1) = 1.835938 + [1.276042 + 2(0.817057) + 2(0.737305) + 0.255941](0.5) = 2.222663 y 2 (1) =1.276042 + [-1.835938 + 2(2.154948) + 2(-2.040202) 2.204590](0.5) = 0.240139
=1.835938+
(0.5)(0.817057)
= 2.040202 + y 2,i hk 2,2 =1.276042+ (0.5)(-2.154948) = 0.737305 k 3,1 = f 1 (0.75, 2.040202, 0.737305) = 0.737305 k 3, 2 = f 2 (0.75, 2.040202, 0.737305) = -2.040202
133
