Apuntes de Métodos Numéricos JFC UPIITA

Ecuaciones diferenciales ordinarias.

Introducción
Los métodos numéricos son adecuados para para resolver ecuaciones diferenciales
ordinarias que no pueden solucionarse analíticamente. La forma general, de una
ecuación diferencial de primer orden es

y ' = f ( x, y) , (1) 1

y las ecuaciones de orden superior se pueden escribir como

y
(n)
(
= f x, y, y '', y '',…, y
(n−1)
). (2)

Las ecs. de la formada dada por la ec. (1) tiene soluciones que se pueden escribir
como

y = ϕ ( x) , (3)

La solución de la ec. (1) queda definida una vez que se conocen las condiciones de
frontera y la relación entre x y y de la ec. (3).

Una forma común de manejar ecuaciones de segundo o mayores ordenes consiste en

reemplazarlas con sistemas equivalentes de ecuaciones de primer grado. Note que
las ecuaciones de mayores órdenes como la ec. (2) se pueden transformar en un
conjunto de n ecuaciones de primer orden. Usando la notación

y0 = y
y1 = y '
y2 = y '' , (4)
⋮
(n−1)
yn−1 = y

entonces las ecuaciones equivalentes de primer orden son

 Apuntes de Métodos Numéricos JFC UPIITA

y '0 = y1
y '1 = y2
y '2 = y3 . (5)
⋮
y 'n = f ( x, y0 , y1 ,…, yn−1 )

Por esta razón, se enfatizarán los métodos para resolver las ecuaciones de primer
orden. 2

Método de Euler
El método más simple para realizar la integración numérica de ecuaciones de
ecuaciones diferenciales ordinarias de primer orden es el método de Euler. Considere
el problema

dy
= f ( x, y) = g ( x) . (6)
dx

El método asume que para una pequeña distancia ∆x a lo largo del eje x desde un
punto inicial x0 , la función g ( x) es una constante igual a g ( x0 ) :

dy
= f ( x0 , y0 ) = g ( x0 ) . (7)
dx x= x
0

Si se asume que g ( x) = g ( x0 ) para todos los valores de x entre x0 y x1 = x + ∆x ,

entonces el cambio en y que corresponde a un pequeño cambio ∆x en x se da
aproximadamente por

∆y
= f ( x0 , y0 ) . (8)
∆x

Sí y1 = y0 + ∆y , entonces la ec. (8) se puede reescribir como:

y1 = y0 + ∆xf ( x0 , y0 ) = y0 + hf ( x0 , y0 ) , (9)

donde h = ∆x . Por medio de la ec. (9) y1 se puede calcular. De forma similar se

puede repetir el proceso que lleve a una solución general de la forma:
 Apuntes de Métodos Numéricos JFC UPIITA

y2 = y1 + hf ( x1 , y1 )
y3 = y2 + hf ( x2 , y2 )
. (10)
⋮
yn+1 = yn + hf ( xn , yn )

Lo cual se conoce como el método de Euler. Este proceso de cálculo, paso a paso, a lo
largo del eje x desde el punto inicial x0 al punto final requerido.
3
El método puede expresarse de forma geométrica, como se muestra en la fig. 1. La
pendiente de la curva de y contra x es dy dx , es decir, f ( x, y) . De aquí que, el punto
de inicio ( xn , yn ) de cada intervalo la pendiente es exactamente f ( xn , yn ) . Este valor
es una buena aproximación, si h es pequeña, el promedio de la pendiente sobre el
intervalo h de xn a xn+1 .

Ejercicio 1. Con el método de Euler resuelva numéricamente la ecuación diferencial

dy
= −2 x3 + 12 x2 − 20 x + 8.5 ,
dx
desde x = 0 hasta x = 4 con un tamaño de paso de h = 0.5 . La condición inicial es
x = 0 es y = 1 . Compare el resultado con la solución exacta que es:

y = −0.5 x4 + 4 x3 − 10 x2 + 8.5 x + 1 .

Ejercicio 2. Utilice el método de Euler para resolver la ecuación diferencial ordinaria

dy
= x+ y
dx

considerando las condiciones iniciales x0 = 0 y y0 = 1 . La ecuación anterior tiene

como solución la función

y = 2e x − x − 1 .

Compare el resultado del cálculo numérico con la solución analítica considerando

que x = 1 y h = 0.001 .
 Apuntes de Métodos Numéricos JFC UPIITA

Método de Runge-Kutta de segundo orden

En este método, de forma semejante al método de Euler se considera una expansión
de la función y ( x) por medio de una serie de Taylor. Lo que se busca es una solución
aproximada al problema de valor inicial siguiente:

(
y ' ( x) = f x, y ( x) , ) (11)

donde y ( x0 ) = y0 . Al realizar la expansión de la serie de Taylor de la función y ( x) 4

alrededor del punto x = x0 se obtiene

y ' ( x0 ) y '' ( x0 ) 2 y '' ( x0 ) 3

y ( x) ≈ y ( x0 ) + ( x − x0 ) + ( x − x0 ) + (x − x ) 0
+⋯ . (12)
1! 2! 3!

Si se trunca la serie después de tres términos y se retienen los demás se obtiene

y ' ( x0 ) y '' ( x0 ) 2
y ( x) ≈ y ( x0 ) + (x − x ) +0 (x − x ) 0
. (13)
1! 2!

Si se propone que h = x − x0 y observando que y ' ( x0 ) = f x0 , y ( x0 ) = f ( x0 , y0 ) y que ( )

y '' ( x0 ) = f ' ( x0 , y0 ) la ec. (13) se puede escribir como

1
yn+1 ≈ yn + f ( xn , yn ) h + f ' ( xn , yn ) h2 . (14)
2!

Adicionalmente, a diferencia del método de Euler que sólo evalúa las derivadas al
inicio del intervalo se puede construir una integral más simétrica al considerar el
punto medio del intervalo. Por este motivo, se introducen los parámetros k1 y k2 de
la siguiente manera:

k1 = f ( xn , yn ) h
 1 1 
k2 = f  xn + h, yn + k1  h , (15)
 2 2 
yn+1 = yn + k2 + O h3 ( )
 Apuntes de Métodos Numéricos JFC UPIITA

donde la elección de la simetría anterior cancela el error de primer orden, haciéndolo

de segundo orden. Entonces, la ec. (15) se puede escribir como
 h h 
yn+1 = yn + hf  xn + , yn + f ( xn , yn ) . (16)
 2 2 

Método Runge-Kutta de cuarto orden

A pesar de que el método de Euler provee una mejor aproximación que el método de
Euler, éste no se usa mucho en la práctica. En su lugar, debido a que se tiene mayor 5
precisión se suele emplear el método de Runge-Kutta de cuarto orden para resolver
ecuaciones diferenciales de primer orden. En esencia éste se obtiene reteniendo más
términos en la serie de Taylor de y ( x) alrededor de x = x0 . La derivación es muy
similar a la del método de Runge-Kutta de segundo orden, y se puede resumir en las
siguientes ecuaciones:

k1 = f ( xn , yn ) h
 1 1 
k2 = f  xn + h, yn + k1  h
 2 2 
 1 1 
k3 = f  xn + h, yn + k2  h . (17)
 2 2 
k4 = f ( xn + h, yn + k3 ) h
1
yn+1 = yn +
6
( ( )
k1 + 2 k2 + 2k3 + k4 ) + O h5

Nótese que el valor yn+1 se determina el valor presente yn más el producto del
tamaño del paso h y un promedio de cuatro pendientes con mayor peso al punto
medio y a las pendientes.

Método Runge-Kutta-Fehlberg
El método de Runge-Kutta óptimo de un orden particular es aquel cuyo error de
truncado es mínimo. Se podrían utilizar métodos de Runge-Kutta de orden superior;
pero se puede producir una sobrecarga en las operaciones que se tienen que realizar
que se va incrementando conforma aumenta el orden del método. El método de
Runge-Kutta-Fehlberg incorpora un procedimiento adaptativo, en el cual se estima
 Apuntes de Métodos Numéricos JFC UPIITA

el error de truncamiento en cada paso y el ancho del paso h se ajusta

automáticamente para mantener el error dentro de unos límites preestablecidos. En
cada iteración, este método compara dos valores estimados de la solución. Sin
embargo, si los resultados no concuerdan con la exactitud preestablecida, el ancho
de paso se reduce; de lo contrario el ancho de pulso se incrementa.

El método de Runge-Kutta-Fehlberg usa un método Runge-Kutta con un error de

truncamiento local del orden de error de cinco, 6

(5) 16 6656 28561 9 2

yn+1 = yn + k1 + k3 + k4 − k5 + k, (18)
135 12825 56430 50 55 6

para estimar el error local en un método de Runge-Kutta de orden cuatro, dado por

( 4) 25 1408 2197 1
yn+1 = yn + k1 + k3 + k4 − k5 , (19)
216 2565 4104 5

donde los coeficientes de las ecuaciones se dan por:

k1 = f ( xn , yn ) h
 1 1 

k2 = f  xn + h, yn + k1  h
 4 4 
 3 3 9 
k3 = f  xn + h, yn + k1 + k h
 8 32 32 2 
 12 1932 7200 7296  (20)

k4 = f  xn + h, yn + k − k + k h
 23 2197 1 2197 2 2197 3 
 439 3680 845 
k5 = f  xn + h, yn + k1 − 8 k2 + k3 − k h
 216 513 4104 4 
 1 8 3544 1859 11 

k6 = f  xn + h, yn − k + 2k2 − k + k − k h
 2 27 1 2565 3 4104 4 40 5 

En este método se hacen 6 evaluaciones por paso, que es menor al número que se
requieren si se usan juntos los métodos de orden cuatro y cinco, los cuales
requerirían de al menos de 10 evaluaciones.

El tamaño de paso óptimo se determina multiplicando el factor de escala s por el

tamaño de paso actual h , donde el factor de escala s se da por
 Apuntes de Métodos Numéricos JFC UPIITA

1
 tolerancia  4
s =   , (21)
 2 × error 

donde la tolerancia es el parámetro de control del error preestablecido y el error es

el error de truncamiento del método dado por:

1 (5) ( 4) 1 1 128 2197 1 2

error = yn+1 − yn+1 = k1 − k3 − k4 + k5 + k . (22)
h h 360 4275 75240 50 55 6
7
Ecuaciones diferenciales de orden superior y sistemas de ecuaciones
diferenciales
Un sistema m − ésimo orden de problemas de valor inicial de primer orden tiene la
forma

du1
= f1 (t, u1 , u2 ,…, um )
dt
du2
= f2 (t, u1 , u2 ,…, um )
dt , (23)
⋮
dum
= fm (t, u1 , u2 ,…, um )
dt

para a ≤ t ≤ b , con condiciones iniciales

u1 ( a) = α1 , u2 ( a) = α2 ,…, um ( a) = αm . (24)

El objetivo consiste en encontrar m funciones

u1 (t), u2 (t),…, um (t) , (25)

las cuales satisfacen cada una de las ecuaciones diferenciales, junto con todas sus
respectivas condiciones iniciales.

Los métodos para resolver sistemas de ecuaciones diferenciales de primer orden son
generalizaciones de los métodos de ecuaciones de primer orden simples que se han
presentado previamente. Por ejemplo, se reescribe con una nomenclatura
ligeramente diferente a le expuesta anteriormente el método de cuarto orden Runge-
Kutta, teniéndose
 Apuntes de Métodos Numéricos JFC UPIITA

y0 = α
k1 = hf (ti , yi )
 h k1 

k2 = hf ti + , yi + 
 2 2 
 h k , (26)
k3 = hf ti + , yi + 2 
 2 2 
k4 = hf (ti + h, yi + k3 )
8
1
yi+1 = yi + (k + 2k2 + 2k3 + k4 )
6 1

para cada i = 0,1,…, N − 1 ; el cual se utiliza para resolver el problema de valor inicial
de primer orden

y ' = f (t, y), a ≤ t ≤ b, y ( a) = α . (27)

Este método se puede extender para resolver ecuaciones diferenciales de orden

superior, si se hacen las siguientes consideraciones:

Se elige un entero N > 0 y se establece que h =

(b − a) . Se parte el intervalo  a, b
N  
en N subintervalos con los puntos de malla

t j = a + jh , (28)

para cada j = 0,1,…, N . Se usa la notación yi, j , para cada j = 0,1,…, N e

i = 1,2,…, m para denotar una aproximación para la i −ésima solución ui (t) de la ec.
(23), en el j − ésimo punto de la malla t j . Para las condiciones iniciales, establezca

y1,0 = α1, y2,0 = α2 ,…, ym,0 = αm . (29)

Suponga que se calculan los valores

y1, j , y2, j ,…, ym, j . (30)

Entonces, los valores

y1, j+1 , y2, j+1 ,…, ym, j+1 , (31)

 Apuntes de Métodos Numéricos JFC UPIITA

se pueden calcular por medio de las siguientes ecuaciones:

k1, j = hfi (t j , y1, j , y2, j ,…, ym, j ) , (32)

para cada i = 1,2,…, m ; i = 1,2,…, m

 h k1,1 k1,2 k1,m 
  ,
k2, j = hfi t j + , y1, j + , y2, j + ,…, ym, j + (33)
 2 2 2 2 
9
para cada i = 1,2,…, m ;
 h k2,1 k2,2 k2,m 
 
k3, j = hfi t j + , y1, j + , y2, j + ,…, ym, j + (34)
 2 2 2 2 

para cada i = 1,2,…, m ;

k4, j = hfi (t j + h, y1, j + k3,1 , y2, j + k3,2 ,…, ym, j + k3,m ) , (35)

para cada i = 1,2,…, m ; y entonces se tiene que

1
yi, j+1 = yi, j + (k + 2k2,i + 2k3,i + k4,i ) ,
6 1,i
(36)

para cada i = 1,2,…, m . Nótese que todos los valores k1,1 , k1,2 ,…, k1,m se tienen que
calcular antes de que se calculen los términos k2,i para cada i = 1,2,…, m .

Ejercicio 3. Por medio de Runge-Kutta de cuarto orden resuelva la ecuación

diferencial

x 2 − y2
y '' = 2
, y (0) = 1, y ' (0) = 0 , (37)
1 + ( y ')

para x = 0.5 y x = 1 .

Ejercicio 4. En la fig. 1 se muestra un circuito eléctrico, por medio del método de

mallas determine las corrientes i1 e i2 por medio del método de Runge-Kutta de
cuarto orden. Considere el intervalo 0,0.5 y h = 0.1 . El interruptor se cierra cuando
t = 0 (s). La solución exacta del sistema es:
 Apuntes de Métodos Numéricos JFC UPIITA

i1 (t) = −3.375e−2 t + 1.875e−0.4 t + 1.5

. (38)
i2 (t) = −2.25e−2 t + 2.25e−0.4 t

10

Fig. 1. Circuito eléctrico.

Ejercicio 5. Transforme el problema de valor inicial de segundo orden

y ''− 2 y '+ 2 y = e2 t sen (t), para 0 ≤ t ≤ 1 (s), con y (0) = −0.4, y y ' (0) = −0.6 ,
(39)

en un sistema de problemas de valor inicial de primer orden y use el método de

Runge-Kutta con h = 0.1 para aproximar la solución.

Referencias

[1] J. Xu, 2008, “Practical Numerical Methods with C#: Numerical Programming and
Math Functions for Real-World .NET Applications with C#”, UniCAD.

[2] W. Dos Passos, 2009, “Numerical Methods, Algorithms and Tools in C#”, CRC
Press.

[3] R. L. Burden, D. J. Faires, A. M. Burden, 2016,“Análisis Numérico”, CENGAGE

Learning.