Revista de Matema´tica: Teor´ıa y Aplicaciones 1998 5(1) : 73–86
cimpa – ucr – ccss issn: 1409-2433
generalizacio´n a Rn de algunos me´todos de
interpolacio´n conocidos en ecuaciones
diferenciales
Vernor Arguedas Troyo*— Roberto Mata Montero**
Recibido: 9 Octubre 1997 / Versio´n revisada: 11 Mayo 1998
Resumen
Presentamos el me´todo de integracio´n Runge-Kutta, de varios niveles de un paso
en Rn, as´ı como cuasi algoritmos para los me´todos impl´ıcito y expl´ıcito. Se estudia el
problema del control del error en Rn y se dan ejemplos de me´todos nume´ricos mediante
su respectiva tabla o esquema parame´trico.
Palabras clave: Me´todos de Runge-Kutta n-dimensionales, expl´ıcitos e impl´ıcitos, inter-
polacio´n, cuasi algoritmos, esquema parame´trico.
Abstract
We present the Runge-Kutta method of several one-step levels in Rn. as well as
algorithms in pseudo-code for the implicit and explicit methods. We study the problem
of error control in Rn and we give numerical examples in tables or parameter schemata.
Keywords: n-dimensional Runge-Kutta methods, explicit and implicit methods, interpo-
lation, pseudo-algorithms, parametric schema.
AMS Subject Classification: 34-04, 65L06
1. Introduccio´n
En los u´ltimos an˜os ha habido una tendencia, a la bu´squeda de me´todos nume´ricos,
que permitan una alta precisio´n o que reduzcan el tiempo de ma´quina, v´ıa simplificacio´n
de los algoritmos.
*Escuela de Matemtica, Universidad de Costa Rica, 2060 San Jose´, Costa Rica
**Sede Regional de Guanacaste, Universidad de Costa Rica, Liberia, Costa Rica
73
74 v. arguedas – r. mata
En este sentido, la teor´ıa de Runge-Kutta ha tenido un desarrollo muy grande; as´ı como
los me´todos de extrapolacio´n, interpolacio´n y mixtos.
Es importante saber diferenciar los tipos de interpolacio´n mencionados en libros y
art´ıculos y que pueden ser clasificados en:
i) Interpolacio´n esta´ndar o finita, como la de Newton, Hermite y otros.
ii) Interpolacio´n finita con polinomios mo´nicos.
iii) La interpolacio´n con polinomios mo´nicos de grado mı´nimo.
2. Fundamentos teo´ricos
Un me´todo Runge-Kutta, en forma matricial, es descrito de la siguiente manera:
Sean f : [a, b]×Rn −→ Rn y y : [a, b] −→ Rn, donde f posee (p+1) derivadas parciales
continuamente diferenciables y:
f =

f1
f2
...
fn
 x ∈ Rn, x =

x1
x2
...
xn
 y =

y1
y2
...
yn
 y′(t) =

y′1(t)
y′2(t)
...
y′n(t)
 (1)
Emplearemos la norma:
‖ x ‖∞= sup | xi | , para i = 1, · · · , n
Dado el problema de valor inicial:
y′(x) = f(x, y(x)),
y(x0) = η0, η0 ∈ Rn,
con:
yn+α := y(xn + αh) ∈ Rn;
xn := a+ hn ≤ b− h;
Definimos, de acuerdo con la literatura usual, el me´todo de integracio´n de un paso y p
niveles, el cual puede ser impl´ıcito o expl´ıcito, de la forma [1]:
η0 = η0(h) ( usualmente y0 = η0 )
ηn+αi = ηn + h
p∑
j=1
βijf(xn + αjh, ηn+αj ) ; i = 1, . . . , p
ηn+αs+1 = ηn + h
p∑
j=1
γjf(xn + αjh, ηn+αj ) ; n = 0, 1, . . .
generalizacio´n a IRn de algunos me´todos de interpolacio´n 75
con αj , βij y las γj constantes del me´todo, y con:
αi =
p∑
j=1
βij , i = 1, . . . , p ;
αp+1 =
p∑
i=1
γi
Donde yn+αi es la solucio´n exacta y ηn+αi es la solucio´n aproximada.
Esta notacio´n es la que usamos en [2].
Observemos que el usar ‖ ‖∞, nos permite trabajar con la funcio´n f : [a, b] −→ Rn,
de una manera muy sencilla, usando la funcio´n componente. De paso digamos que estamos
usando n veces el me´todo Runge-Kutta; y la tolerancia, as´ı como la aceptacio´n del h > 0,
es para todo el sistema de ecuaciones diferenciales.
El me´todo de integracio´n anterior, puede ser representado por el siguiente esquema
parame´trico.
α1 β11
α2 β21 β22
...
...
...
αp βp1 βp2 βp3 · · · βpp
αp+1 γ1 γ2 γ3 · · · γp
3. Cuasi algoritmos
A fin de lograr una mayor comprensio´n, del me´todo de integracio´n de varios niveles
de un paso y de facilitar su empleo, presentamos a continuacio´n, cuasi algoritmos para los
me´todos impl´ıcito y expl´ıcito.
3.1. Cuasi algoritmo para el me´todo impl´ıcito [10]
Calcula la solucio´n aproximada al problema de valor inicial:
y′(xn) = f(xn, y(xn)), con y(xn) ∈ V dado, V espacio vectorial de dimensio´n finita; para
f : [a, b] × V −→ V , con h variable, h > 0 y una tolerancia  > 0.
for i = 1 thru p
ηn+αi ←− 0
αi =
∑p
i=1 βij
ηn+αi ←− ηn + h
∑p
j=1 βijf(xn + αjh, ηn+αj )
{Observemos que el me´todo se llama a s´ı mismo. Lo que estamos
es buscando un “cero inteligente”, lo cual haremos por algu´n
me´todo nume´rico como gradiente descendiente, Newton, punto
fijo, despeje, o alguna modificacio´n de esos me´todos.}
repeat.
76 v. arguedas – r. mata
{ Tomar αp+1 =
∑p
i=1 γi}
for n = 0 thru p
ηn+αp+1 ←− 0
ηn+αp+1 ←− ηn ++h
∑p
j=1 γjf(xn + αjh, ηn+αj )
repeat.
end algoritmo.
En [2] realizamos un trabajo de este tipo, cuando, para el esquema trapezoidal impl´ıcito
en una dimensio´n:
ηi(xi + h) = ηi(xi) +
h
2
[f(xi, ηi(xi)) + f(xi + h, ηi(xi + h)],
tomamos su equivalente:
ηi(xi + h)− h2 · f(xi + h, ηi(xi + h)) = ηi(xi) +
h
2
· f(xi, ηi(xi)) (2)
donde:
ηi(xi) +
h
2
· f(xi, ηi(xi)) := β ∈ R
y definimos:
ϕ(u) := u− h
2
· f(xi + h, u)
Luego, si ϕ(u) = β admite una solucio´n u´nica y el sistema no lineal es computacionalmente
estable y fa´cil de resolver, ponemos u := ηi(xi + h) y retornamos a (2).
Aqu´ı, la idea del me´todo es ejemplarizada con:
f(xi, y(xi)) = xi · sen y(xi)
y′(xi) = xi · sen y(xi) , con y(0) = 1 y h = 0,01
x = (x1, · · · , xn)
Para xi ∈ [0, 1,5], con i = 1, · · · , n. (por satisfacer f las condiciones de Lipshitz en este
intervalo).
Al aplicar la regla del trapecio, obtenemos:
ηi(xi + h)− h2 [(xi + h) · sen(ηi(xi + h))] = ηi(xi) +
h
2
xi · sen(ηi(xi))
ηi(xi + h)− 10
−2
2
· sen(ηi(xi + h)) = 1 + 0,12 · 0 = 1
u− 10
−2
2
· sen(u)− 1 = 0
y, con la ayuda de MATLAB, queda u = 1.
generalizacio´n a IRn de algunos me´todos de interpolacio´n 77
3.2. Cuasi algoritmo para el me´todo expl´ıcito [10]
Calcula la solucio´n aproximada al problema de valor inicial:
y′(xn) = f(xn, y(xn)), con y(xn) ∈ V , para f : [a, b] × V −→ V , con h variable, h > 0 y
una tolerancia  > 0.
{ Tomar yn+1 = 0 }
for i = 1 thru p
xn,i ←− 0
yn,i ←− 0
yi ←− 0
repeat.
{ Tomar f0 = f(xn, yn) }
for j = 1 thru p
xn,j ←− xn + αjh
y0 ←− hγ0f0
yn,j ←− yn + h
∑j−1
k=0 βjkfk
fj ←− f(xn,j, yn,j)
yj ←− yj−1 + hγjfj
repeat.
yn+1 ←− yn + yp
end algoritmo.
Un ejemplo muy ilustrativo de los me´todos expl´ıcitos, es el me´todo Runge-Kutta-
Fehlberg, que fuera presentado en 1970 por Fehlberg [5]. El mismo es generalizado para
Rn, considerando (1), para:
k1 =

k11
k12
...
k1n
 · · · kn =

kn1
kn2
...
knn

con lo que se obtiene la te´cnica de usar el me´todo de Runge-Kutta, con error de trun-
camiento local de orden cinco:
yˆn+1 = yn +
16
135
k1 +
6656
12825
k3 +
28561
56430
k4 − 950k5 +
2
55
k6
para estimar el error local en el me´todo de Runge-Kutta de orden cuatro:
yn+1 = yn +
25
216
k1 +
1408
2565
k3 +
2197
4104
k4 − 15k5
78 v. arguedas – r. mata
donde:
k1 = hf(xn, yn)
k2 = hf(xn +
h
4
, yn +
1
4
k1)
k3 = hf(xn +
3h
8
, yn +
3
32
k1 +
9
32
k2)
k4 = hf(xn +
12h
13
, yn +
1932
2197
k1 − 72002197k2 +
7296
2197
k3)
k5 = hf(xn + h, yn +
439
216
k1 − 8k2 + 3680513 k3 −
845
4104
k4)
k6 = hf(xn +
h
2
, yn − 827k1 + 2k2 −
3544
2565
k3 +
1859
4104
k4 − 1140k5)
teniendo presente que:
hf(xn, yn) =

hf1(xn, yn)
hf2(xn, yn)
...
hfn(xn, yn)

y as´ı sucesivamente para las dema´s expresiones.
Una ventaja clara de este me´todo es que so´lo se requieren seis evaluaciones de f por
paso, mientras que los me´todos Runge-Kutta arbitrarios de orden cuatro y cinco usados
conjuntamente, requerir´ıan diez evaluaciones de f por paso, cuatro evaluaciones para el
me´todo de orden cuatro y seis evaluaciones adicionales para el me´todo de orden quinto.
El siguiente algoritmo usa el me´todo de Runge-Kutta-Fehlberg con control de error,
generalizado para Rn. El paso 9 se an˜ade para eliminar modificaciones grandes en el
taman˜o de paso. Esto se hace para evitar las pe´rdidas de tiempo ocasionadas por los
taman˜os de paso muy pequen˜os, en regiones con irregularidades en la derivada de y, y para
evitar taman˜os de paso grandes, que puedan resultar en la omisio´n de regiones sensibles
cercanas. En algunos casos el procedimiento del incremento del taman˜o de paso se omite
del algoritmo y el procedimiento del decrecimiento del taman˜o de paso se modifica para
que se incorpore solamente cuando sea necesario tener al error bajo control.
Algoritmo de Runge-Kutta-Fehlberg:
{Para aproximar la solucio´n del problema de valor inicial y′(xn) = f(xn, y(xn)), donde
y(xn) ∈ Rn, a ≤ xn ≤ b, y(a) = α, con α = (α1, · · · , αn), y′ : [a, b] −→ Rn, h ∈ Rn y la
norma ‖ ‖∞; con el error de truncamiento local dentro de cierta tolerancia dada.}
ENTRADA puntos extremos a,b; condicio´n inicial α; tolerancia TOL; taman˜o de paso
ma´ximo hma´x; taman˜o de paso mı´nimo hmı´n.
SALIDA x, y, h, donde y aproxima a y(t). Se emplea la funcio´n f por medio de las
componentes f1, · · · , fn y se usa el taman˜o de paso h o un mensaje de que el taman˜o de
paso mı´nimo fue excedido.
generalizacio´n a IRn de algunos me´todos de interpolacio´n 79
Paso 1 Tomar:
x = a
y = α
h = hma´x;
SALIDA (x, y).
Paso 2 Mientras que (x < b) seguir los Pasos 3 - 11.
Paso 3 Tomar:
k1 = hf(x, y) (se calcula para cada fi, con i = 1, . . . , n;)
k2 = hf(x+
1
4
h, y +
1
4
k1);
k3 = hf(x+
3
8
h, y +
3
32
k1 +
9
32
k2);
k4 = hf(x+
12
13
h, y +
1932
2197
k1 − 72002197k2 +
7296
2197
k3);
k5 = hf(x+ h, y +
439
216
k1 − 8k2 + 3680513 k3 −
845
4104
k4);
k6 = hf(x+
1
2
h, y − 8
27
k1 + 2k2 − 35442565k3 +
1859
4104
k4 − 1140k5).
Paso 4 Tomar:
R =‖ 1360k1 − 1284275k3 − 219775240k4 + 150k5 + 255k6 ‖∞ /h (se escoge R = O(h), donde:
O(h) =‖ yˆn+1 − yn+1 ‖∞ /h , el cual
es explicado posteriormente.)
Paso 5 Tomar δ = 0,84(TOL/R)1/4.
Paso 6 Si R ≤ TOL entonces seguir los pasos 7 y 8.
Paso 7 Tomar:
x = x+ h; (aproximacio´n aceptada).
y = y + 25216k1 +
1408
2565k3 +
2197
4104k4 − 15k5.
Paso 8 SALIDA (x, y, h).
Paso 9 Si δ ≤ 0,1 entonces tomar h = 0,1h
si no; si δ ≥ 4 entonces tomar h = 4h
si no tomar h = δh.
Paso 10 Si h > hma´x entonces tomar h = hma´x.
Paso 11 Si h < hmi´n entonces
SALIDA (“h mı´nimo excedido”);
(Procedimiento terminado sin e´xito).
PARAR.
Paso 12 (El procedimiento se completo´). [En la norma ‖ ‖∞ ]
PARAR.
Observemos que en el proceso anterior se esta´n manipulando los vectores de la manera
usual.
80 v. arguedas – r. mata
4. Control del error
A continuacio´n, algunas definiciones que nos van a permitir una mayor comprensio´n,
de las apreciaciones que haremos sobre el control del error.
Si se tiene el problema de valor inicial:
y′(x) = f(x, y(x))
y(x0) = y0 (3)
con f : [a, b]× Rn −→ Rn,
como en (1)
y tenemos el me´todo de un paso, para n dimensiones:
η0 := y0
para i = 1, · · ·
ηi+1 := ηi + hφ(xi, ηi;h)
xi+1 := xi + h
donde φ es la funcio´n generatriz del me´todo, con:
x ∈ Rh := {x0 + ih/i = 1, 2, · · ·}
para el me´todo η(x;h), definido por η(x;h) := ηi , cuando x = x0 + ih
Llamamos error local de truncamiento del me´todo ηi+1, a:
‖ y(xi + h)− y(xi)
h
− ηi(xi) ‖∞
El me´todo η se llama de orden p > 0, si su error de truncamiento local satisface:
‖ y(xi + h)− y(xi)
h
− ηi(xi) ‖∞≤ O(hp) , ∀f ∈ Cp
Observemos que y denota a la solucio´n exacta de (3).
El error global en el punto x, para el me´todo η, se define como: [3]
e(x;h) := η(x;h) − y(x)
Si definimos:
hn ∈ Hx := {x− x0
n
/n = 1, 2, · · ·}
tendremos que el me´todo η es convergente si:
l´ım
n→∞ e(x;hn) = 0 , para todo x ∈ [a, b]
Un teorema cla´sico garantiza que un me´todo de orden p es convergente de orden p,
cuando φ satisface las condiciones de Lipshitz, con respecto a y y y′ [3].
generalizacio´n a IRn de algunos me´todos de interpolacio´n 81
Un ana´lisis del error nos conduce a las siguientes consideraciones:
dado el problema de valor inicial:
y′(xn) = f(xn, y(xn))
con f : [a, b]× Rn −→ Rn
a ≤ xn ≤ b , con y(a) = α,
haremos una estimacio´n del error local para el me´todo de Euler:
y0 = α,
yn+1 = yn + hf(xn, yn),
el cual tiene error de truncamiento local de orden O(h) dado por:
τn+1 =
y(xn+1)− y(xn)
h
− f(xn, y(xn))
Ahora bien, el me´todo de Euler modificado,
ηˆ0 = α
ηˆn+1 = ηˆn +
h
2
[f(xn, ηˆn) + f(xn+1, ηˆn + hf(xn, ηˆn))],
tendra´ error de truncamiento local τˆn+1 de orden O(h2).
Si ηn ≈ y(xn) ≈ ηˆn, entonces:
y(xn+1)− ηn+1 = y(xn+1)− ηn − hf(xn, ηn)
≈ y(xn+1)− y(xn)− hf(xn, y(xn)))
= hτn+1
As´ı que:
τn+1 ≈ 1
h
[y(xn+1)− ηn+1]
=
1
h
[y(xn+1)− ηˆn+1] + 1
h
(ηˆn+1 − ηn+1)
≈ τˆn+1 + 1
h
(ηˆn+1 − ηn+1)
Puesto que τn+1 es O(h) y τˆn+1 = O(h2), la parte ma´s importante de τn+1 se debe a
(ηˆn+1 − ηn+1)/h, con lo que concluimos que :
τn+1 ≈ 1
h
(ηˆn+1 − ηn+1)
se puede usar para aproximar el error de truncamiento local del me´todo de Euler. Esta es
la idea subyacente, en la bu´squeda de pares de me´todos φ, τ , en donde uno es de orden
n y el otro es de orden n + 1.( Casualmente, esto es lo que se utiliza en el algoritmo de
Runge-Kutta-Fehlberg. )
82 v. arguedas – r. mata
La estimacio´n del error de truncamiento local de un me´todo de diferencia puede usarse
ventajosamente para aproximar el taman˜o de paso o´ptimo “n”, para controlar el error
global, como veremos a continuacio´n.
Supongamos que hay dos me´todos cano´nicos de diferencia disponibles, para aproximar
la solucio´n del problema de valor inicial:
y′(xn) = f(xn, y(xn))
con f : [a, b]× Rn −→ Rn
a ≤ xn ≤ b , con y(a) = α
y que uno de estos me´todos,
η0 = α
ηn+1 = ηn + hnφ(xn, hn, ηn),
tiene un error de truncamiento local τn+1 de orden O(hn), mientras que el otro me´todo:
ηˆ0 = α
ηˆn+1 = ηˆn + hnφˆ(xn, hn, ηˆn)
tiene un error de truncamiento local τˆn+1 de orden O(hn+1).
El mismo ana´lisis que se uso´ para el me´todo de Euler, nos permite obtener:
τn+1 ≈ 1
h
(ηˆn+1 − ηn+1)
Sin embargo, τn+1 es de orden O(hn); as´ı que existe una constante k tal que:
τn+1 ≈ khn.
Luego, de las relaciones anteriores se tiene que:
khn ≈ 1
h
(ηˆn+1 − ηn+1).
Ahora consideremos el error de truncamiento, reemplazando a h por qh, donde q es positivo
pero acotado superiormente y no cerca de cero. Si denotamos este error de truncamiento
como τn+1(qh), entonces:
τn+1(qh) ≈ k(qh)n = qn(khn) ≈ q
n
h
(ηˆn+1 − ηn+1).
Para acotar τn+1(qh) por , escojamos q tal que:
qn
h
| ηˆn+1 − ηn+1 |≈| τn+1(qh) |≤ ;
es decir, tal que:
q ≤
(
h
|ηˆn+1−ηn+1|
)1/n
.
generalizacio´n a IRn de algunos me´todos de interpolacio´n 83
Puesto que el usuario especifica una norma ‖ ‖ y una tolerancia de error τ , la bu´sque-
da de h termina cuando formemos una aproximacio´n que denominamos “est”, para el error
local, y que cumpla que:
‖ est ‖≤ τ
Si el “est”no es aprobado, reducimos el valor de h. Si es aprobado, entonces yn+1 es
aceptado y se escoge un h conveniente para el pro´ximo paso.
5. Aplicaciones
Como ejemplos de aplicaciones, presentamos a continuacio´n el esquema parame´trico
de algunos me´todos de uno y de varios pasos; todos en n dimensiones.
5.1. Me´todo Runge-Kutta de cuarto orden, de un paso
Este me´todo queda definido por la siguiente tabla:
1/2 1/2 0 0
1/2 0 1/2 0
1 0 0 1
5/6 1/3 1/3 1/6
con c0 = 1/6
5.2. Me´todo Runge-Kutta de orden tres, de un paso
Su esquema parame´trico es el siguiente:
1/2 1/2 0
1 -1 2
-1/2 -2/3 1/6
con c0 = 1/6
5.3. Me´todo de Heun
Este es un me´todo Runge-Kutta de orden dos, dado mediante el arreglo:
2/3 2/3
3/4 3/4
donde c0 = 1/4.
84 v. arguedas – r. mata
5.4. Me´todo Adams-Bashforth de tres pasos
Este me´todo tiene un error de truncamiento local τn+1 = 38y
(4)(µn)h3, y es dado por
la tabla:
1 1
2 2 0
7
12
−16
12
23
12
con c0 = 512 .
5.5. Me´todo Runge-Kutta-Ralston
Este es un me´todo Runge-Kutta expl´ıcito de orden cuatro, con error de truncamiento
de orden cinco, dado por:
0.4 0.4
0.45573726 0.29697760 0.15875966
1 0.21810038 -3.05096470 3.83286432
0.82523972 -0.55148053 1.20553547 0.17118478
con c0 = 0,17476028.
Podr´ıamos seguir ampliando la lista de esquemas; sin embargo, queremos terminar
brindando los siguientes ejemplos obtenidos de [6], y que fueron dados a conocer inicial-
mente, en [7] y [8].
5.6. El esquema regular
1/2 -
√
3/6 1/4 1/4 -
√
3/6
1/2 +
√
3/6 1/4 +
√
3/6 1/4
1/2 1/2
Otro esquema regular es dado por:
0 0 0 0
1/2 5/24 1/3 -1/24
1 1/6 2/3 1/6
1/6 2/3 1/6
5.7. El esquema irregular
1/2 -
√
15/10 5/36 2/9 -
√
15/15 5/36 -
√
15/30
1/2 5/36 +
√
15/24 2/9 5/36 -
√
15/24
1/2 +
√
15/10 5/36 +
√
15/30 2/9 +
√
15/15 5/36
5/18 4/9 5/18
generalizacio´n a IRn de algunos me´todos de interpolacio´n 85
con:
0 0 0
1/2 1/2 0
0 1
6. Conclusiones
La generalizacio´n de los me´todos Runge-Kutta a Rn fue cano´nica, as´ı como el ana´lisis
de la propagacio´n del error. En el estudio de los me´todos impl´ıcitos, debemos continuar
investigando, pues queda mucho por hacer al respecto.
Referencias
[1] Albrecht, P. (1987) “A new theoretical approach to Runge-Kutta methods”, SIAM
Journal Numerical Analysis 24(2).
[2] Arguedas, V.; Mata, R. (1992) “Un teorema sobre el orden de los me´todos de Runge-
Kutta”, Ciencias Matema´ticas, U.C.R. 3(1).
[3] Burden, R.; Faires, D. (1985) Ana´lisis Nume´rico. Editorial Iberoame´rica, Me´xico.
[4] Collatz, L. (1960) The Numerical Treatment of Differential Equations. Springer-
Verlag, Heidelberg.
[5] Fehlberg, E. (1970) “Klassische Runge-Kutta Formeln vierter und niedrigerer Ord-
nung mit Schrittweiten-Kontrolle und ihre Anwendung auf Wa¨rmeleitungs-probleme”,
Computing 6: 61–71.
[6] Hairer, E.; Iserles, A.; Sanz-Serna, J.M. (1990) “Equilibria of Runge-Kutta methods”,
Numerische Mathematik 58: 243–254.
[7] Kansy, K. (1973) “Elementare Fehlerdarstellung fu¨r Ableitungen bei der Hermite-
Interpolation”, Numer. Math. 21: 350–354.
[8] Mata, R. (1990) Algunos Aspectos Teo´ricos de los Me´todos Runge-Kutta. Tesis de
Licenciatura, Universidad de Costa Rica.
[9] Press, W.; Flannery, B.; Teukolsky, S.; Vetterling, W. (1987) Numerical Recipes, USA.
[10] Press, W.; Flannery, B.; Teukolsky, S.; Vetterling, W. (1987) Numerical Recipes,
Example Book (Fortran), USA.
[11] Shampine, L.F.; Gordon, M.K. (1975) Computer Solution of Ordinary Differential
Equations: the Initial Value Problem. W. H. Freeman, San Francisco.
[12] Shampine, L.F.; Watts, H.A. (1976) “Practical solution of ordinary differential equa-
tions by Runge-Kutta methods”, Rept. SAND 76-0585, Sand´ıa National Laboratories,
Albuquerque, N. M.
86 v. arguedas – r. mata
[13] Shampine, L. F. (1985) “Interpolation for Runge-Kutta methods”, SIAM Journal
Numerical Analysis 22(5).
[14] Stoer, J.; Burlisch, R. (1990) Numerische Mathematik 2. Tercera edicio´n, Springer-
Verlag, Berlin.
[15] Werner, H.; Schabach, R. (1972) Praktische Mathematik II. Springer-Verlag, Berlin.