FUNDACIÓN UNIVERSITARIA KONRAD LORENZ

FACULTAD DE MATEMÁTICAS E INGENIERÍAS

MÉTODOS NUMÉRICOS

SOLUCIÓN NUMÉRICA DE ECUACIONES DIFERENCIALES ORDINARIAS

Por favor revise también el enlace donde se encuentra el archivo de Excel con las Hojas de
cálculo correspondientes, este archivo puede ser consultado a través del siguiente enlace
EJERCICIOS RESUELTOS EN EXCEL

METODO DE TAYLOR DE TRES TÉRMINOS

El desarrollo en series de Taylor de tres términos alrededor del punto x=a es

y( x ) = y( a ) + y' ( a ) ⋅
x−a
+ y' ' ( a ) ⋅
( x − a)
2

1! 2!

donde la función y(x) es dos veces derivable. Si en particular se hace a=xn y x= xn+h, entonces
se transforma en
h2
y ( x n + h ) = y ( x n ) + y ' ( x n ) ⋅ h + y ' ' ( x) ⋅
2
o, en forma abreviada
h2
y n +1 = y n + y n' ⋅ h + y n'' ⋅
2
Esta última expresión se conoce como método de Taylor de orden tres para aproximar
soluciones de la ecuación diferencial
y’=f(x,y), y(xo)=yo

Los métodos que se presentarán a continuación tiene la estructura general:

Nuevo valor= valor anterior + pendiente x tamaño del paso, o en términos matemáticos,
yn+1= yi + φh.
En la cual la pendiente φ se usa para extrapolar, desde un valor anterior y i , un nuevo valor
yi+1 en una distancia h. Esta fórmula se puede aplicar paso a paso para calcular el valor en el
futuro, y, por lo tanto, trazar la trayectoria de la solución.

METODO DE EULER O MÉTODO DE LAS TANGENTES

Se desea aproximar la solución de la ecuación diferencia:

y’=f(x,y), y(xo)=yo

Si h>0, entonces sobre la recta tangente en el punto (xo,yo) a la curva solución desconocida,
se encuentra el punto (x1,y1)=(xo+h, y1). Pues bien, la recta que pasa por los punto (xo,yo) y
(x1, y1) será
y1 − y 0
y 0' = o, despejando y1 = y 0 + h ⋅ y 0
'
( x0 + h ) − x0

donde y 0 = f ( x 0 , y 0 )
'
Si se considera que h toma un valor constante, se puede obtener una sucesión de puntos
(x1,y1), (x2, y2), ..., (xn, yn) que sería aproximaciones de los puntos (x1, y(x1)), (x2, y(x2)), ...,
(xn, y(xn)).

Utilizando (x1,y1) se puede obtener el valor de y2 sobre la nueva recta tangente. Así,
y − y1
y1' = 2 , y despejando y 2 = y1 + h.y1'
h
y, en general.
y n +1 = y n + h ⋅ y n' = y n + f ( x n , y n ) ⋅ h , siendo xn=xo+nh o xn=xn-1+h
Si hacemos comparación con la fórmula
yn+1= yi+φh., entonces se puede ver que la pendiente estimada para extrapolar el valor de yn+1
a partir del valor anterior yn ,es φ=f(xn, yn)
La expresión: y n +1 = y n + f ( x n , y n ) ⋅ h , es conocida como el método de Euler o método de las
tangentes.

METODO DE EULER MODIFICADO O MÉTODO DE HEUN

Una fuente fundamental del error en el método de Euler es que la derivada al inicio del
intervalo supuestamente se aplica a todo el intervalo. Para mejorar la estimación de la
pendiente se involucra la determinación de dos derivadas para el intervalo( una en el punto de
inicio y otra en el punto final). Las dos derivadas se promedian con el fin de obtener una
estimación mejorada de la pendiente para todo el intervalo.
Recuerde que para el método de Euler, la pendiente al inicio de un intervalo
y i' = f ( x i , y i )
se utiliza para interpolar linealmente a yn+1

y n*+1 = y n + f ( x n , y n ) ⋅ h
Esta sería la estimación de yn+1 en el caso del método de Euler estándar. Sin embargo, en el
*
método de Euler modificado( o método de Heun), la y n +1 , calculada anteriormente no es la
repuesta final, sino una predicción intermedia. Es por esto que se la distingue con un asterisco
como superíndice. A la ecuación y n +1 = y n + f ( x n , y n ) ⋅ h , se le denomina ecuación predictor.
*

*
Ahora, con y n +1 y xn+1=xn+h, se estima la pendiente al final del intervalo ( xn, xn+1), como
(
y n' +1 = f x n +1 , y n*+1 )
Así las dos pendientes se pueden combinar para obtener la pendiente promedio en el
intervalo.

y =
' y n' + y n' +1
=
( )
f ( x n , y n ) + f x n +1 , y n*+1
2 2
Esta pendiente promedio es la que se utiliza para extrapolar linealmente desde yn hasta yn+1
usando el método de Euler:
y n +1 = y n +
( )
f ( x n , y n ) + f x n +1 , y n*+1
⋅h
2
La cual se conoce como ecuación corrector.
Por lo tanto, el método de Euler modificado o Método de Heun es un procedimiento
predictor–corrector.
Este método se podría resumir este método como la expresión:
y n +1 = y n +
( )
f ( x n , y n ) + f x n +1 , y n*+1
⋅h
2
en donde
y n*+1 = y n + f ( x n , y n ) ⋅ h .

METODOS DE RUNGE KUTTA

Los métodos de Runge-Kutta(RK) logran la exactitud del procedimiento de una serie de

Taylor sin requerir el cálculo de derivadas superiores. Existen muchas variaciones , pero todas
se pueden denotar en la forma generalizada de la ecuación.
y n +1 = y n + φ ( x n , y n , h ) ⋅ h

donde φ ( x n , y n , h ) se conoce como función incremento, la cual puede interpretarse como una
pendiente representativa del intervalo. La función incremento se escribe por lo general como:

φ ( x n , y n , h ) = a1 k1 + a 2 k 2 + ⋅ ⋅ ⋅ + a n k n
donde las “a” son constantes y las k con
k1=f(xi,yi)
k2=f(xi+p1h. yi+q11k1h)
k3=f(xi+p2h, yi+q21k1h+q22k2h)
.
.
.
kn=f(xi+pn-1h, yi+qn-12k2h+…+qn-1 n-1kn-1h)

En las cuales se puede apreciar que las k son relaciones de recurrencia. Esto es, cada k es
función del k anterior.
Es posible concebir varios tipos de Métodos de Runge Kutta al emplear diferente número de
términos en la función incremento, especificados por n. Obsérvese que el método de Runge-
Kutta de orden 1 es el método de Euler. Una vez que se elige n, se evalúan las a, p y q al
igualr la ecuación φ ( x n , y n , h ) = a1 k1 + a 2 k 2 + ⋅ ⋅ ⋅ + a n k n a los términos de un expansión en
series de Taylor. De esta manera, al menos para las versiones de orden inferior, el número de
términos n, representa el orden de la aproximación.

METODOS DE RUNGE-KUTTA DE CUARTO ORDEN

El más popular de todos los métodos de Runge Kutta es el de cuarto orden, de los cuales hay
infinitas versiones. Sin embargo hay una forma en la que convencionalmente se usa que se
conoce como METODO DE RUNGE KUTTA CLASICO DE CUARTO ORDEN, cuya
ecuación de aproximación es:
y n +1 = y n + ( k1,n + 2k 2, n + 2k 3,n + k 4,n ) ⋅ h
1
6
en donde
k1, n = f ( x n , y n )
  1 1 
k 2,n = f  x n + h, y n + k1 ⋅ h 
 2 2 
 1 1 
k 3, n = f  x n + h, y n + k 2 ⋅ h 
 2 2 
k 4,n = f ( x n + h, y n + k 3 ⋅ h )

Nota: en algunos textos, estas ecuaciones aparecen en la forma equivalente:

k1, n = h ⋅ f ( x n , y n )
 1 1 
k 2,n = h ⋅ f  x n + h, y n + k1 
 2 2 
 1 1 
k 3, n = h ⋅ f  x n + h, y n + k 2 
 2 2 
k 4,n = h ⋅ f ( x n + h, y n + k 3 )

Si se compara la expresión del método de Runge Kutta de cuarto orden :

y n +1 = y n + ( k1,n + 2k 2, n + 2k 3,n + k 4,n ) ⋅ h , con la fórmula:
1
6
yn+1= yi+φ h.
Se puede observar que el método de Runge Kutta de cuarto orden utiliza la pendiente
1
( k1,n + 2k 2,n + 2k 3,n + k 4,n ) , para extrapolar linealmente yn+1, desde yn, un tamaño de paso h.
6
Pero si se observa más detenidamente expresión para esta pendiente, se nota claramente que
cada k es una pendiente, que se ponderan para hallar la pendiente promedio en el intervalo.

EJEMPLO No. 1

1) Utilizar el método de Euler para obtener una aproximación de y(0.5) en la ecuación

diferencial.
y’=(x-y)2 , y(0)=0.5

Siendo h=0.1 y realizando los cálculo con cuatro decimales redondeados

Solución
La fórmula de iteración para el método de Euler es:
y n +1 = y n + f ( x n , y n ) ⋅ h
Para este ejemplo: f ( x, y ) = y ' = ( x − y ) 2
Reemplazando, resulta:

y n +1 = y n + ( x n − y n ) ⋅ h = y n + ( x n − y n ) ∗ 0.1
2 2

Comenzando el proceso iterativo con la condición inicial

xo=0, yo=0.5
x1=xo+h=0+0.1=0.1
y1=y(x1)=y(0.1)=yo+(xo -yo)2*0.1
=0.5+(0-0.5)2*0.1=0.525
Resultado: x1=0.1, y(0.1)=0.525
x2=x1+h=0.1+0.1=0.2
y2=y(x2)=y(0.2)= y1+(x1-y1)2*0.1
=0.525+(0.1-0.525)2*0.1=0.5431
Resultado: x2=0.2, y(0.2)=0.5431
x3=x2+h=0.2+0.1=0.3
y3=y(x3)=y(0.3)= y2+(x2-y2)2*0.1
=0.5431+(0.2-0.5431)2*0.1=0.5548
Resultado: x3=0.3, y(0.3)=0.5548
x4=x3+h=0.3+0.1=0.4
y4=y(x4)=y(0.4)= y3+(x3-y3)2*0.1
=0.5548+(0.3-0.5548)2*0.1=0.5613
Resultado: x4=0.4, y(0.4)=0.5613
x5=x4+h=0.4+0.1=0.5
y5=y(x5)=y(0.5)= y4+(x4-y4)2*0.1
=0.5613+(0.4-0.5613)2*0.1=0.5639
Resultado x5=0.5, y(0.5)=0.5639.

Por lo tanto, la aproximación de la solución de la ecuación diferencial en x=0.5 será

y(0.5)=0.5639
Se podría construir una tabla en EXCEL, con el fin de facilitar estos cálculos:

h: 0.1

n xn yn=y(xn) f(xn,yn)
0 0 0.5 0.25
1 0.1 0.525 0.180625
2 0.2 0.5430625 0.11769188
3 0.3 0.55483169 0.06493919
4 0.4 0.56132561 0.02602595
5 0.5 0.5639282 0.00408682

Se podría graficar la aproximación a la función solución de esta ecuación en el intervalo

[0,0.5], graficando los valores obtenidos

x y
0 0.5
0.1 0.525
0.2 0.5431
0.3 0.5548
0.4 0.5613
0.5 0.5639
 Solución de la ecuación diferencial y'=(x-y)^2, en el intervalo [0,0.5]
utilizando en método de Euler

0.57

0.56

0.55

Aproximaciones 0.54

0.53

0.52

0.51

0.5

0.49
0 0.1 0.2 0.3 0.4 0.5
X

EJEMPLO No. 2

Dada la ecuación diferencial

y
y' = x ⋅ y 2 −
x
con la condición inicial y(1)=1.0. Encontrar una aproximación de la solución en x=1.5 con
cuatro cifras decimales mediante el método de Euler y haciendo h=0.05

Solución:
La expresión para el método de Euler es:
y n +1 = y n + f ( x n , y n ) ⋅ h
y
Para este ejemplo: f ( x, y ) = y ' = x ⋅ y −
2
. Reemplazando, resulta
x
 y   y 
y n +1 = y n +  x n ⋅ y n2 − n  ⋅ h = y n +  x n ⋅ y n2 − n  × 0.05
 xn   xn 
Comenzando el proceso iterativo con la condición inicial
xo=1, yo=1
Se obtiene
x1=xo+h=1.0+0.05=1.05
 y   1.0 
y1 = y ( x1 ) = y (1.05) = y 0 +  x 0 ⋅ y 02 − 0  ⋅ 0.05 = 1.0 + 1.0 ∗ 1.0 2 −  × 0.05
 x 0   1 . 0 
y1=1.0000
Resultado: x1=1.05, y(1.05)=1.0000
x2=x1+h=1.05+0.05=1.10
 y   1.0000 
y 2 = y ( x 2 ) = y (1.1) = y1 +  x1 ⋅ y12 − 1  ⋅ 0.05 = 1.0000 + 1.05 ∗ 1.0000 2 −  × 0.05
 x1   1.05 
y2=1.0049.
Resultado: x2=1.1, y(1.1)=1.0049
x3=x2+h=1.10+0.05=1.15
 y2   1.0049 
y 3 = y ( x3 ) = y (1.15) = y 2 +  x 2 ⋅ y 22 −  ⋅ 0.05 = 1.0049 + 1.10 ∗ 1.0049 2 −  × 0.05
 x2   1.10 
y3=1.0147
Resultado: x3=1.15, y(1.15)=1.0147
x4=x3+h=1.15+0.05=1.20
 y3 
 ⋅ 0.05 = 1.0147 + 1.15 ∗ 1.0147 2 −
1.0147 
y 4 = y ( x 4 ) = y (1.20) = y 3 +  x 3 ⋅ y 32 −  × 0.05
 x3   1.15 
y4=1.0298.
Resultado: x4=1.20, y(1.20)=1.0298.
x5=x4+h=1.20+0.05=1.25
 y4   1.0298 
y 5 = y ( x5 ) = y (1.25) = y 4 +  x 4 ⋅ y 42 −  ⋅ 0.05 = 1.0298 + 1.20 ∗ 1.0298 2 −  × 0.05
 x4   1.20 
y5=1.0506
Resultado: x5=1.25, y(1.25)=1.0506
x6=x5+h=1.25+0.05=1.30
 
 ⋅ 0.05 = 1.0506 + 1.25 ∗ 1.0506 2 −
y5 1.0506 
y 6 = y ( x 6 ) = y (1.30) = y 5 +  x5 ⋅ y 52 −  × 0.05
 x5   1.25 
y6=1.0775
Resultado:x6=1.30, y(1.30)=1.0775
x7=x6+h=1.30+0.05=1.35

 
 ⋅ 0.05 = 1.0775 + 1.30 ∗ 1.0775 2 −
y6 1.0775 
y 7 = y ( x 7 ) = y (1.35) = y 6 +  x 6 ⋅ y 62 −  × 0.05
 x6   1.30 
y7=1.1115
Resultado: x7=1.35, y(1.35)=1.1115
x8=x7+h=1.35+0.05=1.40.
 
 ⋅ 0.05 = 1.1115 + 1.35 ∗ 1.1115 2 −
y7 1.1115 
y 8 = y ( x8 ) = y (1.40) = y 7 +  x 7 ⋅ y 72 −  × 0.05
 x7   1.35 
y8=1.1538
Resultado: x8 =1.40, y(1.40)=1.1538
x9=x8+h=1.40+0.05=1.45

 y   1.1538 
y 9 = y ( x9 ) = y (1.45) = y 8 +  x8 ⋅ y 82 − 8  ⋅ 0.05 = 1.1538 + 1.40 ∗ 1.1538 2 −  × 0.05
 x8   1.40 
y9=1.2057

Resultado: x9=1.45, y(1.45)=1.2057

x10=x9+h=1.45+0.05=1.50
 
 ⋅ 0.05 = 1.2057 + 1.45 ∗ 1.2057 2 −
y 1.2057 
y10 = y ( x10 ) = y (1.50) = y 9 +  x 9 ⋅ y 92 − 9  × 0.05
 x9   1.45 
y10=1.2696
RESULTADO: x10=1.50, Y(1.50)=1.2696.
Por lo tanto, la aproximación de la solución de la ecuación diferencial dad, en x=1.5 es
y(1.5)=1.2696.

Igual que para el ejemplo anterior se podría construir una tabla en EXCEL, con el fin de
facilitar estos cálculos:

h: 0.05

n xn yn=y(xn) f(xn,yn)
0 1 1 0
1 1.05 1 0.09761905
2 1.1 1.00488095 0.19723616
3 1.15 1.01474276 0.30177329
4 1.2 1.02983143 0.41447046
5 1.25 1.05055495 0.53913817
6 1.3 1.07751186 0.68048607
7 1.35 1.11153616 0.84458194
8 1.4 1.15376526 1.03952593
9 1.45 1.20574155 1.27648251
10 1.5 1.26956568 1.5713184

Se podría graficar la aproximación a la función solución de esta ecuación en el intervalo

[0,0.5], graficando los valores obtenidos

x y
1.0 1.0000
1.05 1.0000
1.10 1.0049
1.15 1.0147
1.20 1.0298
1.25 1.0506
1.30 1.0775
1.35 1.1115
1.40 1.1538
1.45 1.2057
1.50 1.2696
 Gráfica de la aproximación de la solución a la
ecuación diferencial y'=xy2+y/x, en el intervalo
[1,1.5], utilizando el método de Euler
1.4
1.2
1

Aproximación 0.8

0.6
0.4
0.2
0
1 1.1 1.2 X 1.3 1.4 1.5

EJEMPLO No. 3

Aplicar el método de Euler para obtener una aproximación con cuatro cifras decimales de
redondeo de y(2.0), en la solución de la ecuación diferencial
y ' =| y | −2 x , y(1.0)=1.0
utilizando h=0.25

Solución:
La expresión del método de Euler es
Solución:
La expresión para el método de Euler es:
y n +1 = y n + f ( x n , y n ) ⋅ h
Para este ejemplo: f ( x, y ) = y ' =| y | −2 x . Reemplazando, resulta
y n +1 = y n + ( | y n | −2 x n ) ⋅ h = y n + ( | y n | −2 x n ) × 0.05
Comenzando el proceso iterativo con la condición inicial
xo=1.0, yo=1.0
x1=x0+h=1.0+0.25=1.25
y1=y(x1)=y(1.25)= y 0 + ( | y 0 | −2 x 0 ) × 0.05 =1.0+(|1.0|-2*1.0)*0.05=0.75
Resultado: x1=1.25, y(1.25)=0.75
x2=x1+h=1.25+0.25=1.5
y2=y(x2)=y(1.5)= y1 + ( | y1 | −2 x1 ) × 0.05 =0.75 +(|0.75|-2*1.25)*0.05=0.3125
Resultado: x2=1.5, y(1.5)=0.3125
x3=x2+h=1.5+0.25=1.75
y3=y(x3)=y(1.75)= y 2 + ( | y 2 | −2 x 2 ) × 0.05 =0.3125 +(|0.3125 |-2*1.5)*0.05= -0.3594
Resultado: x3=1.75, y(1.75)=-0.3594
y4=x3+h=1.75+0.25=2.0
y4=y(x4)=y(2.0)= y 3 + ( | y 3 | −2 x 3 ) × 0.05 =-0.3594 +(|-0.3594|-2*1.75)*0.05= -1.1446
Resultado: x4=1.25, y(2.0)=-1.1446

Por lo tanto, la aproximación para y(0.2) en la ecuación diferencial y ' =| y | −2 x , con

condición inicial y(1.0)=1.0, utilizando el método de Euler con un tamaño de paso de
h=0.25 es, y(2.0)= -1.1446.
 La tabla en Excel sería:

h: 0.25

n xn yn=y(xn) f(xn,yn)
0 1 1 -1
1 1.25 0.75 -1.75
2 1.5 0.3125 -2.6875
3 1.75 -0.359375 -3.140625
4 2 -1.14453125 -2.85546875

La gráfica de la aproximación será

Aproximación de la solución de la ecuación diferencial y'=|y|-2x, en
el intervalo [1,2], mediante el método de Euler, con h=0.25

1.5

0.5
Aproximación

0
1 1.2 1.4 1.6 1.8 2
-0.5

-1

-1.5
X

EJEMPLO No. 4

Mediante el método de Euler modificado encontrar una aproximación de y(1.3) en la solución

de la ecuación diferencial
y’=2x-3y+1, y(1.0)=5.0

Utilizando 4 decimales redondeados en los cálculos y h=01.

Solución.

La fórmula iterativa del método de Euler modificado está dada por:

y n +1 = y n +
(
f ( x n , y n ) + f x n +1 , y n*+1
⋅h
)
2
en donde
y n*+1 = y n + f ( x n , y n ) ⋅ h .
Para este ejemplo: f ( x, y ) = y ' = 2 ⋅ x − 3 ⋅ y + 1
Por lo tanto, reemplazando se encuentran las siguientes fórmulas de iteración:
Predictor:
y n*+1 = y n + f ( x n , y n ) ⋅ h = y n + ( 2 x n − 3 y n + 1) ⋅ h
Corrector:
y n +1 = y n +
( )
f ( x n , y n ) + f x n +1 , y n*+1
⋅ h = yn +
( )
( 2 x n − 3 y n + 1) + 2 x n +1 − 3 y n*+1 + 1
× 0.1
2 2
xo=1, yo=5
x1=xo+h=1.0+0.1=1.1
Predictor:
y1* = y 0 + ( 2 x 0 − 3 y 0 + 1) × 0.1 = 5.0 + ( 2 ∗ 1.0 − 3 ∗ 5.0 + 1) × 0.1 =3.8
Corrector:

y1 = y ( x1 ) = y (1.1) = y 0 +
(
( 2 x0 − 3 y 0 + 1) + 2 x1 − 3 y1* + 1 )
× 0.1
2
= 5.0 +
( 2 ∗ 1.0 − 3 ∗ 5.0 + 1) + ( 2 ∗ 1.1 − 3 ∗ 3.8 + 1) × 0.1
2
y1=3.99
Resultado: x1=1.1, y(1.1)=3.99
X2=x1+h=1.1+0.1=1.2

Predictor:
y 2* = y1 + ( 2 x1 − 3 y1 + 1) × 0.1 = 3.99 + ( 2 ∗ 1.1 − 3 ∗ 3.99 + 1) × 0.1 =3.113
Corrector:

y 2 = y ( x 2 ) = y (1.2) = y1 +
(
( 2 x1 − 3 y1 + 1) + 2 x 2 − 3 y 2* + 1) × 0.1
2
= 5.0 +
( 2 ∗ 1.1 − 3 ∗ 3.99 + 1) + ( 2 ∗ 1.2 − 3 ∗ 3.113 + 1) × 0.1
2
y2=3.2546
Resultado: x2=1.2, y(1.2)=3.2546
X3=x2+h=1.2+0.1=1.3
Predictor:
y 3* = y 2 + ( 2 x 2 − 3 y 2 + 1) × 0.1 = 3.113 + ( 2 ∗ 1.2 − 3 ∗ 3.113 + 1) × 0.1 =3.113
Corrector:

y 3 = y ( x 3 ) = y (1.3) = y 2 +
(
( 2 x 2 − 3 y 2 + 1) + 2 x3 − 3 y 3* + 1 )
× 0.1
2
= 5.0 +
( 2 ∗ 1.1 − 3 ∗ 3.99 + 1) + ( 2 ∗ 1.2 − 3 ∗ 3.113 + 1) × 0.1
2

y3=2.7237

Por lo tanto la aproximación de la solución de la ecuación diferencial: y’=2x-3y+1,

y(1.0)=5.0, en x=1.3, utilizando el método de Euler modificado con un tamaño de paso h=0.1
es y(1.3)=2.7237

Al igual que para los casos anteriores se puede elaborar una tabla en Excel con le fin de
facilitar los cálculos. En la siguiente tabla se muestran los resultados hasta x=2.0
 h 0.1

n xn yn f(xn, yn) xn+1 yn+1 f(xn+1,yn+1)

0 1 5.0000 -12.0000 1.1000 3.8000 -8.2000
1 1.1 3.9900 -8.7700 1.2000 3.1130 -5.9390
2 1.2 3.2546 -6.3637 1.3000 2.6182 -4.2546
3 1.3 2.7236 -4.5709 1.4000 2.2665 -2.9996
4 1.4 2.3451 -3.2353 1.5000 2.0216 -2.0647
5 1.5 2.0801 -2.2403 1.6000 1.8561 -1.3682
6 1.6 1.8997 -1.4990 1.7000 1.7498 -0.8493
7 1.7 1.7823 -0.9468 1.8000 1.6876 -0.4628
8 1.8 1.7118 -0.5354 1.9000 1.6582 -0.1747
9 1.9 1.6763 -0.2288 2.0000 1.6534 0.0398
10 2 1.6668 -0.0005 2.1000 1.6668 0.1997

La gráfica de la aproximación en el intervalo [1.0, 2.0], se muestra a continuación:

Aproxim ación de la solución de la ecuación diferencial y'=2x-

3y+1, y(1.0)=5.0, en el intervalo[1.0,2.0], m ediante el m étodo
de Euler Modificado, con h=0.1
6.0000

5.0000
Aproximación

4.0000

3.0000

2.0000

1.0000

0.0000
1 1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8 1.9 2
X

EJEMPLO No. 5

Aplicar El método de Euler Modificado para encontrar un aproximación de y(0.5) en la

solución de la ecuación diferencial
dy
= y ' = sen( x) + y , y(0)=2.0
dx
Utilícense dos decimales redondeados en los cálculos y h=0.1

Solución:
La fórmula iterativa del método de Euler modificado está dada por:
y n +1 = y n +
(
f ( x n , y n ) + f x n +1 , y n*+1
⋅h
)
2
en donde
y n*+1 = y n + f ( x n , y n ) ⋅ h .
Para este ejemplo: f ( x, y ) = y ' = sen( x) + y
Por lo tanto, reemplazando se encuentran las siguientes fórmulas de iteración:
Predictor:
y n*+1 = y n + f ( x n , y n ) ⋅ h = y n + ( sen( x n ) + y n ) ⋅ h
Corrector:
y n +1 = y n +
( )
f ( x n , y n ) + f x n +1 , y n*+1
⋅ h = yn +
( )
( sen( x n ) + y n ) + sen( x n +1 ) + y n*+1
× 0.1
2 2
Se inicia el proceso de aproximación desde la condición inicial: xo=0, y0=2.0
x1=x0+h=0+0.1=0.1
Predictor
y1* = y 0 + ( sen( x 0 ) + y 0 ) ⋅ h =2.0+(sen(0)+2.0)*0.1=2.2
Corrector:

y1 = y 0 +
(
( sen( x0 ) + y 0 ) + sen( x1 ) + y1* ) × 0.1
2
2.0 +
( sen(0) + 2.0) + ( sen(0.1) + 2.2) × 0.1 = 2.21
2
Resultado: x1=0.1, y(0.1)=2.21
x2=x1+h=0.1+0.1=0.2
Predictor
y 2* = y1 + ( sen( x1 ) + y1 ) ⋅ h =2.21+(sen(0.1)+2.21)*0.1=2.44
Corrector:

y 2 = y1 +
(
( sen( x1 ) + y1 ) + sen( x 2 ) + y 2* ) × 0.1
2
2.0 +
( sen(0.1) + 2.21) + ( sen(0.2) + 2.44) × 0.1 = 2.46
2
Resultado: x2=0.2, y(0.2)=2.46
x3=x2+h=0.2+0.1=0.3
Predictor
y 3* = y 2 + ( sen( x 2 ) + y 2 ) ⋅ h =2.46+(sen(0.2)+2.46)*0.1=3.04
Corrector:

y 2 = y1 +
(
( sen( x1 ) + y1 ) + sen( x 2 ) + y 2* ) × 0.1
2
2.0 +
( sen(0.1) + 2.21) + ( sen(0.2) + 2.44) × 0.1 = 3.06
2

Resultado: x3=0.3, y(0.3)=3.06

x4=x3+h=0.3+0.1=0.4
Predictor
y 4* = y 3 + ( sen( x 3 ) + y 3 ) ⋅ h =3.06+(sen(0.3)+3.06)*0.1=3.40
Corrector:

y 2 = y1 +
(
( sen( x1 ) + y1 ) + sen( x 2 ) + y 2* )
× 0.1
2
2.0 +
( sen(0.1) + 2.21) + ( sen(0.2) + 2.44) × 0.1 = 3.43
2
Por lo tanto la aproximación de la solución de la ecuación diferencial: y ' = sen( x) + y ,
y(1.0)=5.0, en x=0.3, utilizando el método de Euler modificado con un tamaño de paso h=0.1
es y(0.3)=3.43

Al igual que para los casos anteriores se puede elaborar una tabla en Excel con el fin de
facilitar los cálculos. En la siguiente tabla se muestran los resultados hasta x=1.0

h 0.1

n xn yn f(xn, yn) xn+1 yn+1 f(xn+1,yn+1)

0 0 2.0000 2.0000 0.1000 2.2000 2.2998
1 0.1 2.2150 2.3148 0.2000 2.4465 2.6451
2 0.2 2.4630 2.6617 0.3000 2.7292 3.0247
3 0.3 2.7473 3.0428 0.4000 3.0516 3.4410
4 0.4 3.0715 3.4609 0.5000 3.4176 3.8970
5 0.5 3.4394 3.9188 0.6000 3.8313 4.3959
6 0.6 3.8551 4.4198 0.7000 4.2971 4.9413
7 0.7 4.3232 4.9674 0.8000 4.8199 5.5373
8 0.8 4.8484 5.5658 0.9000 5.4050 6.1883
9 0.9 5.4361 6.2195 1.0000 6.0581 6.8995
10 1 6.0921 6.9335 1.1000 6.7854 7.6766

La gráfica de la aproximación en el intervalo [0.0, 1.0], se muestra a continuación:

Aproximación de la solución de la ecuación diferencial
y'=sen(x)+y, y(0)=2.0, mediante el método de Euler
Modificado, con h=0.1

7.0

6.0

5.0
Aproximación

4.0

3.0

2.0

1.0

0.0
0 0.2 0.4 0.6 0.8 1
X
 EJEMPLO No. 6

Se considera la ecuación diferencial

y’=x+y-1, y(1.0)=5.0
Con h=0.1 y trabajando con cuatro cifras decimales, halle una estimación de y(1.5),
utilizando:
a) El método de Euler
b) El método de Euler Modificado
c) El método de Taylor de tres términos.
d) Resolver la ecuación diferencial para obtener el valor exacto de y(0.5)
e) Compárense los resultados.

La expresión para el método de Euler es:

y n +1 = y n + f ( x n , y n ) ⋅ h
Para este ejemplo: f ( x, y ) = y ' = x + y − 1 . Reemplazando, resulta
y n +1 = y n + ( x n + y n − 1) ⋅ h = y n + ( x n + y n − 1) × 0.1
Comenzando con la condición inicial: x0=1.0, yo=5.0
x1=xo+h=1.0+0.1=1.1
y1 = y ( x1 ) = y (1.1) = y 0 + ( x 0 + y 0 − 1) ⋅ 0.1 = 5.0 + (1.0 + 5.0 − 1) × 0.1 = 5.5
Resultado: x1=1.1, y(1.1)=5.5
x2=x1+h=1.1+0.1=1.2
y 2 = y ( x 2 ) = y (1.2) = y1 + ( x1 + y1 − 1) ⋅ 0.1 = 5.5 + (1.1 + 5.5 − 1) × 0.1 = 6.06
Resultado: x2=1.2, y(1.2)=6.06
x3=x2+h=1.2+0.1=1.3
y 3 = y ( x 3 ) = y (1.3) = y 2 + ( x 2 + y 2 − 1) ⋅ 0.1 = 6.06 + (1.2 + 6.06 − 1) × 0.1 = 6.686
Resultado: x3=1.3, y(1.3)=6.686
x4=x3+h=1.3+0.1=1.4
y 4 = y ( x 4 ) = y (1.4) = y 3 + ( x 3 + y 3 − 1) ⋅ 0.1 = 6.686 + (1.3 + 6.686 − 1) × 0.1 = 7.3846
Resultado: x4=1.4, y(1.4)=7.3846
x5=x4+h=1.14+0.1=1.5
y 5 = y ( x 5 ) = y (1.5) = y 4 + ( x 4 + y 4 − 1) ⋅ 0.1 = 7.3846 + (1.4 + 7.3846 − 1) × 0.1 = 8.1631
Resultado: x5=1.5, y(1.5)=8.1631
 h 0.1

n xn yn f(xn, yn)
0 1 5.0000 5.0000
1 1.1 5.5000 5.6000
2 1.2 6.0600 6.2600
3 1.3 6.6860 6.9860
4 1.4 7.3846 7.7846
5 1.5 8.1631 8.6631
6 1.6 9.0294 9.6294
7 1.7 9.9923 10.6923
8 1.8 11.0615 11.8615
9 1.9 12.2477 13.1477
10 2 13.5625 14.5625

Aproximación de la solución de la ecuación diferencial y'=(x+y-1),

y(1.0)=5.0, en el intervalo [1.0,2.0], mediante el método de Euler con
h=0.1
16.0
14.0
12.0
10.0
8.0
6.0
4.0
2.0
0.0
1.0 1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8 1.9 2.0

b) La fórmula de iteración de Euler Modificado

y n +1 = y n +
( )
f ( x n , y n ) + f x n +1 , y n*+1
⋅h
2
en donde
y n*+1 = y n + f ( x n , y n ) ⋅ h .
Para este ejemplo: f ( x, y ) = y ' = x + y − 1
Por lo tanto, reemplazando se encuentran las siguientes fórmulas de iteración:
Predictor:
y n*+1 = y n + f ( x n , y n ) ⋅ h = y n + ( x n + y n − 1) ⋅ h
Corrector:
y n +1 = y n +
( )
f ( x n , y n ) + f x n +1 , y n*+1
⋅ h = yn + n
( )
( x + y n − 1) + x n +1 + y n*+1 − 1
× 0.1
2 2
Iniciando el proceso de iteración desde xo=1.0, yo=5.0
x1=xo+h=1.0+0.1=1.1
Predictor
y1* = y 0 + ( x 0 + y 0 − 1) ⋅ h =5.0+(1.0+5.0-1)*0.1=5.5
Corrector:

y1 = y 0 +
( )
( x0 + y 0 − 1) + x1 + y1* − 1
× 0.1
2
5.0 +
(1.0 + 5.0 − 1) + (1.1 + 5.5 − 1) × 0.1 = 5.53
2
Resultado: x1=1.1, y(1.1)=5.53
x2=x1+h=1.1+0.1=1.2
Predictor
y 2* = y1 + ( x1 + y1 − 1) ⋅ h =5.53+(1.1+5.53-1)*0.1=6.0930
Corrector:

y 2 = y1 +
( )
( x1 + y1 − 1) + x 2 + y 2* − 1
× 0.1
2
5.53 +
(1.1 + 5.53 − 1) + (1.2 + 6.0930 − 1) × 0.1 = 6.1262
2

Resultado: x2=1.2, y(1.2)=6.1262

x3=x2+h=1.2+0.1=1.3
Predictor
y 3* = y 2 + ( x 2 + y 2 − 1) ⋅ h =6.1262+(1.2+6.1262-1)*0.1=6.7588
Corrector:

y3 = y 2 + 2
(
( x + y 2 − 1) + x3 + y 3* − 1 )
× 0.1
2
6.1262 +
(1.2 + 6.1262 − 1) + (1.3 + 6.7588 − 1) × 0.1 = 6.7954
2
Resultado: x3=1.3, y(1.3)=6.7954
x4=x3+h=1.3+0.1=1.4
Predictor
y 4* = y 3 + ( x 3 + y 3 − 1) ⋅ h =6.7954+(1.3+6.7954-1)*0.1=7.5049
Corrector:

y 4 = y3 + 3
(
( x + y 3 − 1) + x 4 + y 4* − 1 )
× 0.1
2
6.7588 +
(1.3 + 6.7588 − 1) + (1.4 + 7.5049 − 1) × 0.1 = 7.5454
2

Resultado: x4=1.4, y(1.4)=7.5454

x5=x4+h=1.4+0.1=1.5
Predictor
y 5* = y 4 + ( x 4 + y 4 − 1) ⋅ h =7.5454+(1.2+6.1262-1)*0.1=8.3400
Corrector:
 ( x 4 + y 4 − 1) + ( x5 + y 5* − 1)
y5 = y 4 + × 0.1
2
7.5454 +
(1.4 + 7.5454 − 1) + (1.5 + 8.3400 − 1) × 0.1 = 8.3847
2

Resultado: x5=1.5, y(1.5)=8.3847

h 0.1

n xn yn f(xn, yn) xn+1 yn+1 f(xn+1,yn+1)

0 1 5.0000 5.0000 1.1000 5.5000 5.6000
1 1.1 5.5300 5.6300 1.2000 6.0930 6.2930
2 1.2 6.1262 6.3262 1.3000 6.7588 7.0588
3 1.3 6.7954 7.0954 1.4000 7.5049 7.9049
4 1.4 7.5454 7.9454 1.5000 8.3400 8.8400
5 1.5 8.3847 8.8847 1.6000 9.2731 9.8731
6 1.6 9.3226 9.9226 1.7000 10.3148 11.0148
7 1.7 10.3694 11.0694 1.8000 11.4764 12.2764
8 1.8 11.5367 12.3367 1.9000 12.7704 13.6704
9 1.9 12.8371 13.7371 2.0000 14.2108 15.2108
10 2 14.2845 15.2845 2.1000 15.8129 16.9129

Aproximación de la solución de la ecuación diferencial y'=(x+y-1), en

el intervalo [1,2], utilizando el método de Euler Modificado, con h=0.1
16.0
14.0
12.0
Aproximación

10.0
8.0
6.0
4.0
2.0
0.0
1 1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8 1.9 2
x

c) La expresión del método de Taylor de tres términos, es:

2
'' h
y n +1 = y n + y n ⋅ h + y n ⋅
'

2
Para este caso: y n = f ( x n , y n ) = ( x n + y n − 1)
'

Derivando esta expresión implícitamente con respecto a xn, se obtiene

( )
y n'' = 1 + y n' = (1 + ( x n + y n − 1) ) = x n + y n
 h2
Reemplazando en y n +1 = y n + y n' ⋅ h + y n'' ⋅ , se obtiene
2
h2
y n +1 = y n + ( x n + y n − 1) ⋅ h + ( x n + y n ) ⋅
2
0.12
y n +1 = y n + ( x n + y n − 1) ⋅ 0.1 + ( x n + y n ) ⋅
2
Comenzando en xo=1.0, yo=5.0
x1=xo+h=1.0+0.1=1.1
0.12
y1 = y ( x1 ) = y (1.1) = y 0 + ( x 0 + y 0 − 1) ⋅ 0.1 + ( x 0 + y 0 ) ⋅
2
0.12
y1 = 5.0 + (1.0 + 5.0 − 1) ⋅ 0.1 + (1.0 + 5.0 ) ⋅ = 5.53
2
Resultado: x1=1.1, y(1.1)=5.53
x2=x1+h=1.1+0.1=1.2
0.12
y 2 = y ( x 2 ) = y (1.2) = y1 + ( x1 + y1 − 1) ⋅ 0.1 + ( x1 + y1 ) ⋅
2
2
0.1
y 2 = 5.53 + (1.1 + 5.53 − 1) ⋅ 0.1 + (1.1 + 5.53) ⋅ = 6.1262
2
Resultado: x2=1.2, y(1.2)=6.1262

x3=x2+h=1.2+0.1=1.3
0 .1 2
y 3 = y ( x 3 ) = y (1.3) = y 2 + ( x 2 + y 2 − 1) ⋅ 0.1 + ( x 2 + y 2 ) ⋅
2
2
0.1
y 3 = 6.1262 + (1.2 + 6.1262 − 1) ⋅ 0.1 + (1.2 + 6.1262) ⋅ = 6.7954
2
Resultado: x3=1.3, y(1.3)=6.7954
x4=x3+h=1.3+0.1=1.4
0.12
y 4 = y ( x 4 ) = y (1.4) = y 3 + ( x 3 + y 3 − 1) ⋅ 0.1 + ( x 3 + y 3 ) ⋅
2
2
0.1
y 4 = 6.7954 + (1.3 + 6.7954 − 1) ⋅ 0.1 + (1.3 + 6.7954) ⋅ = 7.5454
2
Resultado: x4=1.4, y(1.4)=7.5454

x5=x4+h=1.4+0.1=1.5
0 .1 2
y 5 = y ( x 5 ) = y (1.5) = y 4 + ( x 4 + y 4 − 1) ⋅ 0.1 + ( x 4 + y 4 ) ⋅
2
2
0.1
y 5 = 7.5454 + (1.4 + 7.5454 − 1) ⋅ 0.1 + (1.4 + 7.5454) ⋅ = 8.3847
2
Resultado: x5=1.5, y(1.5)=8.3847
 h 0.1

n xn yn
0.0000 1.0000 5.0000
1.0000 1.1000 5.5300
2.0000 1.2000 6.1262
3.0000 1.3000 6.7954
4.0000 1.4000 7.5454
5.0000 1.5000 8.3847
6.0000 1.6000 9.3226
7.0000 1.7000 10.3694
8.0000 1.8000 11.5367
9.0000 1.9000 12.8371
10.0000 2.0000 14.2845

Aproximación de la solución de la ecuación diferencial y'=(x+y-1),

y(1.0)=5.0, mediante el método de Taylor de tres términos con
h=0.1
16.0
14.0
12.0
10.0
8.0
yn

6.0
4.0
2.0
0.0
1.0 1.2 1.4 xn 1.6 1.8 2.0

d) Reordenando la ecuación inicial resulta:

y’-y=x-1
que se trata de una ecuación diferencial lineal( es decir, de la forma y’+P(x)y=Q(x)) con
P(x)=-1, y
Q(x)=x-1
Cuya solución general es: y=-x+Cex
Reemplazando la condición inicial y(1.0)=5.0, resulta: 5.0=-1.0+Ce1.0, de donde
C=6e-1.0
Con lo que la solución particular, queda como: y=-x+6e-1.0ex = -x+6ex-1.0

Así, y(1.5)= -1.5+6*e1.5-1.0=-1.5+6e0.5=8.3923

 Solución exacta de la ecuación diferencial y'=(x+y-1), y(1.0)=5.0,
en el intervalo{1.0,2.0]
16

14

12

10

8
y

0
1 1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8 1.9 2
x

COMPARACIÓN

Euler Mod Taylor 3 term. Sol. Exacta Error Euler Error E. Mod2 Error Taylor2
2
n xnEuler
0 1 5.0000 5.0000 5.0000 5.0000 0.0000 0.0000 0.0000
1 1.1 5.5000 5.5300 5.5300 5.5310 0.0010 0.0000 0.0000
2 1.2 6.0600 6.1262 6.1262 6.1284 0.0047 0.0000 0.0000
3 1.3 6.6860 6.7954 6.7954 6.7992 0.0128 0.0000 0.0000
4 1.4 7.3846 7.5454 7.5454 7.5509 0.0277 0.0000 0.0000
5 1.5 8.1631 8.3847 8.3847 8.3923 0.0526 0.0001 0.0001
6 1.6 9.0294 9.3226 9.3226 9.3327 0.0920 0.0001 0.0001
7 1.7 9.9923 10.3694 10.3694 10.3825 0.1523 0.0002 0.0002
8 1.8 11.0615 11.5367 11.5367 11.5532 0.2418 0.0003 0.0003
9 1.9 12.2477 12.8371 12.8371 12.8576 0.3720 0.0004 0.0004
10 2 13.5625 14.2845 14.2845 14.3097 0.5584 0.0006 0.0006
SUMATORIA DE LOS ERRORES ELEVADOS AL CUADRADO 1.5151 0.0017 0.0017

GRAFICO COMPARATIVO ENTRE LOS DIFERENTES MÉTODOS PARA RESOLVER LA

ECUACIÓN DIFERENCIAL y'=(x+y-1), y(1.0)=5.0

16.0

14.0

12.0

10.0
yn

8.0 Euler
Euler Mod
6.0 Taylor 3 term.
Sol. Exacta

4.0
1 1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8 1.9 2
xn
Como se puede apreciar en la tabla, el orden de aproximación del método de Euler
modificado y el de Taylor de tres términos es muy similar y para este ejemplo su error es muy
pequeño, aunque va creciendo a medida que se aleja de la condición inicial, loc aul es de
esperarse, ya que lo que se está realizando es un proceso de extrapolación desde la condición
inicial, por lo cual, a medida que se aleja de esta condición, la incertidumbre crece, y se va
acumulando los errores de los pasos anteriores. El método de Euler presenta un
comportamiento desfavorable, con respecto a los dos métodos anteriores, lo cual también era
de esperarse dado el orden de aproximación del método de Euler( orden 1). Todo esto se
puede apreciar tanto en la tabla como en las gráficas.

EJEMPLOS MÉTODO DE RUNGE-KUTTA DE CUARTO ORDEN

EJEMPLO No. 7

Use el método de Runge-Kutta de cuarto orden para estimar y(0.5) en la ecuación diferencial
y’=-2x3+12x2-20x+8.5, y(0.0)=1.0. Utilice un h=0.1
Obsérvese que para este caso f(x,y)= -2x3+12x2-20x+8.5
Las fórmulas de iteración del método de Runge-Kutta son las siguientes:
1
yn +1 = yn + ( k1, n + 2k2, n + 2k3, n + k4, n ) ⋅ h
6
en donde
k1, n = f ( x n , y n )
 1 1 
k 2,n = f  x n + h, y n + k1, n ⋅ h 
 2 2 
 1 1 
k 3, n = f  x n + h, y n + k 2, n ⋅ h 
 2 2 
k 4,n = f ( x n + h, y n + k 3,n ⋅ h )

Inciando el proceso iterativo con la condición inicial xo=0.0, yo=1.0

x1=xo+h=0.0+0.1=0.1

k1, 0 = f ( x 0 , y 0 ) = −2 x 03 + 12 x 02 − 20 x 0 + 8.5
= −2 ∗ 0.0 3 + 12 ∗ 0.0 2 − 20 ∗ 0.0 + 8.5 = 8.5

3 2
 1 1   1   1   1 
k 2,0 = f  x0 + h, y0 + k1,0 ⋅ h  = −2 x0 + h  + 12 x0 + h  − 20 x0 + h  + 8.5
 2 2   2   2   2 
3 2
 1   1   1 
= −2 0.0 + 0.1 + 12 0.0 + 0.1 − 20 0.0 + 0.1 + 8.5
 2   2   2 
= 7.52975
3 2
 1 1   1   1   1 
k 3, 0 = f  x 0 + h, y 0 + k 2, 0 ⋅ h  = −2 x 0 + h  + 12 x 0 + h  − 20 x 0 + h  + 8.5
 2 2   2   2   2 
3 2
 1   1   1 
= −2 0.0 + 0.1 + 12 0.0 + 0.1 − 20 0.0 + 0.1 + 8.5
 2   2   2 
= 7.52975
k 4, 0 = f ( x 0 + h, y 0 + k 3, 0 ⋅ h ) = −2( x 0 + h ) + 12( x 0 + h ) − 20( x 0 + h ) + 8.5
3 2

= −2( 0.0 + 0.1) + 12( 0.0 + 0.1) − 20( 0.0 + 0.1) + 8.5
3 2

= 6.618

luego se aplica la fórmula de recursión:

y1 = y 0 + ( k1,0 + 2k 2,0 + 2k 3,0 + k 4, 0 ) ⋅ 0.1
1
6
1
y1 = 1.0 + ( 8.5 + 2 ∗ 7.52975 + 2 ∗ 7.52795 + 6.618) ⋅ 0.1 = 1.75395
6
Resultado: x1=0.1, y(0.1)=1.75395
x1=x1+h=0.1+0.1=0.2

k1,1 = f ( x1 , y1 ) = −2 x13 + 12 x12 − 20 x1 + 8.5

= −2 ∗ 0.13 + 12 ∗ 0.12 − 20 ∗ 0.1 + 8.5 = 6.618

3 2
 1 1   1   1   1 
k 2,1 = f  x1 + h, y1 + k1,1 ⋅ h  = −2 x1 + h  + 12 x1 + h  − 20 x1 + h  + 8.5
 2 2   2   2   2 
3 2
 1   1   1 
= −2 0.1 + 0.1 + 12 0.1 + 0.1 − 20 0.1 + 0.1 + 8.5
 2   2   2 
= 5.76325
3 2
 1 1   1   1   1 
k 3,1 = f  x1 + h, y1 + k 2,1 ⋅ h  = −2 x1 + h  + 12 x1 + h  − 20 x1 + h  + 8.5
 2 2   2   2   2 
3 2
 1   1   1 
= −2 0.1 + 0.1 + 12 0.1 + 0.1 − 20 0.1 + 0.1 + 8.5
 2   2   2 
= 5.76325
k 4,1 = f ( x1 + h, y1 + k 3,1 ⋅ h ) = −2( x1 + h ) + 12( x1 + h ) − 20( x1 + h ) + 8.5
3 2

= −2( 0.1 + 0.1) + 12( 0.1 + 0.1) − 20( 0.1 + 0.1) + 8.5
3 2

= 4.964

luego se aplica la formula de recursión:

y 2 = y1 + ( k1,1 + 2k 2,1 + 2k 3,1 + k 4,1 ) ⋅ 0.1
1
6
1
y 2 = 1.75395 + ( 6.618 + 2 ∗ 5.76325 + 2 ∗ 5.76325 + 4.964) ⋅ 0.1 = 2.3312
6
Resultado: x2=0.2, y(0.2)=2.3312
x3=x2+h=0.2+0.1=0.3
k1, 2 = f ( x 2 , y 2 ) = −2 x 23 + 12 x 22 − 20 x 2 + 8.5
= −2 ∗ 0.2 3 + 12 ∗ 0.2 2 − 20 ∗ 0.2 + 8.5 = 4.964
 3 2
 1 1   1   1   1 
k 2, 2 = f  x 2 + h, y 2 + k1, 2 ⋅ h  = −2 x 2 + h  + 12 x 2 + h  − 20 x 2 + h  + 8.5
 2 2   2   2   2 
3 2
 1   1   1 
= −2 0.2 + 0.1 + 12 0.2 + 0.1 − 20 0.2 + 0.1 + 8.5
 2   2   2 
= 4.21875
3 2
 1 1   1   1   1 
k 3, 2 = f  x 2 + h, y 2 + k 2, 2 ⋅ h  = −2 x 2 + h  + 12 x 2 + h  − 20 x 2 + h  + 8.5
 2 2   2   2   2 
3 2
 1   1   1 
= −2 0.2 + 0.1 + 12 0.2 + 0.1 − 20 0.2 + 0.1 + 8.5
 2   2   2 
= 4.21875
k 4, 2 = f ( x 2 + h, y 2 + k 3, 2 ⋅ h ) = −2( x 2 + h ) + 12( x 2 + h ) − 20( x 2 + h ) + 8.5
3 2

= −2( 0.2 + 0.1) + 12( 0.2 + 0.1) − 20( 0.2 + 0.1) + 8.5
3 2

= 3.526

luego se aplica la fórmula de recursión:

y 3 = y 2 + ( k1, 2 + 2k 2, 2 + 2k 3, 2 + k 4, 2 ) ⋅ 0.1
1
6
1
y 3 = 2.3312 + ( 4.964 + 2 ∗ 4.21875 + 2 ∗ 4.21875 + 3.526 ) ⋅ 0.1 = 2.75395
6

Resultado: x3=0.3, y(0.3)= 2.75395

x4=x3+h=0.3+0.1=0.4

k1,3 = f ( x 3 , y 3 ) = −2 x 33 + 12 x 32 − 20 x3 + 8.5
= −2 ∗ 0.3 3 + 12 ∗ 0.3 2 − 20 ∗ 0.3 + 8.5 = 3.526

3 2
 1 1   1   1   1 
k 2, 3 = f  x 3 + h, y 3 + k1,3 ⋅ h  = −2 x 3 + h  + 12 x 3 + h  − 20 x 3 + h  + 8.5
 2 2   2   2   2 
3 2
 1   1   1 
= −2 0.3 + 0.1 + 12 0.3 + 0.1 − 20 0.3 + 0.1 + 8.5
 2   2   2 
= 2.88425
3 2
 1 1   1   1   1 
k 3,3 = f  x 3 + h, y 3 + k 2, 3 ⋅ h  = −2 x 3 + h  + 12 x 3 + h  − 20 x 3 + h  + 8.5
 2 2   2   2   2 
3 2
 1   1   1 
= −2 0.3 + 0.1 + 12 0.3 + 0.1 − 20 0.3 + 0.1 + 8.5
 2   2   2 
= 2.88425
k 4,3 = f ( x 3 + h, y 3 + k 3,3 ⋅ h ) = −2( x 3 + h ) + 12( x 3 + h ) − 20( x 3 + h ) + 8.5
3 2

= −2( 0.3 + 0.1) + 12( 0.3 + 0.1) − 20( 0.3 + 0.1) + 8.5
3 2

= 2.292

luego se aplica la fórmula de recursión:

y 4 = y3 +
1
( k1,3 + 2k 2,3 + 2k 3,3 + k 4,3 ) ⋅ 0.1
6

1
y 4 = 2.3312 + ( 3.526 + 2 ∗ 2.88425 + 2 ∗ 2.88425 + 2.292) ⋅ 0.1 = 3.0432
6

Resultado: x4=0.4, y(0.4)= 3.0432

x5=x4+h=0.4+0.1=0.5

k1, 4 = f ( x 4 , y 4 ) = −2 x 43 + 12 x 42 − 20 x 4 + 8.5
= −2 ∗ 0.4 3 + 12 ∗ 0.4 2 − 20 ∗ 0.4 + 8.5 = 2.292

3 2
 1 1   1   1   1 
k 2, 4 = f  x 4 + h, y 4 + k1, 4 ⋅ h  = −2 x 4 + h  + 12 x 4 + h  − 20 x 4 + h  + 8.5
 2 2   2   2   2 
3 2
 1   1   1 
= −2 0.4 + 0.1 + 12 0.4 + 0.1 − 20 0.4 + 0.1 + 8.5
 2   2   2 
= 1.74775
3 2
 1 1   1   1   1 
k 3, 4 = f  x 4 + h, y 4 + k 2, 4 ⋅ h  = −2 x 4 + h  + 12 x 4 + h  − 20 x 4 + h  + 8.5
 2 2   2   2   2 
3 2
 1   1   1 
= −2 0.4 + 0.1 + 12 0.4 + 0.1 − 20 0.4 + 0.1 + 8.5
 2   2   2 
= 1.74775
k 4, 4 = f ( x 4 + h, y 4 + k 3, 4 ⋅ h ) = −2( x 4 + h ) + 12( x 4 + h ) − 20( x 4 + h ) + 8.5
3 2

= −2( 0.4 + 0.1) + 12( 0.4 + 0.1) − 20( 0.4 + 0.1) + 8.5
3 2

= 1.25

luego se aplica la formula de recursión:

y 5 = y 4 + ( k1, 4 + 2k 2, 4 + 2k 3, 4 + k 4, 4 ) ⋅ 0.1
1
6

1
y 5 = 3.0432 + ( 2.292 + 2 ∗ 1.74775 + 2 ∗ 1.74775 + 1.25) ⋅ 0.1 = 3.21875
6
Resultado: x5=0.5, y(0.5)=3.21875

Por lo tanto la estimación de y(0.5) en la ecuación diferencial: y’=-2x3+12x2-20x+8.5,

y(0.0)=1.0, con h=0.1, utilizando el método de Runge-Kutta de cuarto orden estándar, es:
y(0.5)= 3.21875

Esta integral es de variables separables, luego realizando el despeje correspondiente se

obtiene
dy
dx
(
= −2 x 3 + 12 x 2 − 20 x + 8.5 ∴ dy = − 2 x 3 + 12 x 2 − 20 x + 8.5 dx )
 y x

∫ dy = ∫ ( − 2 x )
+ 12 x 2 − 20 x + 8.5 dx
3

1.0 0.0

1 4
y − 1 .0 = − x + 4 x 3 − 10 x 2 + 8.5 x
2
1 4
y = 1.0 − x + 4 x 3 − 10 x 2 + 8.5 x , es la solución exacta de la ecuación diferencial.
2
1
y ( 0.5) = 1.0 − 0.5 4 + 4 ⋅ ( 0.5) − 10 ⋅ ( 0.5) + 8.5( 0.5) = 3.21875
3 2

2
Así, como la solución verdadera es de orden cuatro, el método de Runge Kutta de cuarto
orden también da una solución exacta.
La tabla de EXCEL correspondiente se muestra a continuación

h 0.1

n xn yn k1n k2n k3n k4n solución exacta

0 0 1 8.5 7.52975 7.52975 6.618 1
1 0.1 1.75395 6.618 5.76325 5.76325 4.964 1.75395
2 0.2 2.3312 4.964 4.21875 4.21875 3.526 2.3312
3 0.3 2.75395 3.526 2.88425 2.88425 2.292 2.75395
4 0.4 3.0432 2.292 1.74775 1.74775 1.25 3.0432
5 0.5 3.21875 1.25 0.79725 0.79725 0.388 3.21875
6 0.6 3.2992 0.388 0.02075 0.02075 -0.306 3.2992
7 0.7 3.30195 -0.306 -0.59375 -0.59375 -0.844 3.30195
8 0.8 3.2432 -0.844 -1.05825 -1.05825 -1.238 3.2432
9 0.9 3.13795 -1.238 -1.38475 -1.38475 -1.5 3.13795
10 1 3 -1.5 -1.58525 -1.58525 -1.642 3

La gráfica de la aproximación en el intervalo [0.0, 1.0], se presenta aquí

Aproximación de la solución de la ecuación diferencial:

y'=-2x 3+12x 2-20x+8.5, y(0)=1.0, utilizando el método de Runge
Kutta de cuarto orden con h=0.1, en el intervalo [0.0,1.0]
3.5
3
2.5
2
y

1.5
1
0.5
0
0 0.2 0.4 0.6 0.8 1 1.2
x

EJEMPLO No. 8

Use el método de Runge-Kutta de cuarto orden para estimar y(0.5) en la ecuación diferencial
y’= 4e0.8x -0.5y, y(0.0)=2.0. Utilize un h=0.1
Obsérvese que para este caso f(x,y)= 4e0.8x -0.5y
Inciando el proceso iterativo con la condición inicial xo=0.0, yo=2.0
x1=xo+h=0.0+0.1=0.1

k1, 0 = f ( x 0 , y 0 ) = 4e 0.8 x0 − 0.5 y 0

= 4e 0.8( 0.0 ) − 0.5( 2.0) = 3.0

 1 
 1 1  0.8 x0 + h 
 1 
k 2,0 = f  x 0 + h, y 0 + k1,0 ⋅ h  = 4e  2  − 0.5 y 0 + k1,0 ⋅ h 
 2 2   2 
 1 

0.8 0.0 + ⋅0.1 1 
= 4e 
− 0.5 2.0 + 3.0 ⋅ ( 0.1) 
2 

 2 
= 3.088243097
 1 
 1 1  0.8 x0 + h 
 1 
k 3, 0 = f  x 0 + h, y 0 + k 2,0 ⋅ h  = 4e  2  − 0.5 y 0 + k 2, 0 ⋅ h 
 2 2   2 
 1 
0.8 0.0 + ⋅0.1 
 1 
= 4e  2 
− 0.5 2.0 + 3.08824397 ⋅ ( 0.1) 
 2 
= 3.08603702
k 4, 0 = f ( x 0 + h, y 0 + k 3, 0 ⋅ h ) = 4e 0.8( x0 + h ) − 0.5( y 0 + k 3,0 ⋅ h )
= 4e 0.8( 0.0+⋅0.1) − 0.5( 2.0 + 3.08603702 ⋅ ( 0.1) )
= 3.17884642

luego se aplica la formula de recursión:

y1 = y 0 + ( k1,0 + 2k 2, 0 + 2k 3, 0 + k 4,0 ) ⋅ 0.1
1
6
1
y1 = 2.0 + ( 3.0 + 2 ⋅ ( 3.088243097 ) + 2( 3.08603702) + ( 3.17884642) ) ⋅ 0.1 = 2.30879011
6
Resultado: x1=0.1, y(0.1)= 2.30879011
x2=x1+h=0.1+0.1=0.2

k1,1 = f ( x1 , y1 ) = 4e 0.8 x1 − 0.5 y1

= 4e 0.8( 0.1) − 0.5( 2.30879011) = 3.17875322

 1 
 1 1  0.8 x1 + h 
 1 
k 2,1 = f  x1 + h, y1 + k1,1 ⋅ h  = 4e  2  − 0.5 y1 + k1,1 ⋅ h 
 2 2   2 
 1 

0.8 0.1+ ⋅0.1  1 
= 4e 
− 0.5 2.30879011 + 3.17875322 ⋅ ( 0.1) 
2 

 2 
= 3.276123521
  1 
 1 1  0.8 x1 + h 
 1 
k 3,1 = f  x1 + h, y1 + k 2,1 ⋅ h  = 4e  2  − 0.5 y1 + k 2,1 
 2 2   2 
 1 
0.8 0.1+ ⋅0.1 
 1 
= 4e  2 
− 0.5 2.30879011 + 3.276123521 ⋅ ( 0.1) 
 2 
= 3.27368926
k 4,1 = f ( x 0 + h, y 0 + k 3, 0 ⋅ h ) = 4e − 0.5( y1 + k 3,1 )
0.8 ( x1 + h )

= 4e 0.8( 0.1+⋅0.1) − 0.5( 2.30879011 + 3.08603702 ⋅ ( 0.1) )

= 3.37596397

luego se aplica la formula de recursión:

y 2 = y1 + ( k1,1 + 2k 2,1 + 2k 3,1 + k 4,1 ) ⋅ 0.1
1
6
1
y 2 = 2.30879011 + ( 3.17875322 + 2 ⋅ ( 3.276123521) + 2( 3.27368926 ) + ( 3.37596397 ) ) ⋅ 0.1 = 2.63636249
6
Resultado: x2=0.2, y(0.2)= 2.63636249

x3=x2+h=0.2+0.1=0.3

k1, 2 = f ( x 2 , y 2 ) = 4e 0.8 x2 − 0.5 y 2

= 4e 0.8( 0.2 ) − 0.5( 2.63636249) = 3.37586224

 1 
 1 1  0.8 x2 + h 
 1 
k 2, 2 = f  x 2 + h, y 2 + k1, 2 ⋅ h  = 4e  2  − 0.5 y 2 + k1, 2 ⋅ h 
 2 2   2 
 1 

0.8 0.2 + ⋅0.1  1 
= 4e 
− 0.5 2.63636249 + 3.37586224 ⋅ ( 0.1) 
2 

 2 
= 3.483033232
 1 
 1 1  0.8 x2 + h 
 1 
k 3, 2 = f  x 2 + h, y 2 + k 2, 2 ⋅ h  = 4e  2  − 0.5 y 2 + k 2, 2 
 2 2   2 
 1 
0.8 0.2 + ⋅0.1 
 1 
= 4e  2 
− 0.5 2.63636249 + 3.483033232 ⋅ ( 0.1) 
 2 
= 3.48035396
k 4 , 2 = f ( x 2 + h , y 2 + k 3, 2 ⋅ h ) = 4e 0.8( x2 + h ) − 0.5( y 2 + k 3, 2 )
= 4e 0.8( 0.2 +⋅0.1) − 0.5( 2.63636249 + 3.48035396 ⋅ ( 0.1) )
= 3.59279766

luego se aplica la formula de recursión:

y 3 = y 2 + ( k1, 2 + 2k 2, 2 + 2k 3, 2 + k 4, 2 ) ⋅ 0.1
1
6
 1
y 3 = 2.63636249 + ( 3.37586224 + 2 ⋅ ( 3.483033232) + 2( 3.48035396) + ( 3.59279766) ) ⋅ 0.1 = 2.98461973
6
Resultado: x3=0.3, y(0.3)= 2.98461973

x4=x3+h=0.3+0.1=0.4

k1,3 = f ( x3 , y 3 ) = 4e 0.8 x3 − 0.5 y 3

= 4e 0.8( 0.3 ) − 0.5( 2.98461973) = 3.59268674

 1 
 1 1  0.8 x3 + h 
 1 
k 2,3 = f  x 3 + h, y 3 + k1,3 ⋅ h  = 4e  2  − 0.5 y 3 + k1,3 ⋅ h 
 2 2   2 
 1 

0.8 0.3+ ⋅0.1  1 
= 4e − 0.5 2.98461973 + 3.59268674 ⋅ ( 0.1) 
 2 

 2 
= 3.710392217
 1 
 1 1  0.8 x3 + h 
 1 
k 3,3 = f  x3 + h, y 3 + k 2,3 ⋅ h  = 4e  2  − 0.5 y 3 + k 2,3 
 2 2   2 
 1 
0.8 0.3+ ⋅0.1 
 1 
= 4e  2 
− 0.5 2.98461973 + 3.710392217 ⋅ ( 0.1) 
 2 
= 3.70744958
k 4,3 = f ( x 3 + h, y 3 + k 3,3 ⋅ h ) = 4e − 0.5( y 3 + k 3,3 )
0.8 ( x3 + h )

= 4e 0.8( 0.3+⋅0.1) − 0.5( 2.98461973 + 3.70744958 ⋅ ( 0.1) )

= 3.83082871

luego se aplica la formula de recursión:

y 4 = y 3 + ( k1,3 + 2k 2,3 + 2k 3,3 + k 4,3 ) ⋅ 0.1
1
6
1
y 4 = 2.98461973 + ( 3.59268674 + 2 ⋅ ( 3.710392217 ) + 2( 3.70744958) + ( 3.83082871) ) ⋅ 0.1 = 3.35560638
6
Resultado: x4=0.4, y(0.4)= 3.35560638
x5=x4+h=0.4+0.1=0.5

k1, 4 = f ( x 4 , y 4 ) = 4e 0.8 x4 − 0.5 y 4

= 4e 0.8( 0.4 ) − 0.5( 3.35560638) = 3.83070787
  1 
 1 1  0.8 x4 + h 
 1 
k 2, 4 = f  x 4 + h, y 4 3 + k1, 4 ⋅ h  = 4e  2  − 0.5 y 4 + k1, 4 ⋅ h 
 2 2   2 
 1 

0.8 0.4 + ⋅0.1  1 
= 4e 
− 0.5 3.35560638 + 3.83070787 ⋅ ( 0.1) 
2 

 2 
= 3.959746772
 1 
 1 1  0.8 x4 + h 
 1 
k 3, 4 = f  x 4 + h, y 4 + k 2, 4 ⋅ h  = 4e  2  − 0.5 y 4 + k 2, 4 
 2 2   2 
 1 
0.8 0.4 + ⋅0.1 
 1 
= 4e  2 
− 0.5 3.35560638 + 3.959746772 ⋅ ( 0.1) 
 2 
= 3.9565208
k 4 , 4 = f ( x 4 + h , y 4 + k 3, 4 ⋅ h ) = 4e 0.8( x4 + h ) − 0.5( y 4 + k 3, 4 )
= 4e 0.8( 0.4 +⋅0.1) − 0.5( 3.35560638 + 3.9565208 ⋅ ( 0.1) )
= 4.09166956

luego se aplica la formula de recursión:

y 5 = y 4 + ( k1, 4 + 2k 2, 4 + 2k 3, 4 + k 4, 4 ) ⋅ 0.1
1
6
1
y 5 = 3.35560638 + ( 3.83070787 + 2 ⋅ ( 3.959746772 ) + 2( 3.9565208) + ( 4.09166956 ) ) ⋅ 0.1 = 3.75152159
6
Resultado: x5=0.5, y(0.5)= 3.75152159

Solución exacta:
La ecuación diferencial
y’= 4e0.8x -0.5y, y(0.0)=2.0. Utilice un h=0.1
puede reordenarse como:
y’+0.5y=4e0.8x
dy
La cual es una ecuación diferencial con lineal, es decir de la forma + P ( x) y = Q( x)
dx
La cual se resuelve multiplicando toda la ecuación por un factor integrante µ ( x ) = e ∫ P ( x )dx
En este caso el factor integrante sería µ ( x ) = e ∫ 0.5dx = e 0.5 x . Multiplicando la ecuación
diferencial por e0.5x, resulta:
e0.5x y’+ e0.5x 0.5y=4e1.3x. La expresión del miembro de la izquierda es la derivada del producto
e0.5xy
(
d e 0.5 x y ) ( )
= 4e1.3 x ∴ d e 0.5 x y = 4e1.3 x dx , integrando, resulta:
dx
( ) 4 1.3 x
∫ d e y = ∫ 4e dx ∴ e y = 1.3 e + C ∴ y = 1.3 e + Ce
0.5 x 1.3 x 0.5 x 4 0.8 x − 0.5 x

Reemplazando la condición inicial(y(0.0)=2.0, para hallar la solución particular, queda

4 0.8( 0.0 ) x 1.4
0.2 = e + Ce −0.5( 0.0 ) ∴ C = −
1.3 1.3
Reemplazando, se obtiene la solución exacta de la ecuación diferencial, como:
 4 0.8 x 1.4 −0.5 x
y= e − e
1.3 1.3

Ahora, se reemplaza x=0.5, para obtener el valor exacto de y(0.5):

4 0.8( 0.5 ) 1.4 −0.5( 0.5 )

y ( 0 .5 ) = e − e = 3.751521303
1 .3 1 .3
Entonces el error cometido en la estimación es:
e= 3.751521303 -3.75152159=-2.84925E-07, lo que indica la gran exactitud de la aproximación
dad por el método de Runge-Kutta.

A continuación se muestra la Tabla de EXCEL correspondiente, utilizada para facilitar los

cálculos:
h 0.1

n xn yn k1n k2n k3n k4n solución exacta

0 0 2 3 3.088243097 3.08603702 3.17884642 2
1 0.1 2.30879011 3.17875322 3.276123521 3.27368926 3.37596397 2.308790059
2 0.2 2.63636249 3.37586224 3.483033232 3.48035396 3.59279766 2.636362384
3 0.3 2.98461973 3.59268674 3.710392217 3.70744958 3.83082871 2.984619565
4 0.4 3.35560638 3.83070787 3.959746772 3.9565208 4.09166956 3.355606156
5 0.5 3.75152159 4.091538 4.23277963 4.22924859 4.37707439 3.751521303
6 0.6 4.17473274 4.37693124 4.53132095 4.52746121 4.68895057 4.174732384
7 0.7 4.62779017 4.68879492 4.857360243 4.85314611 5.02937113 4.62778975
8 0.8 5.11344315 5.02920194 5.213059305 5.20846287 5.40058812 5.113442656
9 0.9 5.63465706 5.40040431 5.600766246 5.5957572 5.80504733 5.634656485
10 1 6.19463203 5.8048477 6.023030699 6.01757612 6.245404 6.194631377

Como se observa en la tabla, la solución aproximada por medio del método de Runge-Kutta,
es bastante exacta(compárese con la columna: solución exacta).
A continuación se muestra la gráfica de la aproximación en el intervalo [0.0,1.0]

Gráfica de la aproximación de la ecuación diferencial: y'=4e0.8x-0.5y,

y(0.0)=2.0, en el intervalo [0.0,1.0], mediante e método de Runge
Kutta de cuarto orden, con h=0.1
7

4
y

0
0 0.2 0.4 x 0.6 0.8 1