ODE NUMERICAL METHODS 1

COMPUTACION II
GRADO DE FISICAS
Universidad Autónoma de Madrid

Dr. Pablo Luis Pernas

DEFINICIONES BÁSICAS:
ECUACIÓN DIFERENCIAL CONTIENE LAS DERIVADAS DE UNA O MÁS
VARIABLES DEPENDIENTES RESPECTO A UNA O MÁS VARIABLES
INDEPENDIENTES.

ECUACIÓN DIFERENCIAL ORDINARIA : DERIVADAS ORDINARIAS DE

UNA O MÁS VARIABLES DEPENDIENTES RESPECTO A UNA ÚNICA
VARIABLE INDEPENDIENTE.

ECUACIÓN DIFERENCIAL EN DERIVADAS PARCIALES DE UNA O MÁS

VARIABLES DEPENDIENTES RESPECTO A MÁS DE UNA VARIABLE
INDEPENDIENTE.
 Recordar que una ecuación ordinaria de segundo orden lineal:

2
d y dy
2
+ p + Qy = f ( x )
dx dx

Y es no lineal si las constantes son funciones:

d2y dy
2
+ p( x, y ) + Q ( x, y )  y = f ( x )
dx dx

✓Lineal: Ningún producto de variables, sin funciones no lineales en ella.

✓No lineales son “todas las demás”.
 Órdenes superiores

✓ Puede transformarse en un conjunto de n ecuaciones de 1er.

orden, por ejemplo:
d 3x d 2x dx
a 3 + b 2 + c = d (t )  dz 1
 dt = a d ( t ) − bz − cy
dt dt dt

dx d 2x  dy
✓ Haciendo: y , z 2  =z
dt dt  dt
 dx
 dt = y

✓ Son importantes los métodos para soluciones de ecuaciones
de primer orden
 Solución numérica
Ecuaciones diferenciales de 1er orden

Una ecuación diferencial ordinaria de primer orden es una expresión del tipo
siguiente:

y' = f (x, y)

El problema que se suele presentar es el de calcular una función y = f(x) tal que verifique
la ecuación anterior con una condición inicial o de contorno:

y(x0) = y0
El siguiente Teorema de Cauchy sólo garantiza la existencia y unicidad de la solución,
pero no nos dice como encontrarla.
 Problema de Cauchy
Nos limitamos aquí solo al primer orden. Si f(x,y) es analítica en un dominio que contiene al punto (x0,y0), existe
una, y sólo una, función analítica y(x) que verifique la ecuación:

dy
y' = = f (x, y)
dx

con la condición inicial: y(x 0 ) = y0

¡Pero la solución ya la conocemos! Si tomamos suficientes términos en la expansión de la serie de Taylor de f(x,y),
tendríamos un algoritmo de solución para este problema con la precisión que necesitemos.

Sin embargo, desde el punto del código es difícil calcular numéricamente derivadas de orden superior. Surge la idea
de aproximar el valor de las derivadas mediante la evaluación de f entre i e (i + 1): de ahí nacen los métodos que
vemos aquí.
 Método de Euler

y ( x0 + h ) − y ( x 0 )
y  ( x0 ) 
h

hy ( x0 )  y ( x0 + h) − y( x0 )

dy
y ( x=+ hf )(xy,(yx )) + hy ( x )
0 0 0
dx
Siempre que h sea pequeño
Solution :
yi +1 = yi + f ( xi , yi )  h
 Errores: local y global

Error local Error global

y y
yi+1 yi+1
yi yi
o xi xi+1 x o xi xi+1 xi+2 x
 MÉTODO DE EULER
Aunque el cómputo capta la tendencia general de la solución, el error de Euler es considerable.
Este error puede reducirse mediante el uso de un paso h mas pequeño a cambio de mas tiempo de cómputo.
• UNA FUENTE DE ERROR EN EL MÉTODO DE EULER ES QUE LA DERIVADA SE APLICA EN
TODO EL INTERVALO: EXISTEN MODIFICACIONES SIMPLES CON NOMBRES COMO:

✓ HEUN’S METHOD
✓ THE MIDPOINT METHOD
✓ RALSTON’S METHOD
Ejemplo:

y = y − x , y (0) = 0.5
Nos pedía encontrar y(0.8) y h = 0.1
 1
 ( )
y1  y0 + h f ( x0 , y0 ) + f x1 , y1* 
2 
1 
y1  0.5 + 0.1 0.5 + 0.45 = 0.5475
2 
Método de Euler modificado (Heun’s Method)
El método de Euler consiste en quedarse con el primer término del desarrollo de Taylor de la
solución.

Es fácilmente generalizable a ordenes mayores simplemente tomando más términos en el

desarrollo de Taylor. De esta manera se pueden conseguir métodos cuyo error vaya como
O(h2), O(h3), etc.

El problema de estos métodos es que requieren calcular derivadas de ordenes superiores de

la función f(x,y). Los métodos de Runge-Kutta substituyen el cálculo de las derivadas de f(x,y)
por la evaluación de esta función en puntos intermedios, de manera que sigan manteniendo el
mismo error que el método de Taylor correspondiente.

La idea general para estos métodos es substituir en el método de Euler el valor de la pendiente
f(xn,yn), que es exacta sólo en xn, por un especie de valor promedio para el intervalo.
 MÉTODOS RUNGE-KUTTA
Los métodos de Runge-Kutta logran la precisión de una serie de Taylor (truncada) pero sin requerir el cálculo de las derivadas más
altas. En realidad Runge y Kutta desarrollaron estos métodos haciendo coincidir la solución con los n primeros términos de un
desarrollo en serie de Taylor. Análisis numérico con aplicaciones de Gerald Wheatley, ed. Pearson [2].

y i +1 = y i +  (x i , y i , h )h
 = a1k 1 + a2 k 2 + + an k n Función de incremento que representa la
pendiente sobre el intervalo

a ' s = constants
k 1 = f (x i , y i )
k 2 = f (x i + p1h , y i + q11k 1h )
k 3 = f (x i + p 3h , y i + q 21k 1h + q 22 k 2 h )

k n = f (x i + p n −1h , y i + q n −1k 1h + q n −1,2 k 2 h + + q n −1,n −1k n −1h )

p ' s and q ' s are constants
 MÉTODOS RUNGE-KUTTA
Los métodos RK pueden formularse mediante el empleo de diversos número de términos en la función de
incremento .
LOS RK SON UN CONJUNTO O FAMILIA DE MÉTODOS

• El método RK de primer orden (n=1), es sencillamente el método de Euler, es decir, que Euler pasa a ser
un caso especial del método Runge Kutta.

• Los métodos RK de segundo orden: error o(h2)

yi +1 = yi + ( a1k1 + a2k2 )h
k1 = f ( x i , yi )
k2 = f ( xi + p1h, yi + q11k1h )

Los valores de a1, a2, p1 y q11 se evalúan haciendo que la ecuación coincida lo mejor posible con la
expansión de la serie de Taylor (ver libro de referencia página 459)
 MÉTODOS RUNGE-KUTTA N=2

1 1
a1 + a2 = 1, a 2 p1 = , a2 q11 =
2 2
Como se tienen tres ecuaciones con cuatro incógnitas se tiene que suponer el valor de una de ellas.
suponiendo que se especificó un valor para a2, se puede resolver de manera simultánea el sistema de
ecuaciones obtenido:
a1 = 1 − a2
Como se puede elegir un número infinito de valores para a2, evidentemente hay un número infinito de
métodos runge-kutta de segundo orden
Cada versión produciría exactamente los mismos resultados si la solución fuera lineal, cuadrática, o una
constante. sin embargo, da resultados diferentes si la solución es más complicada (lo que casi siempre se da
en el caso de una ode).
Si se toma a1 = ½ = a2 y p1 = 1 = q11, se obtiene el método de Euler modificado!!,
Es decir heun’s es un caso especial del rk de segundo orden.
 dy
= f ( x, y )
RESUMEN MÉTODOS dx
RUNGE-KUTTA N== f2( x, y )
dy
Solution :
yi +1 = yi + (a1k1 + a2 k2 )h
dx
ASÍ LOS
dy TRES DE LOS MÉTODOS MÁS USADOS
Solution SON:: 1 k1 =2 f ( x i , yi )
= f ( x, y ) yi + 1 = yi + ( k1 + k 2 )  h
• HUEN´S METHOD (a2=1/2)
dx
yki +1dy
== f =
y
( ix+f, (
(y
1
x
3
,
) k
y
k
1)+
2 =
3
1 f
k 2
()x h
i + p1 h, y i + q11k1h)
Solution : 1 i 2i dx 2
kk1 =Solution
f ( x , y 3
) 3
2 = f ( xi +
dy 1 1 : h , y i + k 1h )
yi +1 ==yi f+ (( xk,1y+) k 2 )  h
i i

2 2 k 2 =ydyf (==( xiyf +(+x4h,k),y()yhi + 4k1h ) )

dx i +1
dx i 2
k1 = f ( xi , yi )
• THE
Solution
MIDPOINT : METHOD (a2= 1) 1 = f ( xi ,: yi )
kSolution
k 2 = f (( xi + h ), ( yi + k1h ) ) 11 1
dy
yi +1= =f (yxi, +y )k 2  h k 2 = f ( xi + h, yi + k1h )
yi +1 = yi + (2 k1 + 2k 2 )  h
2
dx 3 3
k1 = f ( :xi , METHOD
• RALTSON’S
Solution yi ) (a2= 2/3) k1 = f ( xi , yi )
1 1 2 1 3 3
y = y + ( k + k ) 
ki2+1= f i( xi 3+ 1 h,3 yi2 + k1h )
h k 2 = f ( xi + h, yi + k1h )
2 2 4 4
k1 = f ( xi , yi )
 y
HEUN’S METHOD: CONSISTE EN LA
DETERMINACIÓN DE DOS DERIVADAS, PARA EL
dy
= f ( x, y )
INTERVALO EN EL PUNTO INICIAL Y PARA EL
dxdy
PUNTO FINAL.

Solution = f (: x, y ) f(xi+h,yi+k1h)
dx f(xi,yi)
x
1 1 xi
y i + 1 = yi + ( k1 + k 2 )  h
Solution : xi+h
2 1 21
k1y=i +1f=( xyi i, +yi () 2 k1 + 2 k 2 )  h y
k 2k1==f (f((xxi i +, yhi )), ( yi + k1h ) ) ea
k 2 = f (( xi + h ), ( yi + k1h ) )
Slope: 0.5(k1+k2)
x
xi xi+h
 y

Midpoint Method: Utiliza el método

de Euler para
dypredecir un valor de y en el punto
= f ( x, y )
medio del intervalo:
dy
dx= f ( x, y ) f(xi,yi) f(xi+h/2,yi+k1h/2)
dx
Solution : x
xi xi+h/2
Solution
yi +1 = yi: + k 2  h
yki +1== fy(i x+,ky2 ) h y
1 i i
k1 = f ( xi , yi )1 1 ea
k 2 = f ( xi +1 h, yi +1 k1h )
k 2 = f ( xi + 2h, yi + 2k1h )
2 2 Slope: k2
x
xi xi+h
dy
dy = f ( x, y ) y
dx = f ( x, y )
Ralston’s
dx Method:
Solution :
Solution :
1 2
yi +1 = yi + ( 1k1 + 2k 2 )  h f(xi+ 3/4 h,
yi +1 = yi + (3 k1 +3 k 2 )  h f(xi,yi)
xi
yi+3/4k1h)

k1 = f ( xi , yi )3 3 xi+3/4h x
k1 = f ( xi , y3i ) 3
k 2 = f ( xi + h, yi + k1h ) y
43 4 3 ea
k = f ( x + h, y + k h )
2 i i 1
4 4
Slope:
(1/3k1+2/3k2)
xi xi+h x
22
 Métodos Runge-Kutta n= 3
Métodos de Runge-Kutta de tercer orden

Para n = 3 el resultado son seis ecuaciones con ocho incógnitas, por lo tanto se deben suponer dos valores con
antelación para poder desarrollar el sistema de ecuaciones. Una versión ampliamente usada es:

donde

k1 = f(xi,yi)

k3 = f(xi + h,yi − k1h + 2k2h)

 Métodos Runge-Kutta n = 4 o RK4
El método de Runge-Kutta de cuarto orden

Un miembro de la familia de los métodos Runge-Kutta es usado tan comúnmente que a menudo es referenciado
como “RK4” o como “el método Runge-Kutta”. Definamos un problema de valor inicial como:

Entonces el método RK4 para este problema está dado por la siguiente ecuación:

Donde
 MÉTODOS RK DE CUARTO ORDEN

1
yi +1 = yi + (k1 + 2k2 + 2k3 + k 4 )h
6
Los RK de cuarto orden son los mas usados (y el de la práctica 15 y 17 concretamente) y se obtienen de manera
semejante a los anteriores, pero son de mayor complejidad al tener que comparar términos de hasta h4.

El conjunto de 11 ecuaciones puede resolverse si se escogen dos arbitrariamente y el conjunto mas utilizado
conduce al algoritmo presentado aquí:

k1 = f ( x i , yi )
1 1
k2 = f ( xi + h, yi + k1h )
2 2
1 1
k3 = f ( xi + h, yi + k2h )
2 2
k3 = f ( xi + h, yi + k3h )
 Resumen
En general recordar que el siguiente valor yn+1 es determinado por el presente valor yn más el producto del
tamaño del intervalo h por una pendiente estimada, la función .
Esta pendiente es un promedio ponderado de pendientes: k1 es la pendiente al principio del intervalo; k2 es la
pendiente en el punto medio del intervalo, usando k1 para determinar el valor de “y” en el punto usando el
método de Euler.
k3 es otra vez la pendiente del punto medio, pero ahora usando k2 para determinar el valor de “y” k4 es la
pendiente al final del intervalo, con el valor de y determinado por k3. Promediando las cuatro pendientes, se le
asigna mayor peso a las pendientes en el punto medio:

References:

1. Wikipedia
2. Análisis Numérico con Aplicaciones de Gerald Wheatley, Ed. Pearson.
3. “Lecture 2: Euler's Numerical Method for y'=f(x,y)” https://ocw.mit.edu/courses/mathematics/18-03-
differential-equations-spring-2010/video-lectures/lecture-2-eulers-numerical-method-for-y-f-x-y/