Libro An
Libro An
Libro An
net/publication/321158466
CITATIONS READS
0 2,662
1 author:
José M. Gallardo
University of Malaga
40 PUBLICATIONS 703 CITATIONS
SEE PROFILE
Some of the authors of this publication are also working on these related projects:
Analysis and Numerical Simulation of fluids governed by non-ideal Equations of State View project
All content following this page was uploaded by José M. Gallardo on 20 December 2019.
1
C
A B
!1
h = 1/4 – p.6/64
ÍNDICE
Prefacio V
El presente texto recopila la materia que he impartido durante los últimos años en la asigna-
tura Análisis Numérico, del Grado en Matemáticas de la Universidad de Málaga. Se trata de un
curso básico sobre métodos numéricos para problemas de valor inicial, así como para problemas
de contorno. El trabajo está orientado principalmente a estudiantes de Matemáticas, ya que se
presta especial atención al rigor y a la justificación de los resultados presentados. Sin embargo,
también es perfectamente válido para estudiantes de Física o Ingeniería, ya que los resultados
teóricos se ilustran con multitud de ejemplos y aplicaciones.
El trabajo está estructurado en tres bloques principales. El primero de ellos se dedica al estudio
de los métodos unipaso para problemas de valor inicial. Partiendo de la construcción particular
de algunos métodos clásicos (Euler, Heun, punto medio, . . . ), se establece la forma general de
un método unipaso y se estudian sus propiedades fundamentales: consistencia, estabilidad, con-
vergencia y orden. Después se aplica la teoría general desarrollada al importante caso de los
métodos de Runge-Kutta. Por último, se introduce el concepto de estabilidad absoluta y se aplica
para el estudio de problemas rígidos.
El segundo capítulo se dedica a la teoría general de los métodos multipaso lineales, prestando
especial atención a las familias de métodos de Adams-Bashforth y Adams-Moulton, así como
a los métodos de tipo predictor-corrector. En particular, se estudia la consistencia, estabilidad,
convergencia y orden de métodos multipaso lineales, así como sus propiedades de estabilidad
absoluta.
El último capítulo se dedica al análisis de métodos numéricos para problemas de contorno,
comenzando con una breve estudio del método del tiro. La parte central del capítulo se dedica
al método de diferencias finitas, tanto para problemas unidimensionales como para problemas
dependientes de varias variables. También se hace una introducción a la formulación variacional
de problemas de contorno unidimensionales, así como al método de elementos finitos.
Al final de cada capítulo se pueden encontrar implementaciones en el lenguaje de progra-
mación Python de los distintos métodos y problemas estudiados. Se supone que el lector posee
conocimientos rudimentarios de Python: existen numerosos libros sobre este lenguaje que el lec-
tor interesado puede consultar para ampliar sus conocimientos; por ejemplo, las referencias [12],
[16] y [17] son muy recomendables. La idea de usar Python surgió de una iniciativa, refleja-
da en un proyecto de innovación docente, basada en sustituir el software comercial específico
(Matlab, Mathematica, etc.) utilizado en las asignaturas de Análisis Numérico por sofware libre.
Se eligió Python debido a su versatilidad, el gran número de bibliotecas matemáticas y gráficas
disponibles, y su amplio uso en el mundo del software científico. En particular, todos los códi-
gos presentados en este trabajo utilizarán la biblioteca numpy de cálculo científico y la biblioteca
matplotlib de representación gráfica; de forma más específica, se usará la biblioteca pylab, que
incluye a las dos anteriores bajo una sintaxis similar a la del software Matlab.
Un libro de este estilo no surge de la nada, sino que se apoya en una extensa bibliografía. Para
el lector interesado en ampliar los temas presentados, recomendamos consultar los textos de Ciar-
let [5], Henrici [11], Lambert [15], LeVeque [19] y Ortega-Poole [21], en los que se ha basado gran
parte del presente trabajo. En la Bibliografía pueden encontrarse otros trabajos complementarios.
Un resultado básico de la teoría de ecuaciones diferenciales es que el problema (1.1) posee una
única solución y ∈ C1 ( I ) en el intervalo I = [t0 , t0 + T ], T > 0. A lo largo del capítulo supondre-
mos siempre que se verifican las hipótesis anteriores sobre la función f (t, y).
Señalemos que los métodos y resultados que estudiaremos en este capítulo pueden aplicarse
también a sistemas de ecuaciones, esto es, al caso en que y ∈ Rn y f : I × Rn → Rn . Sin embargo,
por claridad en la exposición, la mayor parte del tiempo nos restringiremos al caso escalar (n = 1).
Nuestro objetivo consiste en construir una aproximación numérica de la solución del proble-
ma (1.1). Para ello, consideremos una discretización o partición del intervalo,
cuyos elementos denominaremos nodos. Por simplicidad, supondremos que la partición es uni-
forme, es decir, la diferencia entre dos nodos consecutivos es constante: tk+1 − tk = h, para
k = 0, . . . , N − 1; la constante h se denomina paso de malla. Un método numérico para resolver
el problema (1.1) es un algoritmo que proporciona una serie de valores {y0 , y1 , . . . , y N } (la solu-
ción numérica) de modo que cada yk es una aproximación del valor exacto de la solución en el
nodo tk :
yk ≈ y(tk ), k = 0, 1, . . . , N.
Nótese que el primer valor es siempre exacto, por tratarse de la condición inicial del problema:
y0 = y ( t0 ).
El error local cometido en cada nodo se define como
e k = | y ( x k ) − y k |, k = 0, 1, . . . , N.
También podemos considerar el error global, que involucra a toda la solución numérica, definido
como
e = máx ek = máx |y( xk ) − yk |.
k =0,1,...,N k =0,1,...,N
Notemos que el error e depende del paso de malla h; si queremos hacer resaltar esta dependen-
cia, escribiremos e(h) en lugar de e. Una propiedad deseable en todo método numérico es la
convergencia, que puede expresarse como
lı́m e(h) = 0.
h →0
A grandes rasgos, esto significa que mientras menor sea el paso de malla mejor será la aproxi-
mación a la solución exacta.
Teniendo que cuenta que y(t) es solución de la ecuación diferencial, podemos sustituir el valor
y0 (t0 ) por f (t0 , y(t0 )):
y1 = y0 + h f (t0 , y(t0 )).
Por último, haciendo uso de la condición inicial, deducimos que
y1 = y0 + h f ( t0 , y0 ).
1.1 El método de Euler 3
Esta es la primera iteración o paso del método de Euler, cuya idea geométrica se representa en la
figura 1.1. Es importante notar que el valor y1 puede determinarse, ya que t0 , y0 y f (t, y) son
datos del problema y h es un valor dado.
siendo y(t) cualquier solución de la ecuación diferencial y0 = f (t, y). Nótese que la función f (t, y) es conocida, mientras
que en general y(t) no lo es; por tanto, en la práctica es más conveniente estudiar la regularidad de f (t, y) que la de y(t).
4 Métodos unipaso para problemas de valor inicial
Figura 1.2: Publicación original del método de Euler. E342, Inst. Calc. Integralis 1768, §650.
1
deducimos la siguiente cota del error:
C
e1 6 Ch2 .
En la deducción de la expresión para y1 se utilizó el hecho de que y0 = y(t0 ). A la hora de
calcular la siguiente iteración del método de Euler, tenemos que tener en cuenta que, en general,
y1 6= y(t1 ); por tanto, no podemos razonar de manera similar a como hicimos anteriormente.
Consideremos entonces, en lugar de la solución del problema (1.1), la solución del problema con
condición inicial y(t1 ) = y1 : A B
!1
(
y0 (t) = f (t, y(t)), t ∈ I = [t1 , t0 + T ],
y ( t1 ) = y1 .
Aproximando mediante la recta tangente de forma análoga al primer paso, podemos construir la
siguiente aproximación de y(t2 ) a partir de y1 :
h = 1/4 – p.6/64
y2 = y1 + h f ( t1 , y1 ).
La dependencia continua de las soluciones de una ecuación diferencial respecto de las condicio-
nes iniciales nos permite asegurar que, para h suficientemente pequeño, el valor y2 va a ser una
buena aproximación del valor y(t2 ).
El proceso anterior puede repetirse para determinar la aproximación yk+1 a partir de yk ,
obteniéndose así la forma general del método de Euler:
y k +1 = y k + h f ( t k , y k ), k = 0, 1, . . . , N − 1.
Se trata de un método unipaso, en el sentido de que para construir la aproximación yk+1 en el nodo
tk+1 tan solo es necesario conocer la aproximación yk en el nodo anterior tk .
Procedamos ahora a hacer una estimación del error cometido en cada iteración. En primer
lugar, podemos escribir
y(tk+1 ) − yk+1 = y(tk+1 ) − (yk + h f (tk , yk ))
= y(tk+1 ) − (yk + h f (tk , yk )) ± (y(tk ) + h f (tk , y(tk )))
= y(tk+1 ) − (y(tk ) + h f (tk , y(tk ))) + y(tk ) − yk + h( f (tk , y(tk )) − f (tk , yk )).
Aplicando el teorema de Taylor y razonando como en la primera iteración del método, se tiene
que
h2
|y(tk+1 ) − (y(tk ) + h f (tk , y(tk )))| = |y00 (ck+1 )|, ck+1 ∈ (tk , tk+1 ).
2
1.1 El método de Euler 5
de donde
ek+1 6 Ch2 + (1 + hL)ek , k = 0, 1, . . . , N − 1.
Este tipo de cota se denomina recurrente, ya que para estimar el error cometido en un paso
es necesario conocer el error en el paso anterior. Es sencillo obtener una cota no recurrente,
simplemente aplicando la desigualdad anterior hasta llegar a e0 :
1 − (1 + hL)k+1 (1 + hL)k+1 − 1
ek+1 6 Ch2 = Ch .
1 − (1 + hL) L
ehL(k+1) − 1 e LT − 1
ek+1 6 Ch 6 Ch ,
L L
donde se ha tenido en cuenta que (k + 1)h 6 Nh = T para k 6 N − 1.
En resumen, si redefinimos la constante C como
1 e LT − 1
C= máx |y00 (t)| ,
2 t∈ I L
hemos deducido la siguiente cota de error:
ek 6 Ch, k = 0, 1, . . . , N.
Teorema 1.1. Supongamos que f ∈ C1 ( I × R) y sea y(t) la única solución del problema de Cauchy (1.1).
Sea {yk }kN=0 la solución numérica generada por el método de Euler:
(
y k +1 = y k + h f ( t k , y k )
, k = 0, 1, . . . , N − 1. (1.2)
t k +1 = t k + h
6 Métodos unipaso para problemas de valor inicial
Existe entonces una constante C > 0, independiente del paso de malla h, tal que
def
e(h) = máx |y(tk ) − yk | 6 Ch, k = 0, 1, . . . , N.
06 k 6 N
Observación 1.1. En la práctica, la condición e(h) 6 Ch se interpreta diciendo que «el error se
comporta aproximadamente como Ch», esto es, e(h) ≈ Ch. Teniendo esto en cuenta, supongamos
que se aplica el método de Euler con pasos de malla h1 y h2 , siendo h2 = h1 /2. En tal caso, se
tendría: )
e(h1 ) ≈ Ch1 e ( h1 ) Ch1 h e ( h1 )
⇒ ≈ = 1 = 2 ⇒ e ( h2 ) ≈ .
e(h2 ) ≈ Ch2 e ( h2 ) Ch2 h2 2
Es decir, al dividir por dos el paso de malla el error también se divide aproximadamente por dos.
En general, al dividir el paso de malla por un cierto valor M el error obtenido también se dividirá
aproximadamente por M. Aunque no sea esta la definición precisa del orden de un método (que
veremos más adelante), por ahora diremos, intuitivamente, que el método de Euler tiene orden
uno.
∂f
Observación 1.2. Supongamos que la función ∂y existe y está acotada en I × R. Fijemos t ∈ I.
Entonces, dados y, z ∈ R, por el teorema del valor medio podemos escribir
∂f
f (t, y) − f (t, z) = (t, w)(y − z),
∂y
para un cierto valor w. Por tanto,
∂f
| f (t, y) − f (t, z)| = (t, w)|y − z| 6 L|y − z|,
∂y
siendo
∂f
L = máx (t, y)) : t ∈ I, y ∈ R .
∂y
Esto proporciona una forma de determinar el valor de la constante de Lipschitz L en la práctica.
Por otra parte, nótese que
d ∂f ∂f
y0 (t) = f (t, y(t)) ⇒ y00 (t) = f (t, y(t)) ⇒ y00 (t) = (t, y(t)) + y0 (t) (t, y(t))
dt ∂t ∂y
∂f ∂f
⇒ y00 (t) = (t, y(t)) + f (t, y(t)) (t, y(t)),
∂t ∂y
de donde
∂f ∂f
00
|y (t)| 6 (t, y(t)) + | f (t, y(t))| (t, y(t)).
∂t ∂y
∂f ∂f
De esta forma, acotando (cuando ello sea posible) | f |, | ∂t | y | ∂y |, es posible determinar una cota
para el máximo de |y00 (t)|.
1.1 El método de Euler 7
y(0) = 1,
h
yk+1 = yk + (t2k − yk ), k = 0, 1, . . . , N − 1,
2
donde h = 1/N.
En la tabla 1.1 se muestran, con una precisión de seis cifras decimales, los valores obtenidos
con el método de Euler para N = 10 (paso de malla h = 0.1), así como los valores exactos de la
solución y los errores cometidos. En este caso, el error global es igual a e = 0.025697.
k tk yk y(tk ) ek
0 0.0 1.000000 1.000000 0.000000
1 0.1 0.950000 0.951394 0.001394
2 0.2 0.903000 0.906138 0.003138
3 0.3 0.859850 0.865044 0.005194
4 0.4 0.821357 0.828885 0.007527
5 0.5 0.788290 0.798395 0.010105
6 0.6 0.761375 0.774272 0.012897
7 0.7 0.741306 0.757183 0.015877
8 0.8 0.728741 0.747760 0.019019
9 0.9 0.724304 0.746603 0.022299
10 1.0 0.728589 0.754285 0.025697
Tabla 1.1: Método de Euler: solución del ejemplo 1.1 con paso h = 0.1.
Nótese que el error aumenta conforme nos aproximamos al extremo superior del intervalo; esto
es debido a la acumulación de errores en las distintas aproximaciones mediante rectas tangentes.
En la figura 1.3 se representan la solución exacta del problema y la solución numérica obteni-
da.
Apliquemos ahora el método de Euler con paso de malla h = 0.05, es decir, la mitad del paso
de malla anterior. En tal caso, se obtiene que
e(0.05) = 0.012830,
y por tanto,
e(0.1) e(0.1)
= 2.002792 ≈ 2 ⇒ e(0.05) ≈ ,
e(0.05) 2
8 Métodos unipaso para problemas de valor inicial
es decir, al dividir entre dos el paso de malla el error cometido se divide aproximadamente entre
dos. Al realizar el experimento un cierto número de veces se obtienen los resultados de la tabla
1.2, donde en la última columna se han escrito los cocientes entre dos errores consecutivos2 .
Como puede verse, cada vez que dividimos por dos el paso de malla, el error correspondiente se
divide aproximadamente por dos, lo que indica que el método de Euler tiene orden uno.
N h e cocientes
10 0.100000 0.025697 −
20 0.050000 0.012830 2.002792
40 0.025000 0.006410 2.001471
80 0.012500 0.003204 2.000753
160 0.006250 0.001602 2.000381
Por último, en la tabla 1.3 se muestran los errores obtenidos al aplicar el método de Euler con
una serie de pasos de malla que tienden a cero. Como puede observarse, dichos errores se van
aproximando a cero, lo que indica la convergencia del método.
Asimismo, en la figura 1.4 se representan la solución exacta y las soluciones aproximadas con
pasos de malla 10−1 y 10−2 (el resto de las soluciones son indistinguibles de la solución exacta
en la figura). Observamos que mientras menor es el paso de malla, más se aproxima la solución
numérica a la solución exacta.
2 Aunque en la tabla se representan los valores redondeados a seis cifras decimales, todas las operaciones se realizan
N h e
10 10−1 2.5697 · 10−2
102 10−2 2.5630 · 10−3
103 10−3 2.5623 · 10−4
104 10−4 2.5626 · 10−5
105 10−5 2.5622 · 10−6
106 10−6 2.5623 · 10−7
En la sección 1.9.1 se realiza una implementación en el lenguaje Python del método de Euler
para una ecuación escalar, y en el apartado 1.9.3 se hace para un sistema.
tras aplicar la regla de Barrow en el primer miembro. Por tanto, obtenemos así la formulación
integral de la ecuación diferencial:
Z t
k +1
y ( t k +1 ) = y ( t k ) + f (t, y(t)) dt. (1.3)
tk
10 Métodos unipaso para problemas de valor inicial
Aplicando la regla del rectángulo3 para aproximar la integral de la función f (t, y(t)) en [tk , tk+1 ],
resulta
y k +1 = y k + h f ( t k , y k ),
donde yk e yk+1 denotan, respectivamente, aproximaciones de y(tk ) e y(tk+1 ). Como puede ob-
servarse, la expresión obtenida es precisamente la del método de Euler.
El desarrollo anterior nos permite plantear el uso de fórmulas de cuadratura más precisas para
construir diferentes métodos numéricos. El planteamiento general de dicha técnica se realizará
en el apartado 1.4; en esta sección nos limitaremos a introducir una serie de métodos clásicos.
En primer lugar, el método de Heun o de Euler mejorado se construye al aplicar la regla del
trapecio4 para aproximar el término integral en (1.3). En tal caso,
Z t
k +1 h
f (t, y(t)) dt ≈ f (tk , y(tk )) + f (tk+1 , y(tk+1 )) .
tk 2
A continuación, el valor y(tk+1 ) que aparece en la expresión anterior se aproxima mediante una
iteración del método de Euler: y(tk+1 ) ≈ y(tk ) + h f (tk , y(tk )). Por último, y teniendo en cuenta
que tk+1 = tk + h, a partir de la formulación (1.3) se obtiene la forma general del método de
Heun:
h
y k +1 = y k + f (tk , yk ) + f (tk + h, yk + h f (tk , yk )) , k = 0, 1, . . . , N − 1.
2
El método de Heun puede escribirse también como
p k +1 = y k + h f ( t k , y k ),
h
y
k +1 = y k + f (tk , yk ) + f (tk + h, pk+1 ) .
2
Visto de esta forma, se trata de un caso particular de método predictor-corrector, en el que la
«predicción» pk+1 dada por el método de Euler se «corrige» para obtener yk+1 (véase la sección
2.5).
Puede demostrarse que, bajo hipótesis de regularidad adecuadas, existe una constante C > 0,
independiente de h, tal que
e(h) 6 Ch2 .
Esta desigualdad puede deducirse de forma similar a como se hizo para el método de Euler.
Sin embargo, no detallaremos dicho proceso aquí, ya que en la sección 1.3 veremos un resultado
general para determinar el orden de un método unipaso.
3 La regla del rectángulo para aproximar la integral de una función g : [ a, b] → R tiene la siguiente forma:
Z b
g( x ) dx ≈ (b − a) g( a).
a
El grado de exactitud de dicha regla es cero, es decir, tan solo integra de forma exacta las funciones constantes.
4 La regla del trapecio tiene la forma
Z b
g( a) + g(b)
g( x ) dx ≈ (b − a) ,
a 2
siendo su grado de exactitud uno, esto es, la fórmula integra exactamente los polinomios de grado menor o igual que
uno.
1.2 Otros métodos clásicos 11
Ejemplo 1.2. Retomemos el problema del ejemplo 1.1. Al aplicar el método de Heun para resol-
verlo con paso de malla h = 0, 1, se obtienen los resultados de la tabla 1.4 (compárese con los
resultados de la tabla 1.1).
k tk yk y(tk ) ek
0 0.0 1.000000 1.000000 0.000000
1 0.1 0.951500 0.951394 1.0597 · 10−4
2 0.2 0.906352 0.906138 2.1380 · 10−4
3 0.3 0.865367 0.865044 3.2306 · 10−4
4 0.4 0.829318 0.828885 4.3334 · 10−4
5 0.5 0.798939 0.798395 5.4429 · 10−4
6 0.6 0.774928 0.774272 6.5559 · 10−4
7 0.7 0.757950 0.757183 7.6693 · 10−4
8 0.8 0.748638 0.747760 8.7805 · 10−4
9 0.9 0.747591 0.746603 9.8870 · 10−4
10 10.0 0.755384 0.754285 1.0987 · 10−3
Tabla 1.4: Método de Heun: solución del ejemplo 1.1 con paso h = 0.1.
En la figura 1.5 se comparan gráficamente los resultados obtenidos con los métodos de Euler
y Heun.
Los resultados obtenidos al dividir sucesivamente entre dos el paso de malla se muestran en
la tabla 1.5 . Como puede observarse, el cociente entre errores sucesivos tiende a 22 = 4.
12 Métodos unipaso para problemas de valor inicial
N h e cocientes
10 0.100000 1.0987 · 10−3 −
20 0.050000 2.7270 · 10−4 4.028877
40 0.025000 6.7926 · 10−5 4.014620
80 0.012500 1.6950 · 10−5 4.007353
160 0.006250 4.2337 · 10−6 4.003687
Al aplicar la fórmula de cuadratura del punto medio5 en (1.3) se obtiene la siguiente cadena
de aproximaciones:
donde el valor y(tk + 2h ) se ha aproximado mediante una iteración del método de Euler con paso
de malla h/2. La expresión anterior da lugar al denominado método del punto medio:
yk+1 = yk + h f tk + 2h , yk + 2h f (tk , yk ) , k = 0, 1, . . . , N − 1.
En el apartado 1.3 demostraremos la existencia de una constante C > 0 tal que e(h) 6 Ch2 .
Por otra parte, los resultados obtenidos al aplicar el método del punto medio al problema del
ejemplo 1.1 son similares a los del método de Heun (ejemplo 1.2), por lo que los omitiremos aquí.
El último método que vamos a definir en esta sección es el método de Runge-Kutta de cuarto
orden, que abreviaremos como RK4. En el apartado 1.4 estudiaremos la clase general de métodos
Runge-Kutta, de los cuales RK4 es un caso particular. La forma del método es la siguiente:
h
yk+1 = yk + (K1 + 2K2 + 2K3 + K4 ), k = 0, 1, . . . , N − 1,
6
donde
K1 = f ( t k , y k ) ,
K2 = f ( t + h , y + h K1 ) ,
k 2 k 2
K3 = f (tk + 2h , yk + 2h K2 ),
K4 = f (tk + h, yk + hK3 ).
Nótese que los coeficientes Ki dependen de k, aunque dicha dependencia no suele explicitarse.
Es posible interpretar el método RK4 a partir de una aplicación de la regla de Simpson6 para
5 La regla del punto medio, que tiene orden de exactitud uno, tiene la siguiente forma:
Z b
a+b
g( x ) dx ≈ (b − a) g .
2
a
aproximar la integral en el segundo miembro de 1.3. Veremos en la sección 1.3 que el método
RK4 tiene orden cuatro, lo que implica la existencia de una constante C > 0 tal que
e(h) 6 Ch4 .
Ejemplo 1.3. Al aplicar el método de Runge-Kutta con paso de malla h = 0.1 al problema del
ejemplo 1.1, se obtienen los resultados de la tabla 1.6. Es claro que el método de Runge-Kutta
produce resultados mucho más precisos que los métodos de Euler o de Heun para el mismo paso
de malla.
k tk yk y(tk ) ek
0 0.0 1.000000 1.000000 0.000000
1 0.1 0.951394 0.951394 7.9633 · 10−9
2 0.2 0.906138 0.906138 1.6420 · 10−8
3 0.3 0.865044 0.865044 2.5303 · 10−8
4 0.4 0.828885 0.828885 3.4550 · 10−8
5 0.5 0.798395 0.798395 4.4106 · 10−8
6 0.6 0.774273 0.774272 5.3917 · 10−8
7 0.7 0.757183 0.757183 6.3936 · 10−8
8 0.8 0.747760 0.747760 7.4120 · 10−8
9 0.9 0.746603 0.746603 8.4429 · 10−8
10 10.0 0.754285 0.754285 9.4825 · 10−8
Tabla 1.6: Método de Runge-Kutta: solución del ejemplo 1.1 con paso h = 0.1.
N h e cocientes
10 0.100000 9.4825 · 10−8 −
20 0.050000 5.9132 · 10−9 16.036346
40 0.025000 3.6911 · 10−10 16.020039
80 0.012500 2.3053 · 10−11 16.011356
160 0.006250 1.4405 · 10−12 16.003314
Asimismo, al realizar el experimento de divisiones sucesivas del paso de malla, se obtienen los
resultados de la tabla 1.7. Dichos resultados ilustran el hecho de que el método de Runge-Kutta
tiene orden cuatro.
Finalizamos esta sección con un experimento donde se comparan los métodos de Euler, Heun,
punto medio y RK4.
Su grado de exactitud es tres, esto es, integra exactamente los polinomios de grado menor o igual que tres.
14 Métodos unipaso para problemas de valor inicial
En las figuras 1.7 y 1.8 se han representado las soluciones numéricas obtenidas con N = 25 y
N = 50 particiones, respectivamente.
1.2 Otros métodos clásicos 15
h
Yk+1 = Yk + (K1 + 2K2 + 2K3 + K4 ), k = 0, 1, . . . , N − 1,
6
donde
K1 = F (tk , Yk ),
K2 = F ( t + h , Y + h K1 ) ,
k 2 k 2
K3 = F (tk + 2h , Yk + 2h K2 ),
K4 = F (tk + h, Yk + hK3 ),
donde ahora Ki ∈ Rn .
La implementación en Python de los métodos estudiados en esta sección se realiza en los
apartados 1.9.2 y 1.9.3 para ecuaciones escalares y sistemas, respectivamente.
y k +1 = y k + h Φ ( t k , y k ; h ), k = 0, 1, . . . , N − 1,
Observación 1.3. El término unipaso hace referencia a que el valor yk+1 se determina usando
únicamente el valor determinado en el paso anterior, yk .
Observación 1.4. Todos los métodos presentados en el apartado 1.2 son métodos unipaso. En
efecto, las correspondientes funciones incremento son:
Φ(t, y; h) = f (t, y)
1
Φ(t, y; h) = f (t, y) + f (t + h, y + h f (t, y))
2
para el método de Heun, y
Φ(t, y; h) = f (t + 2h , y + 2h f (t, y))
para el método del punto medio. Por último, la función incremento dada por
1
Φ(t, y; h) = (K + 2K2 + 2K3 + K4 ),
6 1
1.3 Teoría general de los métodos unipaso 17
donde
K1 = f (t, y),
K2 = f ( t + h , y + h K1 ) ,
2 2
K3 = f (t + 2h , y + 2h K2 ),
K4 = f (t + h, y + hK3 ),
donde Φ : I × Rn × [0, h∗ ] → Rn .
A continuación daremos una serie de definiciones que son fundamentales a la hora de estudiar
un método unipaso, y estableceremos diversas relaciones entre ellas.
Dada un solución y(t) de la ecuación diferencial y0 = f (t, y), se define el error de discretización
local ε k como
ε k = y ( t k +1 ) − y ( t k ) + h Φ ( t k , y ( t k ); h ) .
Es decir, el error ε k mide la diferencia entre el valor exacto de la solución en el nodo tk+1 y el
resultado de aplicar el método unipaso partiendo del valor exacto y(tk ) en lugar de yk . Intuiti-
vamente, ε k nos da una indicación de «cuánto se parece» el valor proporcionado por el método
unipaso, en el nodo tk+1 , al valor de la solución en dicho nodo.
Definición 1.1 (Consistencia). Diremos que un método unipaso es consistente con la ecuación diferen-
cial y0 = f (t, y) si, para cualquier solución y(t) de dicha ecuación, se verifica:
lı́m ∑
h →0 k < N
|ε k | = 0.
Definición 1.2 (Estabilidad). Un método unipaso es estable si existe una constante M > 0, indepen-
diente de h, tal que para cualesquiera sucesiones {yk }kN=0 , {zk }kN=0 y {δk }kN=−01 que verifiquen
y k +1 = y k + h Φ ( t k , y k ; h ),
zk+1 = zk + h Φ(tk , zk ; h) + δk ,
se tiene que
máx |yk − zk | 6 M |y0 − z0 | +
06 k 6 N
∑ |δk | .
k< N
lı́m e(h) = 0,
h →0
donde e(h) = máx06k6 N |y(tk ) − yk |, siendo {yk }kN=0 la solución numérica generada por el método con
paso de malla h.
18 Métodos unipaso para problemas de valor inicial
Observación 1.6. La condición de consistencia nos permite asegurar que la solución numérica,
en caso de converger, lo hace a una solución de la ecuación diferencial considerada.
La estabilidad de un método permite tener bajo control las posibles perturbaciones o errores
en los datos: si estas perturbaciones δk son pequeñas, el error |yk − zk | cometido en cada paso
también lo será. Dicho de otra forma, un método estable no amplifica los errores. Este hecho es
particularmente importante a la hora de implementar el método en el ordenador, debido a la
existencia inevitable de errores de redondeo.
Por último, la noción de convergencia nos dice que conforme el paso de malla tiende a cero,
las soluciones numéricas convergen, en un cierto sentido, a la solución exacta del problema.
donde y(t) es una solución arbitraria de la ecuación diferencial y0 = f (t, y). Entonces, por una
parte, la condición (1.4) implica que
lı́m
h →0 k < N
∑ |ε k | = 0,
1.3 Teoría general de los métodos unipaso 19
Entonces,
de donde, llamando αk = h( f (ck , y(ck )) − Φ(ck , y(ck ); 0)) y β k = Φ(ck , y(ck ); 0) − Φ(tk , y(tk ); h),
podemos escribir
ε k = αk + hβ k .
Por un lado, se verifica que
Z t0 + T
lı́m ∑
h →0 k < N
|αk | = lı́m ∑
h →0 k < N
h| f (ck , y(ck )) − Φ(ck , y(ck ); 0)| =
t0
| f (t, y(t)) − Φ(t, y(t); 0)| dt,
por ser la función t 7→ | f (t, y(t)) − Φ(t, y(t); 0)| continua, y por tanto, integrable Riemann (nótese
que el segundo término en la expresión anterior es una suma de Riemann asociada a dicha
función). Por otra parte, consideremos la función β : [0, h∗ ] → R dada por
La continuidad uniforme de la función (t, h) 7→ Φ(t, y(t); h) sobre el compacto I × [0, h∗ ] implica
que
lı́m β(h) = 0.
h →0
En consecuencia,
lı́m ∑
h →0 k < N
h| β k | = 0.
20 Métodos unipaso para problemas de valor inicial
y por tanto,
Z t0 + T
lı́m ∑
h →0 k < N
|ε k | = lı́m ∑
h →0 k < N
|αk | =
t0
| f (t, y(t)) − Φ(t, y(t); 0)| dt.
Teorema 1.4. Supongamos que la función Φ(t, y; h) es de Lipschitz respecto de la variable y, esto es, existe
una constante Λ > 0 tal que
|yk+1 − zk+1 | 6 |δk | + ehΛ |δk−1 | + e2hΛ |δk−2 | + · · · + ekhΛ |δ0 | + e(k+1)hΛ |y0 − z0 |.
Teniendo en cuenta que e jhΛ 6 e TΛ para j = 0, 1, . . . , k + 1,
k
TΛ TΛ
| y k +1 − z k +1 | 6 e | y0 − z0 | + ∑ | δj | 6 e | y0 − z0 | + ∑ | δj | ,
j =0 j< N
de donde
máx |yk − zk | 6 M |y0 − z0 | +
06 k 6 N
∑ |δk | ,
k< N
∗
Observación 1.7. Supongamos que la función ∂Φ ∂y existe y está acotada en I × R × [0, h ]. Razo-
nando de forma análoga a como hicimos para determinar la constante de Lipschitz L de f (t, y)
(observación 1.2), podemos tomar
∂Φ ∗
Λ = máx (t, y; h) : t ∈ I, y ∈ R, h ∈ [0, h ] .
∂y
Podemos así determinar en la práctica el valor de la constante de Lipschitz Λ.
1.3 Teoría general de los métodos unipaso 21
La combinación de los teoremas 1.2, 1.3 y 1.4 nos proporciona una herramienta para deter-
minar de forma sencilla la convergencia de un método unipaso, que se resume en el siguiente
resultado.
Corolario 1.1. Si la función incremento Φ(t, y; h) verifica Φ(t, y; 0) = f (t, y) para cada y ∈ R, t ∈ I, y
es de Lipschitz en la variable y, entonces el método unipaso asociado a Φ es convergente.
Ejemplo 1.5. La función incremento del método de Euler es Φ(t, y; h) = f (t, y); en particular, la
igualdad se verifica para h = 0. Por otra parte, Φ(t, y; h) es de Lipschitz en y, por serlo la propia
f . Así pues, el corolario 1.1 asegura la convergencia del método de Euler.
Consideremos ahora la función incremento del método de Heun:
1
Φ(t, y; h) = f (t, y) + f (t + h, y + h f (t, y)) .
2
Al evaluar en h = 0 se recupera la función f (t, y), por lo que el método es consistente. Por otra
parte, dados y, z ∈ R, se tiene que
1
|Φ(t, y; h) − Φ(t, z; h)| 6 | f (t, y) − f (t, z)| + | f (t + h, y + h f (t, y)) − f (t + h, z + h f (t, z))|
2
1
6 L|y − z| + L(|y − z| + h| f (t, y) − f (t, z)|)
2
2L + hL2 2L + TL2
6 |y − z| 6 | y − z |,
2 2
2
lo que prueba que Φ es de Lipschitz en y, con constante Λ = 2L+2TL . (También podemos obtener
este resultado basándonos en la observación 1.77 ). Aplicando el corolario 1.1, podemos deducir
que el método de Heun es convergente.
Razonamientos análogos a los anteriores permiten demostrar la convergencia de los métodos
del punto medio y RK4.
Definición 1.4 (Orden). Un método unipaso tiene orden (de consistencia) p > 1 si existe una constante
C, independiente del paso de malla h, tal que
∑ |ε k | 6 Ch p ,
k< N
∑ | ε k | = O ( h p ).
k< N
Nótese que la igualdad anterior esconde realmente una desigualdad, y que no es necesario expli-
citar, en general, la constante C.
Dados p, q ∈ N, se verifican las siguientes propiedades:
Si f (h) = O(h p ) y g(h) = O(hq ) entonces ( f + g)(h) = O(hmı́n( p,q) ); brevemente, escribire-
mos O(h p ) + O(hq ) = O(hmı́n( p,q) ).
Si f (h) = O(h p ) y g(h) = O(hq ) entonces ( f · g)(h) = O(h p+q ); usualmente, escribiremos
O ( h p ) · O ( h q ) = O ( h p + q ).
|y(tk ) − yk | 6 ch p , k = 0, 1, . . . , N,
donde c = MC.
1.3 Teoría general de los métodos unipaso 23
∂ f ( j) ∂ f ( j)
f (0) (t, y) = f (t, y), f ( j+1) (t, y) = (t, y) + f (t, y) (t, y).
∂t ∂y
Puede comprobarse que se verifica trivialmente
Al ser y( p+1) continua en I, está acotada; por eso, el último sumando del desarrollo de Taylor
será O(h p+1 ). Usando la notación introducida anteriormente, tenemos que
h 2 (1) h p ( p −1)
y(tk+1 ) − y(tk ) = h f (tk , y(tk )) + f (tk , y(tk )) + · · · + f (tk , y(tk )) + O(h p+1 ).
2! p!
Por otra parte, se tiene el siguiente desarrollo de Taylor de la función h 7→ Φ(tk , y(tk ); h) en el
punto h = 0:
∂Φ h p −1 ∂ p −1 Φ
Φ ( t k , y ( t k ); h ) = Φ ( t k , y ( t k ); 0) + h ( t k , y ( t k ); 0) + · · · + ( t , y ( t k ); 0) + O ( h p ).
∂h ( p − 1)! ∂h p−1 k
Combinando las expresiones anteriores, deducimos:
ε k = y ( t k +1 ) − y ( t k ) − h Φ ( t k , y ( t k ); h )
2 1 (1)
∂Φ
= h f (tk , y(tk )) − Φ(tk , y(tk ); 0) + h f (tk , y(tk )) − ( t , y ( t k ); 0) + · · ·
2 ∂h k
1 ( p −1) 1 ∂ p −1 Φ
· · · + hp f (tk , y(tk )) − ( t , y ( t ) ; 0 ) + O ( h p +1 ).
p! ( p − 1)! ∂h p−1 k k
donde
1 1 ∂j Φ
Ψ j (t, y) = f ( j) (t, y) − (t, y; 0).
( j + 1) ! j! ∂h j
Nótese que las igualdades Ψ j ≡ 0 para j = 0, . . . , p − 1, equivalen a las condiciones (1.7).
Supongamos que se verifican las condiciones (1.7), de donde Ψ j ≡ 0 para j = 0, . . . , p − 1.
Entonces, usando (1.8), tenemos que
ε k = hs+1 Ψs (tk , y(tk )) + O(hs+2 ) ⇒ |ε k | = hs+1 |Ψs (tk , y(tk ))| + O(hs+2 ) ⇒
1
∑ |ε k | = hs ∑ h|Ψs (tk , y(tk ))| + O(hs+1 ) ⇒ hs ∑ |ε k | = ∑ h|Ψs (tk , y(tk ))| + O(h),
k< N k< N k< N k< N
1.3 Teoría general de los métodos unipaso 25
donde hemos usado que ∑k< N O(hs+2 ) = TO(hs+1 ) = O(hs+1 ). Al ser el método de orden p, se
tiene que ∑k< N |ε k | = O(h p ), por lo que
1 1
∑ |ε k | = O(h p−s ) ⇒ lı́m ∑ |ε k | = 0,
hs k< N h →0 hs k< N
En consecuencia,
Z t0 + T
|Ψs (t, y(t))| dt = 0 ⇒ Ψs (t, y(t)) = 0, ∀ t ∈ I.
t0
Razonando como en la demostración del teorema 1.3, podemos deducir que Ψs ≡ 0, lo que
supone una contradicción.
Observación 1.9. Las condiciones (1.7) nos dicen realmente que el orden del método es al menos
p. Si además se verifica que
∂pΦ 1
p (t, y; 0) 6= f ( p) (t, y),
∂h p+1
entonces el orden del método es exactamente p.
Observación 1.10. La condición Φ(t, y; 0) = f (t, y) equivale, por una parte, a tener orden uno,
y por otra, a que el método sea consistente (teorema 1.3). Por tanto, la condición de consistencia
equivale a decir que el método tiene orden uno.
Observación 1.11. En la demostración del teorema 1.6 se ha visto que el error de discretización
local verifica
ε k = O ( h p +1 ),
para un método de orden p.
Ejemplo 1.6. La función incremento del método de Euler es Φ(t, y; h) = f (t, y). Entonces,
∂Φ 1
Φ(t, y; 0) = f (t, y), (t, y; 0) = 0 6= f (1) (t, y),
∂h 2
lo que prueba que el orden del método es exactamente uno.
Consideremos ahora la función incremento del método del punto medio:
Trivialmente, Φ(t, y; 0) = f (t, y); por otra parte, es inmediato comprobar que
∂Φ 1
(t, y; 0) = f (1) (t, y).
∂h 2
26 Métodos unipaso para problemas de valor inicial
∂2 Φ 1 1
2
(t, y; 0) = f (2) (t, y) 6= f (2) (t, y),
∂h 4 3
el orden del método del punto medio es exactamente dos.
De forma análoga, puede demostrarse que el orden del método de Heun es exactamente dos.
También puede demostrarse que el método RK4 tiene orden cuatro; sin embargo, en la sección
1.4 veremos una forma más simple de demostrar este hecho.
El teorema 1.6 nos permite construir métodos unipaso de orden arbitrario. Para ello, supon-
gamos que la función f (t, y) es de clase C p ( I × R) y consideremos el desarrollo de Taylor
h2 00 hp
y(tk+1 ) = y(tk ) + hy0 (tk ) + y ( t k ) + · · · + y ( p ) ( t k ) + O ( h p +1 )
2! p!
h 2 (1) h p ( p −1)
= y(tk ) + h f (tk , y(tk )) + f (tk , y(tk )) + · · · + f (tk , y(tk )) + O(h p+1 ).
2! p!
donde ci ∈ [0, 1]. Para aproximar los valores y(tk + ci h) consideramos de nuevo la formulación
integral, esta vez en el intervalo [tk , tk + ci h]:
Z t +c h
k i
y ( t k + ci h ) = y ( t k ) + f (t, y(t)) dt.
tk
De nuevo debemos aproximar el término integral. Para ello, consideramos una fórmula de cua-
dratura basada en los mismos nodos tk + c j h:
Z t +c h q
k i
f (t, y(t)) dt ≈ h ∑ aij f (tk + c j h, y(tk + c j h)).
tk j =1
De este modo,
q
y(tk + ci h) ≈ y(tk ) + h ∑ aij f (tk + c j h, y(tk + c j h)), i = 1, . . . , q.
j =1
Sustituyendo los valores exactos por valores aproximados en el procedimiento anterior, obte-
nemos un método de Runge-Kutta (RK) de q etapas (o pasos):
q
(i ) (i )
y k + 1 = y k + h ∑ bi f ( t k , y k ) , (1.9)
i =1
(i ) (i ) q
donde tk = tk + ci h y los elementos {yk }i=1 se determinan resolviendo el sistema de q ecuacio-
nes dado por
q
(i ) ( j) ( j)
yk = yk + h ∑ aij f (tk , yk ), i = 1, . . . , q. (1.10)
j =1
Nótese que en cada iteración del método hay que resolver un sistema del tipo (1.10).
Observación 1.12. Las fórmulas de cuadratura utilizadas en la definición de los métodos RK
deben integrar de forma exacta, al menos, las funciones constantes (su orden de exactitud debe
ser mayor o igual que cero). Por tanto, aplicando las fórmulas a la función f ≡ 1, deducimos por
una parte que
q
∑ bi = 1,
i =1
y por otra, que
q
ci = ∑ aij , i = 1, . . . , q.
j =1
En particular, esta última igualdad muestra que es innecesario dar de forma explícita los valores
de los coeficientes ci .
Todo método RK es un método unipaso:
y k +1 = y k + h Φ ( t k , y k ; h ), k = 0, 1, . . . , N − 1, (1.11)
28 Métodos unipaso para problemas de valor inicial
..
.
y = y + haq1 f (t(1) , y(1) ) + haq2 f (t(2) , y(2) ) + · · · + haq,q−1 f (t(q−1) , y(q−1) ).
(q)
9 Los huecos en el tablero indican que los elementos correspondientes son cero.
1.4 Métodos de Runge-Kutta 29
Es decir, en cada iteración, la solución del sistema (1.12) se escribe de manera explícita. Esto
hace que los métodos explícitos sean los más sencillos de implementar y los que requieren
un menor coste computacional.
Métodos semi-implícitos: corresponden a un tablero de la forma
c1 a11
c2 a21 a22
.. .. .. ..
. . . .
cq aq1 aq2 ··· aqq
b1 b2 ··· bq
En este caso, las ecuaciones del sistema (1.12) están desacopladas, es decir, en cada ecuación
aparece tan solo una incógnita. Concretamente, la primera ecuación es
corresponde al método RK4. Como es obvio, los cuatro métodos son explícitos.
Por otra parte, el método RK con tablero
1 1
1
y k +1 = y k + h f ( t k +1 , y k +1 ).
Una cuestión fundamental que surge al definir un método RK es si este está bien definido,
es decir, si el sistema (1.12) tiene solución única. Vamos a deducir a continuación una condición
suficiente para ello.
Consideremos la función G : Rq → Rq que a cada elemento Y = (y1 , . . . , yq )t ∈ Rq le asigna
el valor
g1 ( Y )
q
G (Y ) = ... , gi (Y ) = y + h ∑ aij f (t( j) , y j ), i = 1, . . . , q.
j =1
g q (Y )
El sistema (1.12) puede escribirse entonces de forma equivalente como una ecuación de punto
fijo:
G (Y ) = Y.
Para asegurar la existencia de un único punto fijo de la función G, vamos a imponer condiciones
para que G sea contractiva10 . Dados Y, Z ∈ Rq , para cada i = 1, . . . , q se tiene que
q q
( j) ( j)
| gi (Y ) − gi ( Z )| = h ∑ aij f (t , y j ) − f (t , z j ) 6 hL ∑ | aij ||y j − z j |.
j =1 j =1
En consecuencia,
| G (Y ) − G ( Z )| 6 hL| A||Y − Z |, (1.14)
donde hemos introducido la notación |Y | = (|y1 |, . . . , |yq |)t , | A| = (| aij |), y la desigualdad debe
k G (Y ) − G ( Z )k 6 ckY − Z k, ∀ Y, Z ∈ Rq ,
donde k · k es una norma en Rq . Al ser el espacio Rq completo, si G es contractiva entonces posee un único punto fijo
(esto es, Ŷ ∈ Rq tal que G (Ŷ ) = Ŷ), que se construye como límite de la sucesión
siendo Y0 ∈ Rq arbitrario. Tanto la definición como el resultado anteriores se generalizan de forma directa al caso de un
espacio métrico completo.
1.4 Métodos de Runge-Kutta 31
la aplicación G será contractiva y, en consecuencia, el sistema (1.12) tendrá solución única. Sin
embargo, esta condición es demasiado restrictiva, como se muestra en el ejemplo siguiente.
Ejemplo 1.8. Consideremos el método RK4, cuya matriz
0 0 0 0
1/2 0 0 0
A=
0 1/2 0 0
0 0 1 0
verifica k Ak∞ = 1. La condición hLk Ak∞ < 1 se satisface entonces si elegimos h < 1/L. Por
otra parte, al tratarse de un método explícito, está bien definido para cualquier h > 0. Por tanto,
en este caso la condición hLk Ak∞ < 1 introduce una restricción totalmente innecesaria sobre el
paso de malla. Lo mismo puede decirse si en lugar de RK4 consideramos cualquier otro método
explícito.
Notemos que para que la aplicación G posea un único punto fijo no es necesario que sea
contractiva, sino que basta con que lo sea alguna potencia suya12 . Sucesivas aplicaciones de la
desigualdad (1.14) implican que, para n ∈ N,
| G n (Y ) − G n ( Z )| 6 (hL| A|)n |Y − Z |,
de donde
kY k∞ = máx |yi |,
16 i 6 q
12 Supongamos que G m es contractiva para un cierto m ∈ N. Existe entonces un único Ŷ ∈ Rq tal que G m (Ŷ ) = Ŷ. Por
tanto, aplicando G, se tiene que G m+1 (Ŷ ) = G (Ŷ ) ⇒ G m ( G (Ŷ )) = G (Ŷ ). En consecuencia, G (Ŷ ) = Ŷ debido a la unicidad
del punto fijo de G m . Esto prueba que Ŷ es también el único punto fijo de G.
32 Métodos unipaso para problemas de valor inicial
donde ρ(| A|) representa el radio espectral13 de la matriz | A|. En tal caso, existe m ∈ N tal que14
hm Lm k| A|m k∞ < 1. Por consiguiente, la aplicación G m es contractiva, lo que implica que G tiene
un único punto fijo. En resumen, si el paso de malla h verifica la condición (1.15), el método RK
con paso h está bien definido.
Observación 1.13. La matriz A asociada a un método RK explícito verifica que ρ(| A|) = 0 (ya que
los autovalores de una matriz triangular son los elementos de su diagonal, que en este caso son
nulos), por lo que la condición (1.15) se verifica trivialmente para cualquier paso de malla.
A continuación vamos a analizar la estabilidad de un método RK; para ello nos basaremos en
el teorema 1.4. Tomemos pues y, z ∈ R, y sean y(i) , z(i) tales que
q q
y(i) = y + h ∑ aij f (t( j) , y( j) ), z(i) = z + h ∑ aij f (t( j) , z( j) ), i = 1, . . . , q.
j =1 j =1
Se tiene que
q
|y(i) − z(i) | 6 |y − z| + hL ∑ | aij ||y( j) − z( j) |,
j =1
de donde
|Y − Z | 6 |y − z|e + hL| A||Y − Z | 6 |y − z|e + h∗ L| A||Y − Z |,
siendo e = (1, 1, . . . , 1)t ∈ Rq . Por recurrencia, se tiene que
de donde
Supongamos ahora que h∗ Lρ(| A|) < 1. Entonces, tomando límite cuando n → ∞ en la desigual-
dad anterior, se deduce15 que
kY − Z k∞ 6 |y − z|k( I − h∗ L| A|)−1 k∞ .
13 El radio espectral de una matriz B ∈ Mq (R) se define como
ρ( B) = lı́m k Bn k1/n ,
n→∞
Por último,
q q
|Φ(t, y; h) − Φ(t, z; h)| 6 ∑ |bi || f (t(i) , y(i) ) − f (t(i) , z(i) )| 6 L ∑ |bi ||y(i) − z(i) |
i =1 i =1
q
6 Lk( I − h∗ L| A|)−1 k∞ ∑ | bi | | y − z | ≡ Λ | y − z |.
i =1
Esto prueba que Φ(t, y; h) es de Lipschitz en la variable y. Por el teorema 1.4, el método corres-
pondiente es estable.
Por otra parte, es inmediato comprobar que
q
Φ(t, y; 0) = ∑ i f (t, y).
b
i =1
Pero, según comentamos en la observación 1.12, dicha condición se verifica por construcción. Por
tanto, todo método RK es consistente. En consecuencia, aplicando el teorema 1.2, si suponemos
que h∗ Lρ(| A|) < 1 entonces el método RK asociado es convergente.
Un último apunte sobre la regularidad: si f ∈ C k ( I × R), por construcción se tiene de forma
inmediata que Φ ∈ C k ( I × R × [0, h∗ ]).
Podemos recopilar los hechos anteriores en el siguiente resultado.
Teorema 1.7. Sea A ∈ Mq (R) la matriz asociada a un método RK. Si se elige el paso de malla h ∈ [0, h∗ ]
tal que hLρ(| A|) < 1, entonces el método está bien definido. Si además h∗ Lρ(| A|) < 1, el método RK es
convergente. Por último, si f ∈ C k ( I × R) entonces Φ ∈ C k ( I × R × [0, h∗ ]).
Dedicaremos la última parte de esta sección a analizar el orden de un método RK. En primer
lugar, según la observación 1.10, la condición de orden uno equivale a la consistencia, que a su
vez es equivalente a (1.16); por simplicidad, escribiremos esta condición como
q q q
∂y(i) ∂f ∂y(i)
∂h
= ∑ aij f (t + c j h, y( j) ) + h ∑ aij ∂h
(t + c j h, y( j) ) ⇒
∂h h=0
= ∑ ij f (t, y),
a
j =1 j =1 j =1
de donde
q q
∂f ∂f
∂Φ
(t, y; 0) = ∑ bi ci (t, y) + f (t, y) (t, y) ∑ aij .
∂h i =1
∂t ∂y j =1
q
Como ci = ∑ j=1 aij (esto es, Ae = Ce; observación 1.12), deducimos que
q q q q
∂f ∂f
∂Φ
(t, y; 0) = ∑ bi ∑ aij (t, y) + f (t, y) (t, y) = ∑ bi ∑ aij f (1) (t, y).
∂h i =1 j =1
∂t ∂y i =1 j =1
Por tanto,
q q
∂Φ 1 1 1
(t, y; 0) = f (1) (t, y) ⇔ ∑ bi ∑ aij = ⇔ bt Ae = .
∂h 2 i =1 j =1
2 2
bt e = 1
(
(condiciones de orden dos). (1.17)
bt Ae = 21
De forma análoga, podemos determinar las condiciones para orden tres y cuatro:
Ejemplo 1.9. En el ejemplo 1.7 dimos los tableros de los métodos de Euler, Heun, punto medio
y RK4. Para el método de Euler se tiene que A = (0) y b = (1), por lo que la condición (1.16) se
verifica trivialmente; por tanto, el método tiene orden uno. Como además bt Ae = 0 6= 21 , el orden
es exactamente uno.
Por otra parte, para el método de Heun se verifican las condiciones (1.17):
1 1 1 1 1 0 0 1 1
bt e = , = 1, bt Ae = , = .
2 2 1 2 2 1 0 1 2
Como además
1 1 0 0 1 1 1
t 2
bC e= , = 6= ,
2 2 0 1 1 2 3
el orden es exactamente dos.
De forma similar se demuestra que el método del punto medio tiene orden dos, y RK4 orden
cuatro, ambos de forma exacta.
1.5 Adaptación de paso: métodos encajados 35
Sea {yk } la solución numérica generada por el método de orden p. Aplicando el teorema del
valor medio a la función y 7→ Φ∗ (tk , y; h) − Φ(tk , y; h) en los puntos y(tk ) e yk , se tiene que
Φ ∗ ( t k , y ( t k ); h ) − Φ ( t k , y ( t k ); h ) = Φ ∗ ( t k , y k ; h ) − Φ ( t k , y k ; h )
∂(Φ∗ − Φ)
+ (y(tk ) − yk ) ( t k , c k ; h ),
∂y
para un cierto ck entre yk e y(tk ). Por tanto,
∂(Φ∗ − Φ)
ε k − ε∗k = h(Φ∗ (tk , yk ; h) − Φ(tk , yk ; h)) + h(y(tk ) − yk ) ( t k , c k ; h ). (1.20)
∂y
36 Métodos unipaso para problemas de valor inicial
∂(Φ∗ − Φ) ∂(Φ∗ − Φ)
(tk , ck ; h) = ( t k , c k ; 0) + O ( h ).
∂y ∂y
La condición de consistencia (teorema 1.3) nos dice que Φ(t, y; 0) = Φ∗ (t, y; 0) = f (t, y), por lo
que
∂(Φ∗ − Φ) ∂f ∂f
( t k , c k ; 0) = ( t , c ) − ( t k , c k ) + O ( h ) = O ( h ).
∂y ∂y k k ∂y
Por otra parte, al ser de orden p el método que genera el valor yk , por el teorema 1.5 se tiene que
yk − y(tk ) = O(h p ). En consecuencia, de (1.20) se deduce la siguiente igualdad:
Teniendo en cuenta que ε k = O(h p+1 ) (de nuevo por la observación 1.11), deducimos que
def
ε k = h(Φ∗ (tk , yk ; h) − Φ(tk , yk ; h)) = O(h p+1 ),
e
|ε k | ≈ |eε k |.
Observación 1.14. Para un método de orden p, el teorema 1.5 nos proporcionaba una cota del
error de la forma
ek 6 Ch p .
Este es un ejemplo de estimación del error a priori, en el sentido de que dicha cota es conocida
antes de ejecutar el método. Aunque este tipo de acotaciones tiene interés teórico, en la práctica
su utilidad es limitada, ya que la constante C no es, en general, conocida. Además, aun en el
caso de que fuese conocida, la cota podría no ser significativa si la magnitud de C es grande en
relación a h p .
Por otra parte, la aproximación |ε k | ≈ |e
ε k | se denomina estimación del error a posteriori, ya
que se determina una vez ejecutado el método. En la práctica, este tipo de estimaciones resulta
especialmente interesante a la hora de controlar el paso de malla para que el error siempre se
mantenga por debajo de una tolerancia establecida.
Observación 1.15. En el desarrollo anterior hemos supuesto que la solución numérica yk se
determina mediante el método de menor orden. También sería posible calcular yk usando el
método de mayor orden. En este caso, si suponemos p > p∗ + 1 (normalmente, p = p∗ + 1),
razonando a partir de la igualdad (1.21) de forma similar a como hicimos anteriormente, se
llegaría a que
ε∗k = h(Φ(tk , yk ; h) − Φ∗ (tk , yk ; h)) + O(h p+1 ),
1.5 Adaptación de paso: métodos encajados 37
∗ +1
ε k = O(h p
de donde e ) y se tiene la aproximación
|ε∗k | ≈ |eε k |.
En resumen, el valor e
ε k constituye una aproximación del error de discretización local asociado
al método de menor orden.
La pregunta ahora es: ¿cómo elegir las funciones Φ y Φ∗ de manera adecuada? Más preci-
samente, queremos hacer un elección óptima, en el sentido de que el número de operaciones
necesarias para calcular e
ε k sea el menor posible. Un posibilidad es considerar métodos RK enca-
jados, denotados como RKp( p∗ ), que se representan mediante un tablero de la forma
o, brevemente,
c A
bt
(b∗ )t
donde se supone que el método determinado por A y b tiene orden p, y el definido por A y
b∗ tiene orden p∗ ; usualmente, p = p∗ − 1 o p = p∗ + 1. Notemos que la parte más costosa
computacionalmente de un método RK es la resolución del sistema determinado por A. Pero
al tener ambos métodos la misma matriz A, dicha resolución solo ha de realizarse una vez. En
resumen, el método RKp( p∗ ) asociado al tablero anterior se definiría como
q
(i ) (i )
y k + 1 = y k + h ∑ bi f ( t k , y k ) ,
i =1
(i ) (i )
donde tk = tk + ci h y los valores yk forman la solución del sistema
q
(i ) ( j) ( j)
yk = yk + h ∑ aij f (tk , yk ), i = 1, . . . , q
j =1
(esto es, la solución numérica se construye usando el método RK de orden p); la fórmula para la
estimación del error quedaría así:
q
(i ) (i )
ε k = h ∑ (bi∗ − bi ) f (tk , yk )
e
i =1
(un poco más adelante veremos cómo utilizar esta fórmula para elegir automáticamente el paso
(i ) (i )
de malla). Nótese que una vez resuelto el sistema, los valores f (tk , yk ) se calculan solo una vez
y se usan para determinar tanto yk+1 como e εk.
38 Métodos unipaso para problemas de valor inicial
Observación 1.16. En la notación RKp( p∗ ), p es el orden del método con el que se calcula la
solución numérica, mientras que p∗ es el orden del método que se utiliza para hacer la estimación
del error.
Ejemplo 1.11. Uno de los métodos RK encajados más utilizados es el método RK2(3), cuyo tablero
es
0 0 0 0
1/2 1/2 0 0
1 −1 2 0
orden 2 0 1 0
orden 3 1/6 2/3 1/6
0 0 0 0 0 0 0
1 1
4 4 0 0 0 0 0
3 3 9
8 32 32 0 0 0 0
12 1932
13 2197 − 7200
2197
7296
2197 0 0 0
439 3680 845
1 216 −8 513 − 4104 0 0
1 8
2 − 27 2 − 3544
2565
1859
4104
11
− 40 0
25 1408 2197
orden 4 216 0 2565 4104 − 51 0
16 6656 28561 9 2
orden 5 135 0 12825 56430 − 50 55
Como vimos en la observación 1.15, también es posible utilizar el método de mayor orden
para determinar la solución. Un ejemplo sería el método RK5(4) de Dormand-Prince, también
denotado como DOPRI5(4):
0 0 0 0 0 0 0 0
1 1
5 5 0 0 0 0 0 0
3 3 9
10 40 40 0 0 0 0 0
4 44 56 32
5 45 − 15 9 0 0 0 0
8 19372
9 6561 − 25360
2187
64448
6561
212
− 729 0 0 0
9017
1 3168 − 355
33
46732
5247
49
176
5103
− 18656 0 0
35 500 125
1 384 0 1113 192 − 2187
6784
11
84 0
35 500 125
orden 5 384 0 1113 192 − 2187
6784
11
84 0
5179 7571 393 92097 187 1
orden 4 57600 0 16695 640 − 339200 2100 40
1.5 Adaptación de paso: métodos encajados 39
En este caso, la solución se determina mediante el método de orden cinco, mientras que el de
orden cuatro se utiliza para la estimación del error. Otro método de Dormand-Prince de alto
orden muy utilizado en la práctica es DOPRI8(7).
Una vez elegido un método RK encajado, procedemos a realizar el equilibrado del paso de
malla, que consiste en utilizar la estimación del error e
ε k para elegir el paso de malla adecuado en
cada iteración. Concretamente, realizamos los siguientes pasos:
De esta forma, se obtiene un método con paso de malla variable: a partir de los valores iniciales
t0 , y0 , y de un paso de malla inicial dado h0 , la solución numérica se va construyendo variando
el paso de malla en cada iteración. Es importante notar que el número total de iteraciones es
desconocido de antemano. Por otra parte, la elección de la tolerancia ε dependerá del problema
que estemos resolviendo.
Observación 1.17. Expliquemos la definición de hk+1 . Por un lado, notemos que |e εk| > ε ⇒
hk+1 < hk , es decir, si la estimación del error es mayor que la tolerancia, el paso de malla dis-
minuye. Por el contrario, |e ε k | < ε ⇒ hk+1 < hk , lo que significa que el paso de malla aumenta
cuando el error se mantiene por debajo de la tolerencia; esto permite evitar la propagación de
pasos de malla pequeños.
p +1
Por otra parte, si p∗ > p + 1 sabemos que e
ε k = O( hk ), lo que en la práctica significa que
p +1
|eε k | ≈ Chk . Por tanto, teniendo en cuenta que p̂ = p,
p +1
p +1
ε
p +1 Chk
|eε k+1 | ≈ Chk+1 =C h =ε ≈ ε.
|eε k | k |eε k |
Esto es, el error de discretización local esperado en la siguiente iteración estará próximo a la
tolerancia. El caso p > p∗ + 1 se trataría de manera análoga, siendo ahora p̂ = p∗ .
Observación 1.18. En la práctica, para evitar problemas con los errores de redondeo, suele defi-
nirse
ε p̂+1 1
hk+1 = αhk , (1.23)
|eε k |
donde α ∈ (0, 1) es un factor de seguridad; habitualmente, se toma α igual a 0.8 o 0.9.
40 Métodos unipaso para problemas de valor inicial
Observación 1.19. Existen implementaciones más sofisticadas para un método RK encajado. Ha-
bitualmente, si el error calculado es mayor que la tolerancia, el paso de malla hk se rechaza, se
redefine usando la fórmula (1.23), y se vuelve a calcular e
ε k . Por otra parte, en determinadas apli-
caciones es conveniente establecer una tolerancia que sea proporcional al paso de malla, en cuyo
caso es necesario modificar la fórmula (1.23) adecuadamente.
Observación 1.20. A la hora de implementar un método con paso variable, es necesario modificar
el último paso de malla para que el último nodo coincida con t0 + T, el extremo superior del
intervalo. Concretamente, si en la iteración k-ésima se verifica que tk + hk > t0 + T, tomamos
hk = t0 + T − tk (o, equivalentemente, definimos tk+1 = t0 + T). Una forma automática para
conseguir esto es redefinir, en cada iteración hk como mı́n(hk , b − tk ).
y0 = −y + e−t cos(t),
(
t ∈ [0, 20],
y(0) = 0,
cuya solución exacta es y(t) = e−t sen(t). En la figura 1.9 se han representado las soluciones
obtenidas con los métodos de Heun y RK2(3). Para el método de Heun se ha tomado paso de
malla h = 0.2, con lo que el número de iteraciones es N = 100; el error global obtenido es
4.18 · 10−3 . Respecto al método RK2(3), se ha tomado paso de malla inicial h0 = 0.2 y como valor
de la tolerancia ε = 10−5 . Los valores máximo y mínimo del paso de malla en el proceso de
adaptación son, respectivamente, hmax ≈ 2.4 y hmin ≈ 0.028. El número de iteraciones realizadas
es 91, mientras que el error global es 2.59 · 10−3 .
Si ahora tomamos paso de malla h = 0.02, el método de Heun realiza N = 1000 iteraciones,
siendo el error global obtenido igual a 3.57 · 10−5 . Tomando paso de malla inicial h0 = 0.02 y
tolerancia ε = 10−6 , el método RK2(3) necesita 198 iteraciones, con un error global de 5.58 ·
10−5 .
Por último, en la figura 1.10 se han representado los resultados obtenidos con los métodos
RK4(5) y DOPRI5(4), tomando como paso de malla inicial h0 = 0.2 y tolerancia ε = 10−6 . Para
RK4(5) el número de iteraciones es 29 y el error es 3.41 · 10−6 , mientras que para DOPRI5(4) el
número de iteraciones es 28 y el error es 8.01 · 10−7 .
y(t) = y0 eλt ,
que verifica
lı́m y(t) = 0.
t→∞
De esta forma, la solución y(t) decrece hacia cero de forma exponencial, siendo la caída de la
función tanto más acusada conforme mayor sea |λ|. Si aplicamos ahora el método de Euler para
resolver el problema (1.24), las iteraciones adoptan la forma
de donde
yk = (1 + hλ)k y0 .
Por tanto, se tiene que
lı́m yk = 0 ⇔ |1 + hλ| < 1.
k→∞
De esta manera, la solución numérica se comportará como la solución exacta cuando t → ∞ si y
solo si elegimos el paso de malla de modo que |1 + hλ| < 1.
El razonamiento anterior es asimismo válido si tomamos λ ∈ C con Re(λ) < 0. Así pues,
para que la solución numérica de (1.24) converja hacia cero será necesario que |1 + ĥ| < 1, siendo
ĥ = λh. Basándonos en este hecho, definimos el polinomio de estabilidad absoluta para el método de
Euler como
R(ĥ) = 1 + ĥ, ĥ ∈ C,
y su dominio de estabilidad absoluta como
Geométricamente, D A representa el disco abierto de centro (−1, 0) y radio 1 (figura 1.11, izquier-
da).
Repitamos ahora el razonamiento anterior con el método de Euler implícito, definido en el
ejemplo 1.7. En tal caso, las iteraciones son de la forma
de donde
yk = (1 − hλ)−1 y0 .
Así pues,
lı́m yk = 0 ⇔ |1 − hλ| > 1.
k→∞
1.6 Estabilidad absoluta 43
Figura 1.11: Dominios de estabilidad absoluta para los métodos de Euler explícito (izquierda) e implícito
(derecha).
Por tanto, el dominio de estabilidad absoluta para el método de Euler implícito será
D A = {ĥ ∈ C : |1 − ĥ| > 1},
que corresponde al exterior del círculo cerrado de centro (1, 0) y radio 1 (figura 1.11, derecha).
Ejemplo 1.13. Consideremos el problema de Cauchy (1.24) con λ = −10. En tal caso, para que la
solución numérica obtenida mediante el método de Euler verifique lı́mk→∞ yk = 0, será necesario
que ĥ = λh ∈ D A . Al ser λ y h reales, basta entonces elegir h tal que −10h ∈ (−2, 0), es decir,
h ∈ 0, 51 . De modo análogo, si consideramos el método de Euler implícito, tendremos que elegir
h de modo que −10h ∈ (−∞, 0), es decir, basta con tomar h > 0.
Tomemos ahora λ = −2 + i ∈ C. En tal caso,
ĥ = (−2 + i )h ∈ D A ⇔ |1 + (−2 + i )h| < 1 ⇔ (1 − 2h)2 + h2 < 1 ⇔ h ∈ 0, 45 .
En consecuencia, para que la solución numérica del problema (1.24) dada por el método de Euler
tienda hacia cero, será necesario tomar h ∈ 0, 45 . Si empleamos el método de Euler implícito, se
deduce fácilmente que basta con tomar h > 0.
A continuación vamos a estudiar la estabilidad absoluta de un método RK arbitrario, definido
por (1.9)-(1.10). Consideremos pues el problema (1.24), donde f (t, y) = λy. En tal caso, el sistema
(1.10) se escribe como
q
(i ) ( j)
yk = yk + ĥ ∑ aij yk , i = 1, . . . , q,
j =1
siendo ĥ = λh. Si introducimos el vector unitario e = (1, . . . , 1)t ∈ Rq , se tiene entonces que
(1) (1) (1) (1)
yk yk yk y
. k
. = yk e + ĥA .. ⇒ ( I − ĥA) .. = yk e ⇒ .. = yk ( I − ĥA)−1 e.
. . . .
(1) (1) (1) (1)
yk yk yk yk
44 Métodos unipaso para problemas de valor inicial
Definición 1.5. Dado un método RK definido por la matriz A y el vector b, definimos su función de
estabilidad absoluta como
R(ĥ) = 1 + ĥbt ( I − ĥA)−1 e, ĥ ∈ C, (1.25)
y su dominio de estabilidad absoluta como el conjunto
Las expresiones (1.9)-(1.10) que definen a un método RK, pueden interpretarse como un sis-
(1) (q)
tema lineal de q + 1 ecuaciones cuyas incógnitas son {yk , . . . , yk , yk+1 }; podemos entonces
determinar yk+1 mediante la regla de Cramer. Tras unos sencillos cálculos, se obtiene que
Por tanto, podemos escribir la función de estabilidad absoluta de forma alternativa como
Figura 1.12: Dominios de estabilidad absoluta para los métodos RK explícitos de orden p y q = p etapas,
para p = 1, 2, 3, 4.
La ventaja de la expresión (1.26) frente a (1.25) es que no es necesario calcular ninguna matriz
inversa. Notemos que tanto el numerador como el denominador en (1.26) son polinomios en la
variable ĥ. Esto muestra que la función de estabilidad absoluta es una función racional. En el
caso particular en que el método RK sea explícito, se tiene que det( I − ĥA) = 1, por lo que R(ĥ)
es, de hecho, un polinomio (véase el ejemplo 1.14). Por tanto, la región D A de un método RK
explícito ha de estar necesariamente acotada, ya que la condición | R(ĥ)| < 1 no puede darse si
|ĥ| → ∞ (véase la figura 1.12). Por el contrario, un método RK no explícito puede tener región de
estabilidad absoluta no acotada.
h
y k +1 = y k + f ( t k , y k ) + f ( t k +1 , y k +1 ) ,
2
es un método RK semi-implícito cuyo tablero es
0 0 0
1 1/2 1/2
1/2 1/2
2 + ĥ
R(ĥ) = .
2 − ĥ
46 Métodos unipaso para problemas de valor inicial
Entonces,
| R(ĥ)| < 1 ⇔ |2 + ĥ| < |2 − ĥ| ⇔ Re(ĥ) < 0,
y por tanto,
D A = { ĥ ∈ C : Re(ĥ) < 0}.
Desde el punto de vista numérico, es conveniente que el dominio de estabilidad absoluta sea
lo más grande posible, para permitir una mayor amplitud en la elección del paso de malla. A este
respecto, haremos la siguiente definición.
Ejemplo 1.16. Ningún método RK explícito puede ser A-estable, por ser su función de estabili-
dad absoluta un polinomio. En cambio, tanto el método de Euler implícito como el método del
trapecio son A-estables.
Para hacer la definición de estabilidad absoluta considerábamos una ecuación lineal escalar
de la forma y0 = λy. ¿Qué sentido tiene entonces dicha definición cuando tratamos con una
ecuación no lineal o con un sistema de ecuaciones? Consideremos en primer lugar el caso de un
sistema lineal
Y 0 = MY,
1.7 Problemas rígidos (stiff ) 47
donde Z = P−1 Y. Al ser la matriz D diagonal, el sistema se reduce de esta forma a un conjunto
de n ecuaciones desacopladas:
0
z1 = λ1 z1 ,
..
.
0
zn = λn zn ,
Y 0 = F (t, Y ), F : I × Rn → Rn ,
donde 0
y1 y1 −100 1
Y= , Y0 = , M= 1 .
y2 y02 0 − 10
16 Una matriz M ∈ Mn (R) es diagonalizable si puede descomponerse en la forma
M = PDP−1 ,
donde D es una matriz diagonal y P es una matriz invertible. Los elementos de la diagonal principal de D son los
autovalores λi de M, y cada columna de P es un autovector por la derecha asociado a λi , siguiendo el orden establecido
en D.
48 Métodos unipaso para problemas de valor inicial
En particular,
lı́m Yk = 0,
k→∞
en cualquier norma que elijamos. Es interesante notar que la solución (1.28) tiene dos componen-
tes cuya escala de variación es muy diferente: la exponencial e−100t decae a cero de forma mucho
más rápida que e−t/10 .
Apliquemos ahora el método de Euler para resolver el problema (1.27). Se tiene entonces
Yk+1 = Yk + hMYk ⇒ Yk = ( I + hM)k Y0 = P( I + hD )k P−1 Y0 .
Operando, se obtiene la siguiente expresión para la solución numérica:
k
h
ke
Yk = (1 − 100h) Y1 + 1 − e2 ,
Y (1.29)
10
donde Y e1 e Y
e2 son los vectores definidos anteriormente. Para que la sucesión {Yk } converja a
1 h
cero, es necesario elegir h > 0 de forma que |1 − 100h| < 1 (esto es, h < 50 ) y |1 − 10 | < 1 (es
1
decir, h < 20); por tanto, en general será necesario que h ∈ (0, 50 ). Una forma más sencilla (pero
h
menos ilustrativa) de obtener este resultado consiste en imponer que −100h ∈ D A y − 10 ∈ DA ,
1
teniendo en cuenta que los autovalores de la matriz M son −100 y − 10 .
Tomemos ahora como valor inicial
1
Y0 = 999 ,
10
1
que es un autovector asociado al autovalor − 10 . En tal caso, Y
e1 = 0 e Y
e2 = Y0 ; por tanto, la
solución numérica queda de la forma
k
h
Yk = 1 − Y0 .
10
1.7 Problemas rígidos (stiff ) 49
1
En consecuencia, lı́mk→∞ Yk = 0 si y solo si h ∈ (0, 20). Elijamos como paso de malla h = 10 ;
entonces,
k
99 99
Yk = Y0 ⇒ ln kYk k∞ = k ln + ln kY0 k∞ .
100 100
Así pues, si representamos en una gráfica los puntos de la forma (k, ln kYk k∞ ), estos deberán
99
estar situados sobre una recta de pendiente ln( 100 ) ≈ −0.01005. Sin embargo, al resolver el
1
problema de Cauchy (1.27) mediante el método de Euler con paso de malla h = 10 , se obtienen
los resultados de la figura 1.14.
1
Figura 1.14: Resolución del problema (1.27) mediante el método de Euler con paso de malla h = 10 .
¿Qué está sucediendo? Cuando se implementa el método de Euler, la solución numérica adop-
ta la forma (1.29). Aunque se haya elegido la condición inicial de modo que el vector Y e1 se anula,
en la práctica esto no sucede debido a la influencia de los errores de redondeo. Estos añaden una
contribución no nula al término (1 − 100h)k Y e1 (que teóricamente debería ser cero), que crece de
1
forma geométrica a menos que h < 50 . De esta forma, al cabo de un número suficiente de itera-
e1 termina imponiéndose a (1 − h )k Y
ciones, el término (1 − 100h)k Y e2 , que crece de forma mucho
10
más lenta. De hecho, en la figura 1.14 puede observarse que los puntos a la derecha de la gráfica
están situados, aproximadamente, sobre una recta de pendiente ln(9) ≈ 2.1972, que corresponde
al término (1 − 100h)k Y0 = 9k Y0 .
En conclusión, aunque teóricamente bastaría con elegir h ∈ (0, 20), en la práctica es necesario
1
tomar un paso de malla mucho más restrictivo, h ∈ (0, 50 ), para evitar la acumulación inde-
seada de errores de redondeo. El motivo de fondo es la existencia en la solución numérica de
componentes que varían a escalas muy diferentes. Cuando esto sucede, diremos que el problema
considerado es rígido.
De forma algo más concreta, consideremos un sistema de ecuaciones de la forma Y 0 = MY,
donde se supone que la matriz M ∈ Mn (R) posee autovalores λi con Re(λi ) < 0, i = 1, . . . , n.
50 Métodos unipaso para problemas de valor inicial
máxi | Re(λi )|
rs = 1.
mı́ni | Re(λi )|
Cuando rs 1, la solución del sistema posee componentes que tienden a cero a velocidades muy
dispares. La componente que tiende a cero de forma más rápida (que está asociada al autovalor
de módulo máximo) es la que establece la restricción más estricta sobre el paso de malla. En el
1
ejemplo anterior los autovalores eran λ1 = −100 y λ2 = − 10 , por lo que rs = 1000. El autovalor
de módulo mayor, λ1 , era el que imponía una restricción más fuerte sobre el paso de malla al
utilizar el método de Euler.
Apliquemos ahora el método del trapecio (véase el ejemplo 1.15) al problema (1.27). En tal
caso, se obtiene:
−1
h h h
Yk+1 = Yk + ( MYk + MYk+1 ) ⇒ Yk+1 = I− M I+ M Yk ,
2 2 2
de donde
−1 k
h h
Yk = I− M I+ M Y0 .
2 2
Usando la diagonalización de la matriz M, la expresión anterior puede escribirse de la siguiente
forma:
k k
1 − 50h e 20 − h e
Yk = Y1 + Y2 ,
1 + 50h 20 + h
que se verifican para cualquier h > 0. De esta forma, al aplicar el método del trapecio al problema
considerado, no es necesario imponer ninguna restricción sobre el paso de malla para que el
comportamiento asintótico de la solución numérica sea el adecuado.
Otra forma de llegar a la conclusión anterior consiste en considerar el dominio de estabilidad
absoluta del método del trapecio. Según vimos en el ejemplo 1.15, D A = C− ; por tanto, al ser
1 h
−100 y − 10 los autovalores de M, debemos imponer las condiciones −100h ∈ C− y − 10 ∈ C− ,
que se verifican para cualquier h > 0.
17 El símbolo se lee «mucho mayor que».
1.7 Problemas rígidos (stiff ) 51
También es posible encontrar problemas rígido con una única ecuación. Un ejemplo es el
siguiente problema de Cauchy:
Figura 1.15: Izquierda: Solución exacta del problema (1.30). Derecha: Ampliación de la zona de transición
rápida.
En este caso, la función f (t, y) = −1200y + 2000 − 1500e−t que define la ecuación diferencial es
no lineal. Por tanto, para hacer el análisis de la estabilidad absoluta, consideramos la ecuación
linealizada
∂f
y0 = (t , y ) ⇒ y0 = −1200y.
∂y k k
Para aplicar el método de Euler, hemos de elegir el paso h de modo que
1
−1200h ∈ D A = {ĥ ∈ C : |1 + ĥ| < 1} ⇒ −1200h ∈ (−2, 0) ⇒ h ∈ 0, 600 .
1
Además, habrá que tomar el paso de malla aún más pequeño que 600 ≈ 0.00167, ya que dicho
valor se ha obtenido a partir de una linealización de la ecuación original. A este respecto, en la
figura 1.16 se muestran los resultados obtenidos con el método de Euler en la zona de transición
1 1
rápida, para los pasos h = 700 y h = 1400 .
52 Métodos unipaso para problemas de valor inicial
1 1
Figura 1.16: Resolución del problema (1.30) mediante el método de Euler con pasos h = 700 yh= 1400 .
1
Observación 1.22. El paso de malla h < 600 , necesario para evitar inestabilidades en la zona de
transición rápida, es el que va a usarse para calcular el resto de la solución numérica. Esto implica
que en la zona de transición lenta se va a usar un paso de malla mucho más pequeño de lo que
sería necesario para que el error de discretización se mantuviese por debajo de una tolerancia
establecida. En general, este problema no puede arreglarse usando una técnica de adaptación de
paso, ya que la propiedad de estabilidad absoluta va a exigir tomar pasos de malla muy pequeños
en todo el dominio de cálculo.
Observación 1.23. Podríamos plantearnos la posibilidad de utilizar, en lugar del método de Eu-
ler, un método de mayor orden (por ejemplo, RK4) para capturar adecuadamente la zona de
transición rápida usando un paso de malla menos restrictivo. Sin embargo, como puede obser-
varse en la figura 1.12, el dominio de estabilidad absoluta del método RK4 no es mucho mayor
que el del método de Euler. Por tanto, el rango de elección del paso de malla que permite el
método RK4 no es significativamente mayor que el correspondiente al método de Euler.
Consideremos ahora el método de Euler implícito
y k +1 = y k + h f ( t k +1 , y k +1 ),
que corresponde a un método RK de una etapa con A = (1) y b = (1). A partir de (1.26), se
deduce fácilmente que el dominio de estabilidad absoluta es
como ya dedujimos también en la sección 1.6. En particular, el método es A-estable, por lo que no
es necesaria ninguna restricción sobre el paso de malla para garantizar la estabilidad absoluta. En
1
la figura 1.17 se comparan los resultados obtenidos mediante el método de Euler con h = 700 y
1
el método de Euler implícito con h = 200 . Como puede observarse, el método de Euler implícito
1.8 Métodos implícitos: implementación 53
captura la solución de manera mucho más precisa que el método de Euler, usando un paso de
malla mucho más pequeño.
1 1
Figura 1.17: Comparación entre los métodos de Euler con h = 700 y Euler implícito con h = 200 en la zona
de transición rápida.
Observación 1.24. En general, la implementación del método de Euler implícito requiere la re-
solución de una ecuación no lineal en cada iteración para determinar yk+1 (sección 1.8). Sin
embargo, en el caso considerado es posible despejar yk+1 de forma directa:
yk + 2000h − 1500he−tk+1
yk+1 = yk + h(−1200yk+1 + 2000 − 1500e−tk+1 ) ⇒ yk+1 = .
1 + 1200h
Esto permite escribir el método de Euler implícito de forma explícita para este problema.
numérica yk+1 , partiendo del dato inicial y0 . Sin embargo, en el caso de un método implícito, en
la definición de yk+1 aparece el propio valor a determinar, yk+1 ; por tanto, no podemos definir
este de forma directa, sino que su valor se obtendrá como solución de una ecuación.
Para fijar ideas, consideremos un método unipaso implícito de la forma
y k +1 = Ξ ( y k +1 , y k , t k ; h ), k = 0, 1, . . . , N − 1, (1.31)
con las notaciones habituales. En cada iteración, el valor de yk+1 será solución de la ecuación de
punto fijo
def
w = Ξ(w, yk , tk ; h) = Ξk (w), (1.32)
donde yk , tk y h están fjados (es importante notar que a lo largo de los siguientes razonamientos
el valor de k está fijo). En primer lugar, la ecuación (1.32) debe poseer solución única, lo que
equivale a decir que el método (1.31) está bien definido. En algunos casos particulares (por ejem-
plo, cuando Ξk es lineal) es posible despejar analíticamente el valor de w, con lo cual podemos
resolver directamente la ecuación (1.32); dicho de otra forma, podemos escribir el método implí-
cito (1.31) de forma explícita (véase la observación 1.24). Sin embargo, en general no podremos
resolver (1.32) de forma directa.
Una forma de estudiar la unicidad de solución consiste en analizar la contractividad de la
aplicación Ξk : R → R, lo cual nos proporciona además un método práctico para aproximar
dicha solución mediante iteraciones de punto fijo. De forma más concreta, supongamos que la
aplicación Ξk es contractiva, esto es, existe una constante c ∈ (0, 1) de modo que
|Ξk (w1 ) − Ξk (w2 )| 6 c|w1 − w2 |, ∀ w1 , w2 ∈ R.
En tal caso, la ecuación (1.32) posee una única solución w,
e y la sucesión definida recursivamente
por
wn+1 = Ξk (wn ), n = 0, 1, . . .
e para cualquier valor inicial w0 ∈ R. En particular, w
converge a w, e verifica
e = Ξk (w
w e ),
por lo que we es el valor buscado de yk+1 .
Normalmente, la expresión de la constante de contractividad c involucrará al paso de malla h,
por lo que la condición c < 1 impondrá una restricción sobre la elección del paso de malla. Para
determinar w e ≡ yk+1 podemos hacer iteraciones de punto fijo hasta que el error (absoluto, rela-
tivo o una combinación de ambos) entre dos iteraciones sucesivas sea menor que una tolerancia
establecida; también es conveniente, para evitar bucles demasiado largos o infinitos, establecer
un número máximo de iteraciones de punto fijo, de modo que si dicho número se alcanza, el
código dé un aviso (warning).
Un inconveniente de las iteraciones de punto fijo es su lentitud, ya que su velocidad de con-
vergencia es lineal. Por ello, una vez que sabemos que la ecuación (1.32) tiene solución única,
podríamos plantearnos el utilizar otro tipo de método numérico para aproximarla. En la práctica
es habitual considerar el método de Newton que, bajo condiciones adecuadas, converge cuadrática-
mente. Para ello es necesario reescribir primero la ecuación de punto fijo (1.32) como un problema
de búsqueda de raíces:
Ξ
e k (w) = 0,
1.8 Métodos implícitos: implementación 55
Ξ
e k ( wn )
w n +1 = w n − , n = 0, 1, . . .
Ξ
e 0 ( wn )
k
donde la semilla w0 suele tomarse como la solución yk obtenida en la iteración anterior del
método numérico.
En resumen, a la hora de programar un método implícito de la forma (1.31), tenemos que
hacer un bucle sobre el índice k (que varía entre 0 y N − 1) y, para cada valor de k, resolver el
problema de punto fijo (1.32) para definir el valor de yk+1 .
Para ilustrar los razonamientos anteriores, consideremos el método de Euler implícito:
y k +1 = y k + h f ( t k +1 , y k +1 ), k = 0, 1, . . . , N − 1.
Ξ
e k (w) = w − Ξk (w) = w − yk − h f (tk + h, w),
∂f
e 0 (w) = 1 − h (tk + h, w). De esta forma, el método de Newton se escribiría como
en cuyo caso Ξ k ∂y
wn − yk − h f (tk + h, wn )
w n +1 = w n − ∂f
, n = 0, 1, . . . (1.34)
1 − h ∂y (tk + h, wn )
56 Métodos unipaso para problemas de valor inicial
cuya solución exacta es la función y(t) = e−t sen(t). Para ello iremos paso a paso, explicando las
diferentes partes del código. Este comienza así:
# Resolucion del problema de valor inicial
# y '= f (t , y ) , y ( t0 )= y0 ,
# mediante el metodo de Euler .
import pylab
o bien usando un alias (esto es conveniente cuando el nombre de la biblioteca es largo o engorro-
so):
import pylab as P
(aquí P es un nombre arbitrario). Cuando queramos usar una función, por ejemplo el seno, ten-
dremos que escribir antes el nombre de la biblioteca: pylab.sin, o bien su alias: P.sin. Aunque
esta sea una buena práctica de programación, en los sucesivo nos decantaremos por la opción
from pylab import * para que los códigos sean más legibles.
Observación 1.26. En los comentarios hemos evitado el uso de tildes para evitar problemas de
codificación; es posible, sin embargo, usar tildes y otros caracteres propios de nuestro idioma
incluyendo las siguientes líneas al principio del código:
# !/ usr / bin / env python
# -* - coding : utf -8 -* -
Por simplicidad, en los sucesivos códigos prescindiremos de las tildes.
58 Métodos unipaso para problemas de valor inicial
Para definir la función f (t, y) que determina la ecuación diferencial y0 = f (t, y), hacemos lo
siguiente:
def fun (t , y ):
""" Funcion f (t , y ) que define la ecuacion diferencial """
return -y + exp ( - t )* cos ( t )
También definimos la solución exacta (en general, esta no será conocida):
def exacta ( t ):
""" Solucion exacta del problema de valor inicial """
return exp ( - t )* sin ( t )
La parte principal del programa es la implementación del método de Euler; para ello, vamos a
definir una función cuyos parámetros de entrada sean el intervalo [ a, b], el número de particiones
N y la condición inicial y0 , y cuya salida sean los nodos {tk }kN=0 y la solución numérica {yk }kN=0 .
De esta forma, la definición del método será general y no dependerá del problema particular que
estamos resolviendo. La función es la siguiente:
def euler (a , b , N , y0 ):
""" Metodo de Euler en el intervalo [a , b ]
usando N particiones y condicion inicial y0 """
h = (b - a )/ N # paso de malla
t = zeros ( N +1) # inicializacion del vector de nodos
y = zeros ( N +1) # inicializacion del vector de resultados
t [0] = a # nodo inicial
y [0] = y0 # valor inicial
# Metodo de Euler
for k in range ( N ):
y [ k +1] = y [ k ]+ h * fun ( t [ k ] , y [ k ])
t [ k +1] = t [ k ]+ h
return (t , y )
Detengámonos un poco en esta función. Lo primero que hay que tener en cuenta es que los
bloques de código se indican en Python mediante la tabulación, por lo que esta es de capital
importancia. La función empieza definiendo el paso de malla h e inicializando a cero los vectores
t e y, de longitud N+1, donde almacenaremos los nodos y la solución numérica; después se
definen los valores iniciales a e y0. El método de Euler propiamente dicho se implementa a
continuación, mediante un bucle for que va recorriendo los índices desde 0 hasta N-1 (en efecto,
range(N) representa el conjunto de valores {0, 1,...,N-1}); de esta forma, el último valor que
se define sería y[k+1]=y[N], cuando k=N-1. Como puede observarse, las definiciones de y[k+1]
y t[k+1] son una transcripción casi literal de la expresión (1.2) para el método de Euler. Por
último, la función devuelve los vectores (elementos de la clase array) t e y.
El programa principal comienza definiendo los datos del problema a resolver; en este caso,
a = 0, b = 5, N = 100 e y0 = 0. A continuación se llama a la función euler con los datos
anteriormente definidos.
1.9 Implementación en Python de métodos unipaso 59
Observación 1.27. Es importante hacer hincapié en que las operaciones sobre arrays están vec-
torizadas, esto es, cuando se aplica una función a un array el resultado es un nuevo array cuyas
componentes se forman evaluando la función en cada una de las componentes del array. Así pues,
en nuestro ejemplo t es un array que contiene los valores (t[0],...,t[N]); al aplicar la función
exacta a t, se crea un nuevo array (que almacenamos en la variable ye) cuyas componentes son
(exacta(t[0]),...,exacta(t[N])). En el cálculo de error, estamos haciendo la diferencia de
dos arrays con la misma longitud (que vuelve a ser un array), y luego calculamos el máximo de
sus componentes.
Presentamos en pantalla los resultados obtenidos:
# Resultados
print ( ' ----- ')
print ( ' Error : ' + str ( error ))
print ( ' Paso de malla : ' + str (( b - a )/ N ))
print ( ' ----- ')
El símbolo + sirve para concatenar (unir) dos cadenas de caracteres; como los valores tfin-tini,
error, etc. son números en coma flotante, antes debemos transformarlos en cadenas de caracteres
con el comando str.
Por último, dibujamos en una misma gráfica las soluciones exacta y aproximada:
# Dibujamos las soluciones
plot (t , y , 'r . - ') # dibuja la solucion aproximada
plot (t , ye , 'b ') # dibuja la solucion exacta
xlabel ( 't ') # etiqueta eje de abscisas
ylabel ( 'y ') # etiqueta eje de ordenadas
legend ([ ' Euler ' , ' exacta ' ]) # leyendas
grid ( True ) # malla de fondo
show () # muestra la grafica
18 EnPython 2.7 y anteriores, el resultado de hacer 1/2 es 0, ya que la división se hace en el conjunto de los números
enteros. Para obtener el valor 0.5 basta con que alguno de los factores esté escrito en coma flotante: 1./2, 1/2. o 1./2.
dan como resultado 0.5. Sin embargo, en Python 3.0 y siguientes el resultado de hacer 1/2 es 0.5.
60 Métodos unipaso para problemas de valor inicial
def fun (t , y ):
""" Funcion f (t , y ) que define la ecuacion diferencial """
return -y + exp ( - t )* cos ( t )
def exacta ( t ):
""" Solucion exacta del problema de valor inicial """
return exp ( - t )* sin ( t )
def euler (a , b , N , y0 ):
""" Metodo de Euler en el intervalo [a , b ]
usando N particiones y condicion inicial y0 """
h = (b - a )/ N # paso de malla
t = zeros ( N +1) # inicializacion del vector de nodos
y = zeros ( N +1) # inicializacion del vector de resultados
t [0] = a # nodo inicial
y [0] = y0 # valor inicial
# Metodo de Euler
for k in range ( N ):
y [ k +1] = y [ k ]+ h * fun ( t [ k ] , y [ k ])
t [ k +1] = t [ k ]+ h
return (t , y )
# Resultados
print ( ' ----- ')
print ( ' Tiempo CPU : ' + str ( tfin - tini ))
print ( ' Error : ' + str ( error ))
print ( ' Paso de malla : ' + str (( b - a )/ N ))
print ( ' ----- ')
h = (b - a )/ N # paso de malla
t = zeros ( N +1) # inicializacion del vector de nodos
y = zeros ( N +1) # inicializacion del vector de resultados
t [0] = a # nodo inicial
y [0] = y0 # valor inicial
# Metodo de Heun
for k in range ( N ):
faux = fun ( t [ k ] , y [ k ])
y [ k +1] = y [ k ]+0.5* h *( faux + fun ( t [ k ]+ h , y [ k ]+ h * faux ))
t [ k +1] = t [ k ]+ h
return (t , y )
Como el valor f (tk , yk ) aparece dos veces en el método de Heun, lo almacenamos en la variable
faux para poder reutilizarlo y, de esta forma, calcularlo solo una vez. En el programa princi-
62 Métodos unipaso para problemas de valor inicial
pal simplemente habría que cambiar la llamada (t, y) = euler(a, b, N, y0) por (t, y) =
heun(a, b, N, y0), y cambiar ‘Euler’por ‘Heun‘ donde proceda.
Para el método del punto medio la función correspondiente sería:
def punto_medio (a , b , N , y0 ):
""" Metodo del punto medio en el intervalo [a , b ]
usando N particiones y condicion inicial y0 """
h = (b - a )/ N # paso de malla
t = zeros ( N +1) # inicializacion del vector de nodos
y = zeros ( N +1) # inicializacion del vector de resultados
t [0] = a # nodo inicial
y [0] = y0 # valor inicial
return (t , y )
mientras que para RK4 tendríamos:
def rk4 (a , b , N , y0 ):
""" Metodo RK4 en el intervalo [a , b ]
usando N particiones y condicion inicial y0 """
h = (b - a )/ N # paso de malla
t = zeros ( N +1) # inicializacion del vector de nodos
y = zeros ( N +1) # inicializacion del vector de resultados
t [0] = a # nodo inicial
y [0] = y0 # valor inicial
# Metodo RK4
for k in range ( N ):
K1 = fun ( t [ k ] , y [ k ])
K2 = fun ( t [ k ]+0.5* h , y [ k ]+0.5* h * K1 )
K3 = fun ( t [ k ]+0.5* h , y [ k ]+0.5* h * K2 )
K4 = fun ( t [ k ]+ h , y [ k ]+ h * K3 )
y [ k +1] = y [ k ]+ h /6.*( K1 +2.* K2 +2.* K3 + K4 )
return (t , y )
Por último, vamos a dar un código donde se resuelve, mediante los cuatro métodos conside-
rados, el problema de valor inicial
cuya solución exacta es y(t) = (π + 1)e−t + sen(t) − cos(t). El resultado puede verse en la figura
1.18.
1.9 Implementación en Python de métodos unipaso 63
Figura 1.18: Izquierda: Comparación entre los métodos de Euler, Heun, punto medio y RK4, usando N = 25
particiones. Derecha: ampliación de la solución para x ∈ [0, 2]. Puede observarse que los métodos de Heun
y del punto medio son indistinguibles, mientras que la solución más precisa la proporciona RK4.
def fun (t , y ):
""" Funcion que define la ecuacion diferencial """
return -y +2.* sin ( t )
def exacta ( t ):
""" Solucion exacta del problema de valor inicial """
return ( pi +1)* exp ( - t )+ sin ( t ) - cos ( t )
h = (b - a )/ N # paso de malla
y [ k +1] = y [ k ]+ h * fun ( t [ k ] , y [ k ])
elif ( metodo == ' Heun ' ):
# metodo de Heun
for k in range ( N ):
faux = fun ( t [ k ] , y [ k ])
y [ k +1] = y [ k ]+0.5* h *( faux + fun ( t [ k ]+ h , y [ k ]+ h * faux ))
elif ( metodo == ' Punto medio ' ):
# metodo del punto medio
for k in range ( N ):
y [ k +1] = y [ k ]+ h * fun ( t [ k ]+0.5* h ,
y [ k ]+0.5* h * fun ( t [ k ] , y [ k ]))
elif ( metodo == ' RK4 ' ):
# metodo RK4
for k in range ( N ):
K1 = fun ( t [ k ] , y [ k ])
K2 = fun ( t [ k ]+0.5* h , y [ k ]+0.5* h * K1 )
K3 = fun ( t [ k ]+0.5* h , y [ k ]+0.5* h * K2 )
K4 = fun ( t [ k ]+ h , y [ k ]+ h * K3 )
y [ k +1] = y [ k ]+ h /6.*( K1 +2.* K2 +2.* K3 + K4 )
return y
h = (b - a )/ N # paso de malla
t = arange (a , b +h , h ) # nodos
ye = exacta ( t ) # solucion exacta
metodos = [ ' Euler ' , ' Heun ' , ' Punto medio ' , ' RK4 '] # metodos
estilos = [ '^ - ' , 's - ' , 'o - ' , 'D - '] # estilos de linea
legend ([ ' Euler ' , ' Heun ' , ' Punto medio ' , ' RK4 ' , ' Exacta ' ])
grid ( True )
1.9 Implementación en Python de métodos unipaso 65
def fun (t , y ):
""" Funcion que define la ecuacion diferencial """
return -y +2.* sin ( t )
def exacta ( t ):
""" Solucion exacta del problema de valor inicial """
return ( pi +1)* exp ( - t )+ sin ( t ) - cos ( t )
met = 0
tini = clock ()
if ( met == 1):
(t , y ) = mu . euler (a , b , N , fun , y0 )
name = ' Euler '
elif ( met == 2):
(t , y ) = mu . heun (a , b , N , fun , y0 )
name = ' Heun '
elif ( met == 3):
(t , y ) = mu . punto_medio (a , b , N , fun , y0 )
name = ' Pto . medio '
else :
(t , y ) = mu . rk4 (a , b , N , fun , y0 )
name = ' RK4 '
tfin = clock ()
dx = 0.01
te = arange (a , b + dx , dx ) # nodos
ye = exacta ( te ) # solucion exacta
dx
= 0.4x − 0.01xy,
dt
(1.37)
dy = −0.01y2 + 0.01xy,
dt
que modeliza la interacción entre dos especies, una de las cuales se alimenta de la otra. La
variable x ≡ x (t) representa la población de presas y la variable y ≡ y(t) la de depredadores,
en el instante de tiempo t. Se supone que inicialmente hay la misma cantidad de especímenes de
ambas especies: x (0) = y(0) = 10 (en las unidades adecuadas), y queremos estudiar la evolución
del sistema hasta el instante de tiempo t = 30.
Comencemos con el método de Euler; como ya vimos al final de la sección 1.2, basta con
escribir dicho método en forma vectorial. En este caso, cada elemento Yk de la solución aproxi-
x
mada tendrá dos componentes: Yk = ykk . Guardaremos estos vectores en un array (matriz) z
de tamaño 2 × ( N + 1) (en general, de tamaño n × ( N + 1) para un sistema de n ecuaciones). La
columna k-ésima de z contendrá la solución numérica en el nodo k-ésimo, Yk . La estructura de z
será pues de la forma
x0 x1 x2 · · · x N
y0 y1 y2 · · · y N
h = (b - a )/ N # paso de malla
t = zeros ( N +1) # vector de nodos
z = zeros ([ len ( z0 ) , N +1]) # matriz de resultados
t [0] = a # nodo inicial
z [: , 0] = z0 # valor inicial
# Metodo RK4
for k in range ( N ):
K1 = fun ( t [ k ] , z [: , k ])
K2 = fun ( t [ k ]+0.5* h , z [: , k ]+0.5* h * K1 )
K3 = fun ( t [ k ]+0.5* h , z [: , k ]+0.5* h * K2 )
K4 = fun ( t [ k ]+ h , z [: , k ]+ h * K3 )
z [: , k +1] = z [: , k ]+ h /6.*( K1 +2.* K2 +2.* K3 + K4 )
t [ k +1] = t [ k ]+ h
return (t , z )
11
12 def exacta ( t ):
13 """ Solucion exacta del problema de valor inicial """
14 return exp ( - t )* sin ( t )
15
16 def rk23 (a , b , y0 ):
17 """ Implementacion del metodo encajado RK2 (3)
18 en el intervalo [a , b ] con condicion inicial y0 """
19
20 h0 = (b - a )/100. # paso de malla inicial
21 hmin = (b - a )*1. e -5 # paso de malla minimo
22 hmax = (b - a )/10. # paso de malla maximo
23 tol = 1. e -6 # tolerancia
24
25 # coeficientes RK
26 q = 3 # numero de etapas del metodo
27 A = zeros ([ q , q ]) # matriz del sistema
28 A [1 , 0] = 1./2.
29 A [2 , 0] = -1.
30 A [2 , 1] = 2.
31
32 B = zeros ( q ) # vector b del metodo de orden 2
33 B [1] = 1.
34
35 BB = zeros ( q ) # vector b ^* del metodo de orden 3
36 BB [0] = 1./6.
37 BB [1] = 2./3.
38 BB [2] = 1./6.
39
40 C = zeros ( q ) # matriz diagonal con los valores c_ij
41 for i in range ( q ):
42 C [ i ] = sum ( A [i ,:])
43
44 # inicializacion de variables
45 t = array ([ a ]) # nodos
46 y = array ([ y0 ]) # soluciones
47 h = array ([ h0 ]) # pasos de malla
48 K = zeros ( q )
49 k = 0 # contador de iteraciones
50
51 while ( t [ k ] < b ):
52 h [ k ] = min ( h [ k ] , b - t [ k ]) # ajuste del ultimo paso de malla
53
54 for i in range ( q ):
55 K [ i ] = fun ( t [ k ]+ C [ i ]* h [ k ] , y [ k ]+ h [ k ]* sum ( A [i ,:]* K ))
56
57 incrlow = sum ( B * K ) # metodo de orden 2
58 incrhigh = sum ( BB * K ) # metodo de orden 3
59 error = h [ k ]*( incrhigh - incrlow ) # estimacion del error
60 y = append (y , y [ k ]+ h [ k ]* incrlow ) # y_ ( k +1)
61 t = append (t , t [ k ]+ h [ k ]); # t_ ( k +1)
62 hnew = 0.9* h [ k ]* abs ( tol / error )**(1./3.) # h_ ( k +1)
63 hnew = min ( max ( hnew , hmin ) , hmax ) # hmin <= h_ ( k +1) <= hmax
72 Métodos unipaso para problemas de valor inicial
64 h = append (h , hnew )
65 k += 1
66
67 h = h [: -1] # elimina el ultimo h calculado
68 hn = min ( h ) # minimo valor utilizado del paso de malla
69 hm = max ( h ) # maximo valor utilizado del paso de malla
70
71 return (t , y , hn , hm )
72
73 # Datos del problema
74 a = 0. # extremo inferior del intervalo
75 b = 30. # extremo superior del intervalo
76 y0 = 0. # condicion inicial
77
78 tini = clock ()
79
80 (t , y , hn , hm ) = rk23 (a , b , y0 ) # llamada al metodo RK2 (3)
81
82 tfin = clock ()
83
84 # calculo de la solucion exacta
85 dx = 0.01
86 te = arange (a , b + dx , dx )
87 ye = exacta ( te )
88
89 # Dibujamos las soluciones
90 plot (t , y , 'bo - ') # dibuja la solucion aproximada
91 plot ( te , ye , 'r ') # dibuja la solucion exacta
92 xlabel ( 't ')
93 ylabel ( 'y ')
94 legend ([ ' RK2 (3) ' , ' exacta ' ])
95 grid ( True )
96
97 # Calculo del error cometido
98 error = max ( abs (y - exacta ( t )))
99
100 # Resultados
101 print ( ' ----- ')
102 print ( ' Tiempo CPU : ' + str ( tfin - tini ))
103 print ( ' Error : ' + str ( error ))
104 print ( ' No . de iteraciones : ' + str ( len ( y )))
105 print ( ' Paso de malla minimo : ' + str ( hn ))
106 print ( ' Paso de malla maximo : ' + str ( hm ))
107 print ( ' ----- ')
108
109 show () # muestra la grafica
Veamos a continuación las particularidades del código anterior. En la función rk23 se definen
el paso de malla inicial h0, los pasos mínimo y máximos hmin y hmax, y la tolerancia del error tol.
Estos valores se establecen dependiendo del problema estudiado, por lo que para un problema
más complejo quizás sería necesario cambiarlos si los resultados no son aceptables. A continua-
ción se definen las matrices A, b, b∗ y C que definen el método RK2(3), y se inicializan los arrays
1.9 Implementación en Python de métodos unipaso 73
donde se almacenarán los nodos, la solución numérica y los pasos de malla utilizados (esto últi-
mo no es necesario); como puede observarse, los arrays t, y y h contienen inicialmente un único
valor: el instante inicial a, el valor inicial y0 y el paso de malla inicial h0 , respectivamente.
La parte más compleja es el bucle while, por lo que vamos a comentarlo en detalle. Definimos
inicialmente K como un array nulo con q = 3 componentes, y ponemos a cero el contador de
iteraciones: k=0. El bucle while se ejecuta mientras que el valor t[k]<b (inicialmente, t[0]=a<b,
por lo que entramos en el bucle al menos una vez). La línea 52 sirve para ajustar el último
paso de malla, de modo que el último nodo sea precisamente el extremo final del intervalo
(i ) (i )
(observación 1.20). A continuación se guarda en la componente i del array K el valor f (tk , yk );
nótese que la expresión de la línea 54 es válida porque inicialmente K es el array nulo. Las
variables incrlow e incrhigh contienen, respectivamente, los valores Φ(tk , yk ; hk ) y Φ∗ (tk , yk ; hk ),
mientras que la estimación del error eε k se define en la variable error. En la línea 58 se actualiza
el array y, añadiendole el valor yk + hk Φ(tk , yk ; hk ); esto se realiza mediante el comando append,
que añade un elemento al final de una lista o un array. También se añade al array t el nuevo
nodo tk+1 = tk + hk . El equilibrado del paso de malla definido en (1.22) se realiza en la línea 60,
donde se ha añadido un factor de seguridad de 0.9. En la línea 61 simplemente se comprueba
que el nuevo paso de malla esté comprendido entre los valores mínimo y máximo establecidos
inicialmente. Por último, se añade al array h el último paso de malla calculado y se aumenta en
uno el contador de iteraciones.
Una vez fuera del bucle while, se elimina de h el último paso de malla calculado (que no se
usa en el cálculo) y determinamos los pasos de malla mínimo y máximo que se han utilizado. Por
último, la función rk23 devuelve el array de nodos t, la solución numérica y, el paso de malla
mínimo hn y el paso de malla máximo hm.
El programa principal es análogo a los ya utilizados para los métodos no encajados, salvo
en un par de detalles. Al llamar a la función rk23 en la línea 78, hay que añadir hn y hm a los
datos de salida; estos datos se imprimen luego en las líneas 103-104. En las líneas 83-85 se calcula
la solución exacta usando una partición uniforme con paso 0.01; esto se hace para lograr una
representación gráfica suave de la solución exacta, que no dependa de los nodos determinados
por el método encajado.
Figura 1.19: Comparación entre los métodos encajados RK2(3), RK4(5) y RK5(4).
74 Métodos unipaso para problemas de valor inicial
En el siguiente código se implementan los métodos RK2(3), RK4(5) y RK(5)4, y se aplican para
resolver el problema (1.36):
# Resolucion del problema de valor inicial
# y '= f (t , y ) , y ( t0 )= y0 ,
# mediante los metodos RK2 (3) , RK4 (5) y RK5 (4).
def fun (t , y ):
""" Funcion que define la ecuacion diferencial """
return -y +2.* sin ( t )
def exacta ( t ):
""" Solucion exacta del problema de valor inicial """
return ( pi +1)* exp ( - t )+ sin ( t ) - cos ( t )
A [5 , 1] = 2.
A [5 , 2] = -3544./2565.
A [5 , 3] = 1859./4104.
A [5 , 4] = -11./40.
# B : metodo de orden 4
B [0] = 25./216.
B [2] = 1408./2565.
B [3] = 2197./4104.
B [4] = -1./5.
# BB : metodo de orden 5
BB [0] = 16./135.
BB [2] = 6656./12825.
BB [3] = 28561./56430.
BB [4] = -9./50.
BB [5] = 2./55.
elif ( metodo == ' RK5 (4) ' ): # RK5 (4) Dormand - Prince
p = 4 # menor de los ordenes
q = 7 # etapas del metodo
A = zeros ([ q , q ])
B = zeros ( q )
BB = zeros ( q )
A [1 , 0] = 1./5.
A [2 , 0] = 3./40.
A [2 , 1] = 9./40.
A [3 , 0] = 44./45.
A [3 , 1] = -56./15.
A [3 , 2] = 32./9.
A [4 , 0] = 19372./6561.
A [4 , 1] = -25360./2187.
A [4 , 2] = 64448./6561.
A [4 , 3] = -212./729.
A [5 , 0] = 9017./3168.
A [5 , 1] = -355./33.
A [5 , 2] = 46732./5247.
A [5 , 3] = 49./176.
A [5 , 4] = -5103./18656.
A [6 , 0] = 35./384.
A [6 , 2] = 500./1113.
A [6 , 3] = 125./192.
A [6 , 4] = -2187./6784.
A [6 , 5] = 11./84.
# B : metodo de orden 5
B [0] = 35./384.
B [2] = 500./1113.
B [3] = 125./192.
B [4] = -2187./6784.
B [5] = 11./84.
# BB : metodo de orden 4
BB [0] = 5179./57600.
BB [2] = 7571./16695.
BB [3] = 393./640.
BB [4] = -92097./339200.
BB [5] = 187./2100.
76 Métodos unipaso para problemas de valor inicial
BB [6] = 1./40.
C = zeros ( q )
for i in range ( q ):
C [ i ] = sum ( A [i ,:])
return (A , B , BB , C , p )
# coeficientes RK
(A , B , BB , C , p ) = matrices ( metodo )
q = size ( B ) # numero de etapas
# inicializacion de variables
t = array ([ a ]) # nodos
y = array ([ y0 ]) # soluciones
K = zeros ( q )
while ( t [ k ] < b ):
h = min (h , b - t [ k ]) # ajuste del ultimo paso de malla
for i in range ( q ):
K [ i ] = fun ( t [ k ]+ C [ i ]* h , y [ k ]+ h * sum ( A [i ,:]* K ))
incrlow = sum ( B * K )
incrhigh = sum ( BB * K )
error = h *( incrhigh - incrlow ) # estimacion del error
y = append (y , y [ k ]+ h * incrlow ) # y_ ( k +1)
t = append (t , t [ k ]+ h ); # t_ ( k +1)
h = 0.9* h * abs ( tol / error )**(1./( p +1)) # h_ ( k +1)
h = min ( max (h , hmin ) , hmax ) # hmin <= h_ ( k +1) <= hmax
hn = min ( hn , h )
hm = max ( hm , h )
k += 1
return (t , y , k , hn , hm )
metodos = [ ' RK2 (3) ' , ' RK4 (5) ' , ' RK5 (4) '] # metodos
63 h = append (h , hnew )
64 k += 1
65
66 h = h [: -1] # elimina el ultimo h calculado
67 hn = min ( h ) # minimo valor del paso de malla
68 hm = max ( h ) # maximo valor del paso de malla
69
70 return (t , z , hn , hm )
71
72 # Datos del problema
73 a = 0. # extremo inferior del intervalo
74 b = 30. # extremo superior del intervalo
75 z0 = array ([10. ,10.]) # condicion inicial
76 z0 = z0 . reshape (2 ,1) # vector columna
77
78 tini = clock ()
79
80 (t , z , hn , hm ) = rk23 (a , b , z0 )
81
82 tfin = clock ()
83
84 # Graficas
85 subplot (2 ,2 ,1)
86 plot (t , z [0] , 'b ')
87 grid ( True )
88 title ( ' Evolucion de las presas ')
89 subplot (2 ,2 ,2)
90 plot (t , z [1] , 'r ')
91 grid ( True )
92 title ( ' Evolucion de los depredadores ')
93 subplot (2 ,2 ,3)
94 plot (t , z [0] , 'b ' , t , z [1] , 'r ')
95 grid ( True )
96 title ( ' Comparacion ')
97 subplot (2 ,2 ,4)
98 plot ( z [0] , z [1])
99 grid ( True )
100 title ( ' Orbita ')
101
102 print ( ' Tiempo CPU : ' + str ( tfin - tini ))
103 print ( ' Paso de malla minimo : ' + str ( hn ))
104 print ( ' Paso de malla maximo : ' + str ( hm ))
105
106 show ()
En este ejemplo hay una serie de detalles técnicos que es conveniente comentar. En primer
lugar, K es una matriz de tamaño n × q (línea 47), donde n = 2 es el número de ecuaciones del
sistema y q = 3 el número de etapas del método. Al definir la columna i-ésima de K en la línea
53, hay que tener en cuenta los tamaños de los arrays para que el producto de A por K se pueda
realizar. Por este motivo se toma la transpuesta de K (por cierto, para multiplicar dos matrices se
usa el comando dot, no *); lo mismo sucede al definir incrlow e incrhigh. El comando append
no funciona para matrices, por lo que para añadir una nueva columna al array z que contiene la
80 Métodos unipaso para problemas de valor inicial
def fun (t , z ):
""" Funcion que define el sistema """
# z = (x , y ); x = z [0] , y = z [1]
return array ([0.4* z [0] -0.01* z [0]* z [1] ,
-0.01* z [1]**2+0.01* z [0]* z [1]])
A [4 , 0] = 439./216.
A [4 , 1] = -8.
A [4 , 2] = 3680./513.
A [4 , 3] = -845./4104.
A [5 , 0] = -8./27.
A [5 , 1] = 2.
A [5 , 2] = -3544./2565.
A [5 , 3] = 1859./4104.
A [5 , 4] = -11./40.
# B : metodo de orden 4
B [0] = 25./216.
B [2] = 1408./2565.
B [3] = 2197./4104.
B [4] = -1./5.
# BB : metodo de orden 5
BB [0] = 16./135.
BB [2] = 6656./12825.
BB [3] = 28561./56430.
BB [4] = -9./50.
BB [5] = 2./55.
elif ( metodo == ' RK5 (4) ' ): # RK5 (4) Dormand - Prince
p = 4 # menor de los ordenes
q = 7 # etapas del metodo
A = zeros ([ q , q ])
B = zeros ( q )
BB = zeros ( q )
A [1 , 0] = 1./5.
A [2 , 0] = 3./40.
A [2 , 1] = 9./40.
A [3 , 0] = 44./45.
A [3 , 1] = -56./15.
A [3 , 2] = 32./9.
A [4 , 0] = 19372./6561.
A [4 , 1] = -25360./2187.
A [4 , 2] = 64448./6561.
A [4 , 3] = -212./729.
A [5 , 0] = 9017./3168.
A [5 , 1] = -355./33.
A [5 , 2] = 46732./5247.
A [5 , 3] = 49./176.
A [5 , 4] = -5103./18656.
A [6 , 0] = 35./384.
A [6 , 2] = 500./1113.
A [6 , 3] = 125./192.
A [6 , 4] = -2187./6784.
A [6 , 5] = 11./84.
# B : metodo de orden 5
B [0] = 35./384.
B [2] = 500./1113.
B [3] = 125./192.
B [4] = -2187./6784.
B [5] = 11./84.
# BB : metodo de orden 4
82 Métodos unipaso para problemas de valor inicial
BB [0] = 5179./57600.
BB [2] = 7571./16695.
BB [3] = 393./640.
BB [4] = -92097./339200.
BB [5] = 187./2100.
BB [6] = 1./40.
C = zeros ( q )
for i in range ( q ):
C [ i ] = sum ( A [i ,:])
return (A , B , BB , C , p )
# coeficientes RK
(A , B , BB , C , p ) = matrices ( metodo )
q = size ( B ) # numero de etapas
# inicializacion de variables
t = array ([ a ]) # nodos
z = z0 . reshape ( len ( z0 ) , 1) # condicion inicial ( vector columna )
K = zeros ([ len ( z0 ) , q ])
while ( t [ k ] < b ):
h = min (h , b - t [ k ]) # ajuste del ultimo paso de malla
for i in range ( q ):
K [: , i ] = fun ( t [ k ]+ C [ i ]* h ,
z [: , k ]+ h * dot ( A [i ,:] , transpose ( K )))
return (t , z , k , hn , hm )
1.9 Implementación en Python de métodos unipaso 83
metodos = [ ' RK2 (3) ' , ' RK4 (5) ' , ' RK5 (4) ']
n =1
y k +1 = y k + h f ( t k +1 , y k +1 ), k = 0, 1, . . . , N − 1,
h = (b - a )/ N
t = zeros ( N +1) # inicializacion del vector de nodos
y = zeros ( N +1) # inicializacion del vector de resultados
for k in range ( N ):
( t [ k +1] , y [ k +1]) = pto_fijo ( t [ k ] , y [ k ] , h )
return (t , y )
Como vemos, la estructura de la función es análoga a la de la función euler, excepto en el bucle
for. En este caso, se hace una llamada a la función punto_fijo con los datos de entrada tk , yk y
h, que nos devuelve los valores de tk+1 e yk+1 . La función punto_fijo resuelve el problema (1.32)
mediante iteraciones de punto fijo:
1 def pto_fijo (t , y , h ):
2 """ " Resolucion del problema de punto fijo """
3
4 tol = 1. e -6 # tolerancia
5 M = 1000 # numero maximo de iteraciones
6 n = 0 # contador de iteraciones
7 error = tol +1.
8 t1 = t + h # define t_ ( k +1)
9 w0 = y # semilla
10
11 while ( error > tol ) and ( n <= M ):
12 w1 = y + h * fun ( t1 , w0 )
13 error = abs ( w1 - w0 )
14 n += 1
15 w0 = w1
16
17 if ( n > M ):
18 print ( ' Aviso : numero maximo de iteraciones en t = ' + str ( t1 ))
19
20 return ( t1 , w1 )
Los datos de entrada de la función son: el tiempo t≡ tk , el valor y≡ yk y el paso de malla h≡ h.
Primero se establecen la tolerancia tol y el número máximo de iteraciones M (estos valores pueden
depender del problema a resolver). Después se define el valor de la variable error de modo que
entremos al menos una vez en el bucle while de la línea 11 (basta con asignarle cualquier valor
mayor que tol). A continuación definimos t1≡ tk+1 y asignamos a la semilla w0 el valor yk .
El núcleo de la función está en el bucle while, que se ejecuta mientras la variable error sea
mayor que la tolerancia tol establecida, y el número de iteraciones n sea menor que el número
1.9 Implementación en Python de métodos unipaso 85
def fun (t , y ):
""" Funcion que define la ecuacion diferencial """
return -1200.* y +2000. -1500.* exp ( - t )
def exacta ( t ):
""" Solucion exacta del problema de valor inicial """
return 5./3. -1495./3597.* exp ( -1200.* t ) -1500./1199.* exp ( - t )
def pto_fijo (t , y , h ):
""" " Resolucion del problema de punto fijo """
tol = 1. e -6 # tolerancia
M = 1000 # numero maximo de iteraciones
n = 0 # contador de iteraciones
error = tol +1.
t1 = t + h # t_ ( k +1)
w0 = y # semilla
if ( n > M ):
print ( ' Aviso : numero maximo de iteraciones en t = ' + str ( t1 ))
return ( t1 , w1 )
86 Métodos unipaso para problemas de valor inicial
def euler_implicito (a , b , N , y0 ):
""" Implementacion del metodo de Euler implicito """
h = (b - a )/ N
t = zeros ( N +1) # inicializacion del vector de nodos
y = zeros ( N +1) # inicializacion del vector de resultados
for k in range ( N ):
( t [ k +1] , y [ k +1]) = pto_fijo ( t [ k ] , y [ k ] , h )
return (t , y )
tini = clock ()
tfin = clock ()
time = tfin - tini
# Resultados
print ( ' ----- ')
print ( ' Tiempo CPU : ' + str ( time ))
print ( ' Error : ' + str ( error ))
print ( ' Paso de malla : ' + str (( b - a )/ N ))
print ( ' ----- ')
Hemos introducido la variable zoom que, cuando toma el valor True, nos permite dibujar una
ampliación de la zona de transición rápida entre t = 0 y t = 0.03.
Por último, también podemos considerar el método de Newton para determinar el valor
de yk+1 . Concretamente, la expresión (1.34) para el problema considerado se implementa en la
siguiente función, cuya estructura es análoga a la de la función pto_fijo:
def newton (t , y , h ):
""" Resolucion de la ecuacion mediante el metodo de Newton """
tol = 1. e -6
M = 1000
n = 0
error = tol +1.
t1 = t + h
w0 = y
if ( n > M ):
print ' Aviso : numero maximo de iteraciones en t = ' + str ( t1 )
return ( t1 , w1 )
En la función euler_implicito, bastaría con cambiar la línea
( t [ k +1] , y [ k +1]) = pto_fijo ( t [ k ] , y [ k ] , h )
por
( t [ k +1] , y [ k +1]) = newton ( t [ k ] , y [ k ] , h )
88 Métodos unipaso para problemas de valor inicial
CAPÍTULO
y sea {tk }kN=0 una partición uniforme del intervalo I, con paso de malla h > 0. Veamos en primer
lugar cómo construir distintos métodos multipaso a partir de un proceso de integración numérica.
90 Métodos multipaso para problemas de valor inicial
Otra posibilidad consiste en añadir el valor f k+1 al proceso de interpolación, esto es, utilizar
el polinomio de interpolación Qq (t) asociado a los puntos
{(tk+1 , f k+1 ), (tk , f k ), . . . , (tk−q+1 , f k−q+1 )}.
De esta forma, se obtiene el método de Adams-Moulton (AM) de profundidad q:
Z t
k +1
y k +1 = y k + Qq (t) dt, k = q − 1, . . . , N − 1. (2.5)
tk
h
y k +1 = y k + ( f k +1 + f k ), k = 0, . . . , N − 1.
2
En este caso no se obtiene un verdadero método multipaso, sino un método unipaso implícito: el
método del trapecio o de Crank-Nicolson.
Consideremos ahora los puntos de interpolación {(tk−1 , f k−1 ), (tk , f k ), (tk+1 , f k+1 )} (esto es,
q = 2). El método AM correspondiente, que denotaremos como AM3 (donde 3 es el orden del
método, como ya veremos), tiene la siguiente expresión:
h
y k +1 = y k + (5 f k +1 + 8 f k − f k −1 ), k = 1, . . . , N − 1. (2.6)
12
Como puede observarse, se trata de un método implícito.
Observación 2.1. Otra clase de métodos multipaso basados en integración numérica son los
métodos de Milne-Simpson, que se construyen a partir de la formulación integral de la ecuación
diferencial en el intervalo [tk−1 , tk+1 ], usando el polinomio Qq (t):
Z t
k +1
y k +1 = y k −1 + Qq (t) dt.
t k −1
En general, esta formulación da lugar a un método implícito. Por ejemplo, para Q1 (t) se obtiene
el método multipaso explícito
yk+1 = yk−1 + 2h f k ,
mientras que el método implícito
h
y k +1 = y k −1 + ( f k +1 + 4 f k + f k −1 )
3
se obtiene al considerar el polinomio Q2 (t).
92 Métodos multipaso para problemas de valor inicial
y k +1 = y k + h ( β 0 f k +1 + β 1 f k + · · · + β q f k − q +1 ),
para unos ciertos coeficientes β j ∈ R. Esta forma resulta tras desarrollar las expresiones (2.3) y
(2.5) y unificar la notación.
Otra clase de métodos multipaso puede construirse mediante una técnica de derivación nu-
mérica. La idea consiste en evaluar la ecuación diferencial en el punto tk+1 ,
se define el método de diferenciación regresiva (en inglés, backward differentiation formula, o BDF) de
profundidad q como
R0q (tk+1 ) = f k+1 .
Se trata siempre de un método implícito. Desarrollando la igualdad anterior, un método de
diferenciación regresiva puede escribirse de forma general como
Las expresiones anteriores pueden combinarse para construir un método del tipo
α 0 y k +1 + α 1 y k + · · · + α p y k − p +1 = h ( β 0 f k +1 + β 1 f k + · · · + β p f k − p +1 ),
α r y k +r + α r −1 y k +r −1 + · · · + α 0 y k = h ( β r f k +r + β r −1 f k +r −1 + · · · + β 0 f k ), (2.7)
se denomina método multipaso lineal de profundidad r, donde se supone que se verifican las siguien-
tes condiciones:
αr 6= 0, |α0 | + | β 0 | > 0.
Observación 2.2. La condición αr 6= 0 hace que el término yk+r aparezca de forma explícita
en la expresión del método; normalizando, es usual tomar αr = 1. Por otra parte, la condición
|α0 | + | β 0 | > 0 hace que aparezca el término yk explícita (si α0 6= 0) o implícitamente (si β 0 6= 0).
De esta forma, la profundidad del método (el número de puntos necesarios para determinar yk+r )
es realmente r. Por último, cuando β r = 0 el método es explícito, siendo implícito si β r 6= 0.
lı́m ys = y0 ,
h →0
para cada s = 1, . . . , r − 1.
Observación 2.4. Por construcción, si consideramos un paso de malla variable hk los coeficientes
del método (2.7) también dependerán del índice k, variando en cada iteración. Esto complica
el análisis del método correspondiente, por lo que no consideraremos aquí el caso de métodos
multipaso con paso de malla variable.
¿Cuándo está bien definido un método multipaso lineal? Si β r = 0 el método (2.7) es explícito,
por lo que está bien definido para cualquier h > 0. Supongamos entonces que β r 6= 0 y escribamos
(2.7) como una ecuación de punto fijo cuya incógita es yk+r , esto es,
y k +r = Ξ k ( y k +r ),
donde
βr
Ξk (w) = h f (tk+r , w) + Ck ,
αr
94 Métodos multipaso para problemas de valor inicial
siendo
h 1
Ck = (β f + · · · + β 0 f k ) − ( α r −1 y k +r −1 + · · · + α 0 y k ),
α r r −1 k +r −1 αk
que permanece constante. Dados w1 , w2 ∈ R, se tiene que
r r
β β
|Ξk (w1 ) − Ξk (w2 )| = h | f (tk+r , w1 ) − f (tk+r , w2 )| 6 hL |w1 − w2 |.
αr αr
Por tanto, la aplicación Ξk : R → R será contractiva si tomamos h > 0 de modo que
r
β
hL < 1 (2.8)
αr
(nótese que la condición obtenida es independiente del índice k). En conclusión, bajo la condición
(2.8) existe un único punto fijo yk+r de Ξk (w), lo cual a su vez implica que el método (2.7) está
bien definido.
Observación 2.5. Para un método explícito se tiene que β r = 0, por lo que la condición (2.8)
no impone ninguna restricción sobre el paso de malla. Podemos por tanto suponer que (2.8) se
verifica, independientemente de si el método multipaso es explícito o implícito.
Las consideraciones anteriores son clave a la hora de realizar la implementación práctica de
un método multipaso implícito. De hecho, se siguen las mismas directrices que ya vimos para
métodos unipaso implícitos en la sección 1.8.
En concreto, fijado el índice k, tenemos que determinar yk+r como solución de la ecuación de
punto fijo
Ξk (w) = w.
Para ello, consideramos la semilla w0 = yk+r−1 (el último valor calculado de la solución numérica)
y hacemos iteraciones de punto fijo:
w n +1 = Ξ k ( w n ), n = 0, 1, . . .
Si elegimos el paso de malla h de modo que se verifique la condición (2.8), la aplicación Ξk es
contractiva y, por tanto, la sucesión {wn } convergerá hacia el único punto fijo yk+r .
Otra posibilidad consiste en utilizar el método de Newton para calcular yk+r , con la idea de
acelerar la convergencia. En ese caso, podríamos considerar la función Ξe (w) = w − Ξk (w) (u otra
elección adecuada), de modo que se tiene la equivalencia
Ξk (w) = w ⇔ Ξ
e k (w) = 0.
Ξ
e k ( wn )
w n +1 = w n − ,
Ξ
e 0 ( wn )
k
donde
e k (w) = w − h β r f (tk+r , w) − Ck ,
Ξ e 0 ( w ) = 1 − h β r ∂ f ( t k +r , w ).
Ξ k
αr αr ∂y
En el apartado 2.6.2 pueden encontrarse implementaciones en Python de diversos métodos
multipaso lineales.
2.3 Convergencia, estabilidad y consistencia 95
h
y k +2 = y k +1 + (5 f k +2 + 8 f k +1 − f k ), k = 0, . . . , N − 2.
12
5
Se trata de un método implícito de profundidad r = 2, con α2 = 1 y β 2 = 12 ; por tanto, la
12
condición (2.8) se reduce a h < 5L .
Las iteraciones de punto fijo quedan de la forma
5h
w n +1 = f (t + 2h, wn ) + Ck , n = 0, 1, . . .
12 k
h
con Ck = yk+1 + 12 (8 f k+1 − f k ) y semilla w0 = yk+1 . Si usamos el método de Newton, las itera-
ciones quedarían de la siguiente forma:
5h
wn − 12 f ( tk + 2h, wn ) − Ck
w n +1 = w n − ∂f
, n = 0, 1, . . .
1 − 5h
12 ∂y ( t k + 2h, w n )
Definición 2.1 (Convergencia). Consideremos el problema de Cauchy (2.1) y sea y(t) su solución exacta.
Diremos que un método multipaso lineal es convergente si
para cada sucesión {yk } obtenida mediante el método con valores iniciales {y0 , y1 , . . . , yr−1 } que verifican
lı́m ys = y0 , (2.10)
h →0
para cada s = 1, . . . , r − 1.
Observación 2.6. Vamos a explicar con algo más de detalle la definición de convergencia. Con-
sideremos una sucesión de particiones del intervalo I obtenidas dividiendo sucesivamente por
dos el paso de malla: h0 , h1 , h2 , . . . , donde hn = h0 /2n . Denotemos por {yk [hn ]}k a la solución
numérica dada por el método con paso hn . Fijado un valor t ∈ I (por ejemplo, t = t0 + 5h0 ),
podemos expresarlo en las distintas particiones como
Como la sucesión de pasos de malla hn converge hacia cero, la idea intuitiva de convergencia se
expresaría diciendo que la sucesión anterior converge a y(t). En general, este comportamiento
debería verificarse cuando h → 0 de cualquier forma. Así pues, debemos considerar un límite
doble en el que h → 0 y k → ∞ simultáneamente, pero con la restricción de que t = t0 + kh
permanezca fijo. Además esta propiedad ha de verificarse para cada punto t ∈ I.
Por otra parte, la condición (2.10) es fundamental, ya que asegura que los valores iniciales
están cerca del único valor conocido de la solución, que es el proporcionado por la condición
inicial y(t0 ) = y0 .
Observación 2.7. Podríamos haber considerado una definición de convergencia similar a la que
vimos para métodos unipaso (definición 1.3). Sin embargo, aunque dicha definición es mucho
más simple, no es lo suficientemente ilustrativa de las dificultades de notación que aparecen al
usar la definición en determinados desarrollos teóricos.
Es habitual definir un método multipaso lineal a partir de sus polinomios característicos. Dado
un método multipaso lineal de la forma (2.7), definimos su primer polinomio característico como
ρ ( z ) = α r z r + α r −1 z r −1 + · · · + α 1 z + α 0 ,
σ ( z ) = β r z r + β r −1 z r −1 + · · · + β 1 z + β 0 .
Nótese que ρ(z) siempre tiene grado r, mientras que σ (z) tiene grado r sólo cuando el método
es implícito (β r 6= 0). Las tablas 2.1, 2.2 y 2.3 muestran, respectivamente, los polinomios carac-
terísticos de los métodos de Adams-Bashforth, Adams-Moulton y de diferenciación regresiva de
profundidad r = 1, 2, 3, 4.
r ρ(z) σ(z)
1 z−1 1
3 1
2 z2 − z 2z − 2
23 2 16 5
3 z3 − z2 12 z − 12 z + 12
55 3 59 2 37 9
4 z4 − z3 24 z − 24 z + 24 z − 24
Definición 2.2 (Estabilidad). Un método multipaso lineal es estable si todas las raíces de su primer
polinomio característico ρ(z) tienen módulo menor o igual que uno, y las de módulo uno, si las hubiere,
son simples.
2.3 Convergencia, estabilidad y consistencia 97
r ρ(z) σ(z)
1 1
1 z−1 2z + 2
5 2 8 1
2 z2 − z 12 z + 12 z − 12
9 3 19 2 5 1
3 z3 − z2 24 z + 24 z − 24 z + 24
251 4 646 3 264 2 106 19
4 z4 − z3 720 z + 720 z − 720 z + 720 z − 720
r ρ(z) σ(z)
1 z−1 z
2 z2 − 43 z + 1
3
2 2
3z
18 2 9 2 6 3
3 z3 − 11 z + 11 z − 11 11 z
48 3 36 2 16 3 12 4
4 z4 − 25 z + 25 z − 25 z + 25 25 z
Observación 2.8. La definición de estabilidad es, a priori, un tanto oscura, por lo que vamos
a tratar de aclararla un poco. Para ello tendremos que echar mano de algunas herramientas
algebraicas.
Supongamos que los coeficientes β i son todos nulos, por lo que el método adopta la siguiente
forma:
αr yk+r + αr−1 yk+r−1 + · · · + α0 yk = 0. (2.11)
La igualdad anterior puede interpretarse como una ecuación en recurrencia lineal homogénea. En este
sentido, una solución de (2.11) es una sucesión {yn }n∈N de valores que verifican la ecuación para
cada índice. Es claro que una solución de (2.11) queda determinada de manera única dando los
primeros r valores {y0 , y1 , . . . , yr−1 }, por lo que el espacio vectorial de soluciones tiene dimensión
r. En este contexto, el polinomio ρ(z) se denomina polinomio característico de la ecuación (2.11).
Dado η 6= 0, es inmediato comprobar que la sucesión dada por yn = η n es solución de (2.11)
si y sólo si ρ(η ) = 0. La condición η 6= 0 se verifica si suponemos α0 6= 0. Podemos entonces
determinar la solución general de (2.11) teniendo en cuenta los siguientes hechos:
Si todas las raíces ηi de ρ(z) son reales y distintas, es fácil probar que las soluciones {ηin }n∈N
son linealmente independientes. En particular, en este caso tendríamos que la solución
general de (2.11) es de la forma
r
yn = ∑ Ai ηin , n ∈ N, (2.12)
i =1
Si ρ(z) tiene, por ejemplo, una raíz real doble η, le asociamos dos soluciones linealmente
independientes: {η n }n∈N y {nη n }n∈N ; si fuese triple, le añadiaríamos además la solución
{n(n − 1)η n }n∈N ; y así sucesivamente.
Si η = ρeiθ es una raíz compleja de ρ(z), tenemos dos soluciones linealmente indepen-
dientes: {ρ cos(nθ )}n∈N y {ρ sen(nθ )}n∈N . Si fuese raíz doble, tendríamos que añadir dos
nuevas soluciones, {nρ cos(nθ )}n∈N y {nρ sen(nθ )}n∈N ; y así sucesivamente.
Recapitulando, a cada una de las raíces de ρ(z) se le asocia las soluciones correspondientes
(i )
{yn }n∈N teniendo en cuenta las indicaciones anteriores. La solución general de (2.11) será en-
tonces de la forma
r
(i )
yn = ∑ Ai y n , n ∈ N,
i =1
con Ai ∈ R, para i = 1, . . . , r.
Si consideramos ahora una ecuación en recurrencia lineal no homogénea,
α r y k + r + α r −1 y k + r −1 + · · · + α 0 y k = γk + r , (2.13)
esto es, introducimos una perturbación de tamaño ε en el coeficiente A j . Operando, es claro que
yεn = yn + εη jn , n ∈ N.
Por otra parte, supongamos que η j fuese raíz real doble (el mismo razonamiento puede adap-
tarse al caso de raíces con mayor multiplicidad) y consideremos una perturbación de la solución
asociada {nη jn }n∈N , de la forma ( A j + ε)nη jn . En tal caso, se verifica que
yεn = yn + εnη jn , n ∈ N.
siendo L la constante de Lipschitz de f (t, y) respecto de y. En tal caso, el método puede escribirse
en la forma
(αr − hβ r L)yk+r + · · · + (α0 − hβ 0 L)yk = 0.
Para h pequeño se tiene que αi − hβ i L ≈ αi , por lo que el comportamiento de las raíces del nuevo
esquema será similar al del caso homogéneo.
En resumen, las condiciones de la definición de estabilidad nos permiten asegurar que las
posibles perturbaciones de la solución numérica van a permanecer acotadas.
Ejemplo 2.5. Todos los métodos basados en integración numérica (Adams-Bashforth, Adams-
Moulton, Milne-Simpson, etc.) son estables. En efecto, cualquiera de estos métodos puede escri-
birse como
y k +r − y k +r −1 = h ( β r f k +r + · · · + β 0 f k ),
por lo que su primer polinomio característico es de la forma ρ(z) = zr − zr−1 = zr−1 (z − 1). Por
tanto, z = 0 es raíz de multiplicidad r − 1 y z = 1 es raíz simple.
Ejemplo 2.6. Consideremos el método multipaso dado por
tk yk y(tk ) ek
0.0 1.000000 1.000000 0.000000
0.1 1.105171 1.105171 0.000000
0.2 1.220684 1.221403 0.000719
0.3 1.346189 1.349859 0.003670
0.4 1.478566 1.491825 0.013258
0.5 1.606461 1.648721 0.042260
0.6 1.694433 1.822119 0.127686
0.7 1.637057 2.013753 0.376696
0.8 1.126040 2.225541 1.099501
0.9 −0.734423 2.459603 3.194026
1.0 −6.541017 2.718282 9.259299
αr yk+r + · · · α0 yk = 0. (2.14)
2.3 Convergencia, estabilidad y consistencia 101
Fijemos t > 0 arbitrario y tomemos h = t/n (de esta forma, las condiciones h → 0 y n → ∞
son equivalentes). Sea µ = ηeθi una raíz del primer polinomio característico ρ(z) del método,
con η > 0 y θ ∈ [0, 2π ). Sabemos que la sucesión η n cos(nθ ) es solución de la ecuación en
recurrencia (2.14); por tanto, yn = hη n cos(nθ ) también lo es. Como, trivialmente, lı́mh→0 yi = 0
para i = 0, . . . , r − 1, entonces lı́mn→∞ yn = 0 por ser el método convergente. Hemos de distinguir
ahora entre dos casos:
Si θ = 0 o θ = π (caso real), entonces |yn | = hη n . Debido a la convergencia, se tiene que
ηn
0 = lı́m |yn | = lı́m hη n = t lı́m ,
n→∞ n→∞ n→∞ n
de donde η 6 1 necesariamente.
Supongamos θ 6= 0 y θ 6= π. Se tiene que
Debido a la convergencia,
2
y 2 − y n +1 y n −1 ηn
0 = lı́m n = lı́m h2 η 2n = t2 lı́m ,
h →0 sen2 (θ ) h →0 h →0 n
lo que implica η 6 1.
En cualquier caso, hemos probado que |µ| = η 6 1, es decir, toda raíz de ρ(z) tiene módulo
menor o igual que uno.
Supongamos ahora que µ es raíz múltiple de ρ(z). Sabemos que la sucesión nη n cos(nθ ) es
solución de (2.14), por lo que yn = h1/2 nη n cos(nθ ) también lo será. Se verifica trivialmente que
lı́mh→0 yi = 0 para i = 0, . . . , r − 1, por lo que lı́mn→∞ yn = 0 debido a la convergencia. Tenemos
de nuevo dos casos:
Si θ = 0 o θ = π, se tiene que
1 1 1
0 = lı́m |yn | = lı́m h 2 nη n = t 2 lı́m n 2 η n ,
n→∞ n→∞ n→∞
de donde η < 1.
Si θ 6= 0 y θ 6= π, consideremos la sucesión zn = h−1/2 n−1 yn , que verifica
− 12
t 1
lı́m zn = lı́m n−1 yn = lı́m (tn)− 2 yn = 0.
n→∞ n→∞ n n→∞
Hemos demostrado así que toda raíz múltiple de ρ(z) ha de tener módulo menor que uno o, lo
que es equivalente, que toda raíz de módulo uno ha de ser simple.
Para determinar las constantes Ci , denominadas constantes de error, consideramos los siguientes
desarrollos de Taylor:
hp
α j y(t + jh) = α j y(t) + jhy0 (t) + · · · + j p y( p) (t) + O(h p+1 ) ,
p!
p −1
p −1 h
0 0 00 ( p) p
hβ j y (t + jh) = hβ j y (t) + jhy (t) + · · · + j y (t) + O(h ) ,
( p − 1) !
para j = 1, . . . , r. Agrupando convenientemente, podemos deducir que
C0 = α0 + α1 + · · · + αr ,
C1 = α1 + 2α2 + · · · + rαr − ( β 0 + β 1 + · · · + β r ),
1
C2 = (α1 + 22 α2 + · · · + r2 αr ) − ( β 1 + 2β 2 + · · · + rβ r ),
2 (2.15)
..
.
1 1
C p = ( α1 + 2 p α2 + · · · + r p αr ) − ( β + 2 p −1 β 2 + · · · + r p −1 β r ).
( p − 1) ! 1
p!
Notemos que las constantes de error sólo dependen de los coeficientes αi y β i , por lo que tiene
perfecto sentido hacer la siguiente definición.
Definición 2.3 (Orden). Un método multipaso lineal tiene orden p si las constantes de error del operador
diferencia asociado verifican
C0 = C1 = · · · = C p = 0, C p+1 6= 0.
2.3 Convergencia, estabilidad y consistencia 103
C0 = 0 ⇔ α0 + · · · + αr = 0 ⇔ ρ(1) = 0,
C1 = 0 ⇔ α1 + 2α2 + · · · + rαr − ( β 0 + · · · + β r ) = 0 ⇔ ρ0 (1) − σ (1) = 0,
siendo ρ(z) y σ(z) los polinomios característicos asociados al método. Podemos hacer entonces
la siguiente definición.
Definición 2.4 (Consistencia). Un método multipaso lineal es consistente si se verifican las condiciones
cuya solución exacta es y(t) ≡ 1. Debido a la convergencia del método, si tomemos las condicio-
nes iniciales yi = 1, para i = 0, . . . , r − 1, se tiene que
lı́m yn = 1 ⇒ lı́m yn = 1,
h →0 n→∞
nh=t
2.3 Convergencia, estabilidad y consistencia 105
para cualquier t > 0 fijado. Al aplicar el método al problema considerado, aquel adopta la forma
αr yk+r + · · · + α0 yk = 0.
que tiene a y(t) = t como solución exacta. En este caso, el método tiene la siguiente forma:
α r y k +r + · · · + α0 y k = h ( β r + · · · + β 0 ).
σ (1)
yn = nh , n ∈ N,
ρ 0 (1)
es solución de la anterior ecuación en recurrencia (nótese que ρ0 (1) 6= 0 como consecuencia del
teorema 2.1 y la observación 2.12). Se verifica trivialmente que lı́mn→∞ yi = 0, i = 0, . . . , r − 1;
fijado entonces t > 0 arbitrario, la propiedad de convergencia implica que
σ (1)
lı́m yn = t ⇒ lı́m t = t ⇒ ρ 0 (1) = σ (1).
h →0 n→∞ ρ 0 (1)
nh=t
Observación 2.13. En los teoremas 2.1 y 2.2 se ha establecido que todo método multipaso lineal
que sea convergente, es necesariamente estable y consistente. El recíproco es asimismo cierto:
todo método multipaso lineal consistente y estable es convergente. La demostración de este re-
sultado es algo laboriosa: el lector interesado puede encontrarla en la referencia [11]. Debido a su
importancia, enunciamos este resultado en forma de teorema.
Teorema 2.3. Todo método multipaso lineal consistente y estable es convergente.
Ejemplo 2.9. Vimos en el ejemplo 2.5 que los métodos basados en integración numérica (Adams-
Bashforth, Adams-Moulton, etc.), cuyos polinomios característicos tenían la forma
ρ ( z ) = z r − z r −1 , σ ( z ) = β r zr + · · · + β 1 z + β 0 ,
eran todos estables. Respecto a la estabilidad, se tiene que ρ(1) = 0 y la condición ρ0 (1) = σ (1)
se reduce a:
β r + · · · + β 0 = 1.
106 Métodos multipaso para problemas de valor inicial
Bajo esta condición, por el teorema 2.3, el método considerado es convergente. En particular,
todos los métodos de Adams-Bashforth y Adams-Moulton dados en las tablas 2.1 y 2.2 son
convergentes. Lo mismo puede demostrarse para los demás métodos de diferenciación regresiva
dados en la tabla 2.3.
Ejemplo 2.10. Consideremos el método BDF3, definido por los polinomios característicos
18 2 9 2 6 3
ρ ( z ) = z3 − z + z− , σ(z) = z .
11 11 11 11
√
Es claro que z = 1 es raíz simple de ρ(z); las otras dos raíces son z = 227
± i 2239 , cuyo módulo
6
es menor que 1. En consecuencia, el método es estable. Por otra parte, ρ(1) = 0 y ρ0 (1) = 11 =
σ (1); por tanto, el método es consistente. Por el teorema 2.3, deducimos que el método BDF3 es
convergente.
Definición 2.5. Dado ĥ ∈ C, diremos que un método multipaso lineal es absolutamente estable para ĥ
si todas las raíces del polinomio de estabilidad Π(z; ĥ) tienen módulo menor que uno. Se define el dominio
de estabilidad absoluta D A del método como el conjunto de valores ĥ ∈ C para los que el método es
absolutamente estable.
Ejemplo 2.11. El método de Adams-Moulton de profundidad r = 1 tiene polinomios caracterís-
ticos ρ(z) = z − 1 y σ (z) = 12 (z + 1). La única raíz del polinomio de estabilidad absoluta
ĥ ĥ
Π(z; ĥ) = ρ(z) − ĥσ(z) = 1− z− 1+ ,
2 2
es z = 1+ĥ/2 , que tiene módulo menor que 1 si y sólo si Re(ĥ) < 0. Por tanto, el dominio de
1−ĥ/2
estabilidad absoluta es
D A = { ĥ ∈ C : Re(ĥ) < 1}.
En este caso D A = C− , por lo que el método es A-estable (definición 1.6).
Observación 2.14. El dominio de estabilidad absoluta de un método multipaso lineal puede
ser vacío, tal y como veremos en el ejemplo 2.13. Sin embargo, no puede ser todo C: puede
demostrarse que, para cualquier método, el dominio de estabilidad absoluta no puede contener
al eje real positivo en un entorno del origen. Este hecho puede constatarse en las figuras 2.2, 2.3
y 2.4.
La determinación del dominio de estabilidad absoluta D A de un método a partir de la defini-
ción no es, en general, una tarea sencilla. A continuación vamos a ver una forma alternativa para
hallar D A , de carácter geométrico, basada en la técnica conocida como localización de la frontera.
Se define la frontera ∂D A como el conjunto de valores ĥ ∈ C tales que Π(z; ĥ) tiene al menos
una raíz de módulo uno. Dicho de otra forma, ĥ ∈ ∂D A si existe θ ∈ R tal que eiθ es raíz de
Π(z; ĥ). En tal caso,
ρ(eiθ )
Π(eiθ ; ĥ) = 0 ⇔ ρ(eiθ ) − ĥσ (eiθ ) = 0 ⇔ ĥ = ,
σ(eiθ )
supuesto σ (eiθ ) 6= 0. Entonces, el conjunto ∂D A puede interpretarse como una curva en C para-
metrizada por la función
ρ(eiθ )
ĥ(θ ) = , θ ∈ R.
σ (eiθ )
Una vez calculada ∂D A , para determinar D A se considera cada una de las componentes conexas
abiertas Ci en las que ∂D A divide al plano complejo (es decir, C es la unión de ∂D A y todos
los conjuntos Ci ). A continuación, en cada componente se toma ĥi ∈ Ci arbitrario y se estudian
las raíces del correspondiente polinomio Π(z; ĥi ) (en general, dichas raíces se determinarán nu-
méricamente). Si todas las raíces tienen módulo menor que uno, entonces el conjunto Ci está
contenido en D A ; en caso contrario, dicha componente queda fuera de D A . La razón se basa en el
hecho de que las raíces de un polinomio dependen de forma continua de sus coeficientes. Por fijar
ideas, supongamos que para un cierto ĥi ∈ Ci todas las raíces de Π(z; ĥi ) tienen módulo menor
que uno. Por continuidad, al modificar el valor de ĥi siguiendo una curva arbitraria contenida en
108 Métodos multipaso para problemas de valor inicial
Ci , las raíces del polinomio correspondiente seguirán teniendo módulo menor que uno, ya que si
una raíz tuviese módulo uno significaría que hemos llegado a la frontera ∂D A . Esto prueba que si
ĥ ∈ D A ∩ Ci entonces Ci ⊂ D A . Un razonamiento similar sirve para demostrar el caso contrario.
Además, teniendo en cuenta la observación 2.14, la componente que contenga al eje real positivo
en un entorno del origen no puede estar contenida en D A .
En el apartado 2.6.1 puede encontrarse un código Python simple para representar gráficamen-
te ∂D A a partir de los polinomios característicos de un método multipaso lineal.
Ejemplo 2.12. En el ejemplo 2.11 vimos que el dominio de estabilidad absoluta del método de
Adams-Moulton de profundidad r = 1 era D A = C− . Vamos a demostrar este hecho mediante
la técnica de localización de la frontera. Teniendo en cuenta que ρ(z) = z − 1 y σ (z) = 21 (z + 1),
consideramos la curva
ρ(eiθ ) eiθ − 1
ĥ(θ ) = = 1 iθ
.
σ (eiθ ) 2 ( e + 1)
Multiplicando y dividiendo por el conjugado de eiθ + 1, obtenemos:
Este eje divide al plano complejo en dos semiplanos, C− y C+ . Debido a la observación 2.14, el
semiplano positivo C+ no puede estar contenido en D A . Si tomamos, por ejemplo, ĥ = −1 ∈ C− ,
entonces
3 1
Π(z; −2) = z − ,
2 2
cuya única raíz z = 31 tiene módulo menor que uno; por tanto −2 ∈ D A , lo que a su vez implica
C− ⊂ D A . Resumiendo, hemos probado que D A = C− .
Ejemplo 2.13. Consideremos el método multipaso lineal asociado a los polinomios
En este caso, la curva ĥ(θ ) que parametriza la frontera ∂D A puede escribirse como
3i sen θ
ĥ(θ ) = , θ ∈ R,
2 + cos θ
una vez hechas las simplificaciones√ pertinentes.
√ Es fácil ver que la imagen de la función f : R →
R, f (θ ) = 23+sen θ
cos θ , es el intervalo [− 3, 3 ] , lo que significa que
√ √
∂D A = {bi ∈ C : b ∈ [− 3, 3]}.
Por tanto, en este caso sólo hay una componente conexa de C. En virtud de la observación 2.14,
deducimos que D A = ∅.
2.4 Estabilidad absoluta 109
La técnica de localización de la frontera se ha utilizado para construir las figuras 2.2, 2.3 y 2.4,
donde se muestran los conjuntos D A y ∂D A correspondientes a los métodos de Adams-Bashforth,
Adams-Moulton y diferenciación regresiva de profundidad r = 1, 2, 3, 4. Como puede observar-
se, el dominio de estabilidad absoluta de un método de Adams-Bashforth es menor que el de
un método de Adams-Moulton con la misma profundidad, siendo este además mucho menor
que el del método de diferenciación regresiva correspondiente. Esto hace que los métodos de
diferenciación regresiva sean especialmente adecuados para el tratamiento de problemas rígidos.
Figura 2.2: Dominio D A (en gris) y frontera ∂D A (en azul) para los métodos de Adams-Bashforth de pro-
fundidad r = 1, 2, 3, 4.
Observación 2.15. Nótese que, en general, la frontera ∂D A no tiene por que coincidir con la
frontera topológica del conjunto D A . Este hecho puede observarse, por ejemplo, en el método
de Adams-Bashforth de profundidad r = 3 o r = 4 (figura 2.2), donde los dos lóbulos de ∂D A
contenidos en C+ no están incluidos en D A . Por el contrario, para el método de Adams-Bashforth
con r = 1 o r = 2 la frontera topológica de D A sí que coincide con ∂D A .
En algunos casos2 tan sólo es necesario determinar el intervalo de estabilidad absoluta, que se
define como el conjunto I A = D A ∩ R. Para ĥ ∈ R el polinomio Π(z; ĥ) tiene coeficientes reales,
2 Por ejemplo, cuando todos los autovalores (de la matriz en el caso lineal, o del jacobiano en el caso no lineal) son
reales.
110 Métodos multipaso para problemas de valor inicial
Figura 2.3: Dominio D A (en gris) y frontera ∂D A (en azul) para los métodos de Adams-Moulton de profun-
didad r = 1, 2, 3, 4.
por lo que es posible utilizar el criterio de Routh-Hurwitz para determinar cuándo las raíces de
Π(z; ĥ) tienen módulo menor que uno. Para ello, primero consideramos la transformación
z−1
w= ,
z+1
1+w
def r
R ( w ) = (1 − w ) Π ; ĥ = a 0 w r + a 1 w r −1 + · · · + a r −1 w + a r ,
1−w
donde r es la profundidad del método. Entonces, condiciones suficientes para que R(w) tenga
todas sus raíces en C− (lo que equivale a que Π(z; ĥ) tenga todas sus raíces en ∆) son:
ai > 0 ∀ i, en el caso r = 2.
Figura 2.4: Dominio D A (en gris) y frontera ∂D A (en azul) para los métodos de diferenciación regresiva de
profundidad r = 1, 2, 3, 4.
De esta forma,
1+w
3
R ( w ) = (1 − w ) Π ; ĥ = a0 z3 + a1 z2 + a2 z + a3 ,
1−w
αr y k +r + · · · + α0 y k = h ( β r f k +r + · · · + β 0 f k ),
donde
Ck = h( β r−1 f k+r−1 + · · · + β 0 f k ) − (αr−1 yk+r−1 + · · · + α0 yk ).
La igualdad (2.16) puede interpretarse como un problema de punto fijo cuya solución es yk+r . Al
estar Ck fija, la función de iteración es contractiva si h verifica la condición hβ r L < 1; en tal caso,
la solución yk+r del problema es única.
Observación 2.18. En general, el valor yk+r puede aproximarse resolviendo la ecuación (2.16)
mediante un método iterativo adecuado, tomando como punto inicial el último valor conocido
yk+r−1 : esta es la forma habitual de implementar un método implícito.
2.5 Métodos predictor-corrector 113
Los métodos predictor-corrector combinan las mejores características de los métodos explíci-
tos (su simplicidad) e implícitos (su mayor dominio de estabilidad absoluta). Para construir un
método predictor-corrector se considera, en primer lugar, una predicción ykP+r del valor yk+r de-
terminada mediante un método explícito, y a continuación se corrige mediante una iteración del
método implícito, de la forma
yC P
k +r = hβ r f ( tk + rh, yk +r ) + Ck .
El valor así obtenido es el que tomamos como yk+r . Este tipo de métodos suele denotarse con las
siglas PECE (del inglés Predict-Evaluate-Correct-Evaluate).
Observación 2.19. En general, puede realizarse más de una corrección, obteniéndose así un mé-
todo de la forma
(0) ( i +1) (i )
yk+r = ykP+r , yk+r = hβ r f (tk + rh, yk+r ) + Ck , i = 0, 1, . . . , m − 1.
(m)
Se toma entonces yk+r como el último valor obtenido yk+r . El método resultante se denota como
( i +1) (i )
P(EC)m E. También pueden realizarse iteraciones hasta que la diferencia |yk+r − yk+r | sea me-
nor que una tolerancia establecida; en este caso, el número de iteraciones a realizar suele estar
limitado para evitar que el proceso de convergencia al punto fijo se ralentice en exceso.
Ejemplo 2.15. Consideremos el par predictor-corrector AB2-AM3. La predicción se realiza con el
método de Adams-Bashforth AB2:
h
ykP+2 = yk+1 + (3 f k+1 − f k ).
2
A continuación, se realiza la corrección con el método de Adams-Moulton AM3:
h
y k +2 ≡ y C
k +2 = y k +1 + (5 f (tk+2 , ykP+2 ) + 8 f k+1 − f k ).
12
En la figura 2.5 se muestra el dominio de estabilidad absoluta del método predictor-corrector
AB2-AM3. Compárese con las figuras 2.2 y 2.3.
Ejemplo 2.16. Vimos en el capítulo 1 que el método de Heun podía escribirse en la forma
p k +1 = y k + h f ( t k , y k ),
h
y
k +1 = y k + f (tk , yk ) + f (tk + h, pk+1 ) .
2
Así pues, el método de Heun puede interpretarse como un método predictor-corrector, donde el
método de Euler actúa como predictor y el método del trapecio como corrector.
Dados un método multipaso lineal explícito de orden p − 1 y otro método multipaso lineal
implícito de orden p, puede demostrarse que el orden del par predictor-corrector que forman
tiene orden p. Por esta razón es habitual considerar métodos predictor-corrector de la forma
AB( p − 1)-AM( p).
Ejemplo 2.17. El método predictor-corrector AM2-AB3 considerado en el ejemplo anterior tiene
orden p = 3.
114 Métodos multipaso para problemas de valor inicial
# polinomios caracteristicos
rho = array ([1. , -1. , 0. , 0. , 0.]) # primero
sigma = array ([0. , 55. , -59. , 37. , -9.])/24. # segundo
plot (x , y )
grid ( True )
axis ( ' equal ')
show ()
En los arrays rho y sigma definimos los coeficientes de los polinomios característicos ρ y σ (el
primer elemento es el coeficiente líder). A continuación se evalúan dichos polinomios en eθi me-
diante la instrucción polyval, tomando los valores de θ en una partición uniforme del intervalo
[0, 2π ]; el elemento 1j representa la unidad imaginaria. Por último, se define la parametrización
de ∂D A dada por ĥ ≡htilde, y se dibuja la curva que determina.
12 def exacta ( t ):
13 """ Solucion exacta del problema de valor inicial """
14 return exp ( - t )* sin ( t )
15
16 def AB2 (a , b , N , y0 ):
17 """ Implementacion del metodo AB2 """
18
19 t = zeros ( N +1)
20 y = zeros ( N +1)
21
22 h = (b - a )/ N
23
24 t [0] = a
25 y [0] = y0
26
27 # calculo de y1 mediante el metodo de Heun
28 t [1] = t [0]+ h
29 valf0 = fun ( t [0] , y [0])
30 aux = y [0]+ h * valf0
31 y [1] = y [0]+0.5* h *( valf0 + fun ( t [1] , aux ))
32
33 # metodo AB2
34 for k in range (N -1):
35 t [ k +2] = t [ k +1]+ h
36 y [ k +2] = y [ k +1]+0.5* h *(3.0* fun ( t [ k +1] , y [ k +1])
37 - fun ( t [ k ] , y [ k ]))
38
39 return (t , y )
40
41 # ## Programa principal ###
42
43 # Datos del problema
44 a = 0. # extremo inferior del intervalo
45 b = 5. # extremo superior del intervalo
46 N = 50 # numero de particiones
47 y0 = 0. # condicion inicial
48
49 tini = clock () # leemos tiempo de inicio
50
51 (t , y ) = AB2 (a , b , N , y0 ) # llamada al metodo AB2
52
53 tfin = clock () # leemos tiempo final
54
55 ye = exacta ( t ) # calculo de la solucion exacta
56
57 # Calculo del error cometido
58 error = max ( abs (y - ye ))
59
60 # Resultados
61 print ( ' ----- ')
62 print ( ' Tiempo CPU : ' + str ( tfin - tini ))
63 print ( ' Error : ' + str ( error ))
64 print ( ' Paso de malla : ' + str (( b - a )/ N ))
118 Métodos multipaso para problemas de valor inicial
t = zeros ( r )
y = zeros ( r )
t [0] = t0
y [0] = y0
return (t , y )
La salida de la función son los nodos t0 , . . . , tr−1 y los valores y0 , . . . , yr−1 , para un r > 1 dado.
El método AM3 se define en la siguiente función:
1 def AM3 (a , b , N , y0 ):
2 """ Implementacion del metodo AM3
3 usando iteraciones de punto fijo """
4
5 t = zeros ( N +1)
6 y = zeros ( N +1)
7
8 tol = 1. e -6 # tolerancia punto fijo
9 M = 100 # numero maximo de iteraciones de punto fijo
10
11 h = (b - a )/ N
12
2.6 Implementación en Python de métodos multipaso 119
13 # puntos iniciales
14 r = 2 # profundidad del metodo
15 ( t [0: r ] , y [0: r ]) = RK4_ini (a , y0 , h , r ) # valores iniciales
16
17 # metodo AM3
18 for k in range (N -1):
19 t [ k +2] = t [ k +1]+ h
20 C = y [ k +1]+ h *(8.0* fun ( t [ k +1] , y [ k +1]) - fun ( t [ k ] , y [ k ]))/12.0
21 wold = y [ k +1]
22 error = tol +1
23 n = 0
24
25 while (( error > tol ) and ( n < M )):
26 wnew = 5.0* h * fun ( t [ k +2] , wold )/12.0+ C
27 error = abs ( wnew - wold )
28 wold = wnew
29 n += 1
30
31 if ( n == M ):
32 print ( 'k = '+ str ( k )+ '. Aviso : no . maximo de iteraciones ')
33
34 y [ k + r ] = wnew
35
36 return (t , y )
Tanto la tolerancia tol como el número máximo de iteraciones M son arbitrarios, y pueden de-
pender del problema considerado. En la línea 15 se calculan los dos puntos iniciales necesarios
usando la función RK4_ini definida anteriormente. El método AM3 propiamente dicho comienza
en la línea 18, donde el último valor que toma k es precisamente N − 2. La implementación del
método AM3 se hace siguiendo las directrices de la observación 2.2 (véase también la sección
1.9.5), utilizando iteraciones de punto fijo (bucle while).
El código completo es el siguiente:
# Resolucion del problema de valor inicial
# y '= f (t , y ) , y ( t0 )= y0 ,
# mediante el metodo AM3 .
def fun (t , y ):
""" Funcion f (t , y ) que define la ecuacion diferencial """
return -y + exp ( - t )* cos ( t )
def exacta ( t ):
""" Solucion exacta del problema de valor inicial """
return exp ( - t )* sin ( t )
def RK4_ini ( t0 , y0 , h , r ):
""" Calculo de r puntos iniciales mediante el metodo RK4 """
t = zeros ( r )
120 Métodos multipaso para problemas de valor inicial
y = zeros ( r )
t [0] = t0
y [0] = y0
return (t , y )
def AM3 (a , b , N , y0 ):
""" Implementacion del metodo AM3
usando iteraciones de punto fijo """
t = zeros ( N +1)
y = zeros ( N +1)
h = (b - a )/ N
# puntos iniciales
r = 2 # profundidad del metodo
( t [0: r ] , y [0: r ]) = RK4_ini (a , y0 , h , r ) # valores iniciales
# metodo AM3
for k in range (N -1):
t [ k +2] = t [ k +1]+ h
C = y [ k +1]+ h *(8.0* fun ( t [ k +1] , y [ k +1]) - fun ( t [ k ] , y [ k ]))/12.0
wold = y [ k +1]
error = tol +1
n = 0
if ( n == M ):
print ( 'k = '+ str ( k )+ '. Aviso : no . maximo de iteraciones ')
y [ k + r ] = wnew
return (t , y )
2.6 Implementación en Python de métodos multipaso 121
# Resultados
print ( ' ----- ')
print ( ' Tiempo CPU : ' + str ( tfin - tini ))
print ( ' Error : ' + str ( error ))
print ( ' Paso de malla : ' + str (( b - a )/ N ))
print ( ' ----- ')
Si en lugar de usar iteraciones de punto fijo para calcular yk+2 quisiéramos utilizar el método
de Newton, la función AM3 correspondiente sería la siguiente:
def AM3 (a , b , N , y0 ):
""" Implementacion del metodo AM3
usando el metodo de Newton """
t = zeros ( N +1)
y = zeros ( N +1)
h = (b - a )/ N
# puntos iniciales
r = 2 # profundidad del metodo
( t [0: r ] , y [0: r ]) = RK4_ini (a , y0 , h , r ) # valores iniciales
122 Métodos multipaso para problemas de valor inicial
# metodo AM3
for k in range (N -1):
t [ k +2] = t [ k +1]+ h
C = y [ k +1]+ h *(8.0* fun ( t [ k +1] , y [ k +1]) - fun ( t [ k ] , y [ k ]))/12.0
df = 1.0+5.0* h /12.0;
wold = y [ k +1]
error = tol +1
n = 0
if ( n == M ):
print ( 'k = '+ str ( k )+ '. Aviso : no . maximo de iteraciones ')
y [ k + r ] = wnew
return (t , y )
∂f
En este caso, la derivada df≡ ∂y es constante, por lo que podemos definirla fuera del bucle while,
para calcularla tan solo una vez. Cuando la derivada no sea constante, será necesario definirla
dentro de dicho bucle, a continuación de la definición de f.
Consideremos ahora el método de diferenciación regresiva BDF3:
18 9 2 6
y k +3 − y k +2 + y k +1 − y k = hf , k = 0, 1, . . . , N − 3.
11 11 11 11 k+3
Su implementación en Python, usando iteraciones de punto fijo, es la siguiente:
def BDF3 (a , b , N , y0 ):
""" Implementacion del metodo BDF3
usando iteraciones de punto fijo """
t = zeros ( N +1)
y = zeros ( N +1)
h = (b - a )/ N
# puntos iniciales
r = 3 # profundidad del metodo
( t [0: r ] , y [0: r ]) = RK4 (a , y0 , h , r ) # valores iniciales
# metodo BDF3
for k in range (N - r +1):
t [ k + r ] = t [ k +r -1]+ h
2.6 Implementación en Python de métodos multipaso 123
if ( n == M ):
print ( 'k = '+ str ( k )+ '. Aviso : no . maximo de iteraciones ')
y [ k + r ] = wnew
return (t , y )
h
ykP+2 = yk+1 + (3 f k+1 − f k ),
2
C h
y
k +2 ≡ y k +2 = y k +1 + (5 f (tk+2 , ykP+2 ) + 8 f k+1 − f k ).
12
Su implementación sería la siguiente:
def AB2_AM3 (a , b , N , y0 ):
""" Implementacion del metodo AB2 - AM3 """
t = zeros ( N +1)
y = zeros ( N +1)
h = (b - a )/ N
# puntos iniciales
( t [0:2] , y [0:2]) = RK4_ini (a , y0 , h , 2) # dos valores iniciales
return (t , y )
124 Métodos multipaso para problemas de valor inicial
x 0 = A + x2 y − ( B + 1) x,
(
y0 = Bx − x2 y,
def fun (t , z ):
""" Funcion que define el sistema """
A = 1.
B = 3.
return array ([ A + z [0]**2* z [1] -( B +1)* z [0] , B * z [0] - z [0]**2* z [1]])
def RK4_ini ( t0 , y0 , h , r ):
""" Calculo de r puntos iniciales mediante el metodo RK4 """
t = zeros ( r )
z = zeros ([ len ( z0 ) , r ])
t [0] = t0
z [: , 0] = z0
return (t , z )
def BDF3 (a , b , N , z0 ):
""" Metodo BDF3 con iteraciones de punto fijo """
t = zeros ( N +1)
z = zeros ([ len ( z0 ) , N +1])
h = (b - a )/ N
# puntos iniciales
r = 3 # profundidad del metodo
2.6 Implementación en Python de métodos multipaso 125
( ti , zi ) = RK4_ini (a , z0 , h , r )
for i in range ( r ):
t[i] = ti [ i ]
z [: , i ] = zi [: , i ]
# metodo BDF3
for k in range (N - r +1):
t [ k + r ] = t [ k +r -1]+ h
C = 18.0* z [: , k +2] -9.0* z [: , k +1]+2.0* z [: , k ]
wold = z [: , k +2]
error = tol +1
n = 0
if ( n == M ):
print ( 'k = '+ str ( k )+ '. Aviso : no . maximo de iteraciones ')
z [: , k + r ] = wnew
return (t , z )
tini = clock ()
(t , z ) = BDF3 (a , b , N , z0 )
tfin = clock ()
time = tfin - tini
subplot (2 , 2 , 1)
plot (t , z [0])
xlabel ( 't ')
ylabel ( 'x ')
legend ([ 'x ( t ) ' ])
grid ( True )
126 Métodos multipaso para problemas de valor inicial
subplot (2 , 2 , 2)
plot (t , z [1])
xlabel ( 't ')
ylabel ( 'y ')
legend ([ 'y ( t ) ' ])
grid ( True )
subplot (2 , 2 , 3)
plot (t , z [0])
plot (t , z [1])
xlabel ( 't ')
ylabel ( 'x , y ')
legend ([ 'x ( t ) ' , 'y ( t ) ' ])
grid ( True )
subplot (2 , 2 , 4)
plot ( z [0] , z [1])
xlabel ( 'x ')
ylabel ( 'y ')
legend ([ ' orbita ' ])
grid ( True )
3
MÉTODOS NUMÉRICOS PARA
PROBLEMAS DE CONTORNO
En los capítulos anteriores hemos estudiado métodos numéricos para la resolución de proble-
mas de Cauchy, o de valor inicial, de la forma
(
y0 (t) = f (t, y(t)), t ∈ I = [t0 , t0 + T ],
y ( t0 ) = y0 .
En las aplicaciones es habitual interpretar la variable t como el tiempo, mientras que y representa
el estado del sistema. De esta forma, y0 es el estado inicial del sistema y la ecuación diferencial
nos permite determinar la evolución temporal del mismo. Este hecho se refleja en los métodos
numéricos asociados, en los que la aproximación a la solución en un nodo se determina a partir
de uno (métodos unipaso) o varios (métodos multipaso) datos anteriormente calculados.
Otro tipo de problema que aparece con frecuencia en la práctica es aquel en el que se conocen
los valores de la solución en los extremos del intervalo de definición del problema. En tal caso,
estamos ante un problema de contorno o de valor en la frontera, que puede escribirse en la forma
00 0
u ( x ) = f ( x, u( x ), u ( x )), x ∈ ( a, b),
u( a) = α, (3.1)
u(b) = β.
difieren notablemente de las de problemas de valor inicial. En este capítulo nos centraremos fun-
damentalmente en el estudio del método de diferencias finitas. Asimismo, veremos cómo extender
el método de diferencias finitas para resolver problemas bidimensionales, y lo aplicaremos para
resolver algunos problemas clásicos relativos a la ecuación de Poisson, la ecuación del calor y la
ecuación de ondas. Asimismo, trataremos brevemente otro método clásico: el método del tiro.
Una dificultad teórica que presentan los problemas de contorno es que, a diferencia de los
problemas de valor inicial, no existen resultados generales de existencia y unicidad de soluciones
(aunque sí pueden establecerse resultados de este tipo en determinados casos particulares). Por
ello, a lo largo del capítulo supondremos que los problemas de contorno considerados siempre
poseen solución única.
3.1. Motivación
Vamos a motivar el contenido de este capítulo mediante el estudio de las soluciones estacio-
narias de la ecuación del calor. Consideremos el flujo de calor en una varilla formada por un
material conductor, sujeta a alguna fuente externa de calor, y ciertas condiciones de contorno.
Supongamos que las propiedades del material, la distribución inicial de temperaturas y la fuente
externa sólo varían en la dirección longitudinal x. En tal caso, es plausible pensar que la dis-
tribución de temperaturas en cualquier instante dependerá tan sólo de x, por lo que podemos
modelizar el sistema mediante una ecuación unidimensional.
Llamemos u( x, t) a la temperatura de la barra en el punto x ∈ ( a, b) (donde a y b representan
los extremos de la barra) y en el instante t > 0. Dicha función ha de ser solución de la ecuación
del calor, que es una ecuación en derivadas parciales que tiene la forma
∂u ∂2 u
( x, t) − µ 2 ( x, t) = f ( x, t),
∂t ∂x
donde µ es el coeficiente de conductividad térmica (que depende del material, que supondremos
homogéneo) y f ( x, t) es la fuente de calor externa. La ecuación se complementa con una condición
inicial
u( x, 0) = u0 ( x ),
que nos dice la distribución inicial de temperaturas, y dos condiciones de contorno de tipo Dirichlet,
que especifican los valores de la temperatura en los extremos de la barra, en cualquier instante
de tiempo t.
Observación 3.1. Otro tipo de condición de contorno que puede considerarse es una condición
de tipo Neumann, en la que se especifica el valor de la derivada ∂ x u en un extremo del intervalo.
Por ejemplo, la condición ∂ x u( a, t) = 0 indicaría que no hay flujo de calor a través del extremo
x = a (condición de aislamiento).
En general, la distribución de temperatura variará con el tiempo. Sin embargo, si f ( x, t), α(t)
y β(t) son independientes de t, es de esperar que la solución u( x, t) alcance un estado estacionario
3.2 El método del tiro 129
u( x ), que permanece invariable en tiempos posteriores. Para determinar dicho estado estacionario
podemos hacer ∂t u = 0 en la ecuación, que se reduce entonces a la ecuación de Poisson:
−u00 ( x ) = ϕ( x ),
donde
1
ϕ( x ) = f ( x ).
µ
El problema así obtenido,
00
− u ( x ) = ϕ ( x ), x ∈ ( a, b),
u( a) = α,
u(b) = β,
con β̄ a determinar. Si denotamos por u( x; β̄) la solución del problema de Cauchy, determinamos
entonces β̄ de modo que se verifique la condición
u(b; β̄) = β. (3.2)
De esta forma, la solución u( x; β̄) también lo será del problema (3.1).
La ecuación (3.2) se resuelve mediante algún método iterativo: punto fijo, Newton, secante,
etc. Nótese que en cada iteración es necesario resolver numéricamente un problema de Cauchy
para determinar u(b; β̄). La convergencia del método del tiro dependerá del método iterativo
considerado, así como del método elegido para resolver los problemas de Cauchy.
Ejemplo 3.1. Apliquemos el método del tiro para resolver el problema de contorno
00 0 2
u ( x ) + u ( x ) + u( x ) = −( x + x + 1), x ∈ (0, 1),
u(0) = 0,
u(1) = 0.
u(1; β̄) = 0 ⇒ β̄ = 1.
u( x ) ≡ u( x; 1) = x − x2 .
Nótese que, numéricamente, la ecuación u(1; β̄) = 0 se resolvería mediante un método iterativo
en el cual, en cada paso, sería necesario resolver, también de forma numérica, el problema de
Cauchy correspondiente
Como se muestra en el siguiente ejemplo, el método del tiro puede ser inestable cuando se
aplica a problemas rígidos.
Ejemplo 3.2. Consideremos el problema de contorno dado por
00
u ( x ) = 100u, x ∈ (0, 1),
u(0) = 1,
u(1) = 0,
que tiene a
10 − β̄ −10x 10 + β̄ 10x
u( x; β̄) = e + e
20 20
como solución exacta. Resolviendo la ecuación (3.2), se obtiene:
1 + e−20
u(1; β̄) = 0 ⇒ β̄ = −10 ≈ −10.0000000412.
1 − e−20
La solución u( x; β̄) es muy sensible a la elección de β̄, ya que crece como e10x . Por ejemplo, se
verifica que
u(1; −9.99) ≈ 11.0133.
3.3 Aproximaciones por diferencias finitas 131
Figura 3.1: Método del tiro aplicado a un problema rígido. Comparación de la solución exacta con la solución
aproximada u( x; −9.99).
De esta forma, en la práctica los errores de redondeo van a producir inestabilidades en la solución
cuando se aplica el método del tiro. En la gráfica 3.1 se comparan la solución exacta del problema
con la solución obtenida para β̄ = −9.99. Como puede observarse, en este caso el método del tiro
es completamente inestable.
Observación 3.2. Una variante del método del tiro es el método del tiro múltiple, en el que se
consideran una serie de puntos intermedios xk ∈ ( a, b) y aproximaciones uk en estos puntos,
que se conectan mediante soluciones de problemas de Cauchy adecuados siguiendo la idea del
método del tiro simple.
u( x̄ + h) − u( x̄ )
u0 ( x̄ ) ≈ D+ u( x̄ ) = ,
h
132 Métodos numéricos para problemas de contorno
para h > 0 «pequeña». El error cometido en la aproximación puede medirse usando un desarrollo
de Taylor. Concretamente, suponiendo que u ∈ C3 ([ a, b]), podemos escribir
h2 00
u( x̄ + h) = u( x̄ ) + hu0 ( x̄ ) + u ( x̄ ) + O(h3 ),
2
de donde
h 00
|u ( x̄ )| + O(h2 ) = O(h).
| D+ u( x̄ ) − u0 ( x̄ )| =
2
En este sentido, la aproximación D+ u( x̄ ) es de primer orden. Razonando de manera análoga, otra
aproximación de primer orden viene dada por
u( x̄ ) − u( x̄ − h)
u0 ( x̄ ) ≈ D− u( x̄ ) = .
h
Las aproximaciones D+ u( x̄ ) y D+ u( x̄ ) se denominan descentradas. Geométricamente, D+ u( x̄ ) es
la pendiente de la recta que interpola a u( x ) en los puntos x̄ y x̄ + h, mientras que D− u( x̄ ) es
la pendiente correspondiente a los puntos x̄ − h y x̄ (figura 3.2). Estas pendientes aproximan a
la pendiente de la recta tangente a u( x ) en el punto x̄, que es la interpretación geométrica de la
derivada u0 ( x̄ ).
También es posible considerar una aproximación centrada, de la forma
u( x̄ + h) − u( x̄ − h) D+ u( x̄ ) + D− u( x̄ )
Dc u( x̄ ) = = . (3.3)
2h 2
En la figura 3.2 parece que la aproximación centrada Dc u( x̄ ) proporciona una mejor aproximación
a u0 ( x̄ ) que cualquiera de las dos aproximaciones descentradas, D+ u( x̄ ) y D− u( x̄ ). Para confirmar
3.3 Aproximaciones por diferencias finitas 133
este hecho, supongamos que u ∈ C4 ([ a, b]) y consideremos los siguientes desarrollos de Taylor:
h2 00 h3 000
u( x̄ + h) = u( x̄ ) + hu0 ( x̄ ) + u ( x̄ ) + u ( x̄ ) + O(h4 ),
2 6 (3.4)
h2 h3 000
u( x̄ − h) = u( x̄ ) − hu0 ( x̄ ) + u00 ( x̄ ) − u ( x̄ ) + O(h4 ).
2 6
Restando ambas expresiones, resulta
h3 000
u( x̄ + h) − u( x̄ − h) = 2hu0 ( x̄ ) + u ( x̄ ) + O(h4 ),
3
de donde
h2 000
|u ( x̄ )| + O(h4 ) = O(h2 ).
| Dc u( x̄ ) − u0 ( x̄ )| =
3
De esta forma, Dc u( x̄ ) constituye una aproximación de segundo orden a u0 ( x̄ ).
Observación 3.3. Pueden construirse aproximaciones de orden arbitrario a la derivada u0 ( x̄ ),
considerando un número suficiente de puntos. Por ejemplo, la fórmula
1
D3 u( x̄ ) = (2u( x̄ + h) + 3u( x̄ ) − 6u( x̄ − h) + u( x̄ − 2h))
6h
proporciona una aproximación (descentrada) de tercer orden, a costa de utilizar cuatro puntos
de interpolación.
Las aproximaciones de la primera derivada construidas anterioremente se denominan fórmulas
de diferencias finitas. Una técnica para construir fórmulas generales de diferencias finitas es el
método de coeficientes indeterminados. Para ilustrar su funcionamiento, supongamos que queremos
construir una aproximación de u0 ( x̄ ) del mayor orden posible, usando para ello los puntos u( x̄ ),
u( x̄ − h) y u( x̄ − 2h). La fórmula buscada tendrá la forma
D2 u( x̄ ) = au( x̄ ) + bu( x̄ − h) + cu( x̄ − 2h),
donde a, b, c ∈ R son coeficientes a determinar. Consideremos los desarrollos de Taylor
h2 00 h3
u( x̄ − h) = u( x̄ ) − hu0 ( x̄ ) + u ( x̄ ) − u000 ( x̄ ) + O(h4 ),
2 6
4h2 00 8h3 000
0
u( x̄ − 2h) = u( x̄ ) − 2hu ( x̄ ) + u ( x̄ ) − u ( x̄ ) + O(h4 ),
2 6
donde hemos supuesto u ∈ C4 ([ a, b]). Entonces, agrupando términos convenientemente, pode-
mos escribir
D2 u( x̄ ) = ( a + b + c)u( x̄ ) − (b + 2c)hu0 ( x̄ ) + 12 (b + 4c)h2 u00 ( x̄ ) − 16 (b + 8c)h3 u000 ( x̄ ) + O(h4 ).
Para que D2 u( x̄ ) aproxime a u0 ( x̄ ) con el mayor orden posible, hacemos
a + b + c = 0
1 3 2 1
b + 2c = − ⇒a= , b=− , c= .
h 2h h 2h
b + 4c = 0
134 Métodos numéricos para problemas de contorno
1
D2 u( x̄ ) = (3u( x̄ ) − 4u( x̄ − h) + u( x̄ − 2h)), (3.5)
2h
que resulta ser de segundo orden:
1 h2
| D2 u( x̄ ) − u0 ( x̄ )| = − (b + 8c)h3 u000 ( x̄ ) + O(h4 ) = u000 ( x̄ ) + O(h3 ) = O(h2 ).
6 12
Observación 3.4. Una forma alternativa de construir fórmulas de diferencias finitas consiste en
aproximar la función u( x ) mediante un polinomio de interpolación p( x ), y después aproximar
u0 ( x̄ ) mediante p0 ( x̄ ).
La fórmula (3.6) puede obtenerse de forma sencilla aplicando el método de coeficientes indeter-
minados o mediante interpolación. Otra manera de deducirla sería la siguiente:
d 0 D− u( x̄ + h) − D− u( x̄ )
u00 ( x̄ ) = (u ( x̄ )) ≈ D+ ( D− u( x̄ )) =
dx h
1 u( x̄ + h) − u( x̄ ) u( x̄ ) − u( x̄ − h)
= − = D2 u( x̄ ),
h h h
siendo h la distancia (uniforme) entre dos puntos consecutivos. Nótese que x̄ no tiene por que
coincidir con ningún x j . Para cada i = 1, . . . , n, se verifica que
n
(
1 1 si i − 1 = k,
( i − 1) ! ∑ c j (x j − x̄)i−1 = 0 si no.
j =1
Al ser los puntos x j distintos, se obtiene un sistema lineal n × n de Vandermonde con solución
única {c1 , . . . , cn }. En general, el orden de la fórmula será al menos p > n − k.
3.4 El método de diferencias finitas 135
00
− u ( x ) = f ( x ), x ∈ ( a, b),
u( a) = α, (3.7)
u(b) = β,
donde f ∈ C ([ a, b]) y α, β ∈ R. Recordemos que la solución del problema (3.7) podía interpre-
tarse como el estado estacionario de la ecuación del calor en una varilla cuando se fijaban las
temperaturas en los extremos de la misma (véase la sección 3.1).
Consideremos una partición uniforme del intervalo [ a, b],
con paso de malla h > 0. Queremos obtener aproximaciones u j del valor exacto u( x j ), para cada
j = 0, . . . , m + 1. En particular, debido a las condiciones de contorno, se tiene que
u( x0 ) = u( a) = α, u( xm+1 ) = u(b) = β,
−u00 ( x j ) = f ( x j ), j = 1, . . . , m.
u j+1 − 2u j + u j−1
− = fj, j = 1, . . . , m. (3.8)
h2
donde, por consistencia en la notación, hemos definido f j = f ( x j ). Se trata de un sistema lineal
de m ecuaciones, cuyas incógnitas son los m valores {u1 , . . . , um }. Nótese que, imponiendo las
condiciones de contorno u0 = α y um+1 = β, la primera y la última ecuaciones del sistema pueden
reducirse, repectivamente, a
α
u − 2u1 +
u0 u − 2u
− 2 = f1 ⇒ − 2 2 1 = f1 + 2
α
h2 h h
136 Métodos numéricos para problemas de contorno
y
β
u − 2um + um−1 −2um + um−1
− m +1
β
= fm ⇒ − = fm + 2 .
h2 h2 h
En forma matricial, el sistema puede escribirse como
AU = B, (3.9)
donde
α
f1 + 2 −1
u1 h2
f −1 2 −1
u2
2
−1 2 −1
.. 1
..
U= , B= , A= 2 ∈ M m (R).
. . .. .. ..
h . . .
f m −1
u m −1
−1 2 −1
um β
fm + 2 −1 2
h
Puede verse que la matriz A es simétrica y definida positiva, lo que implica que el sistema (3.9)
posee solución única dada por
U = A−1 B.
Observación 3.5. Alternativamente, es sencillo demostrar que
m+1
det( A) = 6= 0,
h2m
por lo que la matriz A es invertible y el sistema (3.9) tiene solución única.
Observación 3.6. Aunque se ha considerado la ecuación diferencial −u00 ( x ) = f ( x ) para introdu-
cir el método de diferencias finitas, este puede extenderse de forma natural al caso de ecuaciones
lineales generales, así como para tratar problemas no lineales. Más adelante estudiaremos estas
extensiones del método.
Observación 3.7. Es importante notar que, en la práctica, el número de incógnitas m puede
ser muy elevado, por lo que la elección de un método apropiado para resolver el sistema es
fundamental. Dado que la matriz A es tridiagonal, podemos aprovechar que existen técnicas
específicas para la resolución de este tipo de sistemas, mucho más eficientes que una regla de
carácter general.
Vimos en la sección 3.3 que la aproximación D2 u( x̄ ) era de segundo orden. Tiene entonces
sentido plantearse si el método de diferencias finitas es también de segundo orden, en un sentido
aún por especificar. Para aclarar esta cuestión, vamos a realizar un análisis del error cometido.
Dado el vector Ue = (u( x1 ), . . . , u( xm ))t de valores exactos, definimos el vector error
E = U − Ue ,
cuyas componentes contienen los errores puntuales cometidos en cada nodo interior del mallado,
e j = u j − u( x j ). Para cuantificar la magnitud del vector E, tendremos que considerar una norma
vectorial apropiada.
3.4 El método de diferencias finitas 137
Observación 3.8. En general, la elección de una norma dependerá tanto del problema que se esté
tratando como de las técnicas disponibles para obtener acotaciones. Algunas normas habituales
son la norma infinito (o del máximo), la norma-1 y la norma-2, que se definen, respectivamente,
como1 m m 1/2
k Ek∞ = máx |e j |, k E k1 = h ∑ | e j |, k E k2 = h ∑ | e j |2 ,
16 j 6 m j =1 j =1
para cada E ∈ Rm .
Para estimar el error global E usaremos una técnica que es básica en el análisis de métodos de
diferencias finitas. Primero estimaremos el error de discretización local del método, y luego usare-
mos un resultado de estabilidad para acotar el error global en función del error de discretización
local.
El error de discretización local ε j se define a partir de cada ecuación en (3.8), sustituyendo los
valores aproximados por valores exactos:
h2 h 2 (4)
ε j = − f ( x j ) − u ( x j ) + u (4) ( x j ) + O ( h 4 )
00
=− u ( x j ) + O ( h4 ), j = 1, . . . , m,
12 12
ε j = O ( h2 ), j = 1, . . . , m. (3.10)
E = − B + AUe ⇒ AUe = B + E .
Por tanto,
AE = AU − AUe = B − ( B + E ) = −E ⇒ AE = −E .
1 Estas normas se obtienen como discretizaciones de las normas usuales en los espacios de funciones L p ( a, b ), para
p = ∞, 1, 2, respectivamente:
Z b Z b
1/2
kek∞ = máx |e( x )|, e ∈ L∞ ( a, b); k e k1 = |e( x )| dx, e ∈ L1 ( a, b); k e k2 = |e( x )|2 dx , e ∈ L2 ( a, b).
x ∈( a,b) a a
138 Métodos numéricos para problemas de contorno
e0 = 0, em+1 = 0,
Ah Eh = −Eh .
Por tanto, si k · k es una norma vectorial y denotamos de la misma forma a la norma matricial
asociada, se tiene que
1
k Eh k 6 k A −
h kkE h k. (3.11)
Llegados a este punto, es conveniente establecer una serie de definiciones. Nótese que estas
definiciones se aplican a un problema lineal general.
Definición 3.1 (Estabilidad). Supongamos que un método de diferencias finitas para un problema de
contorno lineal general genera un conjunto de ecuaciones Ah Uh = Bh , siendo h > 0 el paso de malla.
Diremos que el método es estable si existen constantes h∗ > 0 y C > 0 tales que la matriz Ah es invertible
para h < h∗ , y
1
k A−h k 6 C, ∀ h < h∗ .
Definición 3.2 (Consistencia y orden). Un método de diferencias finitas es consistente con la ecuación
diferencial y las condiciones de contorno si
lı́m kEh k = 0.
h →0
Si se verifica kEh k = O(h p ) para un cierto p ∈ N, diremos que el método tiene orden p.
Definición 3.3 (Convergencia). Un método de diferencias finitas es convergente si
lı́m k Eh k = 0.
h →0
Observación 3.9. Es importante notar que las definiciones anteriores dependen de la elección de
la norma k · k. Aunque en el espacio Rm todas las normas son equivalentes, hay que tener en
cuenta que las constantes de equivalencia pueden depender de h (observación 3.8).
Observación 3.10. En general, probar la consistencia de un método suele ser la parte más sim-
ple. La estabilidad, por el contrario, es la parte que normalmente presenta mayores dificultades
técnicas.
3.4 El método de diferencias finitas 139
Observación 3.11. Bajo ciertas condiciones, puede probarse que si un método de diferencias
finitas es estable y de orden p, entonces el error global Eh es O(h p ).
Volviendo a la desigualdad (3.11), si suponemos que el método es estable y consistente enton-
ces
1
k Eh k 6 k A −
h kkE h k 6 C kE h k −−→ 0, h →0
kEh k2 = O(h2 ).
k Ah k2 = ρ( Ah ) = máx |λk |,
16 k 6 m
1
siendo λk , k = 1, . . . , m, los autovalores de Ah . La matriz A−
h también es simétrica, y los autova-
1
lores de A−h son los inversos de los autovalores de A h , por lo que
−1
1 1
k A−
h k2 = ρ( A−
h ) = máx |λ−1 | = mı́n |λk |
16 k 6 m k 16 k 6 m
(la expresión tiene perfecto sentido ya que, al ser Ah invertible, 0 no puede ser autovalor). Puede
demostrarse que los autovalores de Ah tienen la forma
2 kπh
λk = 2 1 − cos , k = 1, . . . , m.
h b−a
Por tanto, el autovalor de menor magnitud será λ1 , que verifica
2 2 π 2 h2 π 4 h4 π2
πh 6
λ1 = 2 1 − cos = 2 − + O ( h ) = + O ( h2 ).
h b−a h 2( b − a )2 24(b − a)4 ( b − a )2
1 −1 está acotada por una constante no nula cuando h → 0. Esto
En consecuencia, k A−
h k2 = | λ1 |
prueba que el método es estable en la norma-2.
Recapitulando, hemos demostrado el siguiente teorema.
Teorema 3.1. El método de diferencias finitas (3.8) para el problema (3.7) es estable y consistente en la
norma-2; en consecuencia, es convergente en dicha norma. Además, el método tiene orden dos.
La implementación en Python del método de diferencias finitas para un problema con condi-
ciones de tipo Dirichlet puede encontrarse en la sección 3.10.1.
2 Para problemas lineales, el teorema de equivalencia de Lax-Richtmyer establece que un método de diferencias finitas
donde se combinan una condición de tipo Neumann, u0 ( a) = γ, con otra de tipo Dirichlet, u(b) = β.
Observación 3.12. El problema de contorno con dos condiciones de tipo Neumann es, en gene-
ral, un problema mal planteado. Esto significa que, dependiendo de los datos, puede haber una,
ninguna o infinitas soluciones. Para verlo, consideremos un problema de la forma
00
−u ( x ) = f ( x ), x ∈ ( a, b),
u0 ( a) = γ,
0
u (b) = δ.
Integrando ambos miembros de la ecuación diferencial entre a y b, resulta que
Z b
−u0 (b) + u0 ( a) = f ( x ) dx,
a
de donde, imponiendo las condiciones de contorno, deducimos una condición necesaria para que
el problema tenga solución:
Z b
γ−δ = f ( x ) dx.
a
De esta forma, si tal condición no se verifica entonces el problema no tiene solución. Pero aun en
el caso en que se verifique la condición, el problema puede tener infinitas soluciones. En efecto,
si tomamos γ = δ = 0 y f ≡ 0, es fácil comprobar que cualquier función de la forma u( x ) ≡ C,
C ∈ R, es solución del problema de contorno.
Es obvio que la condición u(b) = β puede discretizarse como um+1 = β. Pero, ¿cómo podemos
incluir la condición u0 ( a) = γ a la hora de escribir el método de diferencias finitas? Vamos a
mostrar diferentes maneras de hacerlo.
Una primera idea consiste en aproximar u0 ( a) mediante una fórmula descentrada, incluyendo
a u0 como nueva incógnita:
u1 − u0 u −u
=γ⇒ 0 2 1 =− .
γ
h h h
De esta forma, el sistema a considerar puede escribirse como
γ
1 −1 u0 −
h
−1 2 −1 u1 f1
1 2 1 u2 f2
1 − −
.. ,
. . . .. =
2
h
. . . . . . . .
− 1 2 − 1 u m −1 f m −1
−1 2 um β
fm + 2
h
3.5 Otras condiciones de contorno 141
junto con um+1 = β. El inconveniente de esta discretización es que el error global se reduce a
primer orden, debido a que
u ( x1 ) − u ( x0 ) h h
ε0 = γ − = γ − u0 ( x0 ) + u00 ( x0 ) + O(h2 ) = − u00 ( x0 ) + O(h2 ) = O(h).
h 2 2
Otra idea sería utilizar una aproximación centrada de segundo orden para u0 ( a). Una forma de
realizar esto consiste en introducir una incógnita adicional u−1 ≈ u( a − h) y una nueva ecuación,
de la forma (3.8), para j = 0; esto es, consideramos las dos ecuaciones siguientes:
u1 − 2u0 + u−1
− = f 0 (ecuación para j = 0),
h2
u 1 − u −1
= γ (aproximación centrada de u0 ( a)).
2h
Combinando ambas expresiones, es posible eliminar la incógnita u−1 :
2u1 − 2u0 2γ
− = f0 − .
h2 h
Tenemos así un sistema de m + 1 ecuaciones con m + 1 incógnitas {u0 , u1 , . . . , um }, junto con
um+1 = β. En este caso, el sistema a resolver adopta la siguiente forma:
2γ
2 −2 u0 f −
0
h
−1 2 −1 u1 f 1
1 2 1 u2 f 2
1 − −
.
. . . .. = ..
2
h
. . . . . . .
.
− 1 2 − 1 u m −1 f m −1
−1 2 um
β
fm + 2
h
Puede probarse que esta técnica, denominada del nodo fantasma, proporciona una aproximación
de segundo orden a la solución. Esta es el método que utiizaremos en la implementación práctica
del método.
Una variante de la idea anterior consiste en utilizar una aproximación de segundo orden de
u0 ( a), pero esta vez descentrada. Una fórmula análoga a (3.5) nos proporciona una aproximación
de este tipo:
−u2 + 4u1 − 3u0
= γ,
2h
donde de nuevo u0 es una incógnita. En este caso, el sistema quedaría de la forma
2γ
3 −4 1
u0
−
h
−1 2 −1 u1
f1
−1 2 −1 u2 f2
1 ,
.. .. .. .. = ..
h2 . . . . .
−1 2 −1 u m −1 f m −1
−1 2 um
β
fm + 2
h
142 Métodos numéricos para problemas de contorno
3.6. Generalizaciones
Consideremos el problema de contorno con condiciones de tipo Dirichlet para una ecuación
lineal de segundo orden, esto es,
00 0
−u ( x ) + p( x )u ( x ) + q( x )u( x ) = f ( x ), x ∈ ( a, b),
u( a) = α,
u(b) = β,
2 + h2 q1 −1 + hp1 /2
−1 − hp2 /2 2 + h2 q2 −1 + hp2 /2
−1 − hp3 /2 2 + h2 q3 −1 + hp3 /2
1
A= 2 .. .. .. .
h . . .
2 + h 2 q m −1
−1 − hpm−1 /2 −1 + hpm−1 /2
−1 − hpm /2 2 + h2 q m
Observación 3.15. En el caso de que una condición de contorno sea de tipo Neumann, podemos
aplicar alguna de las técnicas introducidas en la sección 3.5.
A continuación vamos a estudiar cómo aplicar el método de diferencias finitas a un problema
no lineal de la forma
−u00 ( x ) = F ( x, u( x ), u0 ( x )), x ∈ ( a, b),
u( a) = α,
u(b) = β,
G (U ) = 0, (3.12)
∂G
G (Uk ) + (U )(Uk+1 − Uk ) = 0,
∂U k
3 El método de Newton se deduce a partir del desarrollo de Taylor
∂G
G (Uk+1 ) = G (Uk ) + (U )(Uk+1 − Uk ) + · · ·
∂U k
despreciando los términos de alto orden y haciendo G (Uk+1 ) = 0.
144 Métodos numéricos para problemas de contorno
∂G
(U )V = − G (Uk ), (3.13)
∂U k k
se tiene que
Uk+1 = Uk + Vk .
u i +1 − u i −1 u i +1 − u i −1
∂F ∂F
ai = x ,u , , bi = xi , ui , ,
∂y i i 2h ∂z 2h
2
2
− ai si j = i,
h
− 1 + bi
si j = i − 1,
∂Gi
(U ) = h2 2h
∂u j 1 b
− 2− i
si j = i + 1,
h 2h
0 en otro caso.
2 − h2 a1 −1 − hb1 /2
−1 + hb2 /2 2 − h2 a2 −1 − hb2 /2
−1 + hb3 /2 2 − h2 a3 −1 − hb3 /2
∂G 1
(U ) = 2 .. .. .. .
∂U h . . .
2 − h 2 a m −1
−1 + hbm−1 /2 −1 − hbm−1 /2
−1 + hbm /2 2 − h2 a m
En resumen, en cada iteración del método de Newton es necesario resolver el sistema tridiagonal
(3.13), cuya forma es análoga a la del sistema que aparecía al aplicar el método de diferencias
finitas a un problema lineal.
Observación 3.16. De nuevo, el tratamiento de condiciones de tipo Neumann sigue las directrices
explicadas en la sección 3.5.
Observación 3.17. Si las condiciones de contorno no son homogéneas (esto es, u(0) = α y u(1) =
β), basta con hacer el cambio w( x ) = u( x ) − α(1 − x ) − βx para pasar a un problema equivalente
con condiciones de contorno w(0) = 0, w(1) = 0. Asimismo, el cambio w( x ) = u( a + (b − a) x )
permite considerar el problema en un intervalo arbitrario [ a, b]. Por otra parte, el añadir un
término de la forma p( x )u0 ( x ) a la ecuación diferencial no introduce dificultades adicionales
en la deducción de la formulación variacional, por lo que no lo consideraremos en aras de la
simplicidad.
Observación 3.18. Aunque solo consideraremos problemas con condiciones de contorno de tipo
Dirichlet, es posible adaptar las técnicas que estudiaremos al caso de problemas con condiciones
de tipo Neumann o de tipo mixto.
Sea u ∈ C2 ([0, 1]) una solución del problema (3.14): diremos que u es una solución fuerte.
Multipliquemos la ecuación diferencial por una función arbitraria v ∈ C1 ([0, 1]), e integremos
sobre todo el intervalo [0, 1]:
Z 1 Z 1 Z 1
00
− u ( x )v( x ) dx + q( x )u( x )v( x ) dx = f ( x )v( x ) dx.
0 0 0
Supongamos además que v(0) = v(1) = 0 e integremos por partes la primera integral:
Z 1 Z 1 1 Z 1
u00 ( x )v( x ) dx = u0 ( x )v0 ( x ) dx − u0 ( x )v( x ) 0 = u0 ( x )v0 ( x ) dx.
−
0 0 0
todas las funciones v continuas en [0, 1] que verifican las condiciones de contorno v(0) = v(1) = 0
y son de clase C1 a trozos en [0, 1] . El espacio V admite la norma
Z 1
1/2
0 2 2
k v kV = |v ( x )| + |v( x )| ) dx , v ∈ V.
0
y
Z 1
f (v) = f ( x )v( x ) dx, v ∈ V.
0
Demostración. Sea v ∈ V y sean x1 < x2 < · · · < xd los puntos de discontinuidad de v0 ; definamos
x0 = 0 y xd+1 = 1. Entonces
Z 1 d Z x i +1 d Z x
i +1 x i +1
00
u ( x )v( x ) dx = ∑ 00
u ( x )v( x ) dx = ∑ 0 0 0
u ( x )v ( x ) dx + u ( x )v( x )
− xi
0 i =0 x i i =0 xi
Z 1
=− u0 ( x )v0 ( x ) dx,
0
donde
1
J (v) = a(v, v) − f (v).
2
3.7 Aproximación variacional de problemas de contorno 147
Demostración. En primer lugar, vamos a probar que existe una constante α > 0 tal que
de donde
Z 1 Z 1 Z 1
kvk2V = |v0 ( x )|2 + |v( x )|2 dx 6 2 |v0 ( x )|2 dx 6 2 |v0 ( x )|2 + q( x )|v( x )|2 dx
0 0 0
= 2a(v, v).
1
J (u + v) − J (u) = a(u, v) − f (v) + a(v, v), ∀ u, v ∈ V.
2
Si u ∈ V es solución de (3.16) entonces
1 α
J (u + v) − J (u) = a(v, v) > kvk2V , ∀ v ∈ V.
2 2
Esto prueba que J (u) = ı́nfv∈V J (v). Para ver el recíproco, fijemos v ∈ V y sea θ ∈ R; entonces
θ2
0 6 J (u + θv) − J (u) = θ ( a(u, v) − f (v)) + a(v, v).
2
Tomando θ > 0, tenemos que
θ
a(u, v) + a(v, v) > f (v)
2
y, haciendo θ → 0+ , deducimos que a(u, v) > f (v). Razonando de forma análoga para θ < 0,
deducimos que a(u, v) 6 f (v). En consecuencia, a(u, v) = f (v)Êpara cada v ∈ V.
La principal dificultad que nos encontramos al intentar resolver la ecuación variacional (3.15)
está en el hecho de que el espacio vectorial V es un espacio de dimensión infinita; esto hace que no
sea factible una implementación directa en el ordenador. El método de aproximación variacional o
método de Galerkin se basa en considerar un subespacio Vh de V de dimensión finita, y a continua-
ción aproximar la solución u de (3.15) por la solución uh ∈ Vh de la siguiente ecuación variacional
discreta:
Z 1 Z 1 Z 1
u0h ( x )v0h ( x ) dx + q( x )uh ( x )vh ( x ) dx = f ( x )vh ( x ) dx, ∀ vh ∈ Vh , (3.17)
0 0 0
148 Métodos numéricos para problemas de contorno
o, equivalentemente, de
a ( u h , v h ) = f ( v h ), ∀ vh ∈ Vh . (3.18)
Diremos que Vh es el espacio de funciones test en este caso. La notación Vh hace referencia a que el
espacio va a depender del paso de malla elegido: mientras menor sea éste, mejor queremos que
sea la aproximación de la solución.
Teorema 3.4. Supongamos q( x ) > 0, ∀ x ∈ [0, 1]. Se verifican las siguientes proposiciones:
3. Existe una constante C, independiente de la solución u de (3.16) y del subespacio Vh , tal que
ku − uh kV 6 C ı́nf ku − vh kV .
vh ∈Vh
Demostración.
1. El problema (3.18) equivale a la resolución de un sistema lineal4 con matriz cuadrada, por
lo que la existencia y la unicidad de solución son propiedades equivalentes. Se tiene que
Ahora bien,
)
a(u, v) = f (v), ∀ v ∈ V
⇒ a(u − uh , vh ) = 0, ∀ vh ∈ Vh .
a(uh , vh ) = f (vh ), ∀ vh ∈ Vh ⊂ V
4 Más adelante, después de la demostración del teorema, se detalla la forma de dicho sistema lineal.
3.7 Aproximación variacional de problemas de contorno 149
Por tanto,
Al ser los coeficientes {v1 , . . . , vm } arbitrarios, deducimos que la ecuación (3.18) es equivalente al
conjunto de m ecuaciones
a(uh , ϕi ) = f ( ϕi ), i = 1, . . . , m. (3.19)
Es decir, basta con que la ecuación (3.18) se verifique para los elementos ϕi de la base.
Escribamos ahora la solución buscada uh en la base { ϕ1 , . . . , ϕm }:
m
uh ( x ) = ∑ u j ϕ j ( x ), u1 , . . . , um ∈ R.
j =1
Si definimos
aij = a( ϕi , ϕ j ), f i = f ( ϕ i ), 1 6 i, j 6 m,
podemos escribir la expresión anterior de la siguiente forma:
m
∑ aij u j = fi , i = 1, . . . , m.
j =1
Hemos obtenido así un sistema lineal de m ecuaciones y m incógnitas (a saber, los coeficientes ui
de la solución uh ). En forma matricial, dicho sistema se escribe
AU = F,
150 Métodos numéricos para problemas de contorno
donde
a11 a12 ··· a1m u1 f1
a21 a22 ··· a2m u2 f2
A= . .. .. .. , U = . , F = . .
.. . . . .. ..
am1 am2 · · · amm um fm
La solución del sistema nos permitirá construir la solución del problema discreto (3.17). 5
Observación 3.19. La solución U = (u1 , . . . , um )t del sistema lineal permite construir la solución
numérica de la forma
u h ( x ) = u1 ϕ1 ( x ) + u2 ϕ2 ( x ) + · · · + u m ϕ m ( x ).
En este caso, dicha solución numérica es una función y no un conjunto de valores aproximados,
como sucedía en el método de diferencias finitas. Podemos por tanto usar uh ( x ) para aproximar
la solución del problema de contorno en cualquier punto del intervalo [0, 1].
En resumen, hemos visto cómo la solución del problema variacional (3.15) puede aproximarse
mediante la solución de un problema discreto (3.17), cuyo cálculo se reduce a la resolución de un
sistema lineal. En la siguiente sección veremos cómo elegir los elementos ϕi de la base de forma
adecuada para que el sistema pueda resolverse de forma eficiente.
φi
xi-1 xi xi+1
cuya gráfica se representa en la figura 3.3. Cada ϕi ( x ) es una función continua de clase C1 a
trozos, cuya derivada es
1
si xi−1 < x < xi ,
h
ϕi0 ( x ) = − 1 si x < x < x ,
i i +1
h
0 en otro caso.
Nótese que las funciones ϕi ( x ) y ϕi0 ( x ) se anulan fuera del intervalo ( xi−1 , xi+1 ).
Puede demostrarse que la familia de funciones { ϕ1 , . . . , ϕm } constituyen una base del subes-
pacio vectorial Vh formado por todas las funciones continuas en [0, 1], lineales en cada intervalo
[ xi , xi+1 ], 0 6 i 6 m, y que se anulan en los extremos x = 0 y x = 1. En particular, la dimensión
de Vh es m.
Observación 3.21. Para condiciones de contorno de tipo Neumann es necesario añadir dos fun-
ciones de base adicionales, ϕ0 ( x ) y ϕ N ( x ), para imponer las condiciones de contorno en cada
uno de los extremos del intervalo. Asimismo, el espacio de funciones test Vh debe modificarse de
manera adecuada.
A continuación vamos a calcular los coeficientes aij de la matriz A, así como los términos
independientes f i . Por simplicidad, a partir de ahora supondremos que las funciones q( x ) y f ( x )
son constantes: q( x ) ≡ q, f ( x ) ≡ f .
φi-1 φi φi φi+1 φi φj
x − x i −1 2 x i +1 − x 2
Z 1 Z x Z x Z x
2h
i +1 i i +1
2
ϕi ( x ) ϕi ( x ) dx = ( ϕi ( x )) dx = dx + dx = ,
0 x i −1 x i −1 h xi h 3
teniendo en cuenta que xi−1 = xi − h y xi+1 = xi + h. De forma similar, se obtienen las igualdades
Z 1 Z x
i x −x x − x i −1 h
i
ϕi ( x ) ϕi−1 ( x ) dx = · dx = ,
0 x i −1 h h 6
y
Z 1 Z x
i +1 x − x x i +1 − x h
i
ϕi ( x ) ϕi+1 ( x ) dx = · dx = .
0 xi h h 6
Razonando de manera análoga, se deducen las siguientes expresiones:
Z 1 Z 1 Z 1
2 1 1
ϕi0 ( x ) ϕi0 ( x ) dx = , ϕi0 ( x ) ϕi0−1 ( x ) dx = − , ϕi0 ( x ) ϕi0+1 ( x ) dx = − .
0 h 0 h 0 h
Finalmente, obtenemos las fórmulas para los coeficientes aij :
2 2
+ qh si j = i,
h 3
aij = − 1 + 1 qh si j = i − 1 o j = i + 1,
h 6
0 en otro caso.
f i = f h, i = 1, . . . , m.
3.8 Método de elementos finitos 153
A partir de los resultados obtenidos, podemos escribir el sistema lineal que define el método
de elementos finitos como sigue:
2 + 23 qh2 −1 + 16 qh2 0 0 0 u1 f
···
00
−u ( x ) + u( x ) = 2, x ∈ [0, 1],
u(0) = 0,
u(1) = 0,
2
u( x ) = 2 − ( e x + e 1− x ).
1+e
Para construir la formulación variacional del problema, multiplicamos la ecuación dife-rencial
por una función arbitraria v( x ), que sea de clase C1 a trozos en [0, 1] y verifique v(0) = v(1) = 0,
y a continuación integramos en [0, 1]:
Z 1 Z 1 Z 1
− u00 ( x )v( x ) dx + u( x )v( x ) dx = 2v( x ) dx.
0 0 0
Sea V el espacio vectorial formado por todas las funciones continuas v( x ) de clase C1 a trozos
en [0, 1] que verifican v(0) = v(1) = 0. La formulación variacional del problema es entonces la
siguiente: hallar u ∈ V tal que
Z 1 Z 1 Z 1
u0 ( x )v0 ( x ) dx + u( x )v( x ) dx = 2v( x ) dx, ∀ v ∈ V.
0 0 0
El siguiente paso consiste en discretizar la ecuación variacional. Para ello, consideremos una
partición del intervalo [0, 1] con m = 3; de esta forma, el paso de malla es h = 14 y los nodos
vienen dados por x0 = 0, x1 = 14 , x2 = 12 , x3 = 34 y x4 = 1. Las correspondientes funciones de
base { ϕ1 ( x ), ϕ2 ( x ), ϕ3 ( x )} se representan en la figura 3.5.
φ1 φ2 φ3
mientras que
Z 1 Z 1 Z 1
u0h ( x ) ϕ30 ( x ) dx + uh ( x ) ϕ3 ( x ) dx = 2ϕ3 ( x ) dx
0 0 0
resulta al elegir v1 = v2 = 0 y v3 = 1.
Como la solución buscada uh es un elemento de Vh , podemos escribirla en la forma
u h ( x ) = u1 ϕ1 ( x ) + u2 ϕ2 ( x ) + u3 ϕ3 ( x ),
Si definimos
Z 1 Z 1 Z 1
aij = ϕi0 ( x ) ϕ0j ( x ) dx + ϕi ( x ) ϕ j ( x ) dx, fi = 2ϕi ( x ) dx,
0 0 0
y
a31 u1 + a32 u2 + a33 u3 = f 3 .
Los coeficientes aij y f i se determinan explícitamente:
49 95
a11 = a22 = a33 = 24 , a12 = a21 = a23 = a32 = − 96 , a13 = a31 = 0, f 1 = f 2 = f 3 = 18 .
En resumen, las incógnitas {u1 , u2 , u3 } se obtienen como solución de un sistema lineal de tres
ecuaciones, que en forma matricial se escribe como
49
− 95
1
24 96 0 u1 8
95 49 95 u = 1 .
−
96 24 − 96 2 8
0 − 95
96
49
24
u 3
1
8
156 Métodos numéricos para problemas de contorno
El error Ei = |u( x j ) − u j | cometido en cada nodo interior puede determinarse en este caso, ya
que conocemos la solución exacta:
o con otros tipos de condiciones de contorno. En tal caso, la solución numérica se determina
a partir de la solución de un sistema de ecuaciones no lineales, que debe resolverse mediante
técnicas apropiadas. El estudio de tales cuestiones queda fuera del alcance de este curso, por lo
que no insistiremos más en ellas.
Los problemas que surgen en las aplicaciones prácticas suelen modelizarse mediante ecua-
ciones en derivadas parciales, tanto estacionarias como de evolución (dependientes del tiempo), en
dominios bi- o tridimensionales. Tanto el método de diferencias finitas como el de elementos
finitos pueden extenderse para resolver problemas de esta índole, siendo el método de elementos
3.9 El método de diferencias finitas en dimensión superior 157
6 Nótese ∂2 u ∂2 u
que si en la ecuación de Poisson sustituimos las derivadas ∂x2
y ∂y2
por x2 e y2 respectivamente, y f por −1,
se obtiene la elipse x2 + y2 = 1 en el plano x-y.
158 Métodos numéricos para problemas de contorno
y denotemos
hx = máx ( xk+1 − xk ), hy = máx (yk+1 − yk ).
k =0,...,m x k =0,...,my
y
hx
y4=d
y3
hy
y2
y1
y0=c
x
x0=a x1 x2 x3 x4 x5=b
Figura 3.7: Mallado del rectángulo ( a, b) × (c, d) con m x = 4 y my = 3. Los nodos interiores se han repre-
sentado con círculos y los nodos frontera con rectángulos.
Si denotamos por uij a la aproximación del valor exacto u( xi , y j ), podemos aproximar las
derivadas parciales mediante cocientes incrementales análogos a los del caso unidimensional:
Evaluando la ecuación −∆u = f en cada nodo interior y reemplazando las derivadas parciales
por las fórmulas anteriores, obtenemos el siguiente conjunto de ecuaciones:
uij = gij
Para simplificar aún más las expresiones, supondremos que h x = hy = h. En tal caso, las
ecuaciones (3.20) quedan de la siguiente forma:
1
− (ui+1,j + ui,j+1 − 4uij + ui,j−1 + ui−1,j ) = f ij , i = 1, . . . , m x , j = 1, . . . , my . (3.21)
h2
(i, j+1)
(i, j-1)
Figura 3.8: Plantilla del esquema de cinco puntos para el operador de Laplace.
Al escribir el sistema (3.21) en forma matricial se obtiene una matriz hueca (en inglés, sparse),
esto es, la mayoría de sus elementos son cero. Sin embargo, en dos dimensiones la estructura
de dicha matriz no es tan simple como en el caso unidimensional. De hecho, dicha estructura
dependerá del orden en que se escriban las incógnitas. Una elección habitual consiste en seguir
el orden lexicográfico, según el cual los nodos se numeran de izquierda a derecha y de arriba
hacia abajo (figura 3.9). De esta forma, el sistema puede representarse mediante una matriz A ∈
Mm (R) con estructura tridiagonal por bloques:
T D
D T D
D T D
1
A= 2 .. .. .. , (3.22)
h . . .
D T D
D T
160 Métodos numéricos para problemas de contorno
4 −1
−1 4 −1
−1 4 −1
T= .. .. .. ∈ M m x (R).
. . .
−1 4 −1
−1 4
No es difícil probar que la matrix A es simétrica y definida positiva, por lo que el sistema de
diferencias finitas posee una única solución uh .
1 2 3 4
5 6 7 8
9 10 11 12
Puede demostrarse que si u ∈ C4 (Ω), existe entonces una constante C > 0 tal que
y g( x, 0) = g( x, 1) = sen(2πx ), para 0 < x, y < 1. La solución exacta del problema de Poisson es,
en este caso,
u( x, y) = sen(2πx ) cos(2πy).
Figura 3.10: Solución de la ecuación de Poisson: desplazamiento transversal de una membrana elástica.
m x = my h ku − uh k∞ orden
25 0.04 5.27 · 10−3 −
50 0.02 1.37 · 10−3 1.94
100 0.01 3.50 · 10−4 1.97
200 0.005 8.85 · 10−5 1.98
400 0.0025 2.22 · 10−5 1.99
En la tabla 3.1 se muestran los errores en norma infinito obtenidos con distintas particiones,
así como las correspondientes estimaciones del orden, que muestran que el método converge con
orden dos.
Observación 3.27. Otros métodos más sofisticados y muy utilizados en la práctica son el método
de elementos finitos y el método de volúmenes finitos.
162 Métodos numéricos para problemas de contorno
donde
u1 ( t ) f 1 (t) u0 ( x1 )
um (t) f m (t) u0 ( x m )
7 Si ∂2 u
en la ecuación del calor sustituimos las derivadas ∂u
∂t y ∂x2
por t y x2 respectivamente, y f por 1, se obtiene la
parábola t − µx2 = 1 en el plano x-t.
3.9 El método de diferencias finitas en dimensión superior 163
U k +1 − U k
= µA(θU k+1 + (1 − θ )U k ) + θF k+1 + (1 − θ ) F k , k = 0, 1, . . . , (3.24)
∆t
con U 0 = U0 . De forma equivalente, podemos escribir
( I − µθ∆tA)U k+1 = (1 + µ∆t(1 − θ ) A)U k + G k+1 ,
siendo G k+1 = ∆t(θF k+1 + (1 − θ ) F k ). En el caso particular θ = 0 (método de Euler) podemos
obtener U k+1 de forma explícita:
U k+1 = (1 + µ∆tA)U k + ∆tF k .
Para θ 6= 0, hay que resolver un sistema lineal con matriz constante I − µθ∆tA en cada iteración
de tiempo.
Con respecto a la precisión del método, puede probarse que el error de discretización local es
O(∆t + h2 ) si θ 6= 1/2, mientras que es O(∆t2 + h2 ) para θ = 1/2.
Las ideas y resultados anteriores pueden extenderse al caso de la ecuación del calor en un
dominio bidimensional Ω ⊂ R2 ,
∂u
( x, y; t) − µ∆u( x, y; t) = f ( x, y; t), ( x, y) ∈ Ω, t > 0.
∂t
El esquema resultante es análogo a (3.24), donde la matriz A se sustituye por la matriz (3.22).
Ejemplo 3.5. Consideremos una barra de acero (con coeficiente µ = 0.1515) de 20 cm de longitud.
Partiendo de una distribución inicial de temperaturas dada por 10 − x si x ∈ [0, 10], y 200 − 10x
en caso contrario, queremos estudiar la evolución de la temperatura u durante un periodo de 200
segundos. Supondremos que los extremos de la barra están aislados, esto es, u(0, t) = u(20, t) =
0. En la figura 3.11 se muestran las soluciones obtenidas con el método de Euler implícito (θ =
1), desde t = 0 hasta t = 200 con pasos de 25 segundos, tomando m = 50. Asimismo, se ha
representado el resultado obtenido en t = 25 al aplicar el método de Euler (θ = 0). Como puede
observarse, en este caso el método es inestable.
Figura 3.11: Evolución de la temperatura en una barra de acero. Izquierda: Método de Euler implícito en
tiempo para t = 0 hasta t = 200. Derecha: Método de Euler en tiempo para t = 25.
u( a, t) = u(b, t) = 0, t > 0.
Observación 3.30. Cuando f ≡ 0, la solución general de (3.25) viene dada por la función de ondas
8 Al ∂2 u ∂2 u
sustituir en la ecuación de ondas las derivadas ∂t2
y ∂x2
por t2 y x2 respectivamente, y f por 1, se obtiene la
hipérbola t2 − cx2 = 1 en el plano x-t.
3.9 El método de diferencias finitas en dimensión superior 165
A = PΛP−1 ,
Aui = λi ui , i = 1, . . . , p.
anterior es simplemente
wi ( x, t) = wi ( x − λi t, 0), i = 1, . . . , p.
166 Métodos numéricos para problemas de contorno
Por tanto, la solución U = PW del sistema (3.26) para F ≡ 0 viene dada por
p
U ( x, t) = ∑ wi (x − λi t, 0)ui .
i =1
Se define la i-ésima curva característica del sistema (3.26) como la curva ( xi (t), t) en el plano x-t
que satisface xi0 (t) = λi . Es obvio que wi es constante a lo largo de dicha curva; entonces, dado
un punto ( x̄, t̄), el valor de U ( x̄, t̄) tan sólo depende del dato inicial en x̄ − λi t̄. Por esta razón,
el conjunto de los p puntos que forman el pie de las curvas características que parten del punto
( x̄, t̄),
D (t̄, x̄ ) = { x ∈ R : x = x̄ − λi t̄, i = 1, . . . , p},
se denomina dominio de dependencia de la solución U ( x̄, t̄).
Si el sistema (3.26) está definido en un intervalo acotado ( a, b), el punto de entrada para cada
variable característica wi lo determina el signo de λi . Así pues, el número de autovalores positivos
establece el número de condiciones de contorno que se pueden asignar en x = a, mientras que
en x = b hay que asignar tantas condiciones de contorno como autovalores negativos haya.
Volviendo
√ al sistema (3.26)-(3.27), resulta que la matriz A es diagonalizable con autovalores
reales ± c (estas cantidades representan las velocidades de propagación de la onda), por lo que
el sistema es hiperbólico. Como hay un autovalor positivo y otro negativo, es necesario prescribir
una condición de contorno en cada extremo del intervalo ( a, b).
A continuación, vamos a discretizar el sistema (3.26). Por simplicidad, consideraremos tan
sólo el caso en que ( a, b) = R, F ≡ 0 y p = 1. Es decir, vamos a estudiar la ecuación de advección o
transporte lineal escalar
∂U ∂U
+α = 0, x ∈ R, t > 0, (3.28)
∂t ∂x
donde α es una constante. Si la condición inicial es U ( x, 0) = u0 ( x ), entonces la solución de (3.28)
es simplemente
U ( x, t) = u0 ( x − αt). (3.29)
Geométricamente, U ( x, t) es una onda cuya forma viene dada por u0 ( x ), que se desplaza hacia
la derecha si α > 0, y hacia la izquierda si α < 0.
Observación 3.31. La solución U de la ecuación de advección (3.28) puede interpretarse como
la concentración o densidad de una sustancia que se desplaza a velocidad constante α en un
fluido.
Dados un paso de malla espacial h (que supondremos uniforme) y un paso de tiempo ∆t,
consideramos los puntos
x j = jh, j ∈ Z; tn = n∆t, n ∈ N,
y denotamos por Ujn a la aproximación del valor U ( x j , tn ).
Una primera idea para discretizar (3.28) consiste en adaptar lo que hicimos en la sección 3.9.2
para la ecuación del calor. Se obtiene así un método de la forma
Ujn+1 − Ujn α
=− (U n − Ujn−1 ),
∆t 2h k+1
3.9 El método de diferencias finitas en dimensión superior 167
1 n α∆t n
Ujn+1 = (Uj+1 + Ujn−1 ) − (Uj+1 − Ujn−1 ),
2 2h
que es de primer orden.
Vamos a construir ahora un método de segundo orden en tiempo y en espacio. Para ello,
consideremos el siguiente desarrollo de Taylor en tiempo:
∂U ∆t2 ∂2 U
U ( x, t + ∆t) = U ( x, t) + ∆t ( x, t) + ( x, t) + · · ·
∂t 2 ∂t2
Teniendo en cuenta que
∂U ∂U ∂2 U ∂2 U ∂2 U
= −α ⇒ 2 = −α = α2 2 ,
∂t ∂x ∂t ∂t∂x ∂x
podemos escribir
∂U α2 ∆t2 ∂2 U
U ( x, t + ∆t) = U ( x, t) − α∆t ( x, t) + ( x, t) + · · ·
∂x 2 ∂x2
(de esta forma hemos sustituido las derivadas temporales por derivadas espaciales). Usando
ahora aproximaciones centradas y despreciando los términos de orden superior, obtenemos el
denominado método de Lax-Wendroff :
λ λ2 2 n
Ujn+1 = Ujn − α(Ujn+1 − Ujn−1 ) + α (Uj+1 − 2Ujn + Ujn−1 ),
2 2
donde se ha definido
∆t
λ= .
h
Los métodos anteriores son métodos centrados, en el sentido de que se basan en aproximacio-
nes centradas de las derivadas. Teniendo en cuenta que la solución de la ecuación (3.25) depende
del signo del coeficiente α (esto es, de la dirección de propagación de la onda), vamos a construir
a continuación un método en el que se tenga en cuenta este hecho. Para ello, consideremos las
aproximaciones descentradas
∂U 1 ∂U 1
( x , tn ) ≈ (Ujn − Ujn−1 ), ( x , tn ) ≈ (Ujn+1 − Ujn ).
∂x j h ∂x j h
Aplicando el método de Euler en tiempo y usando estas aproximaciones, podemos discretizar
(3.25) como sigue:
α∆t n
Ujn+1 = Ujn − (Uj − Ujn−1 ), (3.30)
h
168 Métodos numéricos para problemas de contorno
o bien
α∆t n
Ujn+1 = Ujn − (Uj+1 − Ujn ). (3.31)
h
Teniendo en cuenta la forma de la solución (3.29), es fácil ver que
U ( x j , tn+1 ) = U ( x j − α∆t, tn ).
De esta forma, la solución en el punto x j en el instante tn+1 viene determinada por el valor de la
solución en el tiempo tn a la izquierda de x j si α > 0, y a la derecha si α < 0. Esto sugiere elegir
el método (3.30) si α > 0 y (3.31) si α < 0. Ambas opciones pueden combinarse en una única
expresión, dando lugar así al denominado método upwind:
λ λ
Ujn+1 = Ujn − α(Ujn+1 − Ujn−1 ) + |α|(Ujn+1 − 2Ujn + Ujn−1 ),
2 2
que es de primer orden en tiempo y en espacio.
Los métodos anteriores pueden extenderse de forma natural para resolver la ecuación general
(3.26):
Método de Lax-Friedrichs:
1 n λ
Ujn+1 = (U + Ujn−1 ) − A(Ujn+1 − Ujn−1 ).
2 j +1 2
Método de Lax-Wendroff:
λ λ2 2 n
Ujn+1 = Ujn − A(Ujn+1 − Ujn−1 ) + A (Uj+1 − 2Ujn + Ujn−1 ).
2 2
Método upwind:
λ λ
Ujn+1 = Ujn − A(Ujn+1 − Ujn−1 ) + | A|(Ujn+1 − 2Ujn + Ujn−1 ),
2 2
donde se define | A| = P|Λ| P−1 , siendo |Λ| la matriz diagonal formada por los módulos de
los autovalores de A.
Observación 3.32. Para garantizar la estabilidad de los esquemas numéricos introducidos ante-
riormente, es necesario imponer una condición de tipo CFL (de Courant, Friedrichs y Lewy):
h
∆t < ,
ρ( A)
Figura 3.12: Izquierda: ejemplo (3.6); derecha: ejemplo (3.7). De arriba hacia abajo: métodos de Lax-
Friedrichs, Lax-Wendroff y upwind.
170 Métodos numéricos para problemas de contorno
Ejemplo 3.6. Consideremos la ecuación de advección lineal (3.28) con α = 1 y condición inicial
2
u0 ( x ) = e−8x , para x ∈ (−2, 2). En la columna de la izquierda de la figura 3.12 se muestran los
resultados obtenidos con los métodos de Lax-Friedrichs, Lax-Wendroff y upwind para t = 1, con
∆x = 0.04 y ∆t = 0.008. Como puede observarse, el método de Lax-Wendroff (que, recordemos,
es de segundo orden en tiempo y en espacio) proporciona la aproximación más precisa, seguido
del método upwind.
Ejemplo 3.7. Tomemos ahora como condición inicial una onda cuadrada, de la forma
0 si x 6 −0.3,
u0 ( x ) = 1 si − 0.3 6 x 6 0.3,
0 si x > 0.3,
donde x ∈ (−2, 2). En la columna derecha de la figura 3.12 se han representado las soluciones
obtenidas para t = 1, con ∆x = 0.04 y ∆t = 0.008. En este caso, el método de Lax-Friedrichs pro-
porciona un resultado pobre (nótese la agrupación por pares de puntos consecutivos). El método
de Lax-Wendroff tampoco aproxima la solución correctamente: cerca de las discontinuidades apa-
recen fuertes oscilaciones, debidas a que los desarrollos de Taylor en los que se basaba el método
de Lax-Wendroff dejan de ser válidos cerca de las discontinuidades. Por último, el método upwind
sí captura de forma adecuada la solución.
En el apartado 3.10.4 se ha realizado una implementación en Python del método upwind para
resolver la ecuación de advección lineal.
3.10 Implementación en Python del método de diferencias finitas 171
00 0
−u ( x ) + p( x )u ( x ) + q( x )u( x ) = f ( x ), a < x < b,
u( a) = α,
u(b) = β,
00
−u ( x ) + u( x ) = sen(2πx ), 0 < x < 1,
u(0) = 0,
u(1) = 0,
Primero definimos la solución exacta (líneas 10-12) y las funciones p( x ), q( x ) y f ( x ) (líneas 14-
3.10 Implementación en Python del método de diferencias finitas 173
19) que definen la ecuación diferencial. El programa principal comienza en la línea 21, donde se
definen los datos del problema y el número de particiones deseado. La matriz A del sistema lineal
resultante se define en las líneas 35-43: primero se crea una matriz nula del tamaño adecuado, y
después se definen los elementos de las tres diagonales principales mediante un solo bucle; este
bucle deja sin definir el último elemento de la matriz, que se define en la línea 42. A continuación,
en las líneas 46-47 se modifica el término independiente para incluir las condiciones de contorno.
El sistema lineal resultante (3.9) se resuelve en la línea 49, y en las líneas 52-53 se añaden las
condiciones de contorno en los extremos inicial y final del intervalo. El resto del código se dedica
a calcular el error cometido y a representar gráficamente los resultados obtenidos.
cuya solución exacta es u( x ) = 2x − x2 . Para ello, usaremos la técnica del nodo fantasma intro-
ducida en la sección 3.5 para discretizar la segunda condición de contorno.
# Resolucion del problema de contorno
# -u ' '( x )+ p ( x ) u '( x )+ q ( x ) u ( x ) = f ( x )
# u ( a ) = alpha , u '( b ) = gamma ,
# ( condiciones de tipo Dirichlet y Neumann respectivamente )
# mediante el metodo de diferencias finitas ( con nodo fantasma )
def exacta ( x ):
""" Solucion exacta """
return 2.* x - x **2
def fun ( x ):
""" Funcion f ( x ) """
f = 2.* ones ( len ( x ))
return f
tini = clock ()
174 Métodos numéricos para problemas de contorno
h = (b - a )/ N # paso de malla
x = arange (a , b +h , h ) # discretizacion del intervalo [a , b ]
tfin = clock ()
# Resultados
print ( ' ----- ')
print ( ' Tiempo CPU : ' + str ( tfin - tini ))
print ( ' Error : ' + str ( error ))
print ( ' Paso de malla : ' + str ( h ))
print ( ' ----- ')
que corresponde a la solución del problema de valor inicial linealizado dado por
00
− θ ( t ) = θ ( t ), 0 < t < 2π,
θ (0) = 0.7,
0
θ (0) = 0.5.
# Pendulo no lineal
# ------------------
# Resolucion del problema de contorno
# u ' '( t ) = - sen ( u ( t )) , 0 <t <T ,
# u (0) = alpha , u ( T ) = beta ,
# ( condiciones de tipo Dirichlet )
# mediante el metodo de diferencias finitas .
def cond_ini ( t ):
""" Condicion inicial para el metodo de Newton """
return 0.7* cos ( t )+0.5* sin ( t )
return w
h = (b - a )/ N # paso de malla
t = arange (a , b +h , h ) # discretizacion del intervalo [a , b ]
u = zeros ( N +1) # solucion numerica
# Condiciones de contorno
u [0] = alpha
u [ N ] = beta
# Metodo de Newton
tol = 1. e -6 # tolerancia
error = tol +1. # error entre dos iteraciones sucesivas
maxiter = 10 # numero maximo de iteraciones
k = 0 # contador de iteraciones
λ λ
Ujn+1 = Ujn − α(Ujn+1 − Ujn−1 ) + |α|(Ujn+1 − 2Ujn + Ujn−1 ), (3.32)
2 2
donde λ = ∆t/h es el cociente entre los pasos de malla temporal y espacial, para resolver la
ecuación de advección o transporte
∂U ∂U
+α = 0, x ∈ R, t > 0,
∂t ∂x
donde α es una constante. Recordemos que si la condición inicial es U ( x, 0) = u0 ( x ), entonces la
solución de la ecuación es simplemente
U ( x, t) = u0 ( x − αt).
2
En el código que daremos a continuación tomaremos α = 1, u0 ( x ) = e−8x , tiempo final de
simulación T = 1, h = 0.04 y ∆t = 0.008. Realizaremos la simulación en el intervalo espacial
[−2, 2], lo cual implica que debemos imponer necesariamente algún tipo de condiciones de con-
torno. Para este problema vamos a imponer condiciones transparentes o de salida libre; de esta
forma, podemos pensar que el dominio espacial es todo R y nosotros solo estamos viendo lo
que sucede en el intervalo [−2, 2]. Para ello, si {Uin }iN=0 es la solución numérica en el instante tn ,
consideramos los valores fantasma definidos por
n n n n
U− 1 = U0 , UN +1 = U N .
λ λ
U0n+1 = U0n − α(U1n − U0n ) + |α|(U1n − U0n ),
2 2 (3.33)
n +1 n λ n n λ n n
UN = UN − α(UN − UN −1 ) + |α|(−UN + UN −1 ).
2 2
El código Python es el siguiente:
1 # Ecuacion de ondas
2 # Resolucion mediante el metodo upwind
3
4 from pylab import *
5 from time import clock
6
7 alpha = 1. # variable global ( velocidad de adveccion )
8
178 Métodos numéricos para problemas de contorno
9 def cond_ini ( x ):
10 """ Condicion inicial del problema """
11 return exp ( -8.* x **2)
12
13 def exacta (x , t ):
14 """ Solucion exacta """
15 return cond_ini (x - alpha * t )
16
17 # Datos del problema
18 a = -2. # extremo inferior del intervalo
19 b = 2. # extremo superior del intervalo
20 N = 100 # numero de particiones
21 h = (b - a )/ N # paso de malla espacial ( dx )
22 dt = 0.008 # paso temporal
23 T = 1. # tiempo final
24
25 tini = clock ()
26
27 x = arange (a , b +h , h ) # discretizacion del intervalo [a , b ]
28 t = arange (0 , T + dt , dt ) # discretizacion temporal
29
30 U = cond_ini ( x ) # condicion inicial
31 Unew = zeros ( len ( U )) # solucion en el siguiente paso de tiempo
32
33 lamb = dt / h
34
35 for n in range ( len ( t )): # bucle temporal
36 for j in range (1 , N ):
37 Unew [ j ] = ( U [ j ] -0.5* lamb * alpha *( U [ j +1] - U [j -1])
38 +0.5* lamb * abs ( alpha )*( U [ j +1] -2.* U [ j ]+ U [j -1]))
39 Unew [0] = ( U [0] -0.5* lamb * alpha *( U [1] - U [0])
40 +0.5* lamb * abs ( alpha )*( U [1] - U [0]))
41 Unew [ N ] = ( U [ N ] -0.5* lamb * alpha *( U [ N ] - U [N -1])
42 +0.5* lamb * abs ( alpha )*( - U [ N ]+ U [N -1]))
43
44 U = copy ( Unew ) # actualiza la solucion
45
46 tfin = clock ()
47
48 Uexacta = exacta (x , T ) # solucion exacta en el instante T
49
50 error = max ( abs (U - Uexacta )) # error cometido
51
52 # Resultados
53 print ( ' ----- ')
54 print ( ' Tiempo CPU : ' + str ( tfin - tini ))
55 print ( ' Error : ' + str ( error ))
56 print ( ' dx : ' + str ( h ))
57 print ( ' dt : ' + str ( dt ))
58 print ( ' ----- ')
59
60 # Dibujamos las soluciones
61 plot (x , Uexacta , 'r ') # dibuja la solucion exacta
3.10 Implementación en Python del método de diferencias finitas 179
9 El motivo lo entenderán fácilmente los usuarios algo más avanzados de Python: al igualar dos arrays se está haciendo
[1] G. Allaire. Numerical Analysis and Optimization: An Introduction to Mathematical Modelling and
Numerical Simulation. Oxford University Press, 2007.
[4] J. C. Butcher. Numerical Methods for Ordinary Differential Equations. Wiley, 2008.
[5] P. G. Ciarlet. Introduction to Numerical Linear Algebra and Optimisation. Cambridge University
Press, 1989.
[7] M. Cruzeix, A. L. Mignot. Analyse Numérique des Equations Différentielles. Masson, 1983.
[9] C. W. Gear. Numerical Initial Value Problems in Ordinary Differential Equations. Prentice Hall,
1971.
[10] E. Hairer, S. P. Norsett, G. Wanner. Solving Ordinary Differential Equations. vols. I-II, Springer,
1993.
[11] P. Henrici. Discrete Variable Methods in Ordinary Differential Equations. John Wiley & Sons,
1962.
[13] E. Isaacson, H. B. Keller. Analysis of Numerical Methods. John Wiley & Sons, 1966.
[15] J. D. Lambert. Computational Methods in Ordinary Differential Equations. John Wiley, 1973.
[19] R. J. LeVeque. Finite Difference Methods for Ordinary and Partial Differential Equations. SIAM,
2007.
[23] A. Quarteroni, F. Saleri. Cálculo Científico con MATLAB y Octave. Springer, 2006.
[25] A. H. Stroud. Numerical Quadrature and Solutions of Ordinary Differential Equations. Springer,
1974.
[26] E. Süli, D. Mayers. An Introduction to Numerical Analysis. Cambrigde University Press, 2003.