Unidad 5.

ECUACIONES DIFERENCIALES ORDINARIAS

5.1 Teoría elemental de los problemas a valor inicial

Sea el problema a valor inicial (P.V.I.)

𝑑𝑦
𝑦′ = = 𝑓(𝑥, 𝑦)
{ 𝑑𝑥
𝑦(𝑥0 ) = 𝑦0
Una solución explícita es una función 𝑦 = 𝜙(𝑥) (una solución implícita es 𝑔(𝑥, 𝑦) = 0) definida en un intervalo 𝐼 ∋ 𝑥0 que
satisface la ecuación diferencial
𝑦′ = 𝑓(𝑥, 𝑦) y también la condición inicial 𝑦(𝑥0 ) = 𝑦0 .

Teorema (Existencia y unicidad). Sea el P.V.I.

𝑑𝑦
𝑦′ = = 𝑓(𝑥, 𝑦)
{ 𝑑𝑥
𝑦(𝑥0 ) = 𝑦0

𝜕𝑓
Si 𝑓(𝑥, 𝑦) y existe y ambas son funciones continuas en un rectángulo que incluye a (𝑥0 , 𝑦0 ), entonces existe un intervalo
𝜕𝑦
suficientemente pequeño 𝐼 =]𝑥0 − 𝛿, 𝑥0 + 𝛿[ con 𝛿 > 0 tal que el P.V.I. tiene solución única en 𝐼.

Algunos tipos de ecuaciones diferenciales:

a) Homogéneas:
𝑑𝑦 𝑦 𝑦 𝑑𝑦 𝑑𝑧
= 𝑓 ( ) → 𝑧 = → 𝑦 = 𝑥𝑧 → = 𝑧 + 𝑥 ⟹
𝑑𝑥 𝑥 𝑥 𝑑𝑥 𝑑𝑥
𝑑𝑧 𝑑𝑧 𝑓(𝑧) − 𝑧 𝑑𝑧 𝑑𝑥
𝑧+𝑥 = 𝑓(𝑧) separable → = →∫ =∫
𝑑𝑥 𝑑𝑥 𝑥 𝑓(𝑧) − 𝑧 𝑥
Ejemplo. Sea el P.V.I. (𝑥 2 − 𝑥𝑦)𝑑𝑦 + (2𝑥𝑦)𝑑𝑥 = 0, 𝑦(1) = 2
(𝑥 2 − 𝑥𝑦)𝑑𝑦 = −2𝑥𝑦𝑑𝑥

𝑑𝑦 −2𝑥𝑦
=
𝑑𝑥 𝑥 2 − 𝑥𝑦

2𝑥𝑦
𝑑𝑦 − 2
= 𝑥
𝑑𝑥 𝑥 2 − 𝑥𝑦
𝑥2
𝑦
𝑑𝑦 −2 𝑥 𝑦
= 𝑦 = 𝑓 ( ) EDH
𝑑𝑥 1 − 𝑥
𝑥
𝑦 𝑑𝑦 𝑑𝑧
Cambio de variable: 𝑧 = → 𝑦 = 𝑥𝑧 → =𝑧+𝑥
𝑥 𝑑𝑥 𝑑𝑥

𝑑𝑧 −2𝑧 𝑑𝑧 2𝑧 −2𝑧−𝑧+𝑧 2 𝑑𝑧 −3𝑧+𝑧 2 1−𝑧 𝑑𝑥

Reemp. 𝑧 + 𝑥 = separable → 𝑥 =− −𝑧 = →𝑥 = →∫ 𝑑𝑧 = ∫
𝑑𝑥 1−𝑧 𝑑𝑥 1−𝑧 1−𝑧 𝑑𝑥 1−𝑧 −3𝑧+𝑧 2 𝑥

1 2
1−𝑧 1−𝑧 − − 1 2
3 3
La primera integral: ∫ 𝑑𝑧 = ∫ 𝑑𝑧 = ∫ + 𝑑𝑧 = − ln 𝑧 − ln(𝑧 − 3)
−3𝑧+𝑧 2 𝑧(𝑧−3) 𝑧 𝑧−3 3 3
𝑑𝑥
Segunda integral: ∫ = ln 𝑥
𝑥

1 2
⟹ − ln 𝑧 − ln(𝑧 − 3) = ln 𝑥 + 𝐶
𝑦
3 3
Reemplazar 𝑧 = :
𝑥
1 𝑦 2 𝑦
− ln − ln ( − 3) = ln 𝑥 + ln 𝐶 ,∗ (−3)
3 𝑥 3 𝑥
𝑦 𝑦
ln + 2 ln ( − 3) = −3 ln 𝐶𝑥
𝑥 𝑥
𝑦 𝑦 2
ln [ ( − 3) ] = ln 𝐶𝑥 −3
𝑥 𝑥
 𝑦 𝑦 2
( − 3) = 𝐶𝑥 −3
𝑥 𝑥
𝑦 𝑦 − 3𝑥 2 𝐶
( ) = 3
𝑥 𝑥 𝑥

𝑦(𝑦 − 3𝑥)2 = 𝐶
2
Condición inicial: 𝑦(1) = 2: 2(2 − 3(1)) = 𝐶 → 2 = 𝐶
𝑦(𝑦 − 3𝑥)2 = 2
Comprobación:
[𝑦(𝑦 − 3𝑥)2 ]𝑥 = [2]𝑥

𝑦 ′ ⋅ (𝑦 − 3𝑥)2 + 𝑦 ⋅ 2(𝑦 − 3𝑥)(𝑦 ′ − 3) = 0

𝑦 ′ (𝑦 − 3𝑥) + 2𝑦(𝑦 ′ − 3) = 0
𝑦𝑦 ′ − 3𝑥𝑦 ′ + 2𝑦𝑦 ′ − 6𝑦 = 0
3𝑦𝑦 ′ − 3𝑥𝑦 ′ − 6𝑦 = 0
𝑦𝑦 ′ − 𝑥𝑦 ′ − 2𝑦 = 0
(𝑦 − 𝑥)𝑦 ′ − 2𝑦 = 0 ∗ (−𝑥)
(𝑥 2 − 𝑥𝑦)𝑦′ + 2𝑥𝑦 = 0
𝑑𝑦
(𝑥 2 − 𝑥𝑦) + 2𝑥𝑦 = 0
𝑑𝑥
2
(𝑥 − 𝑥𝑦)𝑑𝑦 + 2𝑥𝑦𝑑𝑦 = 0 (coincide con la ED inicial)
Además se satisface 𝑦(1) = 2:
2
𝑦(𝑦 − 3𝑥)2 = 2 → 2(2 − 3(1)) = 2 → 2 = 2

b) Exactas:
𝑀(𝑥, 𝑦)𝑑𝑥 + 𝑁(𝑥, 𝑦)𝑑𝑦 = 0
Siendo 𝑀𝑦 = 𝑁𝑥
Para resolver:
Se conoce que existe 𝑧 = 𝑓(𝑥, 𝑦) tal que 𝑑𝑧 = 𝑓𝑥 𝑑𝑥 + 𝑓𝑦 𝑑𝑦 = 𝑀𝑑𝑥 + 𝑁𝑑𝑦 = 0 → 𝑧 = 𝑘
𝑓𝑥 = 𝑀 𝑓 = ∫ 𝑀𝑑𝑥
{ →{ ⟹ 𝑧 = 𝑓(𝑥, 𝑦) = 𝑘
𝑓𝑦 = 𝑁 𝑓 = ∫ 𝑁𝑑𝑦
𝑥2
Ejemplo. Resolver el P.V.I. (𝑥 2 + 𝑥𝑦)𝑑𝑥 + 𝑑𝑦 = 0, 𝑦(2) = 1
2
𝑥2
Es exacta ya que 𝑀 = 𝑥 2 + 𝑥𝑦, 𝑁 = 𝑀𝑦 = 𝑥 = 𝑁𝑥
2
Existe 𝑧 = 𝑓(𝑥, 𝑦) tal que 𝑑𝑧 = 𝑓𝑥 𝑑𝑥 + 𝑓𝑦 𝑑𝑦 = 𝑀𝑑𝑥 + 𝑁𝑑𝑦 = 0 → 𝑧 = 𝑘
𝑥 3 𝑥 2𝑦
𝑓𝑥 = 𝑀 = 𝑥 2 + 𝑥𝑦 𝑓 = ∫ 𝑥 2 + 𝑥𝑦𝑑𝑥 = + + 𝐶(𝑦)
3 2 𝑥 3 𝑥 2𝑦
{ 𝑥2 → ⟹ 𝑧 = 𝑓(𝑥, 𝑦) = + =𝑘
𝑓𝑦 = 𝑁 = 𝑥2 𝑥2𝑦 3 2
2 { 𝑓=∫ 𝑑𝑦 = + 𝐶(𝑥)
2 2
23 1 8 14
Condión inicial 𝑦(2) = 1 → + 22 ⋅ = 𝑘 → + 2 = = 𝑘
3 2 3 3
𝑥 3 𝑥 2 𝑦 14
+ =
3 2 3

2𝑥 3 + 3𝑥 2 𝑦 = 28
28 − 2𝑥 3
𝑦=
3𝑥 2
c) Lineales.
𝑑𝑦
+ 𝑃(𝑥)𝑦 = 𝑄(𝑥)
𝑑𝑥
Para resolver: factor integrante 𝜇 = 𝑒 ∫ 𝑃(𝑥)𝑑𝑥
1
𝑦= [∫ 𝜇𝑄(𝑥)𝑑𝑥 + 𝐶]
𝜇
𝑑𝑦
Ejemplo. + cos 𝑥 𝑦 = 2cos 𝑥
𝑑𝑥
𝜇 = 𝑒 ∫ 𝑃(𝑥)𝑑𝑥 = 𝑒 ∫ cos 𝑥𝑑𝑥=sin 𝑥 = 𝑒 sin 𝑥
1
𝑦 = [∫ 𝜇𝑄(𝑥)𝑑𝑥 + 𝐶]
𝜇
𝑦 = 𝑒 −sin 𝑥 [∫ 𝑒 sin 𝑥 2cos 𝑥 𝑑𝑥 + 𝐶 ]
𝑦 = 𝑒 −sin 𝑥 [𝑒 sin 𝑥 2 + 𝐶]
𝑦 = 2 + 𝐶𝑒 − sin 𝑥
d) Separables.
𝑑𝑦
= 𝑓(𝑥)𝑔(𝑦)
𝑑𝑥
1
⟹∫ 𝑑𝑥 = ∫ 𝑔(𝑦)𝑑𝑦
𝑓(𝑥)
𝑑𝑦 𝑥𝑦+2𝑦−𝑥−2
Ejemplo.
𝑑𝑥
= 𝑥𝑦−3𝑦+𝑥−3 ; 𝑦(4) = 2

𝑑𝑦 (𝑥 + 2)𝑦 − (𝑥 + 2)
=
𝑑𝑥 (𝑥 − 3)𝑦 + 𝑥 − 3
𝑑𝑦 (𝑥+2)(𝑦−1)
𝑑𝑥
= (𝑥−3)(𝑦+1) ED separable

𝑦+1 𝑥+2
∫ 𝑑𝑦 = ∫ 𝑑𝑥
𝑦−1 𝑥−3

𝑦+1 𝑦−1
−𝑦 + 1 1
(2)

𝑥+2 𝑥−3
−𝑥 + 3 1
(5)

2 5
∫1 + 𝑑𝑦 = ∫ 1 + 𝑑𝑥
𝑦−1 𝑥−3
𝑦 + 2 ln(𝑦 − 1) = 𝑥 + 5 ln(𝑥 − 3) + 𝐶
𝑦 − 𝑥 − 𝐶 = 5 ln(𝑥 − 3) − 2 ln(𝑦 − 1)
𝑦 − 𝑥 − 𝐶 = ln(𝑥 − 3)5 − ln(𝑦 − 1)2
(𝑥 − 3)5
𝑦 − 𝑥 − 𝐶 = ln
(𝑦 − 1)2
(4−3)5
Condición inicial 𝑦(4) = 2: 2 − 4 − 𝐶 = ln (2−1)2 → 𝐶 = −2
(𝑥 − 3)5
𝑦 − 𝑥 + 2 = ln
(𝑦 − 1)2

e) Bernoulli.

𝑑𝑦
+ 𝑃(𝑥)𝑦 = 𝑄(𝑥)𝑦 𝑛 , 𝑛 ∈ ℝ − {−1,0}
𝑑𝑥
𝑑𝑧 𝑑𝑦
Para resolver plantear 𝑧 = 𝑦1−𝑛 → = (1 − 𝑛)𝑦 −𝑛 ⋅
𝑑𝑥 𝑑𝑥
Luego multiplicar la ecuación diferencial por (1 − 𝑛)𝑦 −𝑛 para obtener:

𝑑𝑦
(1 − 𝑛)𝑦 −𝑛 1−𝑛 ( )
+ (1 − 𝑛) 𝑦⏟ 𝑃 𝑥 = (1 − 𝑛)𝑄(𝑥)
⏟ 𝑑𝑥 𝑧
𝑑𝑧
𝑑𝑥
𝑑𝑧
+ (1 − 𝑛)𝑃(𝑥)𝑧 = (1 − 𝑛)𝑄(𝑥)
𝑑𝑥
Esta última ecuación es de tipo lineal que ya se sabe resolver.
 1 𝑑𝑦
Ejemplo. 𝑦(1) = ; 𝑥 − (1 + 𝑥)𝑦 = 𝑥𝑦 2 ÷𝑥
2 𝑑𝑥
𝑛
𝑑𝑦 1+𝑥
−
𝑑𝑥 ⏟ 𝑥
𝑦 1 𝑦 ⏞2 ED de Bernoulli
= ⏟
𝑄(𝑥)
𝑃(𝑥)
Sea 𝑧 = 𝑦1−𝑛 → 𝑧 = 𝑦1−2 → 𝑧 = 𝑦 −1 Integrando:
𝑑𝑧 𝑑𝑦
= −1𝑦 −2 ⋅ 𝑥 𝑒 𝑥 𝑧 = − ∫ 𝑥𝑒 𝑥 𝑑𝑥
𝑑𝑥 𝑑𝑥
Multiplicando la ED por −1𝑦 −2:
𝑥 𝑒 𝑥 𝑧 = − [𝑥𝑒 𝑥 − ∫ 1𝑒 𝑥 𝑑𝑥]
𝑑𝑦 1 + 𝑥 −1
−1𝑦 −2 + 𝑦⏟ = −1 𝑥𝑒 𝑥 𝑧 = −𝑥𝑒 𝑥 + 𝑒 𝑥 + 𝐶
⏟ 𝑑𝑥 𝑥
𝑧
𝑑𝑧/𝑑𝑥 Devolviendo 𝑧 = 𝑦 −1:
𝑑𝑧 1+𝑥 1
𝑑𝑥
+ 𝑥
𝑧 = −1 ED Lineal (EDL) 𝑥𝑒 𝑥 𝑦 −1 = −𝑥𝑒 𝑥 + 𝑒 𝑥 + 𝐶 ⋅ 𝑒 −𝑥
1+𝑥
𝑑𝑥 𝑥
Factor integrante: 𝜇 = 𝑒 ∫ 𝑥 = 1 𝐶𝑒 −𝑥
1
∫𝑥+1𝑑𝑥
𝑦 −1 = −1 + +
𝑒 = 𝑒 ln 𝑥+𝑥 = 𝑒 ln 𝑥 𝑒 𝑥 = 𝑥𝑒 𝑥 , 𝑥 > 𝑥 𝑥
0 𝑦 −1
Multiplicando la EDL por 𝜇: −𝑥 + 1 + 𝐶𝑒 −𝑥
= elevando a la (−1)
𝑑𝑧 1+𝑥 𝑥
𝑥𝑒 𝑥 + 𝑥𝑒 𝑥 𝑧 = −𝑥𝑒 𝑥 𝑥
𝑑𝑥 𝑥 𝑦=
𝑑𝑧 𝑑(𝑥 𝑒 𝑥 ) 𝑥 + 1 + 𝐶𝑒 −𝑥
𝑥𝑒 𝑥 + ⋅ 𝑧 = −𝑥𝑒 𝑥 La condición inicial es:
𝑑𝑥 𝑑𝑥 1 1
𝑦(1) = : =
1 1
→ =
1
→𝐶=0
2 2 1+1+𝐶𝑒−1 2 2+𝐶𝑒−1
𝑥 𝑥
𝑑(𝑥 𝑒 𝑧) 𝑦=
= −𝑥𝑒 𝑥 𝑥+1
𝑑𝑥
5.2 Método de Euler y cálculo del error de truncamiento
Dado el problema a valor inicial
𝑑𝑦
𝑦′ = = 𝑓(𝑥, 𝑦)
{ 𝑑𝑥
𝑦(𝑥0 ) = 𝑦0

Métodos de Runge Kutta:

Para estimar 𝑦(𝑥𝑖+1 ) se usa la regla

𝑦(𝑥𝑖+1 ) ≈ 𝑦𝑖+1 = 𝑦𝑖 + 𝜙ℎ
Donde 𝜙 es la pendiente de una línea que pasa por
(𝑥𝑖 , 𝑦𝑖 )

𝑑𝑦
= 𝑓(𝑥, 𝑦); 𝑦(𝑥0 ) = 𝑦0
𝑑𝑥
Método de Euler:

En el método de Runge-Kutta se usa como pendiente la

derivada 𝜙 = 𝑓(𝑥𝑖 , 𝑦𝑖 )
𝑦(𝑥𝑖+1 ) ≈ 𝑦𝑖+1 = 𝑦𝑖 + 𝑓(𝑥𝑖 , 𝑦𝑖 )ℎ
El error de truncamiento en un paso es 𝑂(ℎ2 ), i.e.
proporcional a ℎ2

Ejemplo Método de Euler aplicado al problema

𝑑𝑦
−2𝑥 3 + 12𝑥 2 − 20𝑥 + 8.5 ; 𝑦(0) = 1
=⏟
𝑑𝑥 𝑓(𝑥,𝑦)

La solución exacta es

𝑦 = −0.5𝑥⁴ + 4𝑥³ − 10𝑥² + 8.5𝑥 + 1

ℎ = 0.5; 𝑥0 = 0; 𝑦0 = 1
𝑦1 = 𝑦0 + ℎ𝑓(𝑥0 , 𝑦0 ) = 1 + 0.5𝑓(0,1) = 5.2500
𝑥1 = 𝑥0 + ℎ = 0 + 0.5 = 0.5
𝑦2 = 𝑦1 + ℎ𝑓(𝑥1 , 𝑦1 ) = 5.25 + 0.5𝑓(0.5,5.25) = 5.875
Algoritmo 1 del método de Euler

Algoritmo del método de Euler

Datos: una función 𝑓(𝑥, 𝑦) continua con 𝑓𝑦 también continua sobre un rectángulo 𝑅, un punto (𝑥0 , 𝑦0 ) en el
interior de 𝑅, y el número de iteraciones adicionales 𝑛, y el tamaño de paso ℎ
Salida: Vector de las abscisas 𝑋 = (𝑋0 , 𝑋1 , 𝑋2 , … , 𝑋𝑛 ) y vector de las soluciones aproximadas 𝑌 = (𝑌0 , 𝑌1 , 𝑌2 , … , 𝑌𝑛 )
1. Hacer 𝑋0 ≔ 𝑥0 y 𝑌0 ≔ 𝑦0
2. Hacer para 𝑖 = 0, 1, … , 𝑛 − 1:
𝑋𝑖+1 ≔ 𝑋𝑖 + ℎ
𝑌𝑖+1 ≔ 𝑌𝑖 + ℎ𝑓(𝑋𝑖 , 𝑌𝑖 )

Implementación del algoritmo 1 de Euler en lenguaje Python con datos del ejemplo siguiente

#euler.py
#Algoritmo 1 del método de Euler

#Datos
#Defina la función f(x,y) de y'=f(x,y) y las condiciones dentro del programa:
def f(x,y):
return (1-y**3)/(x*y**2)

x0 = 1
y0 = 2
h = 0.1
n=3
X = [] #lista de abscisas
Y = [] #lista de ordenadas

X.append(x0)
Y.append(y0)

for i in range(n): #i=0,1,...,n-1

X.append(X[i] + h)
Y.append(Y[i] + h * f(X[i],Y[i]))
print('Abscisas',X)
print('Ordenadas',Y)
Implementación del algoritmo 1 en Matlab, incluidos errores relativos y gráficas

%Método de Euler
%Datos: una función f(x,y) continua con fy también continua
% sobre un rectángulo R, un punto (x0,y0 ) en el interior
% de R, y el número de iteraciones adicionales n, y el
% tamaño de paso h
%
syms x y y(x)
disp("Problema a valor inicial y' = f(x,y); y(x0) = y0")
disp('Ingrese los datos:')
f = input('f(x,y) = ');
x0 = input('x0 = ');
y0 = input('y0 = ');
h = input('h = ');
n = input('iter = ');

%Calculo con Euler

X = zeros(n+1,1);
Y = zeros(n+1,1);
X(1) = x0;
Y(1) = y0;
for i = 1 : n
X(i+1) = X(i) + h;
Y(i+1) = Y(i) + h * subs(f, {x,y}, {X(i),Y(i)});
end

%Solucion verdadera
yt(x) = dsolve(diff(y) == f, y(x0) == y0);
Yt = yt(X);

%Error relativo verdadero porcentual

et = double((Yt - Y) ./ Yt *100);

%Graficas
close all
plot(X,Y,'r') %euler
hold on
XX = X(1):0.01:X(n+1);
plot(XX,yt(XX),'b') %grafica curva exacta
grid on
legend('Euler', 'Exacta')

%Mostrar resultados
fprintf('Solucion verdadera y(x) = %s \n', char(simplify(yt)))
disp('X(i) = ')
fprintf('%.7f\n', X)
disp('Solucion verdadera, Y(X) = ')
fprintf('%.7f\n', Yt)
disp('Solucion aproximada, Yeuler(X) =')
fprintf('%.7f\n', Y)
disp('Error relativo verdadero procentual, et =')
fprintf('%.7f\n',et)
 1
Ejemplo. Resolver el P.V.I. 𝑥𝑦 ′ + 𝑦 = ; 𝑦(1) = 2
𝑦2
(a) Con el método analítico
(b) Con el método de Euler, ℎ = 0.1, 3 iteraciones adicionales, calcular los errores relativos porcentuales
Solución
1
(a) 𝑥𝑦 ′ + 𝑦 = 2 ; 𝑦(1) = 2
𝑦

La EDO es de tipo 𝑦 ′ + 𝑃(𝑥)𝑦 = 𝑄(𝑥)𝑦 𝑛 Bernoulli Factor integrante:

1 3
𝜇 = 𝑒 ∫𝑥3𝑑𝑥 = 𝑒 ln 𝑥 = 𝑥 3
1 1 1 1
𝑦 ′ + 𝑦 = 𝑦 −2 , 𝑃(𝑥) = , 𝑄(𝑥) = , 𝑛 = −2 1
𝑥 𝑥 𝑥 𝑥 𝑧 = [∫ 𝜇𝑄1 (𝑥)𝑑𝑥 + 𝐶]
𝜇
1 3
Sea 𝑧 = 𝑦1−𝑛 = 𝑦 1−(−2) = 𝑦 3 → 𝑧 = 𝑦 3 𝑧 = 3 [∫ 𝑥 3 𝑑𝑥 + 𝐶]
𝑥 𝑥
1
→ derivar respecto de x: 𝑧 ′ = 3𝑦 2 𝑦′ 𝑧 = 3 [∫ 3𝑥 2 𝑑𝑥 + 𝐶]
𝑥
1
Dividir la ED por 3𝑦 2 : 𝑧 = 3 [𝑥 3 + 𝐶]
1 3 𝑥
3𝑦 2 𝑦 ′ + 3𝑦 3 = 𝐶
𝑥 𝑥 𝑧 =1+ 3
Reemplazamos: 𝑥
Regresando 𝑧 = 𝑦 3 :
1 3 1 3 𝐶
𝑧 ′ + 3𝑧 = , 𝑃1 (𝑥) = 3, 𝑄1 (𝑥) = 𝑦3 = 1 + 3
𝑥 𝑥 𝑥 𝑥 𝑥

Esta última es de tipo Lineal Condición inicial: 𝑦(1) = 2:

𝐶
23 = 1 + 3 → 𝐶 = 7
1
7 𝑥3 + 7
𝑦3 = 1 + 3 =
𝑥 𝑥3
3
√𝑥 3 + 7
𝑦(𝑥) =
𝑥

1
(b) 𝑥𝑦 ′ + 𝑦 = ; 𝑦(1) = 2
𝑦2
1 1 − 𝑦3 1 − 𝑦3
𝑥𝑦 ′ = − 𝑦 = → 𝑦 ′
= = 𝑓(𝑥, 𝑦)
𝑦2 𝑦2 𝑥𝑦 2
1 − 𝑦3
𝑦 ′ = 𝑓(𝑥, 𝑦), 𝑓(𝑥, 𝑦) =
𝑥𝑦 2
Método de Euler:
ℎ = 0.1
𝑥𝑖+1 = 𝑥𝑖 + ℎ; 𝑦(𝑥𝑖+1 ) ≈ 𝑦𝑖+1 = 𝑦𝑖 + 𝑓(𝑥𝑖 , 𝑦𝑖 )ℎ
𝑥0 = 1, 𝑦0 = 2

𝑖 = 0: 𝑥1 = 𝑥0 + ℎ = 1 + 0.1 = 1.1;

𝑦1 = 𝑦0 + 𝑓(𝑥0 , 𝑦0 )ℎ = 2 + 𝑓(1, 2)(0.1) = 1.8250

𝑖 = 1: 𝑥2 = 𝑥1 + ℎ = 1.1 + 0.1 = 1.2;
𝑦2 = 𝑦1 + 𝑓(𝑥1 , 𝑦1 )ℎ = 1.8250 + 𝑓(1.1,1.8250 )(0.1) = 1.6864
𝑖 = 2: 𝑥3 = 𝑥2 + ℎ = 1.2 + 0.1 = 1.3
𝑦3 = 𝑦2 + 𝑓(𝑥2 , 𝑦2 )ℎ = 1.6864 + 𝑓(1.2, 1.6864)(0.1) = 1.5752
Resumen:
𝒙𝒊 𝒚𝒊 sol. aprox. con Euler 𝒚(𝒙𝒊 ) sol. exacta 𝒆𝒕 %
3 3
√𝑥 + 7
𝑦(𝑥) =
𝑥
1 2 2 0%
1.1 1.8250 3
√1.13 + 7 1.8429 − 1.8250
= 1.8429 ∗ 100% = 0.9713%
1.1 1.8429

1.2 1.6864 1.7158 1.7158 − 1.6864

∗ 100% = 1.7135%
1.7158
1.3 1.5752 1.6116 1.6116 − 1.5752
∗ 100% = 2.2586%
1.6116
Método de Heun

Método de Heun:

Primero se estima 𝑦𝑖+1 usando Euler:

0
𝑦𝑖+1 = 𝑦𝑖 + ℎ𝑓(𝑥𝑖 , 𝑦𝑖 )
Luego la pendiente de la curva en (𝑥𝑖+1 , 𝑦(𝑥𝑖+1 )) se
estima con
0
𝑓(𝑥𝑖+1 , 𝑦𝑖+1 )

En el método de Runge-Kutta se usa como pendiente el

promedio entre la pendiente en (𝑥𝑖 , 𝑦𝑖 ) y la estimada en
(𝑥𝑖+1 , 𝑦(𝑥𝑖+1 )):
0
𝑓(𝑥𝑖 , 𝑦𝑖 ) + 𝑓(𝑥𝑖+1 , 𝑦𝑖+1 )
𝜙=
2
Luego la aproximación final es:
0
𝑓(𝑥𝑖 , 𝑦𝑖 ) + 𝑓(𝑥𝑖+1 , 𝑦𝑖+1 )
𝑦(𝑥𝑖+1 ) ≈ 𝑦𝑖+1 = 𝑦𝑖 + ℎ
2
4
El error de truncamiento en un paso es 𝑂(ℎ ), i.e.
proporcional a ℎ4

(Runge Kutta: 𝑦𝑖+1 = 𝑦𝑖 + 𝜙ℎ)

Resumen:
0
Predictor: 𝑦𝑖+1 = 𝑦𝑖 + ℎ𝑓(𝑥𝑖 , 𝑦𝑖 )
ℎ 0
Corrector: 𝑦𝑖+1 = 𝑦𝑖 + [𝑓(𝑥𝑖 , 𝑦𝑖 ) + 𝑓(𝑥𝑖+1 , 𝑦𝑖+1 )]
2

Algoritmo 1 del método de Heun

Algoritmo 1 del método de Heun

Datos: una función 𝑓(𝑥, 𝑦) continua con 𝑓𝑦 también continua sobre un rectángulo 𝑅, un punto (𝑥0 , 𝑦0 ) en el
interior de 𝑅, y el número de iteraciones adicionales 𝑛, y el tamaño de paso ℎ
Salida: Vector de las abscisas 𝑋 = (𝑋0 , 𝑋1 , 𝑋2 , … , 𝑋𝑛 ) y vector de las soluciones aproximadas 𝑌 = (𝑌0 , 𝑌1 , 𝑌2 , … , 𝑌𝑛 )
1. Hacer 𝑋0 ≔ 𝑥0 y 𝑌0 ≔ 𝑦0
2. Hacer para 𝑖 = 0, 1, … , 𝑛 − 1:
𝑋𝑖+1 ≔ 𝑋𝑖 + ℎ
0
𝑌𝑖+1 ≔ 𝑌𝑖 + ℎ𝑓(𝑋𝑖 , 𝑌𝑖 )
ℎ 0
𝑌𝑖+1 = 𝑌𝑖 + [𝑓(𝑋𝑖 , 𝑌𝑖 ) + 𝑓(𝑋𝑖+1 , 𝑌𝑖+1 )]
2
Implementación del algoritmo 1 de Heun en Matlab

%Método de Heun
%Datos: una función f(x,y) continua con fy también continua
% sobre un rectángulo R, un punto (x0,y0 ) en el interior
% de R, y el número de iteraciones adicionales n, y el
% tamaño de paso h
%
syms x y y(x)
f = input('f(x,y) = ');
x0 = input('x0 = ');
y0 = input('y0 = ');
h = input('h = ');
n = input('n = ');

X = zeros(n+1,1);
Y = zeros(n+1,1);
X(1) = x0;
Y(1) = y0;
hm = h /2;
for i = 1 : n
X(i+1) = X(i) + h;
yi0 = Y(i) + h * double(subs(f,{x,y},{X(i),Y(i)})); %predictor
Y(i+1) = Y(i) + hm * (double(subs(f,{x,y},{X(i),Y(i)})) + double(subs(f,{x,y},{X(i+1),yi0}))); %corrector
end

%Solucion verdadera
yt(x) = dsolve(diff(y) == f, y(x0) == y0);
Yt = double(yt(X));

%Error verdadero porcentual

et = double((Yt - Y) ./ Yt *100);

%Graficas
close all
plot(X,Y,'r') %Heun
hold on
XX = X(1):0.01:X(n+1);
plot(XX, double(yt(XX)),'b') %grafica curva exacta
legend('Heun', 'Verdadera')
grid on

%Mostrar resultados
fprintf('Solucion verdadera y(x) = %s \n', char(simplify(yt)))
disp('X(i) = ')
fprintf('%.7f\n', X)
disp('Solucion verdadera, Y(X) = ')
fprintf('%.7f\n', Yt)
disp('Solucion aproximada, Yheun(X) =')
fprintf('%.7f\n', Y)
disp('Error relativo verdadero procentual, et =')
fprintf('%.7f\n',et)
Implementación del algoritmo 1 de Heun en lenguaje Python con datos del ejemplo siguiente

#Algoritmo 1 del método de Heun

#Datos
#Defina la función f(x,y) de y'=f(x,y) y las condiciones dentro del programa:
def f(x,y):
return ((2*y+3)/(4*x+5))**2

x0 = 1
y0 = 2
h = 0.2
n=3
X = [] #lista de abscisas
Y = [] #lista de ordenadas

X.append(x0)
Y.append(y0)

hm = h / 2;

for i in range(n): #i=0,1,...,n-1

X.append(X[i] + h)
yi0 = Y[i] + h * f(X[i],Y[i]) #predictor
Y.append(Y[i] + hm *(f(X[i],Y[i]) + f(X[i+1],yi0))) #corrector
print('Abscisas',X)
print('Ordenadas',Y)

Implementación del algoritmo 1 de Heun en lenguaje Java con datos del ejemplo siguiente
public class Edos {

/**
* @param args the command line arguments
*/
public static void main(String[] args) {
//Algoritmo 1 del método de Heun
//Datos
double x0 = 1;
double y0 = 2;
double h = 0.2;
int n = 3;

double[] X = new double[n+1]; //arreglo de abscisas

double[] Y = new double[n+1]; //arreglo de ordenadas
X[0] = x0;
Y[0] = y0;
double yi0;
double hm = h / 2;
for(int i = 0; i < n; i++){ //i=0,1,...,n-1
X[i+1] = X[i] + h;
yi0 = Y[i] + h * f(X[i],Y[i]); //predictor
Y[i+1] = Y[i] + hm *(f(X[i],Y[i]) + f(X[i+1],yi0)); //corrector
}
System.out.println("Abscisas");
for(int i = 0; i <= n; i++){
System.out.println(i + " " + X[i]);
}
System.out.println("Ordenadas");
for(int i = 0; i <= n; i++){
System.out.println(i + " " + Y[i]);
}
}

//Defina la función f(x,y) de y'=f(x,y):

static double f(double x, double y){
return Math.pow((2*y+3)/(4*x+5),2);
}
}
 2𝑦+3 2
Ejemplo. Resolver 𝑦′ = ( ) ; 𝑦(1) = 2
4𝑥+5
(a) Con el método analítico
(b) Con el método de Heun, ℎ = 0.2, 3 iteraciones, calcular los errores relativos porcentuales

Solución. (a)

2𝑦 + 3 2 Reemplazando la condición inicial 𝑦(1) = 2:

𝑦′ = ( ) 2 1
4𝑥 + 5 = +𝐶
2
𝑑𝑦 (2𝑦+3)
→ = (4𝑥+5)2 es separable: 2(2) + 3 4(1) + 5
𝑑𝑥 2 1 11
1 1 = +𝐶 →𝐶 =
∫ 𝑑𝑦 = ∫ 𝑑𝑥 7 9 63
(2𝑦 + 3)2 (4𝑥 + 5)2
Sustitución: 2 1 11
𝑑𝑢 ⟹ = +
𝑢 = 2𝑦 + 3 → 𝑑𝑢 = 2𝑑𝑦 → 𝑑𝑦 = 2𝑦 + 3 4𝑥 + 5 63
2
𝑑𝑣 2 118 + 44𝑥 1 59 + 22𝑥
𝑣 = 4𝑥 + 5 → 𝑑𝑣 = 4 𝑑𝑥 → 𝑑𝑥 = = → =
4 2𝑦 + 3 63(4𝑥 + 5) 2𝑦 + 3 252𝑥 + 315
1 𝑑𝑢 1 𝑑𝑣 𝑢−2 𝑣 −2
∫ 2 =∫ 2 → ∫ 𝑑𝑢 = ∫ 𝑑𝑣
𝑢 2 𝑣 4 2 4 252𝑥 + 315 252𝑥 + 315
2𝑦 + 3 = → 2𝑦 = −3
59 + 22𝑥 59 + 22𝑥
𝑢−1 𝑣 −1 2 1
− =− +𝐶 → = +𝐶
2 4 𝑢 𝑣 186𝑥 + 138 93𝑥 + 69
2 1 2𝑦 = →𝑦=
→ = +𝐶 59 + 22𝑥 59 + 22𝑥
2𝑦 + 3 4𝑥 + 5

2𝑦+3 2
(b) 𝑦′ = ( ) ; 𝑦(1) = 2
4𝑥+5
Heun ℎ = 0.2, 3 iteraciones
𝑥𝑖+1 = 𝑥𝑖 + ℎ
0
Predictor: 𝑦𝑖+1 = 𝑦𝑖 + ℎ𝑓(𝑥𝑖 , 𝑦𝑖 )
ℎ 0
Corrector: 𝑦𝑖+1 = 𝑦𝑖 + [𝑓(𝑥𝑖 , 𝑦𝑖 ) + 𝑓(𝑥𝑖+1 , 𝑦𝑖+1 )]
2

𝑥0 = 1, 𝑦0 = 2
Iteración 1: 𝒊 = 𝟎:
𝑥1 = 𝑥0 + ℎ = 1 + 0.2 = 1.2
Predictor: 𝑦10 = 𝑦0 + ℎ𝑓(𝑥0, 𝑦0 ) = 2 + 0.2𝑓(1, 2) = 2.1210
ℎ
Corrector: 𝑦1 = 𝑦0 + [𝑓(𝑥0 , 𝑦0 ) + 𝑓(𝑥1 , 𝑦10 )]
2
0.2
𝑦1 = 2 + [𝑓(1, 2) + 𝑓(1.2, 2.1210)] = 2.1151
2
Iteración 2: 𝒊 = 𝟏
𝑥2 = 𝑥1 + ℎ = 1.2 + 0.2 = 1.4
Predictor: 𝑦20 = 𝑦1 + ℎ𝑓(𝑥1, 𝑦1 ) = 2.1151 + 0.2𝑓(1.2, 2.1151)
𝑦20 = 2.2240
ℎ
Corrector: 𝑦2 = 𝑦1 + [𝑓(𝑥1 , 𝑦1 ) + 𝑓(𝑥2 , 𝑦20 )]
2
0.2
𝑦2 = 2.1151 + [𝑓(1.2, 2.1151) + 𝑓(1.4, 2.2240)]
2
𝑦2 = 2.2189
Iteración 3: 𝒊 = 𝟐
𝑥3 = 𝑥2 + ℎ = 1.4 + 0.2 = 1.6

Predictor: 𝑦30 = 𝑦2 + ℎ𝑓(𝑥2, 𝑦2 ) = 2.3174

ℎ
Corrector: 𝑦3 = 𝑦2 + [𝑓(𝑥2 , 𝑦2 ) + 𝑓(𝑥3 , 𝑦30 )] = 2.3130
2

(a) Resumen en la tabla

𝑥𝑖 𝑦𝑖 Heun 𝑦𝑡 𝑖 verdadera 𝑒𝑡 %
93𝑥 + 69 𝑦𝑡 − 𝑦
𝑦𝑡 (𝑥) = ∗ 100%
59 + 22𝑥 𝑦𝑡
1 2 𝑦𝑡(1) = 2 0%
1.2 2.1151 𝑦𝑡(1.2) = 2.1148 2.1148 − 2.1151
∗ 100% = −0.0142%
2.1148
1.4 2.2189 𝑦𝑡(1.4) = 2.2183 −0.0289%
1.6 2.3130 𝑦𝑡(1.6) = 2.312 −0.0384%
5.4.2 Método de Runge Kutta de cuarto orden

La idea es conseguir las pendientes aproximadas de la

𝑘1 = 𝑓(𝑥𝑖 , 𝑦𝑖 ) curva en algunos puntos del intervalo [𝑥𝑖 , 𝑥𝑖+1 ]
ℎ ℎ
𝑘2 = 𝑓 (𝑥𝑖 + , 𝑦𝑖 + 𝑘1 )
2 2
ℎ ℎ
𝑘3 = 𝑓 (𝑥𝑖 + , 𝑦𝑖 + 𝑘2 )
2 2
𝑘4 = 𝑓(𝑥𝑖+1 , 𝑦𝑖 + ℎ𝑘3 )

ℎ
𝑦𝑖+1 = 𝑦𝑖 + (𝑘1 + 2(𝑘2 + 𝑘3 ) + 𝑘4 )
6

El error de truncamiento en un paso es 𝑂(ℎ5 ), i.e.

proporcional a ℎ5

De esta forma la pendiente final es

1 1 1 1
𝜙 = 𝑘1 + 𝑘2 + 𝑘3 + 𝑘4
6 3 3 6
𝑦𝑖+1 = 𝑦𝑖 + ℎ𝜙

Algoritmo 1 del método de Runge Kutta de cuarto orden

Algoritmo 1 del método de Runge Kutta de cuarto orden

Datos: una función 𝑓(𝑥, 𝑦) continua con 𝑓𝑦 también continua sobre un rectángulo 𝑅, un punto (𝑥0 , 𝑦0 ) en el
interior de 𝑅, y el número de iteraciones adicionales 𝑛, y el tamaño de paso ℎ
Salida: Vector de las abscisas 𝑋 = (𝑋0 , 𝑋1 , 𝑋2 , … , 𝑋𝑛 ) y vector de las soluciones aproximadas 𝑌 = (𝑌0 , 𝑌1 , 𝑌2 , … , 𝑌𝑛 )
1. Hacer 𝑋0 ≔ 𝑥0 y 𝑌0 ≔ 𝑦0
2. Hacer para 𝑖 = 0, 1, … , 𝑛 − 1:
𝑋𝑖+1 ≔ 𝑋𝑖 + ℎ
𝑘1 ≔ 𝑓(𝑋𝑖 , 𝑌𝑖 )
ℎ ℎ
𝑘2 ≔ 𝑓 (𝑋𝑖 + , 𝑌𝑖 + 𝑘1 )
2 2
ℎ ℎ
𝑘3 ≔ 𝑓 (𝑋𝑖 + , 𝑌𝑖 + 𝑘2 )
2 2
𝑘4 ≔ 𝑓(𝑋𝑖+1 , 𝑌𝑖 + ℎ𝑘3 )
ℎ
𝑌𝑖+1 ≔ 𝑌𝑖 + (𝑘1 + 2(𝑘2 + 𝑘3 ) + 𝑘4 )
6
Implementación del algoritmo 1 de Runge Kutta de cuarto orden en lenguaje C++ con datos del ejemplo siguiente

#include <iostream>
#include <math.h>

using namespace std;

//Función dato
double f(double x,double y){
return pow((2*y+3)/(4*x+5),2);
}
int main()
{
// Algoritmo 1 del método de Runge Kutta de cuarto orden
// Datos:
double x0 = 1;
double y0 = 2;
double h = 0.2;
int n = 3;

double X[n+1], Y[n+1], k1, k2, k3, k4;

X[0] = x0;
Y[0] = y0;
double hm = h / 2;
double hs = h / 6;
for(int i = 0; i < n; i++){
X[i+1] = X[i] + h;
k1 = f(X[i],Y[i]);
k2 = f(X[i]+hm,Y[i]+hm*k1);
k3 = f(X[i]+hm,Y[i]+hm*k2);
k4 = f(X[i+1],Y[i]+h*k3);
Y[i+1] = Y[i]+hs*(k1+2*(k2+k3)+k4);
}
for(int i = 0; i <= n; i++)
cout << X[i] << " " << Y[i] << endl;
return 0;

}
Implementación del algoritmo 1 de Runge Kutta de cuarto orden en Matlab

%runge Kutta de cuarto orden

%Datos: una función f(x,y) continua con fy también continua
% sobre un rectángulo R, un punto (x0,y0 ) en el interior
% de R, y el número de iteraciones adicionales n, y el
% tamaño de paso h
% syms x y y(x)
f = input('f(x,y) = ');
x0 = input('x0 = ');
y0 = input('y0 = ');
h = input('h = ');
n = input('iter = ');

%Calculo con Runge Kutta de cuarto orden

X = zeros(n+1,1);
Y = zeros(n+1,1);
X(1) = x0;
Y(1) = y0;
hm = h / 2;
hs = h / 6;
for i = 1 : n
X(i+1) = X(i) + h;
k1 = double(subs(f, {x, y} , {X(i), Y(i)})); %k1=f(x0,y0)
k2 = double(subs(f, {x, y}, {X(i) + hm, Y(i) + hm * k1})); %k2=f(x0+h/2,y0+h/2*k1)
k3 = double(subs(f, {x, y}, {X(i) + hm, Y(i) + hm * k2})); %k3=f(x0+h/2,y0+h/2*k2)
k4 = double(subs(f, {x, y}, {X(i+1), Y(i) + h * k3})); %k4=f(x0+h,y0+h*k3)
Y(i+1) = double(Y(i) + hs * (k1 + 2 * (k2 + k3) + k4)); %y1=y0+h/6*(k1+2(k2+k3)+k4)
end

%Solucion verdadera
yt(x) = dsolve(diff(y) == f, y(x0) == y0);
Yt = double(yt(X));

%Error verdadero porcentual

et = (Yt - Y) ./ Yt * 100;

%Graficas
close all
plot(X,Y,'r')
hold on
XX = X(1):0.01:X(n+1);
plot(XX,double(yt(XX)),'b')
legend('Runge Kutta 4', 'Verdadera')
grid on

%Mostrar resultados
fprintf('Solucion verdadera y(x) = %s \n', char(simplify(yt)))
disp('X(i) = ')
fprintf('%.7f\n', X)
disp('Solucion verdadera, Y(X) = ')
fprintf('%.7f\n', Yt)
disp('Solucion aproximada, Yrk4(X) =')
fprintf('%.7f\n', Y)
disp('Error relativo verdadero procentual, et =')
fprintf('%.7f\n',et)
Implementación del algoritmo 1 de Runge Kutta de cuarto orden en lenguaje Haskell para el ejemplo siguiente

rk4 x0 y0 h (-1) = []
rk4 x0 y0 h n = (x0,y0):(rk4 x1 y1 h (n-1))
where
hm = h / 2
hs = h / 6
x1 = x0 + h
k1 = f(x0,y0)
k2 = f(x0+hm,y0+hm*k1)
k3 = f(x0+hm,y0+hm*k2)
k4 = f(x1,y0+h*k3)
y1 = y0 + hs*(k1+2*(k2+k3)+k4)
f(x,y) = ((2*y+3)/(4*x+5))**2

2𝑦+3 2
Ejemplo. Resolver 𝑦 ′ = ( ) = 𝑓(𝑥, 𝑦); 𝑦(1) = 2
4𝑥+5
(a) Con el método analítico
(b) Con el método de Runge-Kutta de cuarto orden, ℎ = 0.2, 2 iteraciones, calcular los errores relativos porcentuales

Solución

𝑑𝑦 2𝑦+3 2 𝑑𝑦 93𝑥+69
(a) =( ) → = (2𝑦 + 3)2 → 𝑦 = (resuelto más arriba)
𝑑𝑥 4𝑥+5 𝑑𝑥 59+22𝑥
(b) 𝑥𝑖+1 = 𝑥𝑖 + ℎ
𝑘1 = 𝑓(𝑥𝑖 , 𝑦𝑖 )
ℎ ℎ
𝑘2 = 𝑓 (𝑥𝑖 + , 𝑦𝑖 + 𝑘1 )
2 2
ℎ ℎ
𝑘3 = 𝑓 (𝑥𝑖 + , 𝑦𝑖 + 𝑘2 )
2 2
𝑘4 = 𝑓(𝑥𝑖+1 , 𝑦𝑖 + ℎ 𝑘3 )
ℎ
𝑦𝑖+1 = 𝑦𝑖 + (𝑘1 + 2(𝑘2 + 𝑘3 ) + 𝑘4 )
6

Primera iteración:

(𝑥0 , 𝑦0 ) = (1,2)

𝑖 = 0:

𝑥1 = 𝑥0 + ℎ = 1 + 0.2 = 1.2

𝑘1 = 𝑓(𝑥0 , 𝑦0 ) = 𝑓(1, 2) = 0.6049

ℎ ℎ 0.2 0.2
𝑘2 = 𝑓 (𝑥0 + , 𝑦0 + 𝑘1 ) = 𝑓 (1 + ,2 + ⋅ 0.6049) = 0.5739
2 2 2 2

ℎ ℎ 0.2 0.2
𝑘3 = 𝑓 (𝑥0 + , 𝑦0 + 𝑘2 ) = 𝑓 (1 + ,2+ ⋅ 0.5739) = 0.5729
2 2 2 2

𝑘4 = 𝑓(𝑥1 , 𝑦0 + ℎ𝑘3 ) = 𝑓(1.2, 2 + 0.2(0.5729)) = 0.5442

ℎ
𝑦1 = 𝑦0 + (𝑘1 + 2(𝑘2 + 𝑘3 ) + 𝑘4 )
6

0.2
𝑦1 = 2 + (0.6049 + 2(0.5739 + 0.5729) + 0.5442)
6

𝑦1 = 2.1148 ≈ 𝑦(1.2)
Segunda iteración: 𝒊 = 𝟏

(𝑥1 , 𝑦1 ) = (1.2,2.1148)

𝑥2 = 𝑥1 + ℎ = 1.2 + 0.2 = 1.4

𝑘1 = 𝑓(𝑥1 , 𝑦1 ) = 𝑓(1.2,2.1148) = 0.5442

ℎ ℎ 0.2 0.2
𝑘2 = 𝑓 (𝑥1 + , 𝑦1 + 𝑘1 ) = 𝑓 (1.2 + , 2.1148 + ⋅ 0.5442) = 0.5176
2 2 2 2

ℎ ℎ 0.2 0.2
𝑘3 = 𝑓 (𝑥1 + , 𝑦1 + 𝑘2 ) = 𝑓 (1.2 + , 2.1148 + ⋅ 0.5176) = 0.5169
2 2 2 2

𝑘4 = 𝑓(𝑥2 , 𝑦1 + ℎ𝑘3 ) = 𝑓(1.4, 2.1148 + 0.2 ⋅ 0.5169) = 0.4921

ℎ
𝑦2 = 𝑦1 + (𝑘1 + 2(𝑘2 + 𝑘3 ) + 𝑘4 )
6

0.2
𝑦2 = 2.1148 + (0.5442 + 2(0.5176 + 0.5169) + 0.4921)
6

𝑦2 = 2.2183 ≈ 𝑦(1.4)

Resumen

𝑥𝑖 𝑦𝑖 93𝑥𝑖 + 69 𝑒𝑡 %
𝑦(𝑥𝑖 ) =
59 + 22𝑥𝑖
1 2 𝑦(1) = 2 0%
1.2 2.1148 𝑦(1.2) = 2.1148 0%
1.4 2.2183 𝑦(1.4) = 2.2183 0%
5.4 Métodos de Runge Kutta
Dado el problema a valor inicial
𝑑𝑦
𝑦′ = = 𝑓(𝑥, 𝑦)
{ 𝑑𝑥
𝑦(𝑥0 ) = 𝑦0

Un método de Runge-Kutta a 𝑠 pisos está dado por La idea es conseguir las pendientes aproximadas de la
curva en algunos puntos del intervalo [𝑥𝑖 , 𝑥𝑖+1 ]
𝑘1 = 𝑓(𝑥𝑖 , 𝑦𝑖 )
𝑘2 = 𝑓(𝑥𝑖 + 𝑐2 ℎ, 𝑦𝑖 + ℎ𝑎21 𝑘1 )
𝑘3 = 𝑓(𝑥𝑖 + 𝑐3 ℎ, 𝑦𝑖 + ℎ(𝑎31 𝑘1 + 𝑎32 𝑘2 ))
⋮
𝑘𝑠 = 𝑓 (𝑥𝑖 + 𝑐𝑠 ℎ, 𝑦𝑖 + ℎ(𝑎𝑠1 𝑘1 + ⋯ + 𝑎𝑠,𝑠−1 𝑘𝑠−1 ))

𝑦𝑖+1 = 𝑦𝑖 + ℎ(𝑏1 𝑘1 + ⋯ + 𝑏𝑠 𝑘𝑠 )

De esta forma la pendiente final es 𝜙 = 𝑏1 𝑘1 + ⋯ + 𝑏𝑠 𝑘𝑠

𝑦𝑖+1 = 𝑦𝑖 + ℎ𝜙
Viene representado por el esquema:

𝑐1 = 0
𝑐2 𝑎21
𝑐3 𝑎31 𝑎32
⋮ ⋮ ⋮ ⋱
𝑐𝑠 𝑎𝑠1 𝑎𝑠2 ⋯ 𝑎𝑠,𝑠−1
𝑏1 𝑏2 ⋯ 𝑏𝑠−1 𝑏𝑠

Convenio: 𝑐1 = 0; 𝑐𝑖 = ∑𝑖−1𝑗=1 𝑎𝑖𝑗 , 2 ≤ 𝑖 ≤ 𝑠

Definición (Orden de un método de runge kutta) Un método de Runge-Kutta tiene el orden 𝑛 si
𝑦(𝑥0 + ℎ) − 𝑦1 = 𝑂(ℎ𝑛+1 )
𝑛+1 )
Donde 𝑂(ℎ es una constante.
Teorema. Bajo el convenio 𝑐1 = 0; 𝑐𝑖 = ∑𝑖−1 𝑗=1 𝑎𝑖𝑗 , 2 ≤ 𝑖 ≤ 𝑠, es suficiente considerar problemas autónomos de la forma
𝑦 ′ = 𝑓(𝑦) para la condición de orden.
5.4.1 Deducción de los métodos de Runge Kutta de orden 1, 2 y 3

Deducción del método de Runge Kutta de orden 1:

Intentemos con el método de 1 piso:

𝑐1 = 0
𝑏1

El método de Runge-Kutta consiste en:

𝑘1 = 𝑓(𝑥0 , 𝑦0 )
𝑦1 = 𝑦0 + ℎ(𝑏1 𝑘1 ) = 𝑦0 + ℎ𝑏1 𝑓(𝑥0 , 𝑦0 )
Por otro lado debido a la serie de Taylor tenemos:
𝑦(𝑥0 + ℎ) = 𝑦(𝑥0 ) + ℎ𝑦 ′ (𝑥0 ) + 𝑂(ℎ2 )

𝑦(𝑥0 + ℎ) = 𝑦0 + ℎ𝑓(𝑥0 , 𝑦0 ) + 𝑂(ℎ2 )

𝑦(𝑥0 + ℎ) − {𝑦0 + ℎ𝑓(𝑥0 , 𝑦0 )} = 𝑂(ℎ2 )
Igualando la expresión de las llaves con 𝑦1 se tiene:

𝑦0 + ℎ𝑏1 𝑓(𝑥0 , 𝑦0 ) = 𝑦0 + ℎ𝑓(𝑥0 , 𝑦0 )

De donde 𝑏1 = 1
Por lo tanto para que el método de un piso tenga orden 1, el método será de la forma:
𝑦1 = 𝑦0 + ℎ𝑓(𝑥0 , 𝑦0 )
Deducción del método de Runge Kutta de orden 2:
Esta vez considere el problema autónomo 𝑦 ′ = 𝑓(𝑦), 𝑦(𝑥0 ) = 𝑦0
Intentemos con el método de 2 pisos:

𝑐1 = 0
𝑐2 𝑐2
𝑏1 𝑏2
El método de Runge-Kutta consiste en:
𝑘1 = 𝑓(𝑦0 )
𝑘2 = 𝑓(𝑦0 + ℎ𝑐2 𝑘1 ) = 𝑓(𝑦0 ) + ℎ𝑐2 𝑓(𝑦0 )𝑓𝑦 (𝑦0 ) + 𝑂(ℎ2 )
𝑦1 = 𝑦0 + ℎ(𝑏1 𝑘1 + 𝑏2 𝑘2 ) = 𝑦0 + ℎ𝑏1 𝑓(𝑦0 ) + ℎ𝑏2 (𝑓(𝑦0 ) + ℎ𝑐2 𝑓(𝑦0 )𝑓𝑦 (𝑦0 )) + 𝑂(ℎ3 )
𝑦1 = 𝑦0 + ℎ𝑏1 𝑓(𝑦0 ) + ℎ𝑏2 𝑓(𝑦0 ) + ℎ2 𝑏2 𝑐2 𝑓(𝑦0 )𝑓𝑦 (𝑦0 ) + 𝑂(ℎ3 )
Por otro lado usando la expansión en serie de Taylor se tiene:
𝑦′′(𝑥0 ) 2
𝑦(𝑥0 + ℎ) = 𝑦(𝑥0 ) + ℎ𝑦 ′ (𝑥0 ) + ℎ + 𝑂(ℎ3 )
2

𝑓𝑦 (𝑦0 )𝑓(𝑦0 ) 2
𝑦(𝑥0 + ℎ) = {𝑦0 + ℎ𝑓(𝑦0 ) + ℎ } + 𝑂(ℎ3 )
2

Igualando la expresión de las llaves con 𝑦1 se tiene:

𝑓𝑦 (𝑦0 )𝑓(𝑦0 ) 2
𝑦0 + ℎ𝑏1 𝑓(𝑦0 ) + ℎ𝑏2 𝑓(𝑦0 ) + ℎ2 𝑏2 𝑐2 𝑓(𝑦0 )𝑓𝑦 (𝑦0 ) = 𝑦0 + ℎ𝑓(𝑦0 ) + ℎ
2
De donde se tiene el sistema:
ℎ𝑏1 + ℎ𝑏2 = ℎ 𝑏1 + 𝑏2 = 1
{ ℎ2 → { 1
ℎ2 𝑏2 𝑐2 = 𝑏2 𝑐2 =
2 2
La solución paramétrica es:
1 1
𝑐2 = 𝑡, 𝑏2 =, 𝑏1 = 1 − , (𝑡 ≠ 0)
2𝑡 2𝑡
El método de segundo orden Runge Kutta es por lo tanto:

0
𝑡 𝑡
1 1
1−
2𝑡 2𝑡

Por ejemplo, si 𝑡 = 1:
0
1 1
1 1
2 2
Al aplicar el método al problema 𝑦 ′ = 𝑓(𝑥, 𝑦), 𝑦(𝑥0 ) = 𝑦0 :
𝑘1 = 𝑓(𝑥0 , 𝑦0 )
𝑘2 = 𝑓(𝑥0 + 1ℎ, 𝑦0 + 1ℎ𝑘1 )
1 1 𝑓(𝑥0 , 𝑦0 ) + 𝑓(𝑥0 + 1ℎ, 𝑦0 + 1ℎ𝑓(𝑥0 , 𝑦0 ))
𝑦1 = 𝑦0 + ℎ ( 𝑘1 + 𝑘2 ) = 𝑦0 + ℎ
2 2 2
Visto de otra forma es
𝑓(𝑥 ,𝑦 )+𝑓(𝑥 ,𝑦 0 )
𝑦1 = 𝑦0 + ℎ 0 0 1 1
siendo 𝑦10 = 𝑦0 + ℎ𝑓(𝑥0 , 𝑦0 )
2
Es decir se trata del método de Heun
Deducción del método de Runge Kutta de orden 3:
Considere el problema autónomo 𝑦 ′ = 𝑓(𝑦), 𝑦(𝑥0 ) = 𝑦0
Intentemos con el método de 3 pisos:

𝑐1 = 0
𝑐2 𝑐2
𝑐3 𝑎31 𝑐3 − 𝑎31
𝑏1 𝑏2 𝑏3
El método de Runge-Kutta consiste en:
(Para simplificar en las funciones que tengan por argumento 𝑦0 se omitirá este)
𝑘1 = 𝑓
1
𝑘2 = 𝑓(𝑦0 + ℎ𝑐2 𝑘1 ) = 𝑓 + ℎ𝑐2 𝑓𝑓𝑦 + [𝑓𝑦𝑦 ℎ2 𝑐22 𝑓 2 ] + 𝑂(ℎ3 )
2
𝑘3 = 𝑓(𝑦0 + ℎ(𝑎31 𝑘1 + (𝑐3 − 𝑎31 )𝑘2 ))
1
= 𝑓 (𝑦0 + ℎ (𝑐3 𝑓 + (𝑐3 − 𝑎31 ) (ℎ𝑐2 𝑓𝑓𝑦 + [𝑓𝑦𝑦 ℎ2 𝑐22 𝑓 2 ] + 𝑂(ℎ3 ))))
2
1
= 𝑓 (𝑦0 + 𝑐3 𝑓ℎ + ℎ2 (𝑐3 − 𝑎31 )𝑐2 𝑓𝑓𝑦 + (𝑐3 − 𝑎31 )ℎ3 𝑐22 𝑓𝑦𝑦 𝑓 2 + 𝑂(ℎ4 ))
2
1
= 𝑓 + (𝑐3 𝑓ℎ + ℎ2 (𝑐3 − 𝑎31 )𝑐2 𝑓𝑓𝑦 + (𝑐3 − 𝑎31 )ℎ3 𝑐22 𝑓𝑦𝑦 𝑓 2 + 𝑂(ℎ4 )) 𝑓𝑦
2
2
1 1
+ (𝑐3 𝑓ℎ + ℎ 3 − 𝑎31 )𝑐2 𝑓𝑓𝑦 + (𝑐3 − 𝑎31 )ℎ 𝑐2 𝑓𝑦𝑦 𝑓 + 𝑂(ℎ )) 𝑓𝑦𝑦 + 𝑂(ℎ3 )
2 (𝑐 3 2 2 4
2 2
𝑦1 = 𝑦0 + ℎ(𝑏1 𝑘1 + 𝑏2 𝑘2 + 𝑏3 𝑘3 )
1
= 𝑦0 + ℎ𝑏1 𝑓 + ℎ𝑏2 (𝑓 + ℎ𝑐2 𝑓𝑓𝑦 + [𝑓𝑦𝑦 ℎ2 𝑐22 𝑓 2 ]) + 𝑂(ℎ4 ) +
2
1
+ℎ𝑏3 (𝑓 + (𝑐3 𝑓ℎ + ℎ2 (𝑐3 − 𝑎31 )𝑐2 𝑓𝑓𝑦 + (𝑐3 − 𝑎31 )ℎ3 𝑐22 𝑓𝑦𝑦 𝑓 2 + 𝑂(ℎ4 )) 𝑓𝑦
2
2
1 1
+ (𝑐3 𝑓ℎ + ℎ2 (𝑐3 − 𝑎31 )𝑐2 𝑓𝑓𝑦 + (𝑐3 − 𝑎31 )ℎ3 𝑐22 𝑓𝑦𝑦 𝑓 2 + 𝑂(ℎ4 )) 𝑓𝑦𝑦 ) + 𝑂(ℎ4 )
2 2
1 1
= 𝑦0 + ℎ𝑏1 𝑓 + ℎ𝑏2 𝑓 + ℎ2 𝑏2 𝑐2 𝑓𝑓𝑦 + ℎ3 𝑏2 𝑐22 𝑓 2 𝑓𝑦𝑦 + ℎ𝑏3 𝑓 + ℎ2 𝑏3 𝑐3 𝑓𝑓𝑦 + ℎ3 𝑏3 (𝑐3 − 𝑎31 )𝑐2 𝑓𝑓𝑦2 + ℎ3 𝑏3 𝑐32 𝑓 2 𝑓𝑦𝑦 + 𝑂(ℎ4 )
2 2

Por otro lado debido a la serie de Taylor tenemos:

𝑦′′(𝑥0 ) 2 𝑦′′′(𝑥0 ) 3
𝑦(𝑥0 + ℎ) = 𝑦(𝑥0 ) + ℎ𝑦 ′ (𝑥0 ) + ℎ + ℎ + 𝑂(ℎ4 )
2 6

𝑓𝑦 𝑓 2 𝑓𝑦𝑦 𝑓 2 + 𝑓𝑦 𝑓𝑦 𝑓 3
𝑦(𝑥0 + ℎ) = {𝑦0 + ℎ𝑓 + ℎ + ℎ } + 𝑂(ℎ4 )
2 6

Igualando la expresión de las llaves con 𝑦1 se tiene:

1 1
𝑦0 + ℎ𝑏1 𝑓 + ℎ𝑏2 𝑓 + ℎ2 𝑏2 𝑐2 𝑓𝑓𝑦 + ℎ3 𝑏2 𝑐22 𝑓 2 𝑓𝑦𝑦 + ℎ𝑏3 𝑓 + ℎ2 𝑏3 𝑐3 𝑓𝑓𝑦 + ℎ3 𝑏3 (𝑐3 − 𝑎31 )𝑐2 𝑓𝑓𝑦2 + ℎ3 𝑏3 𝑐32 𝑓 2 𝑓𝑦𝑦
2 2
𝑓𝑦 𝑓 2 𝑓𝑦𝑦 𝑓 2 + 𝑓𝑦2 𝑓 3
= 𝑦0 + ℎ𝑓 + ℎ + ℎ
2 6
De donde se tiene el sistema:
𝑏1 + 𝑏2 + 𝑏3 = 1 𝑏1 + 𝑏2 + 𝑏3 = 1
1 1
𝑏2 𝑐2 + 𝑏3 𝑐3 = 𝑏2 𝑐2 + 𝑏3 𝑐3 =
2 2
2 2
1 → 2 2
1
𝑏2 𝑐2 + 𝑏3 𝑐3 = 𝑏2 𝑐2 + 𝑏3 𝑐3 =
3 3
1 1
{𝑏3 (𝑐3 − 𝑎31 )𝑐2 = 6 {𝑏3 (𝑐3 − 𝑎31 )𝑐2 = 6
La solución paramétrica es:
2 − 3𝑣 2 − 3𝑢 𝑣 (3 𝑢2 − 3 𝑢 + 𝑣) 6𝑢𝑣 − 3𝑢 − 3𝑣 + 2
𝑐2 = 𝑢, 𝑐3 = 𝑣, 𝑏2 = , 𝑏3 = ; 𝑎31 = ; 𝑏1 =
6 𝑢 (𝑢 − 𝑣) 6 𝑣 (𝑣 − 𝑢) 𝑢 (3 𝑢 − 2) 6𝑢𝑣
(𝑢 ≠ 0, 𝑣 ≠ 0, 𝑢 ≠ 𝑣)
Otras soluciones son:
2 1 1 3 2
𝑏3 = 𝑢, 𝑎31 = − , 𝑏 = , 𝑏 = − 𝑢, 𝑐2 = 𝑐3 = ; 𝑢 ≠ 0
3 4𝑢 1 4 2 4 3
El método de tercer orden de Runge Kutta es por lo tanto:
𝑐1 = 0
𝑢 𝑢
𝑣 𝑣 (3 𝑢2 − 3 𝑢 + 𝑣) 𝑣 (𝑢 − 𝑣)
𝑢 (3 𝑢 − 2) 𝑢 (3 𝑢 − 2)
6𝑢𝑣 − 3𝑢 − 3𝑣 + 2 2 − 3𝑣 2 − 3𝑢
6𝑢𝑣 6 𝑢 (𝑢 − 𝑣) 6 𝑣 (𝑣 − 𝑢)

O bien:
𝑐1 = 0
2 2
3 3
2 2 1 1
−
3 3 4𝑢 4𝑢
1 3 𝑢
−𝑢
4 4
1
Ejemplo. Resolver el P.V.I. 𝑥𝑦 ′ + 𝑦 = ; 𝑦(1) = 2
𝑦2

Con el método de Runge-Kutta de cuarto orden, ℎ = 0.2, 2 iteraciones, calcular los errores relativos porcentuales
3
√𝑥 3 +7
Solución exacta es 𝑦(𝑥) = .
𝑥