Me To Dos

CAPÍTULO Solución Numérica de Ecuaciones Diferenciales

2.1 Introducción Las leyes fundamentales de la física, mecánica, electricidad y termodinámica están basadas con frecuencia en observaciones experimentales que explican variaciones en las propiedades físicas y estados de los sistemas. Estas leyes, mas que describir directamente el estado de los sistemas físicos se expresan en términos de los cambios espaciales y temporales de las variables intervinientes. Es por ello que las ecuaciones diferenciales tienen importancia fundamental en las aplicaciones de ingeniería ya que numerosos procesos físicos son idealizados (o modelizados) matemáticamente por estas ecuaciones. Tales ecuaciones son a veces conocidas como ecuaciones de razón ya que expresan la razón de cambio de una variable como una función de las variables y parámetros del problema. Podemos citar como ejemplos de las leyes fundamentales que se escriben en términos de la razón de cambio de las variables a: •

Segunda ley de Newton del movimiento dv F = dt m

donde v es la velocidad, F es la fuerza, m es la masa y t el tiempo. •

Ley del calor de Fourier q = − k' ⋅

dT dx

donde q es el flujo de calor, k’ es la conductividad térmica, T es la temperatura y x la variable espacial. •

Ley del difusión de Fick J = − D⋅

dc dx

donde q es el flujo másico, D es el coeficiente de difusión, c es la concentración y x la variable espacial. CAPITULO – SOLUCIÓN NUMÉRICA DE ECUACIONES DIFERENCIALES

CÁTEDRA MÉTODOS NUMERICOS

Pág.1

•

Ley de Faraday δVL = L⋅

di dt

donde VL es la caida de voltaje, L es la inductancia, i es la corriente y t la variable temporal. La mayoría de las ecuaciones diferenciales de importancia práctica no se pueden resolver mediante métodos analíticos de calculo, y es debido a esto que los métodos numéricos han adquirido una importancia extraordinaria en todos los campos de la ingeniería sobre todo a partir de la disponibilidad de computadoras que soportan grandes volúmenes de cálculo. Solo a forma de síntesis de lo estudiado en los cursos de matemática superior, se enumeran a continuación un conjunto de definiciones acerca de las ecuaciones diferenciales, su caracterización y solución. •

Variables independientes y dependientes: Se dice que una variable de una ecuación diferencial es independiente si existen una o más derivadas con respecto a esa variable. Una variable es dependiente cuando existen derivadas de esa variable. El número de variables dependientes se denomina a menudo en los problemas de ingeniería como “grados de libertad” del problema.

•

Ecuaciones Diferenciales ordinarias y parciales: Si en una ecuación diferencial hay una sola variable independiente, las derivadas son totales y a la ecuación se la denomina ordinaria. Por el contrario, si aparecen dos o más variables independientes las derivadas serán parciales y la ecuación será diferencial parcial. Orden de una Ecuación diferencial: es la derivada de mayor orden que aparece en la ecuación.

•

•

Ecuación diferencial lineal: Una ecuación diferencial es lineal si en ella no aparecen potencias de la variable dependiente ni de sus derivadas ni productos de la variable dependiente por sus derivadas o productos entre derivadas. Una ecuación diferencial ordinaria lineal es aquella que se ajusta a la forma general: an (x) ⋅ y (n ) +....... +a1 ( x ) ⋅ y ′+a0 ( x ) ⋅ y = f ( x )

donde y (n) es la n-ésima derivada de y con respecto a x, y a (x) y f (x) son funciones específicas de x. CAPITULO – SOLUCIÓN NUMÉRICA DE ECUACIONES DIFERENCIALES

CÁTEDRA MÉTODOS NUMERICOS

Pág.2

•

Solución de una ecuación diferencial: Es cualquier relación funcional que no incluya derivadas o integrales de funciones desconocidas y que la verifique idénticamente por sustitución directa.

2.2 Solución de una Ecuación Diferencial Dada un a ecuación diferencial ordinaria de orden n, cuya forma general puede escribirse como: (n) F(x,y, y’,y’’.........y ) = 0 (2.1) su solución general resulta dependiente de n constantes arbitrarias y puede escribirse en la forma: G (x,y,c1, c2, ...., cn) = 0 (2.2) Se requiere de n condiciones para obtener una solución única. Gráficamente, la ecuación (2.2) representa a una familia de curvas planas, cada una de ellas obtenidas para valores particulares de las n constantes c1, c2, ...., cn como se muestra en la figura siguiente.

y

G

=

4

G

3

0 =

0 G

G

=

2

=

1

0 0

x

Cada una de estas curvas corresponde a una solución particular de la ecuación (2.1) y analíticamente puede obtenerse “sujetando” la solución general (2.2) a n condiciones independientes que permitan valuar las constantes arbitrarias. Dependiendo de cómo se establezcan estas condiciones, se distinguen dos tipos de problemas: los llamados de valores iniciales y los de valores de frontera. a) Problema de valores iniciales: Un problema de valores iniciales está gobernado por una ecuación diferencial de orden n y un conjunto de n condiciones independientes, todas CAPITULO – SOLUCIÓN NUMÉRICA DE ECUACIONES DIFERENCIALES

CÁTEDRA MÉTODOS NUMERICOS

Pág.3

ellas válidas para el mismo punto inicial. Es decir, si se especifican todas las condiciones en el mismo valor de la variable independiente al inicio del dominio al problema se lo denomina como problema de condiciones iniciales o valor inicial. Si en la ecuación (2.1) x = a es el punto inicial, en un problema de valores iniciales, las condiciones de contorno resultarán del tipo: (n)

y(a) = y0 , y’(a) = y’0 , y’’(a)=y’’0................. y (a)= y

(n)

0

y se tratará de encontrar una solución particular para (2.1) que verifique las anteriores, tal como se muestra en la figura siguiente: y

x =

x

a

b) Problemas de frontera: Por el contrario, en los problemas de valores de frontera deben establecerse condiciones en todos y cada uno de los puntos que constituyen la frontera del dominio de definición del problema. Es decir, si las condiciones se especifican sobre distintos valores de la variable independiente en la frontera del dominio al problema se lo denomina problema de condiciones de contorno o valor límite. En particular, en el espacio unidimensional hay dos puntos de frontera, por ejemplo x = a y x = b si el dominio de definición es el intervalo cerrado [a,b].

y

x =

a

x =

b

x

CAPITULO – SOLUCIÓN NUMÉRICA DE ECUACIONES DIFERENCIALES

CÁTEDRA MÉTODOS NUMERICOS

Pág.4

Básicamente, la solución numérica de ecuaciones diferenciales consisten en sustituir el dominio continuo de soluciones por uno discreto, es decir formado por puntos igualmente espaciados entre sí. Así, en un problema de valores iniciales el dominio de definición del problema, x≥ a, se sustituye por el conjunto infinito de puntos: x0 = a, x1 = x0 + ∆ , x2 = x0 + 2∆ , x3 = x0 + 3∆ .. . En los problemas de frontera, el dominio se discretiza en x0 = a, x1 = x0 + ∆ , x2 = x0 + 2∆ , x3 = x0 + 3∆ ....... xn = x0 + n∆ = b, obtenidos al dividir el intervalo [a,b] en n partes iguales. Habiéndose discretizado el problema continuo se tratará de obtener una solución para los puntos considerados, y esto se hará, en general, sustituyendo las derivadas que aparezcan en la ecuación diferencial y en las condiciones de contorno, por expresiones numéricas de derivación que proporcionen una aproximación a las derivadas o tratando de integrar la ecuación y reemplazando al proceso de integración por una fórmula numérica que se aproxima a la integral.

CAPITULO – SOLUCIÓN NUMÉRICA DE ECUACIONES DIFERENCIALES

CÁTEDRA MÉTODOS NUMERICOS

Pág.5

2.3 Ecuaciones Diferenciales Ordinarias 2.3.1 Problemas de Condiciones Iniciales En este tipo de problemas, como ya lo expresamos anteriormente, debemos obtener valores aproximados de la función solución de una ecuación diferencial en m puntos de un intervalo del dominio de la misma, partiendo de condiciones iniciales conocidas de la función a determinar en el extremo inicial del intervalo. En los denominados métodos de un paso, que describiremos a continuación, el valor de la función en el primer punto interior del intervalo se calcula a partir del valor conocido de la función en punto inicial del intervalo. De la misma forma el valor de la función incógnita en el i-ésimo punto del dominio, xi, se calcula a partir del valor de la función en el punto x i-1. A su vez la función en el punto i+1 se calcula a partir del valor de la función en el punto i, . Es decir que se calcula la función en un punto cualquiera del intervalo partiendo de la solución obtenida para el punto anterior, y así sucesivamente. De esta manera, con la aplicación recurrente los algoritmos correspondientes a cualquiera de estos métodos, a lo largo de todo el intervalo de integración, se obtiene la denominada trayectoria de la solución. Los métodos que estudiaremos están limitados a la solución de ecuaciones diferenciales ordinarias de primer orden. Esto, que en principio parece una severa limitación de los mismos, no lo es tanto si recordamos que cualquier ecuación diferencial ordinaria de orden n puede ser transformada en un sistema de n ecuaciones diferenciales ordinarias de primer orden. Por otra parte el procedimiento para obtener la solución de un sistema de EDOs de primer orden es una sencilla extensión del procedimiento para obtener la solución de una sola EDO por lo que estos métodos pueden ser aplicados sin ninguna dificultad a la solución de EDOs de orden n. Dada una ecuación diferencial ordinaria de primer orden de la forma:

dy = f ( x, y ) dx

la solución numérica tendrá la forma:

y i +1 = y i + φ ⋅ h

De acuerdo con esta ecuación, una pendiente estimada  se usa para extrapolar desde un valor anterior yi a un nuevo valor yi+1 en una distancia h. CAPITULO – SOLUCIÓN NUMÉRICA DE ECUACIONES DIFERENCIALES

CÁTEDRA MÉTODOS NUMERICOS

Pág.6

Todos los denominados métodos de un paso, que veremos en este capítulo, se pueden expresar en esta forma general, que solo va a diferir en la manera en la cual se estima la pendiente.

Método de Euler Como sabemos, la primera derivada de una función f(x) valuada en un punto del dominio xi representa la pendiente de la función en xi. Si la ecuación diferencial ordinaria a resolver la expresábamos como: dy = f ( x, y ) dx

entonces la pendiente en xi e yi será:

φ = f ( xi , y i )

En este procedimiento, la pendiente φ = f ( xi , y i ) calculada en el extremo inicial del intervalo, xi, es tomada como una aproximación de la pendiente promedio sobre todo el intervalo. Entonces se predice un nuevo valor de y por medio de la pendiente en xi , que CAPITULO – SOLUCIÓN NUMÉRICA DE ECUACIONES DIFERENCIALES

CÁTEDRA MÉTODOS NUMERICOS

Pág.7

habrá de extrapolarse en forma lineal sobre el paso de longitud h , mediante la aplicación de la expresión: y i +1 = y i + φ ⋅ h = y i + f ( xi , y i ) ⋅ h

Ejemplo 1: Usar el método de Euler para resolver numéricamente la ecuación: dy = - 2 ⋅ x 3 + 12 ⋅ x 2 - 20 ⋅ x + 8.5 dx

desde x = 0 hasta x = 4 con un paso de h = 0.5. La condición inicial en x = 0 es y = 1. Para el primer paso: y ( 0.5 ) = y ( 0 ) + f ( 0 ,1) ⋅ 0.5

donde y (0) = 1 y la pendiente estimada en x = 0 es:  dy  f (0,1) =   ( x =0,y =1) = - 2 ⋅ ( 0 ) 3 + 12 ⋅ ( 0 ) 2 - 20 ⋅ ( 0 ) + 8.5 = 8.5  dx 

reemplazando:

y ( 0.5 ) =1.0 + 8.5 ⋅ ( 0.5 ) = 5.25

La función solución exacta es: y = - 0.5 ⋅ x 4 + 4 ⋅ x 3 - 10 ⋅ x 2 +8.5 ⋅ x +1

y para x = 0.5 la solución exacta es: y = - 0.5 ⋅ ( 0.5 ) 4 + 4 ⋅ ( 0.5 ) 3 - 10 ⋅ ( 0.5 ) 2 + 8.5 ⋅ ( 0.5 ) +1 = 3.21875

Así, el error absoluto es: E t = verdadero – aproximado = 3.21875 – 5.25 = -2.03125 y el error relativo porcentual: εt % =

Et − 2.03125 × 100 = × 100 = − 63.1% verdadero 3.21875

CAPITULO – SOLUCIÓN NUMÉRICA DE ECUACIONES DIFERENCIALES

CÁTEDRA MÉTODOS NUMERICOS

Pág.8

Para el segundo paso: y (1) = y ( 0.5 ) + f ( 0.5 , 5.25) ⋅ 0.5 ⇒

[

]

y (1) = 5.25 + - 2 ⋅ ( 0.5 ) 3 +12 ⋅ ( 0.5 ) 2 - 20 ⋅ ( 0.5 ) + 8.5 ⋅ 0.5 ⇒ y (1) = 5.875

y para x = 1.0 la solución exacta es: y = - 0.5 ⋅ (1.0 )4 + 4 ⋅ (1.0 ) 3 - 10 ⋅ (1.0 ) 2 + 8.5 ⋅ (1.0 ) +1 = 3.0

Así, el error absoluto es: E t = verdadero – aproximado = 3.0 – 5.875 = -2.875 y el error relativo porcentual: εt % =

Et − 2.875 × 100 = × 100 = − 95.8% verdadero 3.0

y si seguimos calculando la solución para el resto de los puntos del dominio correspondientes a un paso h = 0.5, obtenemos la siguiente tabla: Yverdadero

YEuler

r% local

0.0

1.00000

1.00000

-

-

0.5

3.21875

5.25000

-63.1

-63.1

1.0

3.00000

5.87500

-28.0

-95.8

1.5

2.21875

5.12500

-1.41

-131.0

2.0

2.00000

4.50000

20.5

-125.0

2.5

2.71875

4.75000

17.3

-74.7

3.0

4.00000

5.87500

4.0

-46.9

3.5

4.71875

7.12500

-11.3

-51.0

4.0

3.00000

7.00000

-53.0

-133.3

X

r%

global

Si graficamos la solución verdadera y la solución calculada vemos que si bien el error es considerable, tal como se desprende de los cálculos de errores hechos anteriormente, la solución calculada “copia” de manera bastante aproximada la forma de la solución exacta. Intuitivamente es de esperar que si disminuimos el paso h, esto es si la extrapolación de

CAPITULO – SOLUCIÓN NUMÉRICA DE ECUACIONES DIFERENCIALES

CÁTEDRA MÉTODOS NUMERICOS

Pág.9

la pendiente calculada al inicio del intervalo se hace sobre un intervalo menor, entonces la solución calculada se aproxime a la exacta.

Ejemplo 2. Si para la misma ecuación diferencial del ejemplo anterior, con las mismas condiciones iniciales, tomamos ahora un paso h = 0.25. Es claro que al tomar ahora un paso, es decir un incremento de x, que es la mitad del anterior, los puntos donde deberemos valuar a la expresión solución del método de Euler se duplicara. Esto significa que para aproximar a la función solución en el mismo intervalo el esfuerzo computacional se incrementara. Los resultados obtenidos, que compararemos con los obtenidos en el ejemplo anterior, se condensaron en la siguiente tabla: e%r

x

y verdadero

yEuler

0.00

1.00000

1.00000

-

0.25

2.56055

3.12500

-22.0

0.50

3.21875

4.17969

-29.9

0.75

3.27930

4.49219

-37.0

1.00

3.00000

4.34375

-44.8

1.25

2.59180

3.96875

-53.1

1.50

2.21875

3.55469

-60.2

1.75

1.99805

3.24219

-62.3

2.00

2.00000

3.12500

-56.3

2.25

2.24805

3.25000

-44.6

2.50

2.71875

3.61719

-33.0

2.75

3.34180

4.17969

-25.1

3.00

4.00000

4.84375

-21.1

3.25

4.52930

5.46875

-20.7

3.50

4.71875

5.86719

-24.3

global

CAPITULO – SOLUCIÓN NUMÉRICA DE ECUACIONES DIFERENCIALES

CÁTEDRA MÉTODOS NUMERICOS

Pág.10

3.75

4.31055

5.80469

-34.7

4.00

3.00000

5.00000

-66.7

Si ahora graficamos las soluciones calculadas y la solución exacta obtenemos:

Se puede apreciar que, si bien los errores siguen siendo considerables, estos han disminuido apreciablemente respecto de aquellos obtenidos cuando el paso es h = 0.5. Si se sigue disminuyendo el paso h , los errores van disminuyendo según una ley que se grafica a continuación:

Podemos observar que se puede disminuir el error al disminuir el tamaño del paso y que disminuyendo el paso a la mitad se disminuye el error relativo aproximadamente a la mitad. Sin embargo se observa que para h = 0.001 el error relativo porcentual todavía

CAPITULO – SOLUCIÓN NUMÉRICA DE ECUACIONES DIFERENCIALES

CÁTEDRA MÉTODOS NUMERICOS

Pág.11

excede el 0.1 %, es decir que el método de Euler exige un gran esfuerzo computacional para dar niveles de error aceptables. Es claro que una fuente fundamental de error en este método es que la derivada al inicio del intervalo se aplica a través de todo el intervalo. 2.3.1.1.1 Análisis del error para el método de Euler Se puede tener cierto conocimiento acerca de la magnitud y propiedades del error que se comete al aplicar el método de Euler si derivamos a este directamente de la expansión de la serie de Taylor. Es necesario aclarar que esta estimación del error inherente al método en si mismo no contempla, por lo tanto, los errores de redondeo que se producen al almacenar los números provenientes del calculo en la memoria de las computadoras. En el error inherente al método, denominado también error de truncamiento, se distinguen dos componentes: el error de truncamiento local, que resulta de la aplicación del método en cuestión sobre un paso , y el error de truncamiento propagado que resulta de las aproximaciones producidas en los pasos previos. La suma de los dos es el total o el error de truncamiento global. La ecuación diferencial sujeta a integración es de la forma general: y ′ = f ( x, y )

donde y’ = dy / dx, y es la variable dependiente, y x es la variable independiente. Como sabemos, si la función solución y tiene derivadas continuas, puede representarse como una expansión de la serie de Taylor con respecto a los valores de inicio (xi , yi), es decir: y'' y i +1 = y i + y 'i ⋅ h + i ⋅ h2 + 2!

y (n) y 'i'' 3 ⋅ h + .......... + i ⋅ h n + Rn 3! n!

donde h = xi+1 - xi y Rn es el residuo definido como: Rn =

y ( n +1) ( ξ ) n +1 ⋅h ( n + 1) !

donde  esta en algún lugar en el intervalo xi a xi+1. La ecuación anterior también puede expresarse como:

y i +1 = y i + f(x i , y i ) ⋅ h +

f ' (xi , y i ) 2 f ' ' (xi , y i ) 3 f (n −1) (xi , y i ) n ⋅h + ⋅ h + .. + ⋅ h + O(h n +1 ) 2! 3! n!

CAPITULO – SOLUCIÓN NUMÉRICA DE ECUACIONES DIFERENCIALES

CÁTEDRA MÉTODOS NUMERICOS

Pág.12

donde O( h n +1 ) especifica que el error de funcionamiento local, o residuo, es proporcional al tamaño del paso elevado a la potencia (n+1)-ésima. Al comparar esta última ecuación con la del método de Euler, vemos que esta corresponde a la serie de Taylor que incluye hasta el término de la primera derivada:

y i +1 = y i + f(xi , y i ) ⋅ h = y i + y ′⋅ h = y i + φ⋅ h

De la comparación surge claramente que ocurre un error de truncamiento porque aproximamos la solución verdadera mediante un número finito de términos de la serie de Taylor. Entonces el error de truncamiento en el método de Euler es atribuible a los términos remanentes de la expansión de Taylor que no se incluyeron en la ecuación, esto es:

Et =

f ' (xi , y i ) 2 f '' (xi , y i ) y 'i'' 3 f (n −1)(xi , y i ) n ⋅h + ⋅ h + .. + ⋅ h + O(h n +1 ) 2! 3! n!

donde Et es el error de truncamiento local verdadero. Para h lo suficientemente pequeño (menor que 1) los errores en los términos de la ecuación anterior generalmente disminuyen al aumentar el orden por lo que el resultado es a menudo representado mediante la aproximación: f ' (x i , y i ) 2 Ea = ⋅h 2!

o también como E a = O(h 2 ) donde Ea es el error de truncamiento local aproximado. Esto significa que el error aproximado local es proporcional al cuadrado del tamaño del paso y a la primera derivada de la ecuación diferencial. Es decir que si h disminuye el error disminuye con el cuadrado de h. El desarrollo de la serie de Taylor proporciona solo un estimado del error de truncamiento local, pero no proporciona una medida del error de truncamiento propagado y por lo tanto tampoco del error de truncamiento global. Sin embargo se puede demostrar que el error de truncamiento global es O(h), es decir de orden 1 o proporcional al tamaño del paso. Es obvio que el método de Euler proporcionara predicciones libres de error si la función solución de la ecuación diferencial es lineal, debido a que para una función lineal la derivada segunda es nula, o visto de otra forma porque Euler utiliza segmentos de línea CAPITULO – SOLUCIÓN NUMÉRICA DE ECUACIONES DIFERENCIALES

CÁTEDRA MÉTODOS NUMERICOS

Pág.13

recta para aproximar la solución. De allí que el método de Euler se conozca como uno de primer orden. De la misma manera, podemos extrapolar, los métodos de orden n-ésimo darán resultados perfectos si la función solución de la ecuación diferencial a resolver es un polinomio de orden n-ésimo.

Método de Heun Como ya hemos mencionado la fuente fundamental de error en el método de Euler es que la derivada de la función al inicio del intervalo se aplica a través de toda la longitud del mismo. El método de Heun propone mejorar la estimación de la pendiente calculando las derivadas al comienzo y al final del intervalo. Luego, estas derivadas se promedian para obtener una estimación mejorada de la pendiente para todo el intervalo. Este procedimiento se ilustra en la siguiente figura:

En este método primero se calcula la derivada al comienzo del intervalo:

y i′ = f ( xi , y i ) y se usa para extrapolar linealmente a un valor de y i+1, que distinguiremos con el supraíndice 0 por ser la misma una predicción intermedia y no la respuesta final como lo hubiera sido en el método de Euler: y i0+1 = y i + f ( xi , y i ) ⋅ h

CAPITULO – SOLUCIÓN NUMÉRICA DE ECUACIONES DIFERENCIALES

CÁTEDRA MÉTODOS NUMERICOS

Pág.14

Esta ecuación es la llamada ecuación predictor y da una predicción intermedia de yi+1, 0 denominada y i+ 1 , que permite el cálculo de una estimación de la pendiente al final del intervalo:

(

y 'i +1 = f xi +1 , y i0+1

)

Luego se obtiene una pendiente promedio para el intervalo:

(

y ′ + y i′ +1 f ( xi , y i ) + f xi +1 , y i0+1 y′ = i = 2 2

)

Esta pendiente promedio se usa para extrapolar nuevamente desde yi a yi+1 usando el método, ya conocido, de Euler y i +1 = y i +

(

)

f ( xi , y i ) + f xi +1 , y i0+1 ⋅h 2

la cual es conocida como ecuación corrector. El método de Heun es un procedimiento predictor – corrector o multipaso. Este puede expresarse en forma resumida como: PREDICTOR

CORRECTOR

y i0+1 = y i + f ( xi , y i ) ⋅ h y i +1 = y i +

(

)

f ( xi , y i ) + f xi +1 , y i0+1 ⋅h 2

Como la ecuación corrector tiene a yi+1 en ambos miembros, entonces puede aplicarse en forma iterativa, y de esta manera obtener una estimación mejorada de y i+1. El criterio de iteración no necesariamente converge a la solución verdadera sino que lo hará sobre una solución aproximada con un error de truncamiento finito. Como sucede con la mayoría de los métodos iterativos un criterio de terminación para la convergencia del corrector está dado por la ecuación:

εa =

y ij+1 - y ij-+11 y ij+1

⋅100 %

-1 donde yij+1 e y ij+ 1 resultan de dos iteraciones sucesivas, la actual y la anterior respectivamente. Ejemplo. Resolver mediante el método de Heun la siguiente ecuación diferencial:

CAPITULO – SOLUCIÓN NUMÉRICA DE ECUACIONES DIFERENCIALES

CÁTEDRA MÉTODOS NUMERICOS

Pág.15

dy = y ′ = 4 ⋅ e0.8 x - 0.5 ⋅ y dx

desde x = 0 a x = 4, con un tamaño de paso h = 1, y la condición inicial en x = 0 es y = 2. Antes de resolver el problema en forma numérica y a los efectos de comparar error calculamos la solución exacta que es : y=

(

)

4 ⋅ e0.8 x - e-0.5 x + 2 ⋅ e -0.5 x 1.3

Ahora, ya abordando el problema desde el punto de vista numérico, inicialmente calculamos la pendiente en ( x0 , y0 ): ′ = 4 ⋅ e0 - 0.5 ⋅ 2 = 3 y0

Se puede apreciar la falta de precisión de este calculo si tenemos en cuenta que la pendiente promedio real para el intervalo [ 0 , 1 ] es 4.1946. Usando la ecuación predictor obtenemos un estimado de y en x = 1: ′ ⋅ h =2 +3 ⋅1 =5 y10 =y 0 + y 0

Obsérvese que este sería el resultado que se obtendría con Euler, con un error relativo porcentual del 25.3 %. 0 Ahora para mejorar el estimado anterior en Yi+1 usamos el valor y 1 para predecir la pendiente en el final del intervalo:

(

)

y1′ = f x1 , y10 = 4 ⋅ e0.8 ⋅1 - 0.5 ⋅ 5 = 6.402164

Esta pendiente se promedia con la pendiente inicial para obtener una pendiente promedio sobre el intervalo [ 0,1] y′ =

3 + 6.402164 = 4.701082 2

Podemos ver que esta pendiente promedio calculada es mas cercana a la pendiente promedio verdadera de 4.1946. Ahora aplicamos la ecuación corrector, sin iterar: y1 = 2 + 4.701082 ⋅1 = 6.701082

CAPITULO – SOLUCIÓN NUMÉRICA DE ECUACIONES DIFERENCIALES

CÁTEDRA MÉTODOS NUMERICOS

Pág.16

Se observa claramente que esta estimación, que se obtuvo sin iteración sobre la ecuación corrector, tiene un error relativo porcentual del –8.18 % respecto de la solución exacta. Entonces esta es mucho mas aproximada a la solución exacta que la que se obtuvo en el primer paso, que es equivalente a la aplicación simple de Euler. Si refinamos la predicción de y1 al sustituir este valor calculado en el segundo miembro de la ecuación corrector tendremos:

y ′ + yˆ ′ yˆ1 = y0 + 0 1 ⋅ h 2 Ahora calculamos nuevamente la estimación de la pendiente en el punto final del intervalo: yˆ1′ = 4 ⋅ e0.8 ⋅ x1 - 0.5 ⋅ y1 = 4 ⋅ e0.8 ⋅1 - 0.5 ⋅ 6.701082 = 5.551623

Y por último con esta nueva estimación de la pendiente obtenemos una nueva estimación del valor de la función al final del intervalo: yˆ1 = 2 +

3 + 5.551623 ⋅1 = 6.275811 2

Esta nueva estimación tiene un error relativo porcentual es de –1.31 %. Debe tenerse en cuenta de que puede suceder que al efectuar un número mayor de iteraciones el error aumente. Si esto sucediera significaría que el tamaño del paso utilizado no es el adecuado y tendríamos que disminuirlo. Para el caso particular que estamos analizando si iteráramos nuevamente obtendríamos un valor de y1 = 6.382129, lo que significaría un error relativo porcentual de 3.03%. esto nos indica en principio que el paso h = 1 que estamos utilizando es excesivo. Los resultados se sintetizan en la siguiente tabla: X

Y verdadero

1 iteración

15 iteraciones

YHeun

 r

YHeun

r %

0

2.0000000

2.0000000

0.00

2.0000000

0.00

1

6.1946314

6.7010819

8.18

6.3608655

2.68

2

14.8439219

16.3197819

9.94

15.3022367

3.09

3

33.6771718

37.1992489

10.46

34.7432761

3.17

4

75.3389626

83.3377674

10.62

77.7350962

3.18

CAPITULO – SOLUCIÓN NUMÉRICA DE ECUACIONES DIFERENCIALES

CÁTEDRA MÉTODOS NUMERICOS

Pág.17

En el siguiente gráfico podemos ver una comparación, para la ecuación polinómica del ejemplo anterior, entre los métodos de Heun, Euler, calculado con h = 0.5, y la solución verdadera:

2.3.1.2.1 Análisis del error para el método de Heun En las ecuaciones predictor y corrector que vimos precedentemente la derivada es una función tanto de la variable independiente x, como de la variable dependiente y. En casos tales como polinomios, donde la EDO es solo función de la variable independiente la ecuación predictor no es requerida y el corrector se aplica solo una vez en cada paso. Para estos casos la técnica se aplica en forma concisa como:

y i +1 = yi +

f ( xi ) + f ( xi +1 ) ⋅h 2

Por otro lado, si recordamos la expansión de Taylor para una función f(x) alrededor del punto xi, tendremos: f '' (xi ) 2 f '' ' (x i ) 3 f (n) (xi ) n y i +1 = y i + f (xi ) ⋅ h + ⋅h + ⋅ h + .. + ⋅ h + O(h n +1 ) 2! 3! n! '

La derivada segunda de la función f(x) en el punto xi puede expresarse en forma aproximada como:

f ′′(x i ) =

f ′(x i +1 ) − f ′(x i ) h

CAPITULO – SOLUCIÓN NUMÉRICA DE ECUACIONES DIFERENCIALES

CÁTEDRA MÉTODOS NUMERICOS

Pág.18

Si reemplazamos esta en la expansión de Taylor y tomamos el desarrollo hasta el término de orden 3 tendremos que:

f ′(xi +1 ) − f ′(xi ) f ′′′(xi ) 3 f ′(xi +1 ) − f ′(xi ) h y i +1 = y i + f ′(xi ) ⋅ h + ⋅ h2 + ⋅ h = y i + f ′(xi ) ⋅ h + ⋅ h + O(h3 ) 2! 3! 2 Reordenando la ecuación queda:

y i +1 = y i +

f ′(xi +1 ) + f ′(xi ) ⋅ h + O(h3 ) 2

Vemos que los dos primeros términos del segundo miembro se corresponden con la expresión de la ecuación corrector del método de Heun. Por lo tanto, el error de truncamiento local del método de Heun es de orden 3 y esta dado por la expresión: Ea = −

f ′′′( ξ ) 3 ⋅h 12

Queda demostrado de esta manera que el método de Heun es de segundo orden porque se incluyen términos de segundo orden del desarrollo de Taylor. Nótese que si la función solución es un polinomio cuadrático o inferior el método de Heun es exacto puesto que la tercer derivada de la función solución es cero. Puede demostrarse que el error global es de orden O(h2). Esto significa que al disminuir el paso h disminuye el error a una velocidad mayor que para el método de Euler.

Método del punto medio ( o del polígono mejorado) Esta técnica es otra simple modificación del método de Euler. El método consiste en usar el método de Euler para predecir un valor de y en el punto medio del intervalo: y i + 1 / 2 = y i + f ( xi , y i ) ⋅

h 2

Este valor predicho se usa para calcular una pendiente en el punto medio del intervalo:

y i′ + 1 / 2 = f ( xi + 1 / 2 , yi + 1 / 2 )

CAPITULO – SOLUCIÓN NUMÉRICA DE ECUACIONES DIFERENCIALES

CÁTEDRA MÉTODOS NUMERICOS

Pág.19

Esta pendiente se considera como una pendiente promedio en todo el intervalo. Luego con esta pendiente promedio estimada se extrapola linealmente desde xi hasta xi+1: y i +1 = y i + f ( xi + 1 / 2 , y i + 1 / 2 ) ⋅ h

Esta técnica no es iterativa puesto que yi+1 no esta en ambos miembros de la ecuación corrector.

El error en el método del punto medio En un análisis similar que el realizado para el método de Heun se demuestra que el método del polígono mejorado también es un método de segundo orden, que tampoco requiere de la evaluación de derivadas superiores al primer orden. Tambien en este método los errores de truncamiento local y global son de orden O(h3) y O(h2) respectivamente.

Métodos de Runge – Kutta El método de Heun, el del punto medio y la misma técnica de Euler son casos particulares de una clase mas general de procedimientos de un paso denominados métodos de Runge – Kutta. Esta afirmación será demostrada en los desarrollos subsiguientes. CAPITULO – SOLUCIÓN NUMÉRICA DE ECUACIONES DIFERENCIALES

CÁTEDRA MÉTODOS NUMERICOS

Pág.20

Los métodos de Runge – Kutta (RK) tienen la característica de poseer precisiones propias de desarrollos de Taylor que incluyen términos de derivadas de ordenes superiores a uno sin requerir el cálculo de las mismas. Si escribimos la ecuación original del método de Euler en forma generalizada, tendremos:

y i +1 = y i + φ ( xi , y i , h ) ⋅ h donde a φ ( xi , y i , h ) se la denomina función incremento, y puede ser interpretada como una pendiente representativa sobre el intervalo. La función incremento se escribe en general como: φ = a1 ⋅ k1 + a2 ⋅ k2 + .......... .......... + an ⋅ k n

donde las a son constantes y las k son: k1 = f ( xi , y i )

k2 = f ( xi + p1 ⋅ h, y i + q11 ⋅ k1 ⋅ h )

k3 = f ( xi + p2 ⋅ h, y i + q21 ⋅ k1 ⋅ h + q22 ⋅ k2 ⋅ h ) . kn = f ( xi + pn-1 ⋅ h, y i + qn-1,1 ⋅ k1 ⋅ h + qn-1,2 ⋅ k2 ⋅ h + .......... .. + qn-1,n-1 ⋅ kn-1 ⋅ h )

Se observa claramente que las k son relaciones de recurrencia ya que la primera interviene en la segunda, a su vez la primera y la segunda intervienen en la tercera y así sucesivamente. Es posible concebir varios tipos de métodos de Runge – Kutta al emplear diferentes números de términos en la función incremento, que en general tiene n. Si tomamos n = 1, método de Runge – Kutta de primer orden, tendremos que la función incremento es:

φ = a1 ⋅ k1 = a1 ⋅ f ( xi , y i ) Como se observa el método de Runge – Kutta de primer orden es el método de Euler. Una vez que se elige n, se evalúan las a, p y q al igualar la ecuación inicial a los términos correspondientes de la serie de Taylor. Así, al menos para las versiones de orden inferior, el número de términos n representa el orden de la aproximación. Por ejemplo, los métodos de Runge – Kutta de segundo orden, n = 2, para los cuales la función incremento consta de dos términos, serán exactos si la solución de la ecuación diferencial es una cuadrática. Además, como los términos que tienen h3 y mayores no son considerados en el desarrollo el error de truncamiento local es O(h3) y el global es O(h2).

CAPITULO – SOLUCIÓN NUMÉRICA DE ECUACIONES DIFERENCIALES

CÁTEDRA MÉTODOS NUMERICOS

Pág.21

2.3.1.4.1 Métodos de Runge – Kutta de segundo orden. Si la función incremento tiene dos términos, esto es n = 2 entonces la ecuación general del método de Runge – Kutta se escribe como:

y i +1 = y i + ( a1 ⋅ k1 + a2 ⋅ k2 ) ⋅ h Las ai son constantes a determinar y las ki son: k1 = f ( xi , y i )

k2 = f ( xi + p1 ⋅ h, y i + q11 ⋅ k1 ⋅ h )

Las constantes p1 y q11 también deben ser determinadas. En definitiva deberemos determinar los valores de las constantes a1, a2, p1 y q11. Para ello desarrollamos la serie de Taylor de segundo orden para y i+1 en términos de yi y la derivada primera de la función respecto de x valuada en xi, yi , f ( xi , y i ) :

y i +1 = y i + f ( xi , y i ) ⋅ h +

f ′( x i , y i ) 2 ⋅h 2!

Si tenemos en cuenta que : f ′( xi , y i ) =

∂ f ( x, y ) ∂ f ( x, y ) dy + ⋅ ∂x ∂y dx

Sustituyendo esta en la anterior, tendremos que:  ∂ f ( x, y ) ∂ f ( x, y ) dy y i +1 = y i + f ( xi , y i ) ⋅ h +  + ⋅ ∂y dx  ∂x

 h2  ⋅  2!

Por otro lado desarrollamos por serie de Taylor a la función k2, recordando que una serie de Taylor para dos variables se define como: g ( x +r,y +r ) = g ( x, y ) + r ⋅

∂g ∂g + s⋅ + .......... ........ ∂x ∂y

Entonces tendremos que: k2 = f ( xi + p1 ⋅ h, y i + q11 ⋅ k1 ⋅ h ) = f ( xi , y i ) + p1 ⋅ h⋅

∂f ∂f + q11 ⋅ k1 ⋅ h⋅ + ......... ∂x ∂y

CAPITULO – SOLUCIÓN NUMÉRICA DE ECUACIONES DIFERENCIALES

CÁTEDRA MÉTODOS NUMERICOS

Pág.22

Si ahora reemplazamos k1 y k2 desarrollado como serie en la ecuación inicial nos queda: y i +1 = y i + a1 ⋅ h⋅ f ( xi , y i ) + a2 ⋅ h⋅ f ( xi , y i ) + a2 ⋅ p1 ⋅ h2 ⋅

∂f ∂f + a2 ⋅ f ( xi , y i ) ⋅ q11 ⋅ h2 ⋅ + .... ∂x ∂y

Si agrupamos términos:  ∂f ∂f  2 yi+1 = yi + [a1 ⋅ f ( xi, yi ) + a2 ⋅ f ( xi, yi ) ] ⋅ h + a2 ⋅ p1 ⋅ + a2 ⋅ f ( xi, yi ) ⋅ q11 ⋅  ⋅ h + .... ∂x ∂y  

Si comparamos términos de esta última ecuación con:  ∂ f ( x, y ) ∂ f ( x, y ) dy y i +1 = y i + f ( xi , y i ) ⋅ h +  + ⋅ ∂y dx  ∂x

 h2  ⋅  2!

para hacer equivalentes las dos ecuaciones se debe cumplir que: a1 + a2 = 1 , a2 ⋅ p1 =

1 2

y

a2 ⋅ q11 =

1 2

Se observa que estas tres ecuaciones simultáneas contienen las cuatro constantes desconocidas. Al haber una incógnita mas que el número de ecuaciones, no existe un conjunto único de constantes que satisfagan las ecuaciones. Sin embargo, podemos suponer una de las constantes y calcular a las otras tres. En consecuencia, existe una familia de métodos de segundo orden. Si hacemos que todas las constantes queden, por ejemplo, en función de a2: a1 = 1 - a2 , p1 =

1 2 ⋅ a2

y

q11 =

1 2 ⋅ a2

Debido a que podemos elegir infinitos valores para a2, entonces tenemos infinitos métodos de Runge – Kutta de segundo orden. Cada una de estas versiones darían el mismo resultado si la solución de la ecuación diferencial ordinaria fuera cuadrática, lineal o constante, y diferentes resultados si la solución es mas complicada. Método de Heun con un solo corrector : Para a2 = ½ tendremos que, por reemplazo en las ecuaciones anteriores, a1 = ½ , p1 = 1 y q11 = 1. Reemplazando estos valores en la ecuación general de segundo orden del método de Runge – Kutta:

CAPITULO – SOLUCIÓN NUMÉRICA DE ECUACIONES DIFERENCIALES

CÁTEDRA MÉTODOS NUMERICOS

Pág.23

y i +1 = y i + ( a1 ⋅ k1 + a2 ⋅ k2 ) ⋅ h Obtendremos la ecuación: 1 1  y i +1 = y i +  ⋅ k1 + ⋅ k2  ⋅ h 2 2 

Ecuación esta donde: k1 = f ( xi , y i )

k2 = f ( xi + h, y i + k1 ⋅ h )

Se observa que k1 es la pendiente en el inicio del intervalo y k2 la pendiente al final del mismo. Queda claro entonces que se reproduce de esta manera el método de Heun sin iteración. Método del punto medio: Para a2 = 1, tendremos entonces, por reemplazo en las ecuaciones correspondientes, que a1 = 0, p1 = q11 = ½. Si reemplazamos a estas constantes en la ecuación general de segundo orden:

y i +1 = y i + ( a1 ⋅ k1 + a2 ⋅ k2 ) ⋅ h Quedará la ecuación: y i +1 = y i + k2 ⋅ h

Ecuación esta donde: k1 = f ( xi , y i ) 1 1   k2 = f  xi + ⋅ h, y i + ⋅ k1 ⋅ h  2 2  

Se observa que k1 es la pendiente en el inicio del intervalo y k2 la pendiente en el punto medio del mismo. Entonces, con estos valores de constantes se reproduce el método del punto medio. Método de Ralston: Para a2 = 2/3, tendremos que, por reemplazo en las ecuaciones anteriores, a1 = 1/3 , p1 = q11 = ¾. Si reemplazamos estos valores en la ecuación general de segundo orden :

y i +1 = yi + ( a1 ⋅ k1 + a2 ⋅ k2 ) ⋅ h CAPITULO – SOLUCIÓN NUMÉRICA DE ECUACIONES DIFERENCIALES

CÁTEDRA MÉTODOS NUMERICOS

Pág.24

Quedara la ecuación: 2 1  y i +1 = y i +  ⋅ k1 + ⋅ k2  ⋅ h 3 3 

Ecuación donde: k1 = f ( xi , y i ) 3 3   k2 = f  xi + ⋅ h, y i + ⋅ k1 ⋅ h  4 4  

Se observa que k1 es la pendiente en el inicio del intervalo y k2 la pendiente en el punto ubicado a 3/4 del mismo. Este es el denominado método de Ralston.

Ejemplo: Usar el método de punto medio y el método de Ralston para resolver numéricamente la ecuación: dy = - 2 ⋅ x 3 + 12 ⋅ x 2 - 20 ⋅ x + 8.5 dx

desde x = 0 hasta x = 4 con un paso de h = 0.5. La condición inicial en x = 0 es y = 1.

a) Método del punto medio: k1 = f ( xi , y i ) 1 1   k2 = f  xi + ⋅ h, y i + ⋅ k1 ⋅ h  2 2  

reemplazando, y teniendo en cuenta que la ecuación solo depende de x: k1 = f ( xi , y i ) = −2 ⋅ ( 0 ) 3 + 12 ⋅ ( 0 ) 2 − 20 ⋅ ( 0 ) − 8.5 = 8.5 k2 = −2 ⋅ ( 0.25 ) 3 + 12 ⋅ ( 0.25 ) 2 − 20 ⋅ ( 0.25 ) + 8.5 = 4.21875

Luego, reemplazando en: y i +1 = y i + k2 ⋅ h

CAPITULO – SOLUCIÓN NUMÉRICA DE ECUACIONES DIFERENCIALES

CÁTEDRA MÉTODOS NUMERICOS

Pág.25

Obtendremos el valor en x = 0.5: y ( 0.5 ) =1 + 4.21875 ⋅ 0.25 = 3.109375

donde el error relativo porcentual respecto de la solucion exacta es de 3.4%. El cálculo se repite para los otros puntos aplicando el mismo criterio.

b) Método de Ralston: k1 = f ( xi , y i ) 3 3   k2 = f  xi + ⋅ h, y i + ⋅ k1 ⋅ h  4 4  

Reemplazando en las anteriores, y teniendo en cuenta que la ecuación solo depende de x, tendremos: k1 = f ( xi , y i ) = −2 ⋅ ( 0 ) 3 + 12 ⋅ ( 0 ) 2 − 20 ⋅ ( 0 ) − 8.5 = 8.5 k2 = −2 ⋅ ( 0.375 ) 3 + 12 ⋅ ( 0.375 ) 2 − 20 ⋅ ( 0.375 ) + 8.5 = 2.58203125

Reemplazando estos valores en la ecuación del método:

2 1  y i +1 = y i +  ⋅ k1 + ⋅ k2  ⋅ h 3 3 

Y nos queda para x = 0.5 :

2 1  y(0.5) = 1 +  ⋅ 8.5 + ⋅ 2.58203125  ⋅ 0.5 = 1 + 4.5546875 ⋅ 0.5 = 3.27734375 3 3  

El error relativo porcentual respecto de la solución exacta es de es de –1.82%. Luego, el cálculo se repite para los otros puntos aplicando el mismo procedimiento. Heun (sin iteración)

Punto Medio

Ralston

Y

  r %

Y

r %

Y

r %

1.00000

1.00000

0.0

1.000000

0.0

1.000000

0.0

3.21875

3.43750

6.8

3.109375

3.4

3.277344

1.8

X

Y verdadero

0.0 0.5

CAPITULO – SOLUCIÓN NUMÉRICA DE ECUACIONES DIFERENCIALES

CÁTEDRA MÉTODOS NUMERICOS

Pág.26

1.0

3.00000

3.37500

12.5

2.812500

6.3

3.101563

3.4

1.5

2.21875

2.68750

21.1

1.984375

10.6

2.347656

5.8

2.0

2.00000

2.50000

25.0

1.750000

12.5

2.140625

7.0

2.5

2.71875

3.18750

17.2

2.484375

8.6

2.855469

5.0

3.0

4.00000

4.37500

9.4

3.812500

4.7

4.117188

2.9

3.5

4.71875

4.93750

4.6

4.609375

2.3

4.800781

1.7

4.0

3.00000

3.00000

0.0

3.000000

0.0

3.031250

1.0

En la tabla anterior se presentan los valores de la función solución de la ecuación diferencial propuesta calculada según los distintos métodos. En el gráfico inferior se graficaron dichas soluciones para su comparación.

2.3.1.4.2 Métodos de Runge – Kutta de tercer orden. Partiendo de la ecuación general de Runge – Kutta y haciendo n = 3 nos queda: y i +1 = y i + ( a1 ⋅ k1 + a2 ⋅ k2 + a3 ⋅ k 3 ) ⋅ h Se puede hacer un hacer un desarrollo similar al del método de segundo orden. Como

resultado de dicho desarrollo se llegan a 6 ecuaciones con 8 incógnitas por lo que deben especificarse con antelación los valores de 2 de ellas con el fin de establecer todos los parámetros restantes. Una versión común del método de Runge – Kutta de tercer orden que resulta es: y i +1 = y i +

1 ⋅ ( k1 + 4 ⋅ k2 + k3 ) ⋅ h 6

CAPITULO – SOLUCIÓN NUMÉRICA DE ECUACIONES DIFERENCIALES

CÁTEDRA MÉTODOS NUMERICOS

Pág.27

donde: k1 = f ( xi , y i ) 1 1   k2 = f  xi + ⋅ h, y i + ⋅ k1 ⋅ h  2 2   k3 = f ( xi + h, y i − k1 ⋅ h + 2 ⋅ k2 ⋅ h )

Puede observarse que si la derivada de la función solución depende solo de x este método de tercer orden se reduce a la regla de Simpson 1/3. Los métodos de Runge – Kutta de tercer orden tienen errores local y global de O(h4) y O(h3) respectivamente y dan resultados exactos cuando la función es una función cúbica o de menor orden. Si se trata de polinomios la ecuación anterior dará también resultados exactos cuando la función solución de la ecuación diferencial es de cuarto orden debido a que la regla de Simpson 1/3 proporciona estimaciones exactas de la integral de cúbicas.

2.3.1.4.3 Métodos de Runge – Kutta de cuarto orden. De las infinitas versiones de los métodos de Runge – Kutta los de cuarto orden son los más utilizados. Partiendo de la ecuación general y haciendo n = 4 tenemos: y i +1 = y i + ( a1 ⋅ k1 + a2 ⋅ k2 + a3 ⋅ k3 + a4 ⋅ k4 ) ⋅ h

Se puede hacer un hacer un desarrollo similar al del método de segundo orden. Como resultado de dicho desarrollo se llega a un número de ecuaciones inferior a la cantidad de incógnitas por lo que deben especificarse con antelación los valores de algunas de ellas con el fin de establecer todos los parámetros restantes. La forma de uso mas común, de todas las infinitas posibilidades, es la que se denomina método de Runge – Kutta clásico de cuarto orden. La expresión resultante es: y i +1 = y i +

1 ⋅ ( k1 + 2 ⋅ k2 + 2 ⋅ k3 + k4 ) ⋅ h 6

donde:

CAPITULO – SOLUCIÓN NUMÉRICA DE ECUACIONES DIFERENCIALES

CÁTEDRA MÉTODOS NUMERICOS

Pág.28

k1 = f ( xi , y i ) 1 1   k2 = f  xi + ⋅ h, y i + ⋅ k1 ⋅ h  2 2   1 1   k3 = f  xi + ⋅ h, y i + ⋅ k2 ⋅ h  2 2   k4 = f ( xi + h, y i + k3 ⋅ h )

Puede observarse que si la derivada de la función solución depende solo de x el método de Runge – Kutta clásico de cuarto orden es similar a la regla de Simpson 1/3. También presenta alguna similitud con el método de Heun, en el sentido que son desarrolladas estimaciones múltiples de las pendientes en el punto medio, para finalmente, combinadas con las pendientes obtenidas al inicio y final del intervalo, obtener una pendiente promedio mejorada para el intervalo. En esta versión del método, como en las anteriores de distintos ordenes, cada una de las ki representa una pendiente. Luego, reemplazadas estas en la primera expresión, se obtiene una pendiente media mejorada representativa del intervalo. La interpretación gráfica de las pendientes estimadas ki se presenta a continuación:

Ejemplo: Resolver, utilizando el método de Runge – Kutta de cuarto orden, las siguientes ecuaciones: a)

dy = f ( x, y ) = -2 ⋅ x 3 + 12 ⋅ x 2 - 20 ⋅ x + 8.5 dx

mediante un tamaño de paso de h = 0.5 y una condición inicial de y = 1 en x = 0. En principio debemos calcular las pendientes ki, según las ecuaciones dadas para las mismas. Obsérvese que la derivada solo depende de x:

CAPITULO – SOLUCIÓN NUMÉRICA DE ECUACIONES DIFERENCIALES

CÁTEDRA MÉTODOS NUMERICOS

Pág.29

k1 = f ( xi , y i ) = - 2 ⋅ ( xi ) 3 + 12 ⋅ ( xi ) 2 - 20 ⋅ ( xi ) + 8.5 3

2

3

2

1 1 1  1  1       k2 = f  xi + ⋅ h, y i + ⋅ k1 ⋅ h  = - 2 ⋅  xi + ⋅ h  + 12 ⋅  xi + ⋅ h  - 20 ⋅  xi + ⋅ h  + 8.5 2 2 2  2  2       1 1 1  1  1       k3 = f  xi + ⋅ h, y i + ⋅ k2 ⋅ h  = -2 ⋅  xi + ⋅ h  + 12 ⋅  xi + ⋅ h  - 20 ⋅  xi + ⋅ h  + 8.5 2 2 2  2  2       k4 = f ( xi + h, y i + k3 ⋅ h ) = -2 ⋅ ( xi + h ) 3 + 12 ⋅ ( xi + h ) 2 - 20 ⋅ ( xi + h ) + 8.5

Entonces tenemos que xi = 0, xi+1/2 = 0.25 y xi+h = 0.5. Reemplazando en las anteriores: k1 = -2 ⋅ ( 0 ) 3 + 12 ⋅ ( 0 ) 2 - 20 ⋅ ( 0 ) + 8.5 = 8.5 k2 = -2 ⋅ ( 0.25 ) 3 + 12 ⋅ ( 0.25 ) 2 - 20 ⋅ ( 0.25 ) + 8.5 = 4.21875 k3 = -2 ⋅ ( 0.25 ) 3 + 12 ⋅ ( 0.25 ) 2 - 20 ⋅ ( 0.25 ) + 8.5 = 4.21875 k4 = -2 ⋅ ( 0.5 ) 3 + 12 ⋅ ( 0.5 ) 2 - 20 ⋅ ( 0.5 ) + 8.5 = 1.25

Luego estas pendientes son reemplazadas en la expresión del método de Runge – Kutta clásico de cuarto orden: y i +1 = y i +

1 ⋅ ( k1 + 2 ⋅ k2 + 2 ⋅ k3 + k4 ) ⋅ h 6

Aplicada esta ecuación, para obtener el valor de la función en 0.5, y(0.5), partiendo del valor conocido de la función en 0, y(0): y (0.5) = 1 +

1 ⋅ ( 8.5 + 2 ⋅ 4.21875 + 2 ⋅ 4.21875 + 1.25 ) ⋅ 0.5 = 3.21875 6

Esta solución coincide con la exacta. Esto se debe a que la solución verdadera es de cuarto orden porque estamos integrando un polinomio de tercer orden. Entonces este método, que para polinomios reproduce la solución de Simpson 1/3 nos proporciona la solución exacta.

b)

dy = f ( x, y ) = 4 ⋅ e0.8 ⋅x − 0.5 ⋅ y dx

mediante un tamaño de paso de h = 0.5 y una condición inicial de y = 2 en x = 0. De la misma manera procedemos al cálculo de las pendientes partiendo de que en xi = 0, yi = 2:

CAPITULO – SOLUCIÓN NUMÉRICA DE ECUACIONES DIFERENCIALES

CÁTEDRA MÉTODOS NUMERICOS

Pág.30

k1 = f ( 0,2 ) = 4 ⋅ e0.8 ⋅( 0 ) − 0.5 ⋅ ( 2 ) = 3

Este valor de k1 se usa para calcular un valor de y y la pendiente k2 en el punto medio: y ( 0.25 ) = y ( 0 ) + k1 ⋅

h  0.5  = 2 + 3.  = 2.75 2  2 

k2 = f ( 0.25,2.75 ) = 4 ⋅ e0.8 ⋅( 0.25 ) − 0.5 ⋅ ( 2.75 ) = 3.510611

Esta pendiente k2 se usa a su vez para calcular otro valor de y y otra pendiente, k3, en el punto medio del intervalo: y ( 0.25 ) = y ( 0 ) + k2 ⋅

h  0.5  = 2 + 3.510611.   = 2.877653 2  2 

k3 = f ( 0.25,2.877 653 ) = 4 ⋅ e0.8 ⋅( 0.25 ) − 0.5 ⋅ ( 2.877653 ) = 3.446785

Después, esta pendiente k3 se usa a su vez para calcular el valor de y y la pendiente k4 en el punto final del intervalo: y ( 0.5 ) = y ( 0 ) + k3 ⋅ h = 2 + 3.446785. ( 0.5 ) = 3.723392

k4 = f ( 0.5,3.7233 92 ) = 4 ⋅ e0.8 ⋅( 0.5 ) − 0.5 ⋅ ( 3.723392 ) = 4.105603

Por ultimo, las cuatro estimaciones de la pendiente se combinan para obtener una pendiente promedio, que a su vez es utilizada para realizar la predicción al final del intervalo: y(0.5) = 2 +

1 ⋅ ( 3 + 2 ⋅ 3.510611 + 2 ⋅ 3.446785 + 4.105603 ) ⋅ 0.5 = 3.751669 6

La solución verdadera es y(0.5) = 3.751521.

2.3.1.4.4 Métodos de Runge – Kutta de orden superior. Si se requiere mayor exactitud en las estimaciones es recomendable utilizar alguno de los métodos de Runge – Kutta de quinto orden. Entre estos se destaca el método de Butcher: y i +1 = y i +

1 ⋅ (7 ⋅ k1 + 32 ⋅ k3 + 12 ⋅ k4 + 32 ⋅ k5 + 7 ⋅ k6 ) ⋅ h 90

donde: CAPITULO – SOLUCIÓN NUMÉRICA DE ECUACIONES DIFERENCIALES

CÁTEDRA MÉTODOS NUMERICOS

Pág.31

k1 = f ( xi , y i ) 1 1   k2 = f  xi + ⋅ h, y i + ⋅ k1 ⋅ h  4 4   1 1 1   k3 = f  xi + ⋅ h, y i + ⋅ k1 ⋅ h + ⋅ k2 ⋅ h  4 8 8   1 1   k4 = f  xi + ⋅ h, y i − ⋅ k2 ⋅ h + k3 ⋅ h  2 2   3 3 3   k5 = f  xi + ⋅ h, y i + ⋅ k1 ⋅ h + ⋅ k4 ⋅ h  4 16 16   3 2 12 12 8   k6 = f  xi + h, y i − ⋅ k1 ⋅ h + ⋅ k2 ⋅ h + ⋅ k3 ⋅ h − ⋅ k 4 ⋅ h + ⋅ k5 ⋅ h  7 7 7 7 7  

Este método es mas preciso que el de cuarto orden, pero es mucho mayor la complejidad adicional y el esfuerzo computacional.

2.3.1.4.5 Comparación de los métodos de Runge – Kutta La comparación se plantea en términos de resolver la ecuación diferencial: dy = f ( x, y ) = 4 ⋅ e0.8 ⋅ x − 0.5 ⋅ y dx

mediante distintos tamaños de paso. La condición inicial es que la función solución vale : y = 2 en x = 0 y el intervalo de solución es el [0,4]. Se compararon las exactitudes de los distintos métodos al calcular el resultado de la función solución en x = 4. La respuesta exacta para este punto del dominio es: y(4) = 75.33896 Se realizó el calculo mediante los métodos de Euler, Heun sin iteración, Runge – Kutta de tercer orden, Runge – Kutta clásico de cuarto orden y el Runge – Kutta de Butcher de quinto orden. La comparación se realizo calculando el error relativo porcentual de cada método para distintos tamaños de paso. El esfuerzo computacional involucrado en la obtención de cada solución se define como: Esfuerzo = nf ⋅

b- a h

CAPITULO – SOLUCIÓN NUMÉRICA DE ECUACIONES DIFERENCIALES

CÁTEDRA MÉTODOS NUMERICOS

Pág.32

donde nf es el número de evaluaciones de la función involucradas en el cálculo para un determinado método. Para nf menores o iguales que 4, nf es igual al orden del método. Para ordenes superiores nf es mayor que el orden del método (el método de Butcher es de orden 5 y requiere de 6 evaluaciones). Por otro lado, la cantidad ( b – a ) / h (el intervalo dividido por el tamaño del paso) es la cantidad de aplicaciones del método para obtener el resultado. Si consideramos que las evaluaciones de la función son con frecuencia los pasos que consumen la mayor cantidad de tiempo, la ecuación anterior proporciona una cierta medida del tiempo de ejecución requerido para alcanzar la respuesta. Por ultimo se graficó el valor absoluto del error relativo porcentual contra el esfuerzo computacional para cada uno de los métodos utilizados en la comparación.

La simple inspección del l gráfico anterior nos permite concluir que: •

Los métodos de orden superior alcanzan mayor exactitud para el mismo esfuerzo computacional. • La ganancia en exactitud con el esfuerzo adicional (disminuyendo h) tiende a disminuir después de un punto. Las curvas primero caen con rapidez y luego tienden a nivelarse.

CAPITULO – SOLUCIÓN NUMÉRICA DE ECUACIONES DIFERENCIALES

CÁTEDRA MÉTODOS NUMERICOS

Pág.33

Resolución de Ecuaciones Diferenciales Ordinarias de orden n mediante las técnicas de Runge - Kutta. Transformación de una EDO de orden n en un sistema de n EDOs de primer orden Considérese una ecuación diferencial de la siguiente forma:

 dy d 2 y d n −1y = f  x, y, , ,......,  dx dx 2 dx n dx n −1 

d ny

   

para x [a,b]

Esta ecuación diferencial ordinaria involucra a la función y y a sus n primeras derivadas, puesto que la derivada n-ésima depende, según una función conocida f, de x, y y las n – 1 primeras derivadas. Para que la ecuación anterior tenga una solución única son necesarias n condiciones adicionales sobre la función incógnita y. Como sabemos, estas condiciones adicionales se llaman condiciones iniciales si están dadas en un mismo punto del intervalo [a,b], o condiciones de contorno si están dadas en más de un punto del intervalo [a,b]. Un caso habitual de condiciones iniciales es que la función y y sus n - 1 primeras derivadas tengan valores prescritos conocidos  0, 1, 2,.........,n-1 en el extremo a del intervalo:

y ( a ) = α0 ;

2 dy ( a ) = α1 ; d 2y ( a ) = α2 dx dx

;......... ..;

d n-1y dx n-1

( a ) = αn −1

Entonces, si se complementa la ecuación diferencial ordinaria con las condiciones iniciales anteriores, se obtiene un problema de valor inicial. Partiendo de la información que se tiene de la función y en el punto x = a debemos integrar la ecuación diferencial ordinaria para hallar la evolución de la función y en todo el intervalo [a,b]. Por otro lado sabemos que una ecuación diferencial ordinaria de orden n puede ser transformada en un sistema de n ecuaciones diferenciales ordinarias de primer orden con n funciones incógnitas. Es decir que podemos reducir el orden de las derivadas a costa de aumentar el número de incógnitas

CAPITULO – SOLUCIÓN NUMÉRICA DE ECUACIONES DIFERENCIALES

CÁTEDRA MÉTODOS NUMERICOS

Pág.34

En nuestro caso esta transformación es necesaria puesto que las técnicas numéricas que hemos visto hasta este momento están diseñadas para resolver problemas de primer orden. La idea básica de la transformación es tratar explícitamente como funciones incógnita a las n – 1 primeras derivadas de la función y. Esto puede expresarse como:

y1 ≡ y y2 ≡ y3 ≡

dy dy1 = dx dx d2y dx

2

=

dy 2 dx

=

dy n −1 dx

(1)

. yn ≡

dny dx

n

Con la ayuda de las ecuaciones anteriores la ecuación diferencial ordinaria original puede escribirse como:

dy n = f ( x, y1 , y 2 , y 3 ,......, y n ) dx

para x  [a,b]

Queda claro que, en definitiva, solo hemos hecho un cambio de notación. Si tomamos esta última ecuación y la combinamos con las (1) a excepción de la primera y se transforman también las condiciones iniciales se obtiene:

CAPITULO – SOLUCIÓN NUMÉRICA DE ECUACIONES DIFERENCIALES

CÁTEDRA MÉTODOS NUMERICOS

Pág.35

    dy 2 d 2 y y3 = ≡  dx dx 2   :  dy n −1 d n y  yn = ≡  n dx dx   dy n = f ( x, y1 , y 2 , y 3 ,......, y n )  dx  y2 =

dy1 dy ≡ dx dx

n ecuaciones con n incognitas

y1 ( a ) = y ( a ) = α0 dy y2 ( a) = ( a ) = α1 dx

     d 2 y1  y3 ( a) = ( a ) = α2  2 dx  .   n-1 d y . yn ( a) = ( a ) = αn −1  n-1 dx 

n condicione s iniciales

Ejemplo: Expresar la siguiente ecuación diferencial ordinaria de orden 3 en un sistema de tres ecuaciones diferenciales ordinarias de orden 1. Dada la siguiente ecuación: x⋅

d3y( x) dx

3

+2⋅

d2y( x) dx

2

+ 3 ⋅ x2 ⋅

dy ( x ) + 4 ⋅ y ( x ) = 4 ⋅ ex + 2 dx

Donde las condiciones iniciales son en x = 0, y(0) =4, y’(0) = -1, y’’(0) = 2. Transformarla en un sistema de 3 ecuaciones de primer orden. Para ello definimos las siguientes variables dependientes: y' ( x ) = v ( x ) v' ( x ) = u( x )

Entonces puede apreciarse que: y' ( x ) = v ( x ) y '' ( x ) = v ' ( x ) = u ( x ) y '' ' ( x ) = v '' ( x ) = u' ( x ) CAPITULO – SOLUCIÓN NUMÉRICA DE ECUACIONES DIFERENCIALES

CÁTEDRA MÉTODOS NUMERICOS

Pág.36

La ecuación original podrá entonces escribirse en términos del siguiente sistema: y' ( x ) = v ( x ) v' ( x ) = u( x ) u' ( x ) =

4 ⋅ e x + 2 − 2 ⋅ v' ( x ) − 3 ⋅ x 2 ⋅ y' ( x ) − 4 ⋅ y x

Y las condiciones iniciales en x = 0:

y (0) = 4;

v(0) = −1

y

u(0) = 2

Resolución de un sistema de n EDOs de primer orden Ya sea que estemos afrontando un problema de ingeniería cuya solución implique la resolución de una ecuación diferencial de orden n, o uno que implique la resolución de un sistema de ecuaciones diferenciales ordinarias de primer orden, nos enfrentaremos a la necesidad de resolver un sistema que podremos expresar como sigue: dy 1 = f1 ( x , y1 , y 2 ,......... ..., y n ) dx dy 2 = f2 ( x , y1 , y 2 ,......... ..., y n ) dx . . dy n = fn ( x , y1 , y 2 ,......... ..., y n ) dx

Por supuesto, la solución de tal sistema requiere de que se conozcan las n condiciones iniciales en el valor inicial del intervalo correspondiente a x. Todos los métodos vistos anteriormente para simples ecuaciones pueden extenderse a la resolución de sistemas como el anterior. El procedimiento para resolver un sistema de ecuaciones simplemente involucra aplicar las técnicas conocidas para cada ecuación en cada paso, antes de proceder con el siguiente. Esto quedará claramente ilustrado con el siguiente ejemplo donde hemos aplicado el método de Euler.

CAPITULO – SOLUCIÓN NUMÉRICA DE ECUACIONES DIFERENCIALES

CÁTEDRA MÉTODOS NUMERICOS

Pág.37

Ejemplo 1: Resolución de un sistema de EDOs mediante el método de Euler. Resolver el siguiente conjunto de ecuaciones diferenciales ordinarias: dy1 = - 0.5 ⋅ y1 dx dy 2 = 4 - 0.3 ⋅ y 2 - 0.1 ⋅ y1 dx

Resolveremos el sistema en el intervalo entre x = 0 y x = 2, con condiciones iniciales en x = 0, y1 = 4 , y2 = 6. Utilizaremos un paso h = 0.5. Se implementa el método de Euler para cada variable mediante la ya conocida expresión: y i +1 = y i + f ( xi , y i ) ⋅ h

Primero calculamos las pendientes: dy1 = f1 ( 0 , 4, 6) = - 0.5 ⋅ 4 = - 2 dx dy 2 = f2 ( 0 , 4, 6) = 4 - 0.3 ⋅ 6 - 0.1 ⋅ 4 = 1.8 dx

Y luego los valores de la función para el primer paso: y1 ( 0.5 ) = y1 ( 0 ) + f1 ( 0 , 4 , 6) ⋅ h = 4 + ( - 2 ) ⋅ 0.5 = 3 y 2 ( 0.5 ) = y 2 ( 0 ) + f2 ( 0 , 4 , 6) ⋅ h = 6 + (1.8 ) ⋅ 0.5 = 6.9

Para un segundo paso volvemos a calcular las pendientes: dy1 = f1 ( 0.5 , 3, 6.9) = - 0.5 ⋅ 3 = - 1.5 dx dy 2 = f2 ( 0.5 , 3, 6.9) = 4 - 0.3 ⋅ 6.9 - 0.1 ⋅ 3 = 1.63 dx

Y luego los valores de la función para el segundo paso:

CAPITULO – SOLUCIÓN NUMÉRICA DE ECUACIONES DIFERENCIALES

CÁTEDRA MÉTODOS NUMERICOS

Pág.38

y1 (1.0 ) = y1 ( 0.5 ) + f1 ( 0.5 , 3, 6.9) ⋅ h = 3 + ( - 1.5 ) ⋅ 0.5 = 2.25 y 2 (1.0 ) = y 2 ( 0.5 ) + f2 ( 0.5 , 3, 6.9) ⋅ h = 6.9 + (1.63 ) ⋅ 0.5 =7.715

Y así continúa el cálculo hasta el final. Los resultados se resumen en la siguiente tabla: x

y1

y2

0.0 4.000000 6.000000 0.5 3.000000 6.900000 1.0 2.250000 7.715000 1.5 1.687500 8.445250 2.0 1.265625 9.094870

Ejemplo 2: Resolución de un sistema de EDOs mediante el método de Runge – Kutta clásico de cuarto orden. Resolver el mismo conjunto de ecuaciones diferenciales ordinarias anterior: dy1 = - 0.5 ⋅ y1 dx dy 2 = 4 - 0.3 ⋅ y 2 - 0.1 ⋅ y1 dx

Resolveremos el sistema en el intervalo entre x = 0 y x = 2, con condiciones iniciales en x = 0, y1 = 4 , y2 = 6. Utilizaremos un paso h = 0.5. Si bien la aplicación del método de Runge – Kutta a un sistema de ecuaciones diferenciales ordinarias no presenta mayor complicación, desde el punto de vista conceptual, respecto a su aplicación a una sola ecuación, debe tenerse cuidado al determinarse las pendientes y al calcular los valores intermedios necesarios. Por ello antes de abordar la resolución planteada describiremos en detalle los pasos a seguir:

• Primero se calculan las pendientes para todas las variables en el valor inicial (las k1). • Esas pendientes se usarán entonces para hacer predicciones de las variable dependientes yi en el punto medio del intervalo.

CAPITULO – SOLUCIÓN NUMÉRICA DE ECUACIONES DIFERENCIALES

CÁTEDRA MÉTODOS NUMERICOS

Pág.39

• Dichos valores del punto medio se utilizan, a su vez, para calcular un conjunto de pendientes en el punto medio (las k2). • Esas nuevas pendientes se usarán entonces para hacer otro conjunto de predicciones de las variables dependientes yi en el punto medio, que a su vez se utilizarán para una nueva predicción de la pendiente en el punto medio (las k3). • Estas después se emplearán con el fin de hacer las predicciones de las variables dependientes yi al final del intervalo que serán usadas para calcular las pendientes del final del intervalo (las k4). • Por ultimo las ki se combinan para formar un conjunto de funciones incremento, que se utilizan para hacer la predicción final de las variables dependientes. Entonces, comenzando el procedimiento antes descripto primero calculamos las pendientes al inicio del intervalo: k1,1 = f1 ( 0 , 4, 6) = - 0.5 ⋅ 4 = - 2 k1,2 = f2 ( 0 , 4, 6) = 4 - 0.3 ⋅ 6 - 0.1 ⋅ 4 =1.8

Donde la nomenclatura utilizada debe interpretarse de la siguiente manera: ki,j es el iésimo valor de k para la j-ésima variable dependiente. Luego calculamos los primeros valores de y1 e y2 en el punto medio:

( 2 ) = y1 ( 0 ) + k1,1 ⋅ 2h = 4 + ( - 2 ) ⋅ 0.5 = 3.5 2

y1 h

( 2 ) = y 2 ( 0 ) + k1,2 ⋅ 2h = 6 + (1.8 ) ⋅ 0.5 = 6.45 2

y2 h

Estos valores se usarán para calcular el primer conjunto de pendientes de punto medio: k2,1 = f1 ( 0.25 , 3.5, 6.45) = - 0.5 ⋅ 3.5 = - 1.75 k2,2 = f2 ( 0.25 , 3.5, 6.45) = 4 - 0.3 ⋅ 6.45 - 0.1 ⋅ 3.5 =1.715

Estas pendientes, a su vez, se usan para determinar el segundo conjunto de predicciones de las variables dependientes en el punto medio:

CAPITULO – SOLUCIÓN NUMÉRICA DE ECUACIONES DIFERENCIALES

CÁTEDRA MÉTODOS NUMERICOS

Pág.40

( 2 ) = y1 ( 0 ) + k2,1⋅ 2h = 4 + ( - 1.75 ) ⋅ 0.5 = 3.5625 2

y1' h

( 2 ) = y2 ( 0 ) + k2,2 ⋅ 2h = 6 + (1.715 ) ⋅ 0.5 = 6.42875 2

y '2 h

Estos valores se utilizarán para calcular el segundo conjunto de pendientes de punto medio: k3,1 = f1 ( 0.25 , 3.5625, 6.42875) = - 0.5 ⋅ 3.5625 = - 1.78125 k3,2 = f2 ( 0.25 , 3.5625, 6.42875)

= 4 - 0.3 ⋅ 6.42875 - 0.1 ⋅ 3.5625 =1.715125

Estas pendientes, a su vez, se utilizarán para determinar las predicciones de las variables dependientes al final del intervalo: y1 ( h ) = y1 ( 0 ) + k3,1 ⋅ h = 4 + ( - 1.78125 ) ⋅ 0.5 = 3.109375 y 2 ( h ) = y 2 ( 0 ) + k3,2 ⋅ h = 6 + (1.715125 ) ⋅ 0.5 = 6.857563

Estas predicciones de las variables dependientes al final del intervalo serán utilizadas para calcular las pendientes al final del intervalo: k4,1 = f1 ( 0.5 , 3.109375, 6.857563)

= - 0.5 ⋅ 3.109375

= - 1.554688

k4,2 = f2 ( 0.5 , 3.109375, 6.857563)

= 4 - 0.3 ⋅ 6.857563 - 0.1 ⋅ 3.109375

= 1.631794

Los valores de k se pueden utilizar entonces para calcular las predicciones de las variables dependientes , mediante la ecuación: y i +1 = y i +

1 ⋅ ( k1 + 2 ⋅ k2 + 2 ⋅ k3 + k4 ) ⋅ h 6

Reemplazando los valores obtenidos calculamos las variables dependientes y1 e y2 respectivamente: y1 ( 0.5 ) = y1 ( 0 ) + y1 ( 0.5 ) = 4 +

1 ⋅ ( k1,1 + 2 ⋅ k2,1 + 2 ⋅ k3,1 + k4,1 ) ⋅ h 6

1 ⋅ ( - 2 + 2 ( - 1.75 ) + 2 ( - 1.78125 ) + ( −1.554688 6

) ) ⋅ 0.5

= 3.115234

y CAPITULO – SOLUCIÓN NUMÉRICA DE ECUACIONES DIFERENCIALES

CÁTEDRA MÉTODOS NUMERICOS

Pág.41

y 2 ( 0.5 ) = y 2 ( 0 ) + y 2 ( 0.5 ) = 6 +

1 ⋅ ( k1,2 + 2 ⋅ k2,2 + 2 ⋅ k3,2 + k4,2 ) ⋅ h 6

1 ⋅ (1.8 + 2 ⋅ (1.715 ) + 2 ⋅ (1.715125 ) + (1.631794 6

) ) ⋅ 0.5

= 6.857670

Procediendo de la misma manera para los puntos restantes se obtiene los resultados que se condensan en la siguiente tabla: x

y1

y2

0.0 4.000000 6.000000 0.5 3.115234 6.857670 1.0 2.426171 7.632106 1.5 1.889523 8.326886 2.0 1.471577 8.946865

Resolución de Ecuaciones Diferenciales Ordinarias de orden n mediante el método de diferencias finitas. La integración hacia delante paso a paso de ecuaciones diferenciales de orden superior puede también efectuarse sustituyendo en la ecuación diferencial y en sus condiciones iniciales, las derivadas por las expresiones ya conocidas de derivación numérica por diferencias finitas. Como ya sabemos el reemplazo debe hacerse siempre por expresiones de derivación consistentes, es decir, que tengan el mismo orden de interpolación, o dicho de otra manera el mismo orden de error. Los operadores para interpolaciones de distinto orden se obtienen según lo ya visto en cursos anteriores y son validas todas las conclusiones vistas entonces. Un breve resumen de este tema se desarrolla en los Apéndices 2.1 y 2.2 del presente capítulo. Ejemplo: Resolver, por el método de diferencias finitas, la siguiente ecuación diferencial ordinaria de segundo orden: d2y dx

2

−

dy −2 ⋅y = 0 dx

CAPITULO – SOLUCIÓN NUMÉRICA DE ECUACIONES DIFERENCIALES

CÁTEDRA MÉTODOS NUMERICOS

Pág.42

en el intervalo entre x = 0 y x = 0.5, con paso h = 0.1, y las siguientes condiciones ′ = 0.2 . y(0) iniciales: y(0) = 0.1 y Sustituyendo las derivadas primera y segunda que aparecen en la ecuación diferencial por las expresiones numéricas correspondientes a una interpolación limitada de segundo orden, tendremos, para el i-ésimo punto: 1 h

2

⋅ ( y i-1 − 2 ⋅ y i + y i +1 ) −

1 ⋅ ( - y i-1 + 0 ⋅ y i + y i +1 ) − 2 ⋅ y i = 0 2 ⋅h

Reemplazando h = 0.1 nos queda, para el i-ésimo punto, la siguiente ecuación:

100 ⋅ ( y i-1 − 2 ⋅ y i + y i +1 ) − 5 ⋅ ( - y i-1 + y i +1 ) − 2 ⋅ y i = 0 Si operamos y reordenamos los términos nos queda: 105 ⋅ y i-1 − 202 ⋅ y i + 95 ⋅ y i +1 = 0

Entonces, si despejamos el valor de y en el punto (i + 1), este nos queda en función de los valores de y ya conocidos en los puntos (i) e (i – 1): y i +1 = 2.126 ⋅ y i − 1.105 ⋅ y i-1

Las condiciones iniciales también deben se discretizadas. Entonces tenemos: y ( 0 ) = 0.1 y (′0 ) =

1 ⋅ ( - y i-1 + y i +1 ) = 0.2 2 ⋅h

Se puede apreciar que la segunda ecuación nos permite expresar un punto exterior al dominio (yi-1), en función de otro interior al mismo (yi+1). Entonces: y -1 = y1 - 0.04

Si aplicamos el operador obtenido los puntos x = 0.1, x = 0.2, .... tendremos: Para x = 0.1

CAPITULO – SOLUCIÓN NUMÉRICA DE ECUACIONES DIFERENCIALES

CÁTEDRA MÉTODOS NUMERICOS

Pág.43

y1 = 2.126 ⋅ y 0 − 1.105 ⋅ y -1 = 2.126 ⋅ ( 0.1 ) − 1.105 ⋅ ( y1 − 0.04 ) y1 = 0.213 − 1.105 ⋅ y1 + 0.044 2.105 ⋅ y1 = 0.257 y1 = 0.117

Para x = 0.2 y 2 = 2.126 ⋅ y1 − 1.105 ⋅ y 0 = 2.126 ⋅ ( 0.117 ) − 1.105 ⋅ ( 0.1 ) = y 2 = 0.138

Para x = 0.3 y 3 = 2.126 ⋅ y 2 −1.105 ⋅ y1 = 2.126 ⋅ ( 0.138 ) −1.105 ⋅ ( 0.117 ) = y 3 = 0.164

para x = 0.4 y 4 = 2.126 ⋅ y 3 − 1.105 ⋅ y 2 = 2.126 ⋅ ( 0.164 ) − 1.105 ⋅ ( 0.138 ) = y 4 = 0.196

para x = 0.5 y 5 = 2.126 ⋅ y 4 − 1.105 ⋅ y 3 = 2.126 ⋅ ( 0.196 ) − 1.105 ⋅ ( 0.164 ) = y 4 = 0.235

2.3.1

Problemas de valores en la frontera.

Un problema de valores en la frontera se definió como aquel en el cual las condiciones de contorno se imponen sobre ambos extremos del intervalo que constituye el dominio de definición. Para ecuaciones diferenciales ordinarias de orden 2n, en general habrá n condiciones de frontera en cada borde x = a y x = b y estas condiciones contendrán derivada hasta de orden 2n-1. La solución de estos problemas será planteada en este apartado en base al procedimientos de Diferencias Finitas, el cual consiste, básicamente, en reemplazar a cada una de las derivadas que aparecen en la expresión de la ecuación diferencia por su aproximación en Diferencias Finitas. CAPITULO – SOLUCIÓN NUMÉRICA DE ECUACIONES DIFERENCIALES

CÁTEDRA MÉTODOS NUMERICOS

Pág.44

El procedimiento de cálculo se ilustrará con el siguiente ejemplo: Ejemplo: resolver la siguiente ecuación diferencial: d2y

−y =0

dx 2

0 ≤ x ≤1

y sujeta a las siguientes condiciones de contorno: y(0)=0 ; y(1)=1 El problema se resolverá dividiendo el dominio de definición en 4 partes iguales, es decir ∆ =0.25. Paso 1: Discretización del dominio x

0

= x0 1

x

2

x x 3

4

=

1

Paso 2: Modelo matemático de la ecuación de gobierno Si se considera una interpolación limitada de segundo orden, la derivada puede aproximarse en el punto x = xi mediante la siguiente expresión en diferencias finitas: d2y dx

2

≈ X =Xi

1 Δ2

[ y i −1 −2 y i + y i +1 ]

la cual, reemplazada en la ecuación de gobierno, proporciona la expresión de la ecuación diferencial valuada en ese punto: 1 Δ2

[ y i −1 − 2 y i + y i +1 ] − y i

=0

la cual puede representarse por el siguiente operador:  1 2 1    =0 , −1 ,  Δ2  Δ2 Δ2  

; [16 , - 33 , 16 ] = 0

el cual constituye el “operador de la ecuación diferencial”. Paso 3: Valoración de la ecuación discretizada en cada uno de los puntos incógnitas: •

En x = x1 → 16 y0 –33 y1 +16 y2 = 0 ∴ -33 y1 +16 y2 = 0 (considerando que y0=0)

CAPITULO – SOLUCIÓN NUMÉRICA DE ECUACIONES DIFERENCIALES

CÁTEDRA MÉTODOS NUMERICOS

Pág.45

• •

En x = x2 → 16 y1 –33 y2 +16 y3 = 0 En x = x3 → 16 y2 –33 y3 +16 y4 = 0 ∴ +16 y2 -33 y1 = -16 (considerando que y4=1)

Las ecuaciones anteriores pueden escribirse en el siguiente como : -33 16 0

16 -33 16

0 16 -33

y1 y2 y3

=

0 0 -16

∴

y1 y2 y3

0.22 = 0.44 0.70

Ejemplo: resolver la ecuación diferencial : d4y dx 4

−16. y = x

0 ≤ x ≤1

sujeta a las condiciones de contorno: y/0) = y’’(0) = 0 ; y(1) = y’(1) = 0 Paso 1: Discretización del dominio

Caso a) ∆ = 1/4

Caso b) ∆ = 1/8

Paso 2: Modelo matemático de la ecuación de gobierno Para representar a las derivadas que aparecen tanto en la ecuación de gobierno como en las condiciones de contorno se considerará una interpolación limitada de segundo orden. De acuerdo a ello: 1 [ −1 2Δ

0

+ 1]

[+1

-2

+ 1]

[+ 1

- 4 + 6 - 4 + 1]

a)

y' i =

b)

y' ' i =

c)

y IV =

1 Δ2 1 Δ4

CAPITULO – SOLUCIÓN NUMÉRICA DE ECUACIONES DIFERENCIALES

CÁTEDRA MÉTODOS NUMERICOS

Pág.46

+1

+1

-2

-2

+1

+1

-2

+1

+1

-2

+1

+6

-4

+1

+1 1/∆

4

+1

-4

Operador de la Ecuación Diferencial: 1 4

Δ

[+ 1

]

- 4 + 6 − 16 Δ4 - 4 + 1 = x

Paso 3: Valoración del operador diferencial en los puntos del dominio. Caso a) •

En el punto x = 0 la solución es conocida, y(0)=0, por lo cual no se valoriza el operador en este punto. 4 4 • En x = 1 → 1/∆ (+1 y* -4 y0 +(6-16∆ ).y1 – 4.y2 + 1 y3) = x1 . En esta ecuación el valor y* no se conoce ya que x* (x =-∆ ) se ubica fuera del dominio del problema. Sin embargo, de acuerdo a las condiciones de borde enunciadas anteriormente, en x=0 debe cumplirse la condición de y’’(0)=0. Discretizando esta condición, resulta: y' '0 =

1 Δ2

[ + 1. y*

- 2 * y 0 + 1. y1 ] = 0

reemplazando arriba, resulta en x =1 → 1/∆ 4

4

∴ y * = − y1

4

(+(5-16∆ ).y1 – 4.y2 + 1 y3) = x1

4

•

En x =2 → 1/∆ (+1y0 -4y1 +(6-16∆ ).y2–4.y3+1y4) = x2 ∴ 4 4 1/∆ (-4 y1+(6-16∆ ).y2 –4.y3) = x2

•

En x =3 → 1/∆ (+1y1 -4y2 +(6 -16∆ ).y3 – 4.y4 + 1y*) = x3 , nuevamente, en la expresión anterior y* correspondería al valor de la función en un punto x* (x =1+∆ ), ubicado fuera del intervalo. Considerando que en x=1 debe cumplirse la condición de borde y’(1) = 0, discretizando esta última, resulta:

4

4

y' 4 =

1 [ − y 3 + y* ] = 0 ∴ 2Δ

y * = y3

CAPITULO – SOLUCIÓN NUMÉRICA DE ECUACIONES DIFERENCIALES

CÁTEDRA MÉTODOS NUMERICOS

Pág.47

que reemplazado en la ecuación anterior, resulta para x = 3 4

4

1/∆ (+1y1 -4y2 +(7-16∆ ).y3 ) = x3 Agrupando las ecuaciones que resultan de valuar el operador diferencial en los puntos x1, x2, x3, y x4 resulta el siguiente sistema de ecuaciones lineales: 4

+5-16∆ -4 +1 4 -4 +6-16∆ -4 +1 -4 +7-16∆

y1 y2 y3

4

x1.∆ = x2.∆ x3.∆

4 4 4

y1 y2 y3

∴

0.0025 0.0034 0.0020

=

Caso b) • • • • • •

4

4

En x =1 → 1/∆ (+(5-16∆ ).y1 – 4.y2 + 1 y3) = x1 4 4 En x =2 → 1/∆ (-4 y1+(6-16∆ ).y2 –4.y3) = x2 4 4 En x =3 → 1/∆ (+1y1 -4y2 +(6 -16∆ ).y3 – 4.y4 + 1y5) = x3 4 4 En x =4 → 1/∆ (+1y2 -4y3 +(6 -16∆ ).y4 – 4.y5 + 1y6) = x4 4 4 En x =5 → 1/∆ (+1y3 -4y4 +(6 -16∆ ).y5 – 4.y6 + 1y7) = x5 4 4 En x =6 → 1/∆ (+1y4 -4y5 +(6 -16∆ ).y6 – 4.y7 ) = x6 4 4 • En x =7 → 1/∆ (+1y5 -4y6 +(7-16∆ ).y7 ) = x7 4

+5-16∆ -4 +1 0 0 0 0 4 -4 +6-16∆ -4 +1 0 0 0 4 +1 -4 +6-16∆ -4 +1 0 0 4 0 +1 -4 +6-16∆ -4 +1 0 4 0 0 +1 -4 +6-16∆ -4 +1 4 0 0 0 +1 -4 +6-16∆ -4 0 0 0 0 +1 -4 +7-16∆

=

0.0012 0.0021 0.0027 0.0027 0.0023 0.0015 0.0005

4

x1.∆ x2.∆ x3.∆ = x4.∆ x5.∆ x6.∆ x7.∆

Análisis de Convergencia 0.0040 0.0035 0.0030 0.0025 Y

y1 y2 y3 y4 y5 y6 y7

y1 y2 y3 y4 y5 y6 y7

0.0020 0.0015 0.0010 0.0005 0.0000 0

0.2

0.4

CAPITULO – SOLUCIÓN NUMÉRICA DE ECUACIONES DIFERENCIALES

CÁTEDRA MÉTODOS NUMERICOS

0.6

0.8

X D=1/8

D=1/4Pág.48

1

4 4 4 4 4 4 4

Finalmente, en la figura de la derecha se representa la solución obtenida en ambos casos. De los 2 ejemplos resueltos se obtienen las siguientes conclusiones: •

• •

La discretización de las condiciones de contorno permite en todos los casos obtener un sistema de ecuaciones del cual solo participan como incógnitas aquellos puntos en donde la solución no es conocida a priori. El sistema de ecuaciones resultantes es simétrico y con una distribución en banda. La solución debe buscarse por el camino de estudiar su convergencia a medida que ∆→ 0.

APÉNDICE 2.1. : Interpolación Polinomial. Diferencias Finitas Considérese la función y = f(x) definida en forma tabular, para la cual se desconoce su expresión analítica. Supóngase que los valores de x0, x1 = xo + ∆ x, x2 = xo + 2.∆ x, ... xn= xn+n.∆ x, todos ellos igualmente espaciados una magnitud ∆ x, se conocen los correspondientes valores y0, y1, ... yn. Estos valores pueden arreglarse como se muestra en la tabla siguiente: Xi

Yi

x0

y0

x1=x0+∆ x

y1

x2=x0+2.∆ x

y2

CAPITULO – SOLUCIÓN NUMÉRICA DE ECUACIONES DIFERENCIALES

CÁTEDRA MÉTODOS NUMERICOS

Pág.49

.

.

xn=x0+n.∆ x

yn

Tabla 1: Función definida en forma tabular Se llaman primeras diferencias hacia delante a las diferencias entre dos valores consecutivos de y. En la tabla 1, las primeras diferencias hacia delante son: ao = y1 – y0 a1 = y2 – y1 a2 = y3 – y2 ............... yn.1=yn – yn-1 que se representan por ∆ yi. De igual forma, las diferencias de las primeras diferencias se llaman segundas diferencias hacia delante, y valen: bo = a1 – ao = y2 – 2.y1 + y0 b1 = a2 – a1 = y3 – 2.y2 + y1 b2 = a3 – a2 = y4 – 2.y3 + y2 .............. bn-2 = an-1 – an-2 = yn – 2.yn-1 + yn-2 las cuales se representan por ∆ 2yi. Nuevamente, las diferencias de las segundas diferencias son las terceras diferencias hacia delante, ∆ 3yi, las cuales resultan: co = b1 – bo = y3 – 3.y2 + 3.y1 – y0 c1 = b2 – b1 = y4 – 3.y3 + 3.y2 – y1 ............................. cn-3 = bn-2 – bn-3 = yn – 3.yn-1 + 3.yn-2 – yn-3 Siguiendo el proceso se definen las cuartas, quintas, etc. diferencias hacia delante. Todas las diferencias pueden ordenarse en una tabla denominada Tabla de Diferencias en donde cada diferencia se indica entre los dos elementos que la producen, como se muestra a continuación: x0

y0 a0=y1-y0

CAPITULO – SOLUCIÓN NUMÉRICA DE ECUACIONES DIFERENCIALES

CÁTEDRA MÉTODOS NUMERICOS

Pág.50

x1=x0+∆ x

y1

b0 = a1-a0 a1=y2-y1

x2=x0+2.∆ x y2

c0=b1-b0 b1 = a2-a1

a2=y3-y2 x3=x0+3.∆ x y3

...... c1=b2-b1

b2 = a3-a2 a3=y4-y3

x4=x0+4.∆ x y4

...... c0=b3-b2

·······

......

······· ·······

.

······· ·······

······· ·······

.

...... cn-3=bn-2-bn-3

bn-2=an-1-an-2 an-1=yn-yn-1

xn=x0+n.∆ x yn Tabla Nº 2: Tabla de Diferencias Aparentemente se tiene la impresión de que el proceso de determinar diferencias es infinito, pero un ejemplo nos demostrará que no lo es para el caso en que los puntos dados estén ubicados sobre un polinomio. Supongamos que los puntos están ubicados sobre un polinomio de grado 2, entonces se tendrá: xk = x0 + k.∆ x → yk = A(x0 + k.∆ x)2 + B(x0 + k.∆ x) + C xk+1 = x0 + (k+1)..∆ x → yk = A(x0 + (k+1)..∆ x)2 + B(x0 + (k+1)..∆ x) + C xk+2 = x0 + (k+2)..∆ x → yk = A(x0 + (k+2)..∆ x)2 + B(x0 + (k+2)..∆ x) + C En base a ello, las primeras diferencias adelante resultarán: ∆ yk = yk+1 – yk = 2.A. ∆ x.x0 + B. ∆ x + (2k+1).A. ∆ x2 ∆ yk+1 = yk+2 – yk+1 = 2.A. ∆ x.x0 + B. ∆ x + (2k+3).A. ∆ x2 y las segundas diferencias: ∆ 2yk = ∆ yk+1 - ∆ yk = 2.A. ∆ x2 Como las segundas diferencias hacia delante son constantes (no dependen de k) e iguales entre si, las terceras, cuartas, etc. diferencias serán nulas. En este ejemplo, se ve que en un polinomio de grado dos se llega a las diferencias de orden dos. Una consecuencia razonable es pensar que el grado de desarrollo de las CAPITULO – SOLUCIÓN NUMÉRICA DE ECUACIONES DIFERENCIALES

CÁTEDRA MÉTODOS NUMERICOS

Pág.51

diferencias es el mismo que el grado del polinomio. En efecto, es posible demostrar que para un polinomio del tipo y= xn + ... la enésima diferencia se hace constante e igual a : ∆ nyk = n!. ∆ xn Lo anterior puede establecerse como un teorema, que tiene como consecuencia el siguiente corolario: “Si en el proceso de obtención de las diferencias hacia delante sucesivas de una función, una de estas se vuelve constante (o aproximadamente constante), puede afirmarse que el conjunto de valores tabulados queda satisfecho exactamente (o muy aproximadamente) por un polinomio de grado igual al orden de la diferencia constante. El ejemplo siguiente ilustra este caso. X

Y

0

2

∆y

∆ 2y

∆ 3y

6 2

8

48 54

4

62

48 96

150 6

212

48 144

294 8

506

48 192

486 10

992

Como las terceras diferencias son constantes, el conjunto de valores indicados están representados por un polinomio de grado 3. En efecto, los mismos han sido tabulados a partir de : y = x3 – x + 2

Ejemplo 1 : De la tabla siguiente obtener el valor de y para x=0.22

CAPITULO – SOLUCIÓN NUMÉRICA DE ECUACIONES DIFERENCIALES

CÁTEDRA MÉTODOS NUMERICOS

Pág.52

Y

1.20 1.00

Y

0.80 0.60 0.40 0.20 0.00 0

0.05

0.1

0.15

0.2

0.25

0.3

0.35

X

0.00

0.00

0.05

0.48

0.10

0.84

0.15

1.00

0.20

0.91

0.25

0.60

0.30

0.14

a) Tabla de Diferencias 0

0.00 0.48

0.05

0.48

-0.12 0.36

0.1

0.84

-0.09 -0.21

0.16 0.15

1.00

-0.04 -0.24

-0.09 0.2

0.91

0.06

-0.22

0.60

0.01

0.02

-0.31 0.25

0.05 -0.01 -0.01 0.05 0.08

-0.15 -0.46

0.3

0.14

CAPITULO – SOLUCIÓN NUMÉRICA DE ECUACIONES DIFERENCIALES

CÁTEDRA MÉTODOS NUMERICOS

Pág.53

VALORES DE LA FUNCIÓN “Y” INTERPOLADA EN X=0.22 Origen

Orden Interpolación

x

xo

yo

k

1

2

3

4

Exac

0.22

0.00

0.00

4.4

2.109

1.231

0.701

0.807

0.808

0.05

0.48

3.4

1.710

0.870

0.797

0.808

0.10

0.84

2.4

1.216

0.806

0.810

0.809

0.15

1.00

1.4

0.997

0.997

0.997

0.997

0.20

0.91

0.4

0.874

0.901

0.906

0.906

Errores[%] Orden Interpolación xo

yo

1

2

3

4

0.00

0.00

161

52

13

0

0.05

0.48

112

8

1

0

0.10

0.84

50

0

0

0

0.15

1.00

23

23

23

23

0.20

0.91

8

11

12

12

CAPITULO – SOLUCIÓN NUMÉRICA DE ECUACIONES DIFERENCIALES

CÁTEDRA MÉTODOS NUMERICOS

Pág.54

APÉNDICE 2.2. : Derivación Numérica. Operadores de Diferencias Finitas. Dada una función y=f(x) definida en forma tabular se trata de calcular el valor de sus derivadas en alguno de los puntos x=x0, x1, x2, .....xn. Si se acepta aproximar a la función f(x) con el polinomio que pasa por los n+1 puntos, se tiene:

y K = y 0 + k . Δy 0 +

k(k −1) 2 k(k −1)(k − 2) 3 Δ y0 + Δ y 0 + ... 2! 3!

(b)

en donde k = [X – X0]/∆ x Derivando ambos miembros de (b) con respecto a “x”, y teniendo en cuenta que el segundo miembro es una función compuesta de x, se tiene: df(x) d  k(k −1) 2 k(k −1)(k − 2) 3  dk = y 0 + k . Δy 0 + Δ y0 + Δ y 0 + ..  dx dk  2! 3!  dx  df(x) 1  2 k −1 2 3 k 2 −6 k + 2 3 = Δ y0 + Δ y 0 + .....  . Δy 0 + dx Δx  2 6   

Derivando esta última expresión con respecto a “x” se obtiene la derivada segunda: d 2 f(x) dx

2

=

1

[. Δ

2

Δx 2

y 0 + (k −1). Δ3 y 0 + .....

]

y derivando nuevamente, se obtiene la derivada tercera: d 3 f(x) dx

3

=

1 Δx

[Δ

3

3

y 0 + .....

]

Si la derivada fuera limitada de primer orden, es decir si solo se consideraran los términos asociados a la primer diferencia, la derivada primera resultaría:

CAPITULO – SOLUCIÓN NUMÉRICA DE ECUACIONES DIFERENCIALES

CÁTEDRA MÉTODOS NUMERICOS

Pág.55

df(x) 1 [. Δy 0 ] = dx Δx

expresión que indica que la derivada primera resultaría constante entre los puntos y1 e y0, por lo cual, la derivada segunda resultaría nula en ese intervalo. Considerando que ∆ y0 = y1 – y0 , la derivada primera en x = x0 resulta: df(x) 1 = [- y 0 + y1 ] dx X =X 0 Δx

Notacionalmente, la expresión anterior se representa como: y' 0 =

1 { −1 Δx

+1 }

en donde se ha subrayado el coeficiente del pivote en donde se está calculando la derivada. Geométricamente, equivale a tomar como primera derivada a la pendiente de la recta que une los puntos (x0,y0) y (x1,y1). La expresión indicada puede utilizarse corriendo el pivote sobre todos los puntos del dominio en donde quiera evaluarse la primer derivada a excepción del último punto del cual se conoce información, (xn, yn).. Para él será necesario utilizar un “operador corrido a la derecha” que determine y’1. Como el operador dado no depende de k, y la interpolación supuesta entre los puntos es lineal, significa que la derivada primera es constante y el operador para el punto extremo a la derecha del intervalo puede escribirse como: y'0 =

1 { -1 Δx

+1 }

Lo cual, como es obvio, implica que la derivada primera para los dos últimos puntos del intervalo es la misma, solo válido cuando h→0. Ejemplo 2: La función tabulada en el Ejemplo 1 corresponde a y=seno(10.x). Se calculará a continuación la derivada numérica en 15.00 cada uno de los puntos tabulados. 10.00

X

Y

0.00

0.00

Y’ [Exacta] Y’[D.F.] 10.00

9.59

y'(x)

5.00 0.00 -5.00 -10.00 CAPITULO – SOLUCIÓN NUMÉRICA DE ECUACIONES DIFERENCIALES

CÁTEDRA MÉTODOS NUMERICOS

-15.00 0.00

0.05

0.10 Exacta

0.15 x

Pág.56

0.20

0.25

D.F. 1er Orden

0.30

0.05

0.48

8.78

7.24

0.10

0.84

5.40

3.12

0.15

1.00

0.71

-1.76

0.20

0.91

-4.16

-6.22

0.25

0.60

-8.01

-9.15

0.30

0.14

-9.90

-9.15

Si ahora se considera que la interpolación es limitada de segundo orden, las expresiones para la estimación de las derivadas resultarían: df(x) 1 = dx Δx d 2 f(x) dx

2

2 k −1 2   . Δy 0 + 2 Δ y 0   

=

1 Δx

[. Δ y ] 2

2

0

y la derivada tercera resultaría nula dentro del intervalo [y0, y2]. Como: ∆ y0 = y1 – y0 ∆ 2y0 = y2 – 2.y1 + y0 las expresiones para la estimación de las derivadas resultan: df(x) 1  2 k −1 = y1 − y 0 + ( y 2 − 2 y1 + y0 )   dx Δx  2 

d 2 f(x) dx

2

=

1 Δx 2

[ y 2 − 2 y1 + y0 ]

Como se observa, es necesario contar con la información de tres puntos para la estimación de ambas derivadas con este orden de interpolación. Si bien la derivada segunda resulta la misma en todos los puntos del intervalo, no es así para la derivada primera la cual depende de k, es decir, de cual de los tres puntos del intervalo [y 0,y2] estemos calculando la derivada. En efecto, en x =x0 (k =0) resultará: df(x) 1 = [ − 3. y 0 + 4. y1 − y 2 ] dx X =X 0 2. Δx

CAPITULO – SOLUCIÓN NUMÉRICA DE ECUACIONES DIFERENCIALES

CÁTEDRA MÉTODOS NUMERICOS

Pág.57

Por otra parte, si lo que se busca es calcular la derivada primera en el punto x =x1, es decir k =1, resultará: df(x) 1 = [− y0 + y2 ] dx X =X 1 2. Δx

y finalmente, en el punto x =x2 (k =2), resultará: df(x) 1 = [ +1. y0 − 4. y1 + 3 y 2 ] dx X =X 2 2. Δx

expresiones que pueden ser escritas en forma de operador de la siguiente manera: y' 0 =

1 [ −3 2. Δx

+ 4 − 1]

y'1 =

1 [ −1 2. Δx

0 + 1]

y'2 =

1 [ +1 2. Δx

- 4 + 3]

Los operadores de derivada segunda no dependen de k, por lo cual: y' '0 =

1 Δx 2

[ +1

− 2 + 1]

1 y' '1 = 1 [ + 1 − 2 + 1] y' ' 2 = Δx 2 [ + 1 − 2 + 1] Δx 2

Ejemplo 3: Calcular los valores de la derivada primera de la función tabulada en el Ejemplo 1 utilizando operadores corridos a izquierda, centrados y corridos a derecha. X

Y

Y’ [Exac]

CI

CE

CD

0.00

0.00

10.00

10.76

0.05

0.48

8.78

9.30

8.41

0.10

0.84

5.40

5.56

5.18

6.07

0.15

1.00

0.71

0.46

0.68

1.06

CAPITULO – SOLUCIÓN NUMÉRICA DE ECUACIONES DIFERENCIALES

CÁTEDRA MÉTODOS NUMERICOS

Pág.58

0.20

0.91

-4.16

0.25

0.60

-8.01

0.30

0.14

-9.90

-3.99

-4.21

-7.68

-8.44 -10.61

15.00 10.00 5.00 Y`(X)

En el cuadro de la derecha se representaron las diferentes aproximaciones junto a la solución exacta. Como se observa, en todos los casos los resultados son satisfactorios y significativamente mas precisos aquellos estimados en base a interpolaciones limitadas de primer orden

-4.75

0.00 -5.00 -10.00 -15.00 0.00

que 0.05

0.10

0.15

0.20

0.25

0.30

X Exacta

CI

CE

CD

Por otra parte, si se compara los resultados detallados en la tabla se concluye que aquellos obtenidos utilizando operadores centrados generan, en promedio, errores más reducidos que en caso de utilizar errores corridos. Esta conclusión es coherente con lo que se intuye si se considera que los operadores centrados recogen “información” de la forma de la curva hacia ambos lados del punto en donde se está estimando la derivada, mientras que los operadores corridos lo hacen solamente considerando lo que ocurre hacia atrás o hacia delante según los casos.

CAPITULO – SOLUCIÓN NUMÉRICA DE ECUACIONES DIFERENCIALES

CÁTEDRA MÉTODOS NUMERICOS

Pág.59