UNIVERSIDAD NACIONAL DE

LOMAS DE ZAMORA

Facultad de ingenierı́a

Métodos
Numéricos

Autor: Pablo Azcurra

1
 Índice general

Capı́tulo 1. Métodos numéricos 3

1.1. Método de Euler 4
1.2. Método de Taylor de segundo orden 7
1.3. Métodos de Runge-Kutta 8

2
 Capı́tulo 1

Métodos numéricos

Hasta ahora hemos trabajado con métodos analı́ticos para encontrar la fórmula
explı́cita o implı́cita de la solución y(x) de una ED, ahora presentaremos un enfoque
numérico para obtener, no la función solución, sino una estimación de los valores
que toma esa función en ciertos valores discretos de su dominio.

El problema que vamos a plantear con estos métodos numéricos es el siguiente:

Dada una ecuació diferencial de la forma y ′ (x) = f (x, y), con valor inicial y(a) =
y0 buscaremos una aproximación del valor que toma y en un tiempo posterior x = b.

Es decir, si y(x) es la solución analı́tica de la ED y ′ (x) = f (x, y) con y(a) = y0 ,

que podemos representar gráficamente (ver Figura 1.1), y x = b es un punto del
dominio de y, encontraremos una estimación del valor y(b).

Figura 1.1.

Para llegar a la aproximación de y(b) será necesario estimar algunos valores

intermedios. El algoritmo que vamosa estudiar plantea el siguiente patrón:
1. Elegimos n + 1 valores de x entre a y b.
a = x0 < x1 < x2 < ... < xn = b

Figura 1.2.

2. Tomamos y(a) = y(x0 ) = y0 como condición inicial.

3. Buscamos una estimación de y(x1 ) usando algún método numérico, a esta
estimación la llamamos y1 , repetimos el método para buscar una estimación

3
 1.1. MÉTODO DE EULER 4

de y(x2 ), la llamamos y2 y volvemos a repetir este paso hasta llegar a la

estimación de y(xn ) = y(b) ≈ yn .

Notación: a los puntos x0 , x1 , x2 , ..., xn se los llama partición del intervalo[a, b]:

a = x0 < x1 < x2 < ... < xn = b

Los métodos numéricos que estudiaremos en este curso trabajan con particiones
regulares, es decir que la distancia entre dos puntos consecutivos xi y xi+1 será la
misma para todos los puntos de la partición. A esa distancia la llamaremos h, por
lo tanto:
b−a
h= y xi+1 = xi + h
n

1.1. Método de Euler

Este método propone buscar los valores y1 , y2 ,...,yn usando aproximaciones con
rectas tangentes.

Dada la ED y ′ = f (x, y) con condición inicial y(a) = y0 , supongamos que f y

fy son continuas, esto nos garantiza que la ED tiene solución única. Tomamos una
partición del intervalo [a, b] de n intervalitos.

a = x0 < x1 < x2 < ... < xn = b

b−a
Y llamamos h a la longitud del paso: h = .
n
Usamos la recta tangente al gráfico de y(x) en el punto (x0 , y0 ) para buscar la
aproximación de y(x1 ). La recta tangente en el punto (x0 , y0 ) es
y = y ′ (x0 ).(x − x0 ) + y(x0 )
A esta recta la evaluamos en x1 para conseguir la estimación de y(x1 ), a esta
estimación la llamamos y1 . Y como la ED es de la forma y ′ (x) = f (x, y), podemos
asegurar que y ′ (x0 ) = f (x0 , y0 ). Por lo tanto:

y1 = y ′ (x0 ). (x1 − x0 ) +y(x0 )

| {z }
h

y1 = f (x0 , y0 ).h + y0

Ahora nos paramos en el punto (x1 , y1 ) y consideramos la ED y ′ (x) = f (x, y)

con condición inicial y(x1 ) = y1 . Existe también una única solución para esta ED
cuya gráfica pasa por elpunto (x1 , y1 ) y no es muy distinta a la gráfica de la solución
anterior ya que las dos provienen de la misma ED, solo que con distinta condición
inicial. A esta nueva solución también la llamaremos y(x). En la figura 1.3 podemos
ver un gráfico que representa esta situación.
 1.1. MÉTODO DE EULER 5

Hx1,y1L
y0

x0 x1

Figura 1.3.

Para aproximar y(x2 ) usamos la recta tangente que pasa por el punto (x1 .y1 ) y
la valuamos en x2 , de la misma manera que hizimos antes. Nos queda:

y = y ′ (x1 ).(x − x1 ) + y(x1 )

y2 = f (x1 , y1 ). (x2 − x1 ) +y1

| {z }
h

y2 = f (x1 , y1 ).h + y1

Podemos generalizar esta idea afirmando que cada vez que conseguimos un valor
yi , para conseguir yi+1 debemos hacer

yi+1 = f (xi , yi ).h + yi

En definitiva, el método de Euler se puede resumir con los siguientes pasos:

Método de Euler: dada la ED y ′ (x) = f (x, y) con condición inicial y(a) = y0 ,

se pretende encontrar una aproximación del valor y(b) con b > a.
1. Tomamos una partición a = x0 < x1 < x2 < ... < xn = b del intervalo [a, b],
b−a
con n ∈ N, y llamamos h = a la longitud del paso.
n
2. Llamamos y0 = y(x0 ) = y(a)
3. Calculamos de forma recursiva y1 , y2 , ..., yn de esta manera:
yi+1 = f (xi , yi ).h + yi
4. yn es la aproximación de y(b).
Ejemplo 1.1.
Sea y(x) la solución de la ecuación diferencial:
y ′ (x) = x + x.y(x) con condición inicial y(0) = 1.
Estimar y(1) con el método de Euler usando un tamaño de paso h = 14 .
 1.1. MÉTODO DE EULER 6

y ′ = x + x.y =⇒ f (x, y) = x + x.y.

Condición inicial y(0) = 1 =⇒ (x0 , y0 ) = (0, 1).

La relación recursiva es:

yi+1 = h.f (xi , yi ) + yi
1
yi+1 = .(xi + xi .yi ) + yi
4
El intervalo en el cual vamos a trabajar es el que empieza en x = 0 (por la
condición inicial) y termina en x = 1 (porque nos piden encontrar una aproximación
de y(1)) y como el tamaõ del paso es 14 , vamos a partir este intervalo en cuatro
intervalitos.

Empezamos con el valor y0 = 1 y calculamos recursivamente y1 , y2 , y3 y y4 .


 y1 = h.f (x0 , y0 ) + y0





 y1 = 41 .f (0, 1) + 1

1
y1 : =⇒ (x1 , y1 ) = ( , 1)

 y1 = 14 .(0 + 0,1) + 1 4






y1 = 1

 y2 = h.f (x1 , y1 ) + y1





 y2 = 14 .f ( 41 , 1) + 1

1 9
y2 : =⇒ (x2 , y2 ) = ( , )

 y2 = 14 .( 41 + 14 ,1) + 1 2 8






y2 = 98

 y3 = h.f (x2 , y2 ) + y2





 y3 = 14 .f ( 21 , 89 ) + 98

3 89
y3 : =⇒ (x3 , y3 ) = ( , )

 y3 = 41 .( 21 + 12 . 89 ) + 89 4 64






y3 = 89
64

 y4 = h.f (x3 , y3 ) + y3





 y4 = 14 .f ( 34 , 64
 89
) + 64 89
1883
y4 : =⇒ (x4 , y4 ) = (1, )

 y = 1 3
.( + 3 89
. ) + 89 1024

 4 4 4 4 64 64



 1883
y4 = 1024
1883
Por lo tanto, el método de Euler nos da la aproximación y(1) ≈ 1024
.
 1.2. MÉTODO DE TAYLOR DE SEGUNDO ORDEN 7

En este caso podemos ver qué tan buena es la aproximación porque la ecuación
diferencial y ′ = x + x.y se puede resolver por el método de separación de variables.
x2
La solución a esta ecuación es: y(x) = 2.e 2 − 1 (los cálculos quedan a cargo del
lector).
1 √
Podemos afirmar entonces que y(1) = 2.e 2 − 1 = 2. e − 1 ≈ 2, 297, mientras
que nuestra aproximación nos dió: 1883
1024
≈ 1, 839.

El error en la estimación de y(1) resultó ser aproximadamente: 2, 297 − 1, 839 =

0, 458. Si volvemos a aplicar el método tomando tamaños de paso más chicos la
aproximación será mejor.

1.2. Método de Taylor de segundo orden

A diferencia del método de Euler, en este método reemplazamos las rectas tan-
gentes (polinomio de Taylor de orden 1) por polinomios de Taylor de orden 2.

Observación: si bien no es objetivo de este curso analizar el orden del error que se
comete al utilizar cada uno de los métodos propuestos, el lector podrá percibir que
se obtendrá una mejor aproximación si se utilizan polinomios de orden 2 en lugar
de rectas.

Si y(x) es la solución de la ED y ′ = f (x, y) con condición inicial y(x0 ) = y0 , el

polinomio de Taylor de orden 2 en x0 es
y ′′ (x0 )
′
T (x) = y(x0 ) + y (x0 ).(x − x0 ) + .(x − x0 )2
2
Para desarrollarlo tenemos que calcular y (x).
′′

Si y ′ = f (x, y) =⇒ y ′′ = fx (x, y) + fy (x, y).y ′ = fx (x, y) + fy (x, y).f (x, y).

El polinomio de Taylor nos queda entonces:

fx (x0 , y0 ) + fy (x0 , y0 ).f (x0 , y0 )
T (x) = y(x0 ) + f (x0 , y0 ).(x − x0 ) + .(x − x0 )2
2

Ahora podemos obtener la aproximación de y(x1 ) ≈ y1 evaluando este polinomio

en x1 .
fx (x0 , y0 ) + fy (x0 , y0 ).f (x0 , y0 )
y1 = y(x0 ) + f (x0 , y0 ). (x1 − x0 ) + . (x1 − x0 )2
| {z } 2 | {z }
h h2

h2
y1 = y(x0 ) + f (x0 , y0 ).h + (fx (x0 , y0 ) + fy (x0 , y0 ).f (x0 , y0 )).
2

El mismo procedimiento realizamos para aproximar y(x2 ) ≈ y2 , y(x3 ) ≈ y3 , hasta

llegar a y(xn ) ≈ yn
 1.3. MÉTODOS DE RUNGE-KUTTA 8

Método de Taylor de segundo orden: dada la ED y ′ (x) = f (x, y) con con-

dición inicial y(a) = y0 , se pretende encontrar una aproximación del valor y(b) con
b > a.
1. Tomamos una partición a = x0 < x1 < x2 < ... < xn = b del intervalo [a, b],
b−a
con n ∈ N, y llamamos h = a la longitud del paso.
n
2. Llamamos y0 = y(x0 ) = y(a)
3. Calculamos de forma recursiva y1 , y2 , ..., yn de esta manera:
h2
yi+1 = y(xi ) + f (xi , yi ).h + (fx (xi , yi ) + fy (xi , yi ).f (xi , yi )).
2
4. yn es la aproximación de y(b).

Desafortunadamente, el métode de Taylor de orden 2 requiere el cálculo de 2

derivadas parciales y 3 evaluaciones de funciones lo que lo hace algo engorroso al
trabajar, nosotros estamos interesados en trabajar con métodos que no presenten
este grado de dificultad.

De la misma manera que usamos polinomios de Taylor de orden 2 para hallar las
aproximaciones en los puntos de la partición podemos usar polinomios de Taylor de
mayor orden lo cuál mejorará el orden del error pero, al mismo tiempo, aumentará
considerablemente la complejidad del cálculo. Esta idea da origen los Métodos de
Taylor de orden superior que no serán trabajados en este curso.

1.3. Métodos de Runge-Kutta

Los métodos de Runge-Kutta buscan obtener aproximaciones que tengan un error
del mı́smo orden que los métodos de Taylor pero sin tener que calcular las derivadas
de la función f . Podemos ver que la fórmula recursiva del método anterior puede
verse escrita:
h
yi+1 = yi + h.(f (xi , yi ) + (fx (xi , yi ) + fy (xi , yi ).f (xi , yi )). )
2
En los métodos de Runge-Kutta se buscan funciones de la forma ϕ(x, y, h), que
no involucren el cálculo de las deribadas de f , para armar una fórmula recursiva
yi+1 = yi + h.ϕ(xi , yi , h)
El método que estudiaremos nosotros presenta un error del mismo orden que
el método de Taylor de orden 2, es por eso que se lo llama método de Runge -
Kutta de orden 2. Para motivarlo vamos a analizar la integral de y ′ entre dos puntos
consecutivos de la partición:
Z xi+1 Z xi+1
′
y dx = f (x, y) dx
xi xi
Z xi+1
y(xi+1 ) − y(xi ) = f (x, y) dx
xi

Ahora vamos a aproximar la integral de la derecha usando la regla del trapezoide.

 1.3. MÉTODOS DE RUNGE-KUTTA 9

Observación: La regla del trapezoide proviene de aproximar el área bajo la

curva de una función continua y positiva con el área de un trapecio, lo que da origen
a la igualdad:
Z b
b−a
f (x) dx = .(f (a) + f (b))
a 2

De esta manera llegamos a la aproximación:

h
y(xi+1 ) − y(xi ) ≈ .[f (xi , y(xi )) + f (xi+1 , y(xi+1 ))]
2
h
y(xi+1 ) ≈ y(xi ) + .[f (xi , y(xi )) + f (xi+1 , y(xi+1 ))]
2

El problema que observamos en esta aproximación es que para hacer las cuentas
del lado derecho aun no conocemos el valor de y(xi+1 ), ya que este es el valor que se
desea aproximar con yi+1 . Para salvar este problema vamos a estimar el valor y(xi+1 )
del lado derecho usando el método de Euler, es por este motivo que a este método
también se lo conoce como el método de Euler modificado.

Nos queda la siguiente fórmula recursiva:

h
yi+1 = yi + .[f (xi , yi ) + f (xi+1 , yi + h.f (xi , yi ))]
2

Método de Runge-Kutta de orden 2 o método de Euler modificado:

dada la ED y ′ (x) = f (x, y) con condición inicial y(a) = y0 , se pretende encontrar
una aproximación del valor y(b) con b > a.
1. Tomamos una partición a = x0 < x1 < x2 < ... < xn = b del intervalo [a, b],
b−a
con n ∈ N, y llamamos h = a la longitud del paso.
n
2. Llamamos y0 = y(x0 ) = y(a)
3. Calculamos de forma recursiva y1 , y2 , ..., yn de esta manera:
h
yi+1 = yi + .[f (xi , yi ) + f (xi+1 , yi + h.f (xi , yi ))]
2
4. yn es la aproximación de y(b).

Trabajemos ahora el mismo ejemplo que desarrollamos con el método de Euler

para ver como mejora la aproximación.
Ejemplo 1.2.
Sea y(x) la solución de la ecuación diferencial:
y ′ (x) = x + x.y(x) con condición inicial y(0) = 1.
Estimar y(1) con el método de Runge-Kutta de orden 2 y h = 41 .

La condición inicial es el punto (0, 1) y la fórmula recursiva nos queda:

h
yi+1 = yi + .[f (xi , yi ) + f (xi+1 , yi + h.f (xi , yi ))]
2
 1.3. MÉTODOS DE RUNGE-KUTTA 10

1 1
yi+1 = yi + .[xi .yi + xi + f (xi+1 , yi + .(xi .yi + xi ))]
8 4
1 1
yi+1 = yi + .[xi .yi + xi + xi+1 .(yi + .(xi .yi + xi )) + xi+1 ]
8 4

 1 1

 y1 = y0 + .[x0 .y0 + x0 + x1 .(y0 + .(x0 .y0 + x0 )) + x1 ]


 8 4


 1 17
1 1 1 1
y1 : y1 = 1 + .[0 . 1 + 0 + .(1 + .(0 . 1 + 0)) + ] =⇒ (x1 , y1 ) = ( , )

 8 4 4 4 4 16



 y1 = 17



16

 1 1

 y2 = y1 + .[x1 .y1 + x1 + x2 .(y1 + .(x1 .y1 + x1 )) + x2 ]


 8 4


 1 5177
17 1 1 17 1 1 17 1 1 17 1 1
y2 : y2 = + .[ . + + .( + .( . + )) + ] =⇒ (x2 , y2 ) = ( , )

 16 8 4 16 4 2 16 4 4 16 4 2 2 4096



 y2 = 5177



4096
Como podrán observar los resultados que vamos obteniendo son fracciones con
coeficientes muy alto, creemos que carece de sentido seguir mostrando las cuentas
tan detalladamente dada la complejidad de la escritura, ası́ que pasaremos a mostrar
cómo quedan los resultados finales.
3 1724051
(x3 , y3 ) = ( , )
4 1048576
613259859
(x4 , y4 ) = (1, )
268435456
613259859
El método de Runge-Kutta nos da la aproximación y(1) = ≈ 2, 284.
268435456
Lo cual es una mejor aproximación
√ que la que conseguı́mos con Euler ya que recor-
demos, el valor real es y(1) = 2 e − 1 ≈ 2, 297.

A continuación mostramos un cuadro en el cual podemos comparar los resultados

obtenidos al usar los métodos de Euler y Runge-Kutta y como mejora la aproxima-
ción al usar un h más pequeño.

h Euler Runge-Kutta Valor real

0,25 1,838867 2,284571 2,297442
0,1 2,094221 2,295763 2,297442
0,01 2,275641 2,297428 2,297442