Get documents about "6 Una introducci�n a los m�todos gausslanos para sis- temas
Document Sample
CAPÍTULO 2 :
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
dT
q k'
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
dc
J D
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 2 – SOLUCIÓN NUMÉRICA DE ECUACIONES DIFERENCIALES
CÁTEDRA MÉTODOS COMPUTACIONALES 2 Pág.1
Ley de Faraday
di
δVL L
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 2 – SOLUCIÓN NUMÉRICA DE ECUACIONES DIFERENCIALES
CÁTEDRA MÉTODOS COMPUTACIONALES 2 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 =0
4
G =0
3
G =0 2
G =0
1
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
ellas válidas para el mismo punto inicial. Es decir, si se especifican todas las condiciones
CAPITULO 2 – SOLUCIÓN NUMÉRICA DE ECUACIONES DIFERENCIALES
CÁTEDRA MÉTODOS COMPUTACIONALES 2 Pág.3
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) (n)
y(a) = y0 , y’(a) = y’0 , y’’(a)=y’’0................. y (a)= y 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
x= a x= b
CAPITULO 2 – SOLUCIÓN NUMÉRICA DE ECUACIONES DIFERENCIALES
CÁTEDRA MÉTODOS COMPUTACIONALES 2 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, xa, 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.
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.
CAPITULO 2 – SOLUCIÓN NUMÉRICA DE ECUACIONES DIFERENCIALES
CÁTEDRA MÉTODOS COMPUTACIONALES 2 Pág.5
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 xi-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.
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.
CAPITULO 2 – SOLUCIÓN NUMÉRICA DE ECUACIONES DIFERENCIALES
CÁTEDRA MÉTODOS COMPUTACIONALES 2 Pág.6
2.3.1.1 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
habrá de extrapolarse en forma lineal sobre el paso de longitud h , mediante la aplicación
de la expresión:
CAPITULO 2 – SOLUCIÓN NUMÉRICA DE ECUACIONES DIFERENCIALES
CÁTEDRA MÉTODOS COMPUTACIONALES 2 Pág.7
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:
y0.5 y0 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:
y0.5 1.0 8.5 0.5 5.25
La función solución exacta es:
y - 0.5 x4 4 x3 - 10 x2 8.5 x 1
y para x = 0.5 la solución exacta es:
y - 0.5 0.54 4 0.53 - 10 0.52 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:
Et 2.03125
εt % 100 100 63.1%
verdadero 3.21875
Para el segundo paso:
CAPITULO 2 – SOLUCIÓN NUMÉRICA DE ECUACIONES DIFERENCIALES
CÁTEDRA MÉTODOS COMPUTACIONALES 2 Pág.8
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.04 4 1.03 - 10 1.02 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:
Et 2.875
εt % 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:
x yverdadero yEuler r% local r% global
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
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
la pendiente calculada al inicio del intervalo se hace sobre un intervalo menor, entonces
la solución calculada se aproxime a la exacta.
CAPITULO 2 – SOLUCIÓN NUMÉRICA DE ECUACIONES DIFERENCIALES
CÁTEDRA MÉTODOS COMPUTACIONALES 2 Pág.9
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:
x yverdadero yEuler global
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
3.75 4.31055 5.80469 -34.7
4.00 3.00000 5.00000 -66.7
CAPITULO 2 – SOLUCIÓN NUMÉRICA DE ECUACIONES DIFERENCIALES
CÁTEDRA MÉTODOS COMPUTACIONALES 2 Pág.10
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
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.
CAPITULO 2 – SOLUCIÓN NUMÉRICA DE ECUACIONES DIFERENCIALES
CÁTEDRA MÉTODOS COMPUTACIONALES 2 Pág.11
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'i' 2 y'i'' 3 y i(n) n
y i 1 y i y i h h h .......... h Rn
2! 3! n!
donde h = xi+1 - xi y Rn es el residuo definido como:
y n 1ξ n 1
Rn 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:
f ' (x i , y i ) 2 f ' ' (x i , y i ) 3 f (n 1) (x i , y i ) n
y i 1 y i f(x i , y i ) h h h .. h O(h n 1 )
2! 3! n!
donde O( hn1 ) especifica que el error de funcionamiento local, o residuo, es proporcional
al tamaño del paso elevado a la potencia (n+1)-ésima.
CAPITULO 2 – SOLUCIÓN NUMÉRICA DE ECUACIONES DIFERENCIALES
CÁTEDRA MÉTODOS COMPUTACIONALES 2 Pág.12
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:
f ' (x i , y i ) 2 f '' (x i , y i ) y 'i' ' 3 f (n 1) (x i , y i ) n
Et 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 Ea O(h2 ) 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
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
CAPITULO 2 – SOLUCIÓN NUMÉRICA DE ECUACIONES DIFERENCIALES
CÁTEDRA MÉTODOS COMPUTACIONALES 2 Pág.13
darán resultados perfectos si la función solución de la ecuación diferencial a resolver es
un polinomio de orden n-ésimo.
2.3.1.2 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 yi+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:
yi0 1 yi f ( xi , yi ) h
CAPITULO 2 – SOLUCIÓN NUMÉRICA DE ECUACIONES DIFERENCIALES
CÁTEDRA MÉTODOS COMPUTACIONALES 2 Pág.14
Esta ecuación es la llamada ecuación predictor y da una predicción intermedia de yi+1,
0
denominada y i1 , que permite el cálculo de una estimación de la pendiente al final del
intervalo:
y'i 1 f xi 1 , yi0 1
Luego se obtiene una pendiente promedio para el intervalo:
y y i 1 f xi , y i f xi 1 , y i01
y i
2 2
Esta pendiente promedio se usa para extrapolar nuevamente desde y i a yi+1 usando el
método, ya conocido, de Euler
y i 1 y i
f xi , y i f xi 1 , y i01h
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 y i0 1 y i f ( xi , y i ) h
CORRECTOR 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:
y ij1 - yij-1
1 100 %
εa
j
y i 1
-
donde yij1 e yij1 resultan de dos iteraciones sucesivas, la actual y la anterior
1
respectivamente.
CAPITULO 2 – SOLUCIÓN NUMÉRICA DE ECUACIONES DIFERENCIALES
CÁTEDRA MÉTODOS COMPUTACIONALES 2 Pág.15
Ejemplo. Resolver mediante el método de Heun la siguiente ecuación diferencial:
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
1.3
e0.8 x - e-0.5 x 2 e-0.5 x
Ahora, ya abordando el problema desde el punto de vista numérico, inicialmente
calculamos la pendiente en ( x0 , y0 ):
y0 4 e0 - 0.5 2 3
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:
0
y1 y0 y0 h 2 3 1 5
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 y1 para predecir la
pendiente en el final del intervalo:
y1 f x1 , y1 4 e0.8 1 - 0.5 5 6.402164
0
Esta pendiente se promedia con la pendiente inicial para obtener una pendiente promedio
sobre el intervalo [ 0,1 ]
3 6.402164
y 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:
CAPITULO 2 – SOLUCIÓN NUMÉRICA DE ECUACIONES DIFERENCIALES
CÁTEDRA MÉTODOS COMPUTACIONALES 2 Pág.16
y1 2 4.701082 1 6.701082
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 y1
y1 y0 0
ˆ h
2
Ahora calculamos nuevamente la estimación de la pendiente en el punto final del
intervalo:
y1 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:
3 5.551623
y1 2
ˆ 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 yverdadero 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 2 – SOLUCIÓN NUMÉRICA DE ECUACIONES DIFERENCIALES
CÁTEDRA MÉTODOS COMPUTACIONALES 2 Pág.17
En el siguiente gráfico podemos ver una comparación, para la ecuación polinómica del
ejemploanterior, 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:
f xi f xi 1
y i 1 y i h
2
Por otro lado, si recordamos la expansión de Taylor para una función f(x) alrededor del
punto xi, tendremos:
f ' ' (x i ) 2 f ''' (x i ) 3 f (n) (x i ) n
y i 1 y i f ' (x i ) 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 (xi 1 ) f (xi )
f (xi )
h
CAPITULO 2 – SOLUCIÓN NUMÉRICA DE ECUACIONES DIFERENCIALES
CÁTEDRA MÉTODOS COMPUTACIONALES 2 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 )
h f (xi ) 3 f (xi 1 ) f (xi )
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:
f (xi 1 ) f (xi )
y i 1 y i 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:
f ξ 3
Ea 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.
2.3.1.3 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:
h
y i 1 / 2 y i f ( xi , y i )
2
Este valor predicho se usa para calcular una pendiente en el punto medio del intervalo:
y i 1 / 2 f ( xi 1 / 2 , y i 1 / 2 )
CAPITULO 2 – SOLUCIÓN NUMÉRICA DE ECUACIONES DIFERENCIALES
CÁTEDRA MÉTODOS COMPUTACIONALES 2 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.
2.3.1.3.1 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.
2.3.1.4 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 2 – SOLUCIÓN NUMÉRICA DE ECUACIONES DIFERENCIALES
CÁTEDRA MÉTODOS COMPUTACIONALES 2 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 kn
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. Entonces, 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 2 – SOLUCIÓN NUMÉRICA DE ECUACIONES DIFERENCIALES
CÁTEDRA MÉTODOS COMPUTACIONALES 2 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 yi+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 :
f xi , y i 2
y i 1 y i f xi , y i h h
2!
Si tenemos en cuenta que :
f x, y f x, y dy
f xi , y i
x y dx
Sustituyendo esta en la anterior, tendremos que:
f x, y f x, y dy h2
y i 1 y i f xi , y i h
x y dx 2!
Por otro lado desarrollamos por serie de Taylor a la función k 2, recordando que una serie
de Taylor para dos variables se define como:
g g
gx r,y r gx,y r s ..................
x y
Entonces tendremos que:
CAPITULO 2 – SOLUCIÓN NUMÉRICA DE ECUACIONES DIFERENCIALES
CÁTEDRA MÉTODOS COMPUTACIONALES 2 Pág.22
f f
k2 f xi p1 h, y i q11 k1 h f xi , y i p1 h q11 k1 h .........
x y
Si ahora reemplazamos k1 y k2 desarrollado como serie en la ecuación inicial nos queda:
f f
y i 1 y i a1 h f xi , y i a2 h f xi , y i a2 p1 h2 a2 f xi , y i q11 h2 ....
x y
Si agrupamos términos:
f f
yi1 yi a1 f xi, yi a2 f xi, yi h a2 p1 a2 f xi, yi q11 h2 ....
x y
Si comparamos términos de esta última ecuación con:
f x, y f x, y dy h2
y i 1 y i f xi , y i h
x y dx 2!
para hacer equivalentes las dos ecuaciones se debe cumplir que:
1 1
a1 a2 1 , a2 p1 y a2 q11
2 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:
1 1
a1 1 - a2 , p1 y q11
2 a2 2 a2
Debido a que podemos elegir infinitos valores para a 2, 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.
CAPITULO 2 – SOLUCIÓN NUMÉRICA DE ECUACIONES DIFERENCIALES
CÁTEDRA MÉTODOS COMPUTACIONALES 2 Pág.23
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 :
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. Quda 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.
CAPITULO 2 – SOLUCIÓN NUMÉRICA DE ECUACIONES DIFERENCIALES
CÁTEDRA MÉTODOS COMPUTACIONALES 2 Pág.24
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 y i a1 k1 a2 k2 h
Quedara la ecuación:
1 2
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:
CAPITULO 2 – SOLUCIÓN NUMÉRICA DE ECUACIONES DIFERENCIALES
CÁTEDRA MÉTODOS COMPUTACIONALES 2 Pág.25
k1 f xi , y i 2 0 3 12 0 2 20 0 8.5 8.5
k2 2 0.253 12 0.252 20 0.25 8.5 4.21875
Luego, reemplazando en:
y i 1 y i k2 h
Obtendremos el valor en x = 0.5:
y0.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.3753 12 0.3752 20 0.375 8.5 2.58203125
Reemplazando estos valores en la ecuación del método:
1 2
y i 1 y i k1 k2 h
3 3
Y nos queda para x = 0.5 :
1 2
y(0.5) 1 8.5 2.58203125 0.5 1 4.5546875 0.5 3.27734375
3 3
CAPITULO 2 – SOLUCIÓN NUMÉRICA DE ECUACIONES DIFERENCIALES
CÁTEDRA MÉTODOS COMPUTACIONALES 2 Pág.26
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
x yverdadero
y r % y r % y r %
0.0 1.00000 1.00000 0.0 1.000000 0.0 1.000000 0.0
0.5 3.21875 3.43750 6.8 3.109375 3.4 3.277344 1.8
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 k3 h
CAPITULO 2 – SOLUCIÓN NUMÉRICA DE ECUACIONES DIFERENCIALES
CÁTEDRA MÉTODOS COMPUTACIONALES 2 Pág.27
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:
1
y i 1 y i k1 4 k2 k3 h
6
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:
CAPITULO 2 – SOLUCIÓN NUMÉRICA DE ECUACIONES DIFERENCIALES
CÁTEDRA MÉTODOS COMPUTACIONALES 2 Pág.28
1
y i 1 y i k1 2 k2 2 k3 k4 h
6
donde:
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:
dy
a) f x, y -2 x 3 12 x 2 - 20 x 8.5
dx
CAPITULO 2 – SOLUCIÓN NUMÉRICA DE ECUACIONES DIFERENCIALES
CÁTEDRA MÉTODOS COMPUTACIONALES 2 Pág.29
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:
k1 f xi , y i - 2 xi 3 12 xi 2 - 20 xi 8.5
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
3 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:
1
y i 1 y i 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):
1
y(0.5) 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.
CAPITULO 2 – SOLUCIÓN NUMÉRICA DE ECUACIONES DIFERENCIALES
CÁTEDRA MÉTODOS COMPUTACIONALES 2 Pág.30
dy
b) 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:
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:
h 0.5
y 0.25 y 0 k1 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:
h 0.5
y 0.25 y 0 k2 2 3.510611.
2.877653
2 2
k3 f 0.25,2.877653 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.723392 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:
1
y(0.5) 2 3 2 3.510611 2 3.446785 4.105603 0.5 3.751669
6
La solución verdadera es y(0.5) = 3.751521.
CAPITULO 2 – SOLUCIÓN NUMÉRICA DE ECUACIONES DIFERENCIALES
CÁTEDRA MÉTODOS COMPUTACIONALES 2 Pág.31
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:
1
y i 1 y i 7 k1 32 k3 12 k4 32 k5 7 k6 h
90
donde:
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 k4 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
CAPITULO 2 – SOLUCIÓN NUMÉRICA DE ECUACIONES DIFERENCIALES
CÁTEDRA MÉTODOS COMPUTACIONALES 2 Pág.32
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:
b- a
Esfuerzo nf
h
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:
CAPITULO 2 – SOLUCIÓN NUMÉRICA DE ECUACIONES DIFERENCIALES
CÁTEDRA MÉTODOS COMPUTACIONALES 2 Pág.33
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.
2.3.1.5 Resolución de Ecuaciones Diferenciales Ordinarias de orden n mediante
las técnicas de Runge - Kutta.
2.3.1.5.1 Transformación de una EDO de orden n en un sistema de n EDOs de orden 1
Considérese una ecuación diferencial de la siguiente forma:
d ny dy d 2 y d n 1y
f x, y, , ,......, para x [a,b]
dx n dx dx 2 dx n 1
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:
2 n-1
dy
y a α0 ;
y
a α1 ; d 2 a α2 ;...........; d n-y a αn 1
dx dx dx 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
CAPITULO 2 – SOLUCIÓN NUMÉRICA DE ECUACIONES DIFERENCIALES
CÁTEDRA MÉTODOS COMPUTACIONALES 2 Pág.34
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
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
dy dy1
y2
dx dx
d2y dy2
y3 (1)
2 dx
dx
.
d ny dy n 1
yn
n dx
dx
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 para x [a,b]
dx
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 2 – SOLUCIÓN NUMÉRICA DE ECUACIONES DIFERENCIALES
CÁTEDRA MÉTODOS COMPUTACIONALES 2 Pág.35
dy1 dy
y2
dx dx
dy2 d 2 y
y3
dx dx2
: n ecuaciones con n incognitas
dyn 1 d n y
yn
dx dx n
dy n
f x, y1 , y 2 , y 3 ,......, y n
dx
y1 a y a α0
y 2 a
dy
a α1
dx
d 2 y1
y 3 a a α2 n condiciones iniciales
2
dx
.
n-1
. y n a
d y
n-1
a αn 1
dx
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:
d 3 y x d 2 y x dy x
x 2 3 x2 4 y x 4 e x 2
3 2 dx
dx 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:
CAPITULO 2 – SOLUCIÓN NUMÉRICA DE ECUACIONES DIFERENCIALES
CÁTEDRA MÉTODOS COMPUTACIONALES 2 Pág.36
y ' x v x
y ' ' x v ' x u x
y ''' x v '' x u' x
La ecuación original podrá entonces escribirse en términos del siguiente sistema:
y ' x v x
v ' x u' x
4 e x 2 2 v ' x 3 x 2 y ' x 4 y
u' x
x
Y las condiciones iniciales en x = 0:
y(0) 4; v(0) 1 y u(0) 2
2.3.1.5.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:
dy1
f1 x , y1 , y 2 ,............, y n
dx
dy2
f2 x , y1 , y 2 ,............, yn
dx
.
.
dy n
fn x , y1 , y 2 ,............, yn
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.
CAPITULO 2 – SOLUCIÓN NUMÉRICA DE ECUACIONES DIFERENCIALES
CÁTEDRA MÉTODOS COMPUTACIONALES 2 Pág.37
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.
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
dy2
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
dy2
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 y2 0 f2 ( 0 , 4 , 6) h 6 1.8 0.5 6.9
Para un segundo paso volvemos a calcular las pendientes:
CAPITULO 2 – SOLUCIÓN NUMÉRICA DE ECUACIONES DIFERENCIALES
CÁTEDRA MÉTODOS COMPUTACIONALES 2 Pág.38
dy1
f1 ( 0.5 , 3, 6.9) - 0.5 3 - 1.5
dx
dy2
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:
y1 1.0 y1 0.5 f1 ( 0.5 , 3, 6.9) h 3 - 1.5 0.5 2.25
y 2 1.0 y2 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
dy2
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
CAPITULO 2 – SOLUCIÓN NUMÉRICA DE ECUACIONES DIFERENCIALES
CÁTEDRA MÉTODOS COMPUTACIONALES 2 Pág.39
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.
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 k 3).
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 2 4 - 2 0.5 3.5
y1 h
h
2
2 y 2 0 k1,2 2 6 1.8 0.5 6.45
y2 h
h
2
Estos valores se usarán para calcular el primer conjunto de pendientes de punto medio:
CAPITULO 2 – SOLUCIÓN NUMÉRICA DE ECUACIONES DIFERENCIALES
CÁTEDRA MÉTODOS COMPUTACIONALES 2 Pág.40
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:
'
2 y1 0 k2,1 2 4 - 1.75 0.5 3.5625
y1 h
h
2
2 y2 0 k2,2 2 6 1.715 0.5 6.42875
y '2 h
h
2
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:
1
y i 1 y i k1 2 k2 2 k3 k4 h
6
CAPITULO 2 – SOLUCIÓN NUMÉRICA DE ECUACIONES DIFERENCIALES
CÁTEDRA MÉTODOS COMPUTACIONALES 2 Pág.41
Reemplazando los valores obtenidos calculamos las variables dependientes y1 e y2
respectivamente:
k1,1 2 k2,1 2 k3,1 k4,1 h
1
y1 0.5 y1 0
6
1
y1 0.5 4 - 2 2 - 1.75 2 - 1.78125 1.554688 0.5 3.115234
6
y
k1,2 2 k2,2 2 k3,2 k4,2 h
1
y 2 0.5 y2 0
6
1
y 2 0.5 6 1.8 2 1.715 2 1.715125 1.631794 0.5 6.857670
6
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
2.3.1.6 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.
CAPITULO 2 – SOLUCIÓN NUMÉRICA DE ECUACIONES DIFERENCIALES
CÁTEDRA MÉTODOS COMPUTACIONALES 2 Pág.42
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 dy
2y 0
dx2 dx
en el intervalo entre x = 0 y x = 0.5, con paso h = 0.1, y las siguientes condiciones
iniciales: y(0) 0.1 y
y(0) 0.2 .
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 1
y i-1 2 y i y i 1 - y i-1 0 y i y i 1 2 y i 0
h2 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
1
y 0 - y i-1 y i 1 0.2
2 h
CAPITULO 2 – SOLUCIÓN NUMÉRICA DE ECUACIONES DIFERENCIALES
CÁTEDRA MÉTODOS COMPUTACIONALES 2 Pág.43
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
y1 2.126 y0 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 y0 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
y5 2.126 y 4 1.105 y 3 2.126 0.196 1.105 0.164
y 4 0.235
CAPITULO 2 – SOLUCIÓN NUMÉRICA DE ECUACIONES DIFERENCIALES
CÁTEDRA MÉTODOS COMPUTACIONALES 2 Pág.44
2.3.2 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.
El procedimiento de cálculo se ilustrará con el siguiente ejemplo:
Ejemplo: resolver la siguiente ecuación diferencial:
d2y
y 0 0 x 1
dx2
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 x
x x
0
0 x 1
= 1 2 3 4
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
1
y i 1 2 y i y i 1
dx2 X Xi Δ2
la cual, reemplazada en la ecuación de gobierno, proporciona la expresión de la ecuación
diferencial valuada en ese punto:
1
y i 1 2 y i y i 1 y i 0
Δ2
CAPITULO 2 – SOLUCIÓN NUMÉRICA DE ECUACIONES DIFERENCIALES
CÁTEDRA MÉTODOS COMPUTACIONALES 2 Pág.45
la cual puede representarse por el siguiente operador:
1 1
; , - 33 , 16 0
2
, - 1 , 0 16
Δ2
Δ2 Δ2
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)
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 y1 0 y1 0.22
16 -33 16 y2 = 0 y2 = 0.44
0 16 -33 y3 -16 y3 0.70
Ejemplo: resolver la ecuación diferencial :
d4y
16. y x 0 x 1
dx4
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
CAPITULO 2 – SOLUCIÓN NUMÉRICA DE ECUACIONES DIFERENCIALES
CÁTEDRA MÉTODOS COMPUTACIONALES 2 Pág.46
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:
a) y'i
1
1 0 1
2Δ
b) y' ' i
1
1 -2 1
Δ2
c) y IV
1
1 - 4 6 - 4 1
Δ4
+1 +1 -2 +1
-2 +1 -2 +1
+1 +1 -2 +1
4
1/ +1 -4 +6 -4 +1
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
1. y* - 2 * y 0 1. y1 0 y * y1
Δ2
4 4
reemplazando arriba, resulta en x =1 1/ (+(5-16 ).y1 – 4.y2 + 1 y3) = x1
4 4
En x =2 1/ (+1y0 -4y1 +(6-16 ).y2–4.y3+1y4) = x2
CAPITULO 2 – SOLUCIÓN NUMÉRICA DE ECUACIONES DIFERENCIALES
CÁTEDRA MÉTODOS COMPUTACIONALES 2 Pág.47
4 4
1/ (-4 y1+(6-16 ).y2 –4.y3) = x2
4 4
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:
y'4
1
y3 y* 0 y * y3
2Δ
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 x 1,
x2, x3, y x4 resulta el siguiente sistema de ecuaciones lineales:
4 4
+5-16 -4 +1 y1 x1. y1 0.0025
4 4
-4 +6-16 -4 y2 = x2. y2 = 0.0034
4 4
+1 -4 +7-16 y3 x3. y3 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 4
+5-16 -4 +1 0 0 0 0 y1 x1.
4 4
-4 +6-16 -4 +1 0 0 0 y2 x2.
4 4
+1 -4 +6-16 -4 +1 0 0 y3 x3.
4 4
0 +1 -4 +6-16 -4 +1 0 y4 = x4.
4 4
0 0 +1 -4 +6-16 -4 +1 y5 x5.
4 4
0 0 0 +1 -4 +6-16 -4 y6 x6.
4 4
0 0 0 0 +1 -4 +7-16 y7 x7.
CAPITULO 2 – SOLUCIÓN NUMÉRICA DE ECUACIONES DIFERENCIALES
CÁTEDRA MÉTODOS COMPUTACIONALES 2 Pág.48
Análisis de Convergencia
y1 0.0012
0.0040
y2 0.0021
0.0035
y3 0.0027 0.0030
y4 = 0.0027 0.0025
y5 0.0023
Y
0.0020
y6 0.0015 0.0015
0.0010
y7 0.0005
0.0005
0.0000
Finalmente, en la figura de la 0 0.2 0.4 0.6 0.8 1
derecha se representa la X
solución obtenida en ambos D=1/8 D=1/4
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.
CAPITULO 2 – SOLUCIÓN NUMÉRICA DE ECUACIONES DIFERENCIALES
CÁTEDRA MÉTODOS COMPUTACIONALES 2 Pág.49
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
. .
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:
CAPITULO 2 – SOLUCIÓN NUMÉRICA DE ECUACIONES DIFERENCIALES
CÁTEDRA MÉTODOS COMPUTACIONALES 2 Pág.50
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
x1=x0+x y1 b0 = a1-a0
a1=y2-y1 c0=b1-b0
x2=x0+2.x y2 b1 = a2-a1 ......
a2=y3-y2 c1=b2-b1
x3=x0+3.x y3 b2 = a3-a2 ......
a3=y4-y3 c0=b3-b2
x4=x0+4.x y4 ······· ......
······· ·······
······· . ······· ......
······· 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:
CAPITULO 2 – SOLUCIÓN NUMÉRICA DE ECUACIONES DIFERENCIALES
CÁTEDRA MÉTODOS COMPUTACIONALES 2 Pág.51
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
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 y 2y 3y
0 2
6
2 8 48
54 48
4 62 96
150 48
6 212 144
294 48
8 506 192
486
10 992
CAPITULO 2 – SOLUCIÓN NUMÉRICA DE ECUACIONES DIFERENCIALES
CÁTEDRA MÉTODOS COMPUTACIONALES 2 Pág.52
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
X Y 1.20
0.00 0.00 1.00
0.05 0.48 0.80
0.10 0.84
Y
0.60
0.15 1.00 0.40
0.20 0.91 0.20
0.25 0.60 0.00
0.30 0.14 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35
X
a) Tabla de Diferencias
0 0.00
0.48
0.05 0.48 -0.12
0.36 -0.09
0.1 0.84 -0.21 0.05
0.16 -0.04 0.01
0.15 1.00 -0.24 0.06 -0.01
-0.09 0.02 -0.01
0.2 0.91 -0.22 0.05
-0.31 0.08
0.25 0.60 -0.15
-0.46
0.3 0.14
CAPITULO 2 – SOLUCIÓN NUMÉRICA DE ECUACIONES DIFERENCIALES
CÁTEDRA MÉTODOS COMPUTACIONALES 2 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 2 – SOLUCIÓN NUMÉRICA DE ECUACIONES DIFERENCIALES
CÁTEDRA MÉTODOS COMPUTACIONALES 2 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:
k(k 1) 2 k(k 1)(k 2) 3
y K y0 k . Δy0 Δ y0 Δ y0 ... (b)
2! 3!
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 2! Δ y 0 Δ y 0 ..
dx dk 3! dx
df(x) 1 2 k 1 2 3 k2 6 k 2 3
. Δy 0 Δ y0 Δ y 0 .....
dx Δx
2 6
Derivando esta última expresión con respecto a “x” se obtiene la derivada segunda:
d 2f(x)
dx 2
1
Δx 2
. Δ y 2 3
0 (k 1). Δ y0 .....
y derivando nuevamente, se obtiene la derivada tercera:
d 3f(x)
dx 3
Δx
1
3
Δ y3
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:
df(x)
1
. Δy0
dx Δx
CAPITULO 2 – SOLUCIÓN NUMÉRICA DE ECUACIONES DIFERENCIALES
CÁTEDRA MÉTODOS COMPUTACIONALES 2 Pág.55
expresión que indica que la derivada primera resultaría constante entre los puntos y 1 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
- y0 y1
dx X X 0 Δx
Notacionalmente, la expresión anterior se representa como:
1
y'0 { 1 1 }
Δx
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:
1
{ -1 1 }
y'0
Δx
Lo cual, como es obvio, implica que la derivada primera para los dos últimos puntos del
intervalo es la misma, solo válido cuando h0.
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 cada uno de los puntos tabulados.
15.00
X Y Y’ [Exacta] Y’[D.F.] 10.00
0.00 0.00 10.00 9.59 5.00
0.05 0.48 8.78 7.24
y'(x)
0.00
0.10 0.84 5.40 3.12
-5.00
0.15 1.00 0.71 -1.76
-10.00
0.20 0.91 -4.16 -6.22
-15.00
0.25 0.60 -8.01 -9.15
0.00 0.05 0.10 0.15 0.20 0.25 0.30
0.30 0.14 -9.90 -9.15 x
CAPITULO 2 – SOLUCIÓN NUMÉRICA DE ECUACIONES DIFERENCIALES Exacta D.F. 1er Orden
CÁTEDRA MÉTODOS COMPUTACIONALES 2 Pág.56
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 2 k 1 2
. Δy 0 2 Δ y 0
dx Δx
d 2f(x)
dx 2
1
Δx 2
. Δ y
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:
2 k 1
df(x)
1
y1 y 0 y 2 2 y1 y0
dx Δx 2
d 2 f(x)
2
1
y2 2 y1 y0
dx Δx 2
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 [y0,y2]
estemos calculando la derivada.
En efecto, en x =x0 (k =0) resultará:
df(x)
1
3. y0 4. y1 y2
dx X X 0 2. Δx
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
CAPITULO 2 – SOLUCIÓN NUMÉRICA DE ECUACIONES DIFERENCIALES
CÁTEDRA MÉTODOS COMPUTACIONALES 2 Pág.57
y finalmente, en el punto x =x2 (k =2), resultará:
df(x)
1
1. y0 4. y1 3 y2
dx X X 2 2. Δx
expresiones que pueden ser escritas en forma de operador de la siguiente manera:
y'0
1
3 4 1
2. Δx
y'1
1
1 0 1
2. Δx
y'2
1
1 - 4 3
2. Δx
Los operadores de derivada segunda no dependen de k, por lo cual:
y' '0
1
1 2 1
Δx 2
y' '1
1
1 2 1
Δx 2
y' '2
1
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
0.20 0.91 -4.16 -4.75 -3.99 -4.21
0.25 0.60 -8.01 -7.68 -8.44
0.30 0.14 -9.90 -10.61
CAPITULO 2 – SOLUCIÓN NUMÉRICA DE ECUACIONES DIFERENCIALES
CÁTEDRA MÉTODOS COMPUTACIONALES 2 Pág.58
En el cuadro de la derecha se 15.00
representaron las diferentes 10.00
aproximaciones junto a la 5.00
solución exacta. Como se
Y`(X)
0.00
observa, en todos los casos los
-5.00
resultados son satisfactorios y
-10.00
significativamente mas precisos que
aquellos estimados en base a -15.00
0.00 0.05 0.10 0.15 0.20 0.25 0.30
interpolaciones limitadas de
X
primer orden
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 2 – SOLUCIÓN NUMÉRICA DE ECUACIONES DIFERENCIALES
CÁTEDRA MÉTODOS COMPUTACIONALES 2 Pág.59
