Calculo Numerico II

Fernando Quiros Gracian

ii

Cap
tulo 1

Metodos de un paso

1.1 Introduccion
El objetivo de esta primera parte del curso es encontrar soluciones aproximadas
para el problema de valor inicial
(

y 0(x) = f (x; y (x));

y (a) = ;

(1.1)

donde f : D  R  R d 7! R d es una funcion continua y (a;  ) 2 D, siendo

D un dominio. Estamos tratando con sistemas, y por tanto y , f y  tienen d
componentes cada uno,

y = (y 1 ; y 2; : : : ; y d)T ;

 = ( 1 ;  2 ; : : : ;  d )T ;

f = (f 1 ; f 2 ; : : : ; f d )T :

Este problema tiene una solucion unica bajo condiciones bastante generales,
que supondremos a partir de ahora.

Teorema 1.1 (Euler, 1768). Sea f : D  R R d 7! R d una funcion continua

y (a;  ) 2 D. Supongamos que f satisface la condicion de Lipschitz

kf (x; y) f (x; y^)k  Lky y^k

1

8(x; y); (x; y^) 2 D

para algun L  0. Entonces hay un intervalo I = [a

una unica solucion de (1:1) en I .

; a + ] tal que existe

Aunque el problema (1.1) tenga solucion, en la mayor parte de los casos no

es posible encontrar una forma cerrada para la misma por metodos analticos.
Habra que contentarse con obtener una aproximacion numerica.
Buscaremos una solucion de (1.1) en el intervalo [a; b], donde a y b son
nitos, con f satisfaciendo las hipotesis del teorema. Cuando haga falta, o
convenga para simplicar razonamientos, haremos hipotesis adicionales de regularidad sobre f .
Todos los metodos numericos para resolver (1.1) que veremos involucran la
idea de discretizacion; esto es, el intervalo continuo [a; b] se reemplaza por el
conjunto discreto de puntos fxn g denido por

xn = a + nh;

n = 0; 1; 2 : : : ; N = (b a)=h:

El parametro h es la longitud del paso. De momento lo consideraremos constante, aunque conviene mencionar que gran parte de la potencia de los algoritmos modernos proviene de su capacidad de cambiar h automaticamente a
medida que se realizan los calculos.
Denotamos por yn a la aproximacion de la solucion y (xn ) de (1.1) en xn ,

yn  y (xn ):
Nuestro objetivo es encontrar una forma de producir una sucesion de valores
fyng que aproxime a la solucion de (1.1) en el conjunto discreto de puntos fxng;
tal sucesion constituye una solucion numerica del problema (1.1). Los valores
yn se pueden interpolar para obtener una aproximacion de y (x) en puntos x
que no pertenezcan al conjunto xn .
Un metodo numerico sera un procedimiento para producir la sucesion de
valores aproximados fyng.
2

1.2 Metodo de Euler

Nuestro primer metodo es el metodo de Euler,

yn+1 = yn + hf (xn ; yn):

(1.2)

Vamos a llegar a el desde cuatro puntos de vista diferentes.

Punto de vista geometrico (d = 1). Utilizamos la recta tangente en xn

como aproximacion de la curva solucion en xn  x  xn+1 .

y6

Recta tangente
- ....
...
r
...
6
...
.....................
...
...
...
...
... y (x)
...
...
..
-x
.
xn xn+1

Entonces, llamando
r al incremento de la recta tangente en [xn ; xn+1 ],

se tiene

r
= y 0(xn ) = f (xn ; y (xn));
h
y (xn+1 ) y (xn) 
r = hy 0 (xn );

y (xn+1 )  y (xn ) + hf (xn ; y (xn)):

El metodo de Euler se obtiene imponiendo que la solucion numerica

satisfaga esta relacion de forma exacta.
Serie de Taylor. Desarrollamos y (xn+1 ) alrededor de xn y usamos la
ecuacion diferencial para obtener

h2 00
y ( n )
2
h2
= y (xn) + hf (xn ; y (xn)) + y 00 ( n ):
2

y (xn+1 ) = y (xn) + hy 0(xn ) +

La notacion  n indica que el punto intermedio  n 2 [xn ; xn+1 ] puede

variar de componente a componente.
Si despreciamos el termino de error obtenemos el metodo de Euler. El
residuo o cantidad despreciada,

h2 00
Rn = y (xn+1 ) y (xn) hf (xn ; y (xn)) = y ( n );
2

recibe el nombre de error de truncacion, o error de discretizacion del

metodo. Es lo que le falta a la solucion teorica para satisfacer la ecuacion
del metodo. Tambien se puede ver como el error que se comete al dar
un paso con el metodo, suponiendo que el valor en xn era exacto.
Diferenciacion numerica. De la denicion de derivada,

y (xn+1 ) y (xn)
h

 y0(xn ) = f (xn; y(xn));

de donde, despejando, de nuevo obtenemos el metodo de Euler.

Cuadratura numerica. Integramos y 0(x) = f (x; y (x)) en [xn ; xn+1 ], y
obtenemos
Z xn
y (xn+1 ) = y (xn) +
f (t; y (t)) dt:
+1

xn

Para calcular la integral utilizamos la formula de cuadratura numerica

Z c+h
c

g  hg (c);

llamado regla rectangular izquierda, y llegamos al metodo de Euler.

Todos estos puntos de vista se utilizan para obtener otros metodos mas
precisos.
Para generar una solucion numerica de (1.1) usando (1.2) necesitamos un
valor de arranque y0 . Uno deseara tomar y0 = y (x0 ), pero en general esto
no es posible, debido a la precision nita del ordenador. As pues, habra que
contentarse con tomar una aproximacion y0  y (x0 ).
4

>Da el metodo de Euler una buena aproximacion a la verdadera solucion

del problema (1.1)? Una medida posible de la bondad de la aproximacion es
el mayor error cometido
max ky (xn ) yn k:
1nN

Denicion 1.2. Decimos que un metodo numerico es convergente si para todo

problema (1:1) con solucion sucientemente regular se tiene que
lim max ky (xn) ynk = 0

si lim ky (x0)

h!0+ 1nN

h!0+

y0 k = 0:

Estamos pidiendo que el valor de arranque y0 aproxime bien al dato inicial

y (x0 ); si esto no es as no hay por que esperar que la solucion numerica se
parezca a la teorica.
Observacion. El producto Nh permanece constante, Nh = b

a.

Se tiene el siguiente resultado de convergencia.

Teorema 1.3. El metodo de Euler es convergente. Es mas, se tiene la siguiente cota para el error

ky(xn) ynk  e(b

a)L

ky(x0) y0k + 2CL e(b

a)L

1 h;

donde C = maxx2[a;b] ky 00(x)k y L > 0 es una constante de Lipschitz para

f (x; y ) con respecto a y .
Prueba. Para simplicar la prueba supondremos que f 2 C 1 ; derivando la
ecuacion se ve inmediatamente que esto implica que y 2 C 2 .

Tenemos que estimar en = y (xn)

la ecuacion diferencial tenemos que

yn. Utilizando la formula de Taylor y

h2 00
y ( n )
2
h2
= y (xn) + hf (xn ; y (xn)) + y 00 ( n );
2

y (xn+1 ) = y (xn) + hy 0(xn ) +

con  n 2 [xn ; xn+1 ]. Restandole a esta ecuacion la ecuacion (1.2) queda

en+1 = en + h(f (xn ; y (xn)) f (xn ; yn)) +

h2 00
y ( n ):
2

Usando la desigualdad triangular y la condicion de Lipschitz para f obtenemos

ken+1 k  kenk + hkf (xn; y(xn)) f (xn; yn)k + Ch2 =2

 kenk + hLkenk + Ch2 =2;
donde C = maxx2[a;b] ky 00 (x)k; es decir,

ken+1k  (1 + hL)kenk + Ch2 =2:

De aqu se tiene que

ke1 k 
ke2 k 


(1 + hL)ke0 k + Ch2 =2;

(1 + hL)ke1 k + Ch2 =2


(1 + hL)2 ke0 k + (1 + hL) + 1 Ch2 =2;

kenk 


(1 + hL)ken 1 k + Ch2 =2


(1 + hL)n ke0 k + (1 + hL)n 1 +


+ (1 + hL) + 1 Ch2 =2:

..
.

Usando la formula de la suma de los n primeros terminos de una progresion

geometrica tenemos que

kenk  (1 + hL)nke0k + Ch
((1 + hL)n
2L

1) :

Por otra parte, como hL es positivo, 1 + hL < ehL , y as

enhL
kenk  enhLke0k + Ch
2L

1 = e(b

a)L

ke0k + Ch
e(b
2L

si 0  n  N .
Ejemplo. Aplicamos el metodo de Euler al problema

y 0 (x) = g (x);

y (a) = 0;
6

a)L

para calcular de forma aproximada y (b) =

Rb
a g.

Tomando y0 = 0 tenemos que

y (b)  yN = yN 1 + h g (xN 1 ) = yN 2 + h g (xN 2) + h g (xN 1 )

=

=

As pues,

jy(xN )

NP1
j =0

h g (xj ):

Z
 b
yN =  g
 a

N
X1

j =0



h g (xj ) ;


y la convergencia del metodo de Euler nos garantiza la convergencia de la suma

de Riemann
N
X1

a la integral

j =0

Rb
a g.

h g (xj )

Ejemplo. Consideramos el problema lineal

y 0 (x) = y (x);

y (0) = :

La solucion exacta es y (x) = ex . Aplicando el metodo de Euler con y0 = 

obtenemos
yn+1 = yn + hyn = (1 + h)yn;
e iterando se llega a

yn = (1 + h)n :

As pues,

b N
= eb ;
lim yN = lim  1 +
N !1
h!0
N
lo que prueba la convergencia del metodo de Euler directamente. Sin embargo
esta prueba no da ninguna estimacion del error.
+

La cota proporcionada por el Teorema 1.3 en general sobreestima el error

en muchos ordenes de magnitud, y no se debe por tanto utilizar en la practica.
Sin embargo da una informacion importante, y es que el error es una O(h),
siempre y cuando ky (x0 ) y0 k = O(h). Eso motiva la siguiente denicion.
7

Denicion 1.4. Decimos que un metodo numerico es convergente de orden p

si para todo problema (1:1) con solucion sucientemente regular se tiene que
max ky (xn ) yn k = O(hp) si ky (x0 )

1nN

y0 k = O(hp) cuando h ! 0+ :

Que un metodo sea de orden p signica que para todos los problemas razonables el error que se comete con el metodo es de ese orden. Pero puede
suceder que para algun problema concreto el metodo sea mejor.
Ejemplo. Consideramos el problema

y 0(x) = 1;

y (0) = 0;

cuya solucion exacta es y (x) = x. Aplicando Euler con valor de arranque

y0 = 0 tenemos yn+1 = yn + h, y por tanto yn = nh = xn . Euler nos da por
tanto la solucion para este problema sin ningun error.
Si consideramos el problema

y 0(x) = x;

y (0) = 0;

la solucion exacta es y (x) = x2 =2. Al aplicar Euler con valor de arranque

y0 = 0 tenemos que

y1 = y0 + hx0 = 0;
y2 = y1 + hx1 = h2 ;
y3 = y2 + hx2 = h2 (1 + 2)
..
.
yn = h2 (1 +


+ (n 1)):
Es decir,

yn =
y por tanto

max ky (xn)

1nN

x2n
2

xn h
;
2

ynk = 1max
kx h=2k = bh=2;
nN n
8

lo que demuestra que el metodo no es de orden 2.

Para este problema la funcion f satisface una condicion de Lipschitz en y
para cualquier L  0; por tanto el Teorema 1.3 nos proporciona una estimacion
del error
max ky (xn) ynk  (ebL 1)h=2L;
1nN

para cada L > 0. Pasando al lmite L ! 0+ llegamos a

max ky (xn ) yn k  bh=2;

1nN

y la cota es en este caso optima.

1.3 Un metodo implcito: la regla del trapecio

El metodo de Euler aproxima la derivada en [xn ; xn+1 ] por una constante,
concretamente por su valor en xn . >Por que privilegiar al punto xn ? >No
sera mejor tomar, por ejemplo, como aproximacion constante de la derivada el
promedio de sus valores en los extremos del intervalo? En ese caso

f (xn ; y (xn)) + f (xn+1 ; y (xn+1))

:
2
Esta aproximacion conduce a la llamada regla del trapecio,


h
yn+1 = yn + f (xn ; yn ) + f (xn+1 ; yn+1) :
2
Veamos que converge con un orden mejor que el del metodo de Euler.
y (xn+1)  y (xn) + h

(1.3)

Teorema 1.5. La regla del trapecio es convergente de orden 2.

Prueba. Como en el metodo de Euler, denimos el error de truncacion Rn
como lo que le falta a la verdadera solucion para satisfacer la ecuacion del
metodo, es decir


h
y (xn+1) = y (xn) + f (xn ; y (xn )) + f (xn+1 ; y (xn+1)) + Rn :
(1.4)
2

Restando la ecuacion del metodo (1.3) de la ecuacion (1.4) obtenemos que

el error en = y (xn ) yn satisface la siguiente recurrencia

en+1 = en + h2 f[f (xn ; y (xn)) f (xn ; yn)]

+[f (xn+1 ; y (xn+1)) f (xn+1 ; yn+1 )]g + Rn :
Usando la desigualdad triangular y la condicion de Lipschitz para f obtenemos
que
ken+1k  kenk + h2 (Lkenk + Lken+1 k) + kRnk:
Suponiendo que hL=2 < 1 llegamos a

ken+1k 

1 + hL
2
1 hL
2

kenk + 1kRnhLk :
2

Para proseguir necesitamos una estimacion de kRn k que no dependa de n.

Desarrollando por Taylor alrededor de x = xn y usando la ecuacion diferencial
tenemos que

Rn = y (xn+1 ) fy (xn) + h2 [f (xn ; y (xn)) + f (xn+1 ; y (xn+1 ))]g

= y (xn ) + hy 0(xn ) + h2 y 00 (xn ) + h3! y 000 ( n )
f(y(xn) + h2 [y0(xn ) + y0(xn) + hy00(xn ) + h2 y000(n)]g
 000
000
= y 3!(n ) y (4n ) h3 :
3

As pues,

ken+1k 

1 + hL
2
hL
1 2

kenk + 1ChhL
2

con C =

5
max ky 000 (x)k:
12 x2[a;b]

Iterando, y usando la formula de la suma de los n primeros terminos de una

progresion geometrica tenemos que

kenk 

1 + hL
2
hL
1 2

!n

Ch2

ke0k + L

1 + hL
2
hL
1 2

!n

Usando ahora que

1 + hL
2 = 1 + hL < ehL=(1
hL
1 2
1 hL
2
10

hL )
2

1 :

concluimos que

kenk  e(b

a)L=(1

hL )
2

2

ke0 k + Ch
e(b
L

a)L=(1

hL )
2

si 0  n  N;

de donde se deduce inmediatamente la convergencia de orden 2.

La regla del trapecio tiene un orden de convergencia mayor que el metodo
de Euler. Eso no signica que sea un metodo mejor, pues hay otra diferencia
importante entre ambos metodos. El metodo de Euler es explcito; el valor
yn+1 viene dado explcitamente en terminos del valor anterior yn y se puede
calcular facilmente mediante una evaluacion de f y unas pocas operaciones
aritmeticas. Por el contrario, la regla del trapecio es un metodo implcito;
para calcular yn+1 hay que resolver un sistema de ecuaciones no lineales, lo
que en general es computacionalmente costoso.
En resumen, el orden no lo es todo. Aunque en un metodo de orden alto en
principio hay que dar menos pasos, estos pueden ser muy costosos; puede ser
preferible un metodo de orden mas bajo, en el que, a pesar de dar mas pasos,
cada uno de ellos lleve menos trabajo computacional.
Observacion. A pesar de todo, los metodos implcitos son utiles a veces, por
ejemplo para tratar problemas sti (ver el captulo 2).

1.4 Metodos de Taylor

El metodo de Euler es tan solo de primer orden. Eso supone un inconveniente. Si queremos un error muy peque~no nos veremos forzados a tomar h
muy peque~no, lo que supondra tener que dar muchos pasos, con un coste computacional que puede ser inaceptable. La regla del trapecio es de segundo
orden, pero a cambio es implcita. >Podemos construir metodos \mejores", de
mayor orden que el de Euler, y que a la vez sean explcitos?
La idea para ello sera tomar mas terminos en el desarrollo de Taylor, lo
11

que da lugar a los llamados metodos de Taylor. Como ejemplo construimos el

metodo de Taylor de orden 2, en el caso d = 1.
El desarrollo de Taylor produce

y (xn+1 ) = y (xn) + hy 0 (xn ) +

h3
h2 00
y (xn ) + y 000 (n ):
2
3!

El termino de error, de orden h3 , se desprecia. Para calcular y 0 (xn ) e y 00 (xn )

usamos la ecuacion diferencial,

y 0 (x) = f (x; y (x));

y 00 (x) = fx (x; y (x)) + fy (x; y (x))y 0(x) = fx (x; y (x)) + fy (x; y (x))f (x; y (x)):
Obtenemos el metodo

yn+1 = yn + hf (xn ; yn ) +

h2
(f (x ; y ) + fy (xn ; yn)f (xn ; yn)):
2 x n n

En este caso el error de truncacion Rn = h3! y 000 (n ) satisface que

3

kRnk  Ch3;
y esperamos que el metodo sea de orden 2 (lo veremos mas adelante).
Por el mismo procedimiento podemos generar metodos de Taylor del orden
que queramos. El inconveniente es que involucran cada vez mas derivadas.
Comparemos con el metodo de Euler. Los metodos de Taylor tienen mayor
orden, luego habra que dar menos pasos. A cambio, dar un paso es computacionalmente mas caro. No solo hay que evaluar f , sino tambien sus derivadas,
lo que puede ser muy costoso, sobre todo si estas son difciles. As, si no queremos mucha precision y la funcion es complicada, Euler es preferible. Si, por
el contrario, queremos bastante precision y la funcion y sus derivadas no son
complicadas, es preferible usar un metodo de Taylor de orden mas alto. En
resumen: no hay un metodo que sea \el mejor". Cual sea el mas conveniente
depende de la ecuacion diferencial y de la precision que queramos conseguir.
12

1.5 Metodos de un paso

El metodo de Euler, los metodos de Taylor y muchos otros, entre ellos los
metodos de Runge-Kutta, que veremos mas adelante, se pueden escribir como

yn+1 = yn + h(xn ; yn ; yn+1; h);

(1.5)

donde  es una funcion que depende de f , y que se denomina funcion de

incremento. Los metodos de esta forma reciben el nombre de metodos de un
paso, pues para calcular el paso yn+1 , solo necesitamos conocer un paso, el
anterior, yn .
Si (xn ; yn ; yn+1; h) no depende de yn+1 el metodo es explcito; en caso
contrario sera implcito y en cada paso tendremos que resolver un sistema no
lineal para obtener yn+1 .
Ejemplos.





El metodo de Euler es un metodo de un paso explcito con funcion de

incremento (xn ; yn ; yn+1; h) = f (xn ; yn).
La regla del trapecio es un metodo de un paso implcito con funcion de
incremento (xn ; yn ; yn+1; h) = 21 (f (xn ; yn ) + f (xn + h; yn+1)).
El metodo de Taylor de orden 2 es un metodo de un paso explcito
con funcion de incremento (para d = 1) dada por (xn ; yn; yn+1; h) =
f (xn ; yn) + h2 (fx (xn ; yn ) + fy (xn ; yn )f (xn ; yn)).

Queremos un teorema de convergencia para metodos de un paso. Para ello

tenemos que estimar en = y (xn) yn . Con el n de imitar las pruebas de
convergencia del metodo de Euler y de la regla del trapecio, denimos el error
de truncacion como lo que le falta a la verdadera solucion, y (x), para satisfacer
la ecuacion del metodo.

Denicion 1.6. El error de truncacion, Rn , del metodo (1:5) se dene por

y (xn+1 ) = y (xn) + h(xn ; y (xn); y (xn+1); h) + Rn :
13

(1.6)

Por conveniencia separamos un factor h y escribimos Rn = hn . Restando

a la ecuacion (1.6) la ecuacion del metodo queda

en+1 = en + h((xn ; y (xn); y (xn+1); h) (xn ; yn; yn+1 ; h)) + hn :

(1.7)

Usando la desigualdad triangular obtenemos

ken+1k  kenk + hk(xn; y(xn); y(xn+1); h) (xn; yn; yn+1; h)k + hkn k:
Para poder seguir pedimos que  satisfaga la condicion de Lipschitz

k(x; y; z; h) (x; y^; z^; h)k  L(ky y^k + kz z^k):

(1.8)

Notese que la constante L puede depender de f . Tenemos ahora que

ken+1k  kenk + hL(kenk + ken+1k) + hkn k:

Sea  = max0nN kn k. Entonces, si hL < 1,


h
k
en k +
:
ken+1k  11 + hL
hL
1 hL
A partir de esta desigualdad se tiene, repitiendo lo que hicimos para la regla
del trapecio el siguiente resultado.
Teorema 1.7. Si  satisface la condicion de Lipschitz (1:8), entonces, para
hL < 1 se tiene que


ky(xn) ynk  e2(b a)L=(1 hL) ky(x0) y0k + 2L e2(b a)L=(1 hL) 1 ;
donde  = max0nN kn k. Por tanto, si limh!0  = 0, el metodo converge.
Si ademas  = O(hp), entonces el metodo converge con orden p.
+

Este teorema motiva la siguiente denicion.

Denicion 1.8. Un metodo se dice consistente (respectivamente consistente

de orden p) si para todo problema (1:1) con solucion sucientemente regular
lim  = 0

h!0+

(respectivamente  = O(hp)):
14

Hemos probado que un metodo de un paso con  satisfaciendo la condicion de Lipschitz (1.8) y consistente (consistente de orden p) es convergente
(convergente de orden p).
>Que quiere decir la consistencia? Tenemos que
y (xn + h) y (xn)
= (xn ; y (xn ); y (xn + h); h) + n :
h
Tomamos xn = a + nh = x jo, y pasamos al lmite h ! 0+ . La consistencia
resulta ser equivalente a que y 0 (x) = (x; y (x); y (x); 0), es decir, a

(x; y; y ; 0) = f (x; y ):
Esta condicion garantiza que estamos resolviendo la ecuacion diferencial correcta, que estamos siendo consistentes con ella.
Ejemplos.

Para el metodo de Euler (x; y; z ; h) = f (x; y ); usando la condicion de

Lipschitz en y para f se tiene que

k(x; y; z; h) (x; y^; z^; h)k = kf (x; y) f (x; y^)k  Lky y^k;

de donde se sigue inmediatamente la condicion (1.8). Por otra parte

n = h2 y 00( n); por tanto  = O(h) si y 2 C 2 ([a; b]). Combinando ambos
hechos tenemos que el metodo de Euler es convergente de orden 1.
Para la regla del trapecio (x; y; z ; h) = 21 (f (x; y ) + f (x + h; z )); usando
la condicion de Lipschitz en y para f se tiene que

k(x; y; z; h) (x; y^; z^; h)k  21 kf (x; y) f (x; y^)k

+ 12 kf (x + h; z ) f (x + h; z^)k  L2 (ky y^k + kz z^k);
es decir, se cumple (1.8). Por otra parte
 000
y (

y 000 (n ) 2
n =
h;
3!
4
as  = O(h2 ) si y 2 C 3 ([a; b]), y por tanto la regla del trapecio es
convergente de orden 2.
n)

15

Para el metodo de Taylor de orden k, n = (kh+1)! y (k+1) ( n ), y si la solucion

es sucientemente regular se tiene que  = O(hk ).
k

1.6 Metodos de Runge-Kutta

Los metodos de Taylor son conceptualmente sencillos. Sin embargo:




Exigen calcular derivadas de orden superior, lo que puede ser engorroso.

Exigen evaluar dichas derivadas, lo que puede ser computacionalmente
costoso.

Para evitar estas dicultades, los metodos de Runge-Kutta evaluan f (x; y )

en mas puntos, intentando igualar la precision de los metodos de Taylor. Concretamente toman

(xn ; yn; yn+1; h) =

s
X
i=1

Ai f (i ; i )

de forma que (xn ; yn ; yn+1; h) sea una media ponderada de las pendientes
cerca de (xn ; yn ). La idea es elegir los parametros Ai , i , i de manera que la
funcion de incremento coincida lo mas posible con la funcion de incremento de
un metodo de Taylor.
Veamos por ejemplo como hacer esto con dos evaluaciones de la funcion y
d = 1. Parece razonable tomar (1 ; 1 ) = (xn ; yn ). Para 2 tomamos 2 = xn +
h. Una eleccion sensata para 2 consiste en estimar el valor y (2 ) avanzando
con el metodo de Euler, es decir

2 = yn + hf (xn ; yn):
Queda por tanto

yn+1 = yn + h A1 f (xn ; yn) + A2 f (xn + h; yn + hf (xn ; yn)) :

16

>Como elegir A1 , A2 y ? Desarrollamos por Taylor en potencias de h y

obtenemos

yn+1 = yn + (A1 + A2 )hf (xn ; yn)

+ A2 h2 (fx (xn ; yn ) + fy (xn ; yn)f (xn ; yn)) + O(h3 ):
Si hacemos

A1 + A2 = 1;

A2 = 1=2;

entonces la funcion de incremento  coincide hasta orden uno con la funcion

de incremento T 2 del metodo de Taylor de orden 2,

(xn ; yn; yn+1; h) = T 2 (xn ; yn; yn+1; h) + O(h2 ):

Veamos cual es el orden de consistencia del metodo:

hn = y (xn+1 ) y (xn) h(xn ; y (xn); y (xn+1 ); h)

= y (xn+1 ) y (xn) hT 2 (xn ; y (xn ); y (xn+1); h) + O(h3 )
= O(h3 ):
Concluimos que  = O(h2), es decir, el metodo es consistente de orden 2.
Por tanto, para ver que el metodo es convergente de orden 2 basta con
comprobar que la funcion de incremento

(x; y; z ; h) = A1 f (x; y ) + A2 f (x + h; y + hf (x; y )):

satisface la condicion de Lipschitz (1.8). Esto se prueba facilmente usando que
f es Lipschitz en y .
Hemos obtenido A2 = 21 , A1 = 1
unico metodo, sino toda una familia.

1
2

y  6= 0 arbitrario. No tenemos un

Ejemplos.

Tomando  = 1=2 obtenemos el llamado metodo de Euler modicado,

h
h
yn+1 = yn + hf (xn + ; yn + f (xn ; yn)):
2
2
17

Tomando  = 1 se tiene el metodo de Euler mejorado, tambien conocido

como metodo de Heun de orden 2,
h
yn+1 = yn + (f (xn ; yn ) + f (xn + h; yn + hf (xn ; yn)):
2
En ambos casos se realizan dos evaluaciones de funcion, frente a las tres
necesarias en el metodo de Taylor de orden 2 (a saber, f (xn ; yn), fx(xn ; yn) y
fy (xn ; yn)) cuando d = 1, siendo la cosa aun peor para problemas de dimension
mayor.
El metodo de Euler modicado se puede escribir alternativamente como

k1 = f (xn ; yn);
k2 = f (xn + h2 ; yn + h2 k1 );
yn+1 = yn + hk2 :
En cuanto a Euler mejorado, se puede escribir como

k1 = f (xn ; yn);
k2 = f (xn + h; yn + hk1 );
yn+1 = yn + h2 (k1 + k2 ):
Siguiendo estos ejemplos, denimos un metodo de Runge-Kutta de s etapas
para el problema (1.1) por
8
>
>
< ki

= f (xn + ci h; yn + h

>
>
: yn+1

= yn + h

s
P
j =1

s
P

j =1

aij kj );

i = 1; 2; : : : ; s

bj kj :

Cada una de las evaluaciones de funcion ki es una etapa.

El metodo (1.9) se representa por medio de su tablero de Butcher:

c1 a11 a12
c2 a21 a22
..
.. ..
. .
.
cs as1 as2
b1 b2
18

:::
:::
...
:::
...

a1s
a2s
..
.
ass
bs .

(1.9)

Ejemplos.

El tablero del metodo de Euler modicado es

0 0
0
1=2 1=2 0
0 1.

El tablero del metodo de Euler mejorado es

0
1

0
0
1
0
1=2 1=2.

Si escribimos

c = (c1 ; c2 ; : : : ; cs)T ;

b = (b1 ; b2 ; : : : ; bs )T ;

A = (aij );

el tablero se puede resumir como

A
bT .

Si aij = 0 para j  i, i = 1; 2; : : : ; s, es decir, si la matriz A es triangular

inferior estricta, entonces cada uno de los ki viene dado explcitamente en
terminos de los anteriormente calculados, kj , j = 1; 2; : : : ; i 1. En este caso
se dice que el metodo es explcito. Al escribir su tablero se suelen omitir los
ceros sobre y por encima de la diagonal principal.
Si el metodo no es explcito, es implcito. En general es necesario entonces
resolver en cada paso un sistema no lineal para calcular los ki . Este sistema
tiene dimension ds. Hay una situacion intermedia; si aij = 0 para j > i,
i = 1; 2; : : : ; s, entonces cada ki esta denido individualmente por

ki = f (xn + ci h; yn + h

i
X
j =1

aij kj );

19

i = 1; 2; : : : ; s:

En vez de resolver en cada paso un sistema no lineal de dimension ds, tendremos que resolver s sistemas no acoplados de dimension d. Estos metodos se
llaman semi-implcitos. Al escribir sus tableros se omiten los ceros por encima
de la diagonal principal.

1.7 Condiciones de orden

Los metodos de Runge-Kutta son metodos de un paso con funcion de incremento
s
X
(xn ; yn; yn+1; h) =
bj kj :
j =1

Si hacemos h = 0, entonces ki = f (xn ; yn) para todo i = 1; 2; : : : ; s. As,

(x; y; y ; 0) =

s
X
j =1

bj f (x; y ):

Por tanto, un metodo de Runge-Kutta es consistente si y solo si

s
X
j =1

bj = 1:

(1.10)

Por otra parte, puesto que las etapas ki son evaluaciones de la funcion f , no
es difcil convencerse de que  satisface la condicion de Lipschitz (1.8) si f
satisface una condicion de Lipschitz en y . As pues, por el Teorema 1.7, la
condicion de consistencia (1.10) es suciente para garantizar la convergencia.
Veamos que tambien es necesaria.

Teorema 1.9. Un metodo de Runge-Kutta es convergente si y solo si es consistente.

Prueba. Si es convergente, en particular lo es para el problema

y 0(x) = 1;

y (0) = 0;
20

cuya solucion es y (x) = x. El metodo, aplicado a este problema, se reduce a

yn+1 = yn + h

s
X
j =1

bj :

Tomando como valor de arranque y0 = 0 tenemos que la solucion numerica

viene dada por

yn = nh
Por tanto

s
X
j =1

bj = xn

y (xn) yn = xn 1

s
X
j =1

s
X
j =1

bj :
!

bj ;

y la solucion numerica converge a la teorica si y solo si se cumple la condicion

de consistencia (1.10).
>Que condiciones sobre los coecientes garantizan que el metodo es convergente de orden p?
Supongamos que y (xn ) = yn , la llamada hipotesis de localizacion. Usando (1.7) y el Teorema del Valor Medio tenemos que


@
(1.11)
Rn = I h (xn ; yn ;  ; h) (y (xn+1) yn+1 ):
@y
Esta igualdad nos dice que el termino principal del desarrollo de y (xn+1 ) yn+1
en potencias de h coincide con el de Rn . Por tanto y (xn+1 ) yn+1 = O(hp+1)
si y solo si Rn = O(hp+1 ).
As pues, lo que hay que hacer es obtener desarrollos en potencias de h para
y (xn+1 ), la solucion teorica, e yn+1 , la solucion numerica. Si estos desarrollos
coinciden hasta orden p, entonces Rn = O(hp+1) y tendremos convergencia de
orden p.
Para que los desarrollos de la solucion numerica y la teorica sean iguales
hasta orden p, es necesario que sus derivadas hasta orden p en h = 0 coincidan.
Eso dara lugar a unas ciertas condiciones sobre los coecientes del metodo, las
llamadas condiciones de orden. Obtenerlas sera nuestro siguiente objetivo.
21

Notacion. Los superndices en mayusculas denotaran componentes, y variaran

por tanto entre 1 y d. Los subndices en minusculas denotan etapas; varan
entre 1 y s. Finalmente, los subndices en mayusculas signican derivacion
parcial, y tambien varan entre 1 y d. Ejemplo:
J =
fKL

@f J
:
@y K @y L

Si el sistema es autonomo, y 0 = f (y ), podemos obtener mas simetra en

las formulas sustituyendo ki por gi tal que ki = f (gi). De esta forma el metodo (1.9) se transforma (para sistemas autonomos) en
8
>
J
>
< gi

= ynJ + h

>
J
>
: yn+1

= ynJ

s
P

aij f J (gj );

j =1
s
P

+h

j =1

bj

f J (g

j );

i = 1; : : : ; s;
J = 1; : : : ; d:

(1.12)

As, la unica diferencia entre las deniciones de giJ e ynJ+1 es que las constantes
aij en la primera formula pasan a ser las constantes bj en la segunda. Esta
semejanza nos permitira simplicar notablemente algunos calculos.
Si el sistema no es autonomo, pero el metodo de Runge-Kutta es consistente
y satisface la condicion de suma por las

ci =

s
X
j =1

aij ;

i = 1; 2 : : : ; s;

(1.13)

el metodo tambien se puede escribir en la forma (1.12). Para ello basta con
denir

y^n1 = xn ;
y^nJ +1 = ynJ ;

f^1 = 1;
f^J +1 = f J ;

g^i1 = xn + ci h;
s
P
g^iJ +1 = ynJ + h aij kjJ ; J = 1; : : : ; d:
j =1

En lo que sigue escribiremos siempre los metodos en la forma (1.12). Por

ello, para no dejar fuera los sistemas no autonomos, de aqu en adelante nos
restringiremos a metodos que cumplan la condicion de suma por las.
22

Obtengamos las condiciones de orden 3. Tenemos que calcular las derivadas

con respecto a h en h = 0 hasta orden 3 de la solucion numerica, ynJ+1 , y
compararlas con las de la teorica, y J (xn+1 ) = y J (xn + h), que tambien tenemos
que calcular.
Empezamos calculando las derivadas de la solucion numerica (siempre bajo
la hipotesis de localizacion). Aplicando la regla de Leibniz para la derivacion
de un producto se tiene que

(ynJ+1)(q) h=0

s
X

bj q (f J (gj ))(q

j =1


1) 

h=0 :

As, la derivada primera de la solucion numerica es


ynJ+1 (1) 

h=0

s
X
j =1


bj (f J (gj ))h=0

La derivada segunda viene dada por



s
X

(ynJ+1 )(2) h=0 = 2

Pero

f J (gj ) (1)

j =1
d
X

K =1

s
X

j =1

(1.14)

bj f J (yn):

(1.15)

bj (f J (gj ))(1) h=0 :

fKJ (gj ) gjK

(1)

Usando la semejanza de las deniciones de giJ e ynJ+1 junto con (1.15) se tiene
que
!
s
X

(giJ )(1) h=0 =

j =1

por lo que concluimos que

(2) 
J
yn+1 
h=0

=2

s
X
j;k=1

aij f J (yn);

bj ajk

d
X
K =1

(1.16)

fKJ (yn)f K (yn ):

(1.17)

Para calcular la derivada tercera necesitamos la derivada segunda de f J (gj ),

que viene dada por


f J (gj ) (2)

d
X
K;L=1

J (g )
fKL
j

gjK (1)

23

gjL (1)

d
X
K =1

fKJ (gj ) gjK

(2)

La derivada segunda de gjK se obtiene utilizando de nuevo la semejanza entre

giJ e ynJ+1 . Finalmente llegamos a

ynJ+1

(3) 

h=0

=3

s
P
j;k;l=1
s
P

bj ajk ajl

d
P

J (y )f K (y )f L (y )
fKL
n
n
n

K;L=1
d
P
+6
bj ajk akl
fKJ (yn )fLK (yn)f L (yn ):
K;L
=1
j;k;l=1

(1.18)

En cuanto a la solucion teorica, se tiene que

yJ

(1)

(xn+1 )

y J (2) (x



n+1 )

y J (3) (x



n+1 )

h=0

h=0

h=0

= f J (y (xn + h))h=0 = f J (yn );

=
=
+

d
P
K;L=1

d
P
K =1

fKJ (yn)(y K )(1) (xn ) =

d
P
K;L=1
d
P
K;L=1

d
P
K =1

fKJ (yn)f K (yn );

J (y )f K (y )(y L )(1) (x )
fKL
n
n
n

fKJ (yn)fLK (yn )(y L)(1) (yn ):

J (y )f K (y )f L (y ) +
fKL
n
n
n

d
P
K;L=1

fKJ (yn )fLK (yn)f L (yn):

Comparando las derivadas de la solucion numerica y de la teorica obtenemos las condiciones de orden 3:

 Orden 1:
 Orden 2:
 Orden 3:

s
P

bj
j =1
s
P
2
3

j;k=1
s
P

= 1;

bj ajk = 1;

j;k;l=1

bj ajk ajl = 1;

s
P
j;k;l=1

bj ajk akl = 1:

Un metodo de Runge-Kutta tendra orden de consistencia 3 si y solo si cumple

estas cuatro condiciones.
Observacion. La condicion de consistencia coincide con la condicion de consistencia de orden 1.

24

Observacion. Algunas de estas condiciones se pueden simplicar usando la

condicion de suma por las, quedando

 Orden 1:
 Orden 2:
 Orden 3:

s
P

bj = 1;

j =1
s
P

2
3

j =1
s
P
j =1

bj cj = 1;
bj c2j = 1;

(1.19)
6

s
P
j;k=1

bj ajk ck = 1:

Observacion. Las sumas que aparecen en las condiciones de orden se pueden

interpretar como productos de matrices, lo que puede ser comodo a la hora de
s
P
hacer calculos. Por ejemplo, la condicion 6
bj ajk ck = 1 se puede reescribir
j;k=1

como 6 bT Ac = 1.

Ejemplo. El metodo de Heun de orden 3, de tablero

0
0
1=3 1=3
2=3 0 2=3
1=4 0 3=4,
cumple la condicion de suma por las y las cuatro condiciones (1.19). Es por
tanto de orden 3.

 rboles de Butcher
1.8 A
La continuacion del proceso anterior para obtener condiciones de orden para
ordenes mas altos esta clara desde el punto de vista teorico. Pero conduce
rapidamente a formulas muy complicadas, cuya obtencion es larga y engorrosa.
Conviene utilizar una representacion graca que simplique las cuentas: los
arboles de Butcher (1964).
Empezamos observando que en (1.18) los ndices j; k; l aparecen emparejados en ajk , ajl , . . . exactamente igual que los subndices y superndices J; K; L
25

J , f J , . . . Concretamente, en el primer sumando,

en las expresiones fKL
K

s
X
j;k;l=1

bj ajk ajl

d
X
K;L=1

J (y )f K (y )f L (y );
fKL
n
n
n

J , yendo
encontramos la asociacion k 7! j , l 7! j , a traves de ajk , ajl y de fKL
de derecha a izquierda en los primeros y de abajo hacia arriba en el ultimo
caso. Esta asociacion se puede representar gracamente como
k

t31 =
j

En el segundo sumando
6

s
X
j;k;l=1

(1.20)

bj ajk akl

d
X
K;L=1

fKJ (yn )fLK (yn)f L (yn)

encontramos la asociacion k 7! j , l 7! k, que se puede representar como

l
k

t32 =

(1.21)

Observacion. En la representacion graca usamos siempre letras minusculas.

No obstante, si la asociacion que estamos representando es la dada por los
J , . . . , habr
subndices y superndices en las expresiones fKL
a que interpretar
los ndices del graco como mayusculas.

Los grafos t31 y t32 son ejemplos de arboles etiquetados con raz.

Denicion 1.10. Sea A un conjunto ordenado de ndices, A = fj1 < j2 <

: : : g. Denotamos por Aq al subconjunto de los q primeros ndices, Aq = fj1 <
j2 <


< jq g. Un arbol etiquetado con raz de orden q , q  1, es una
aplicacion
t : Aq n fj1 g 7! Aq

tal que t(ji ) < ji ;

26

i = 2; : : : ; q:

Notacion. Casi siempre, salvo que q sea un numero indeterminado o grande,

usaremos el conjunto ordenado de ndices A = fj < k < l < m < : : : g.

La aplicacion t se representa gracamente de la siguiente manera: cada

ndice ji 2 Aq se representa por un vertice, y la correspondencia t(jk ) = ji
mediante una arista que une los vertices ji y jk , situandose el vertice ji por
debajo del jk .
Ejemplos. El arbol de orden 3, t(k) = j , t(l) = j es el arbol t31 representado
en (1.20), y el arbol de orden 3, t(l) = k, t(k) = j es el arbol t32 representado
en (1.21).

El conjunto de todos los arboles etiquetados con raz de orden q se denota

por LTq (por labelled trees). Decimos que ji es el hijo de t(ji ) y que t(ji ) es
el padre de ji . Con este lenguaje la condicion t(ji ) < ji dice que el padre esta
antes que el hijo en el arbol genealogico. El vertice j1 , antecesor de toda la
dinasta, recibe el nombre de raz de t. El orden q de un arbol etiquetado con
raz es igual al numero de sus vertices.

Denicion 1.11. Se llama diferencial elemental asociado al arbol etiquetado

t 2 LTq a la expresion
F J (t)(y ) =
1

d
X

q
Y

J2 ;:::;Jq =1 i=1

ftJi (Ji ) (y ):
1

Es decir, el diferencial elemental de un arbol se obtiene sumando sobre

todos los ndices del arbol menos el de la raz, siendo cada sumando el producto
de q factores f , donde los superndices recorren todos los vertices de t y los
subndices son los hijos correspondientes.
Ejemplos.

 F J (t31 )(y) =
 F J (t32 )(y) =

d
P
K;L=1
d
P
K;L=1

J (y )f K (y )f L (y );
fKL

fKJ (y )fLK (y )f L(y ).

27

Observacion. Los tres arboles etiquetados

m

m
k

tienen el mismo aspecto topologico. Es mas, sus diferenciales elementales

d
X
K;L;M =1

d
X

J f K f M f L;
fKM
L

K;L;M =1

d
X

J f K f Lf M ;
fKL
L

K;L;M =1

J f L f Kf M;
fLK
M

son todos el mismo (basta con cambiar los nombres de losndices de sumacion).
Esto motiva la siguiente denicion.

Denicion 1.12. Dos arboles t; t0 2 Tq se dicen equivalentes si existe una

permutacion  : Aq 7! Aq tal que  (j1 ) = j1 y t  =  t0 en Aq n fj1 g.
Ejemplo. Los arboles

t=

t0 =

son equivalentes. La permutacion que permite pasar de un arbol a otro es

=

j k l m
j k m l

Observese que la imagen por  de una etiqueta de un nodo en el arbol t es

la etiqueta que tiene ese mismo nodo en el arbol t0 ; con la permutacion  lo
unico que estamos haciendo es cambiar los nombres de los ndices.

Denicion 1.13. Cada clase de equivalencia de arboles etiquetados con raz

de orden q constituye un arbol con raz de orden q .
28

El conjunto de todos los arboles de orden q se denota por Tq . El orden de

un arbol se dene como el orden de un representante. Dado t 2 Tq , denotamos
por (t) al numero de elementos en la clase de equivalencia de t.

Denicion 1.14. Sea t 2 LTq , (Aq = fj1 <


< jq g). Denimos j (t) por
1

j (t) =
1

s
X

q
Y

j1 ;:::;jq =1 i=2

at(ji )ji :

Es decir, j (t) es la suma a todos los ndices del arbol menos el de la raz
de los productos de las a@ (aristas) , donde @ ( j ) = ij:
1

Ejemplo.

t=

l
m

j (t) =

Ps
k;l;m=1 ajk akl ajm .

Observacion. Si t y t0 son equivalentes entonces j (t) = j (t0 ). Notese que

unica diferencia son los nombres de los ndices de sumacion, que son mudos.

Denicion 1.15. Denotamos por t = [t1 ; : : : ; tm ] al arbol que produce los

arboles t1 ; : : : ; tm cuando se le corta la raz. Los arboles t1 ; : : : ; tm son las
ramas del arbol t.
Ejemplo.

t1 =

t2 =

t = [t1 ; t2 ]

Denicion 1.16. Dado t 2 LTq se dene

(t) = 1
r
Q
(t) = q (ti );
i=1

si q = 1;
si t = [t1 ; : : : ; tr ]; q > 1:

29

Es decir, (t) es el producto de (t) y de los ordenes de los arboles que van
apareciendo al ir quitando las races.
Ejemplo.

7!

7!

= 9

7!

2
6

1
4
1

1
1
1

Observacion. Si t y t0 son equivalentes entonces (t) = (t0 ).

Las derivadas de las soluciones teorica y numerica se pueden escribir de

forma compacta usando el lenguaje de arboles, como se hace en los siguientes
teoremas, debidos a Butcher. Posponemos las pruebas, que se pueden omitir
en una primera lectura, hasta la seccion siguiente.

Teorema 1.17. La solucion teorica satisface

(y J )(q) =

X
t2LTq

F J (t)(y ):

Teorema 1.18. La solucion numerica satisface


(ynJ+1 )(q) h=0

X
t2LTq

(t)

s
X
j =1

bj j (t)F J (t)(yn):

(1.22)

Combinando los teoremas 1.18 y 1.17 obtenemos las condiciones de orden.

Teorema 1.19. Un metodo de Runge-Kutta es de orden p si y solo si

(t)

s
X
j =1

bj j (t) = 1

para todos los arboles t 2 Tq con q  p.

30

(1.23)

Prueba. Comparando las derivadas de la solucion numerica y de la teorica,

vemos que si se cumple (1.23) para cada arbol t 2 LTq , el metodo es de
orden p. Ahora bien, puesto que (t) y j (t) coinciden para todos los arboles
de la misma clase, basta con comprobar las condiciones para los arboles sin
etiquetar t 2 Tq .

Solo queda probar la implicacion recproca. Para ello basta con comprobar
que los diferenciales elementales de los diferentes arboles son independientes.
En efecto, dado t 2 Tq , el sistema

yj0 i =

Y
jr 2 t 1 (ji )

yjr ;

ji = 1; : : : ; q;

verica que F 1 (t)(yn ) = 1 y F 1 (u)(yn) = 0 para todos los demas arboles

u 2 Tq .
Ejemplo. Consideramos el metodo de Runge-Kutta explcito de cuatro etapas
de tablero
0
1=2 1=2
(1.24)
1 1=2 3=2
1
0
4=3
1=3
1/6 2/3
0
1/6.

Notese que el tablero cumple la condicion de suma por las.

Hay un unico arbol de orden 1,

=

que da lugar a la condicion de orden 1

s
X
j =1

bj = 1;

que evidentemente se cumple para nuestro tablero.

De orden 2 tenemos el arbol
31

t21 =
que produce la condicion de orden 2

j
s
P
j;k=1

bj ajk = 1. Usando la condicion de

suma por las la condicion se transforma en

2

s
X
j =1

bj cj = 1:

Nuestro tablero la satisface.

Observacion. El efecto de la condicion de suma por las es que desaparece la
s
P
suma en el ndice del nodo terminal (sin hijos) k, ajk , y se sustituye por la
k=1
componente de c correspondiente al padre de k, en este caso cj .

Hay dos arboles de orden 3,

t31 =

t32 =

j
j
que dan lugar (usando la condicion de suma por las) a las condiciones de
orden
3

s
X
j =1

bj c2j

= 1;

s
X

j;k=1

bj ajk ck = 1:

Nuestro metodo cumple ambas.

Hay cuatro arboles de orden 4,
k l m
t41 =

m
t42 =

j
t43 =

k
j
m

t44 =

k
j

k
j

32

Las condiciones correspondientes son

4

s
P

bj c3j = 1;

j =1
s
P

12

j;k=1

bj ajk c2k = 1;

s
P

bj cj ajl cl = 1;

j;l=1
s
P

24

j;k;l=1

bj ajk akl cl = 1:

Nuestro tablero cumple las condiciones correspondientes a t41 y t42 , pero no

las otras dos, pues
s
X
j;l=1

s
X

bj cj ajl cl = 1=6;

j;k=1

bj ajk c2k = 0:

El metodo es por tanto de orden 3, pero no de orden 4.

1.9 Demostracion de los teoremas de Butcher

En esta seccion probamos los teoremas 1.17 y 1.18, que proporcionan respectivamente las derivadas de las soluciones teorica y numerica.
Prueba del teorema 1.17. La prueba es por induccion en q .

Empezamos viendo que el resultado es cierto para q = 1. En efecto, solo

hay un arbol de orden q = 1, que se suele denotar por  ,

=

j1

para el que se tiene F J ( ) = f J .

1

Supongamos ahora que el resultado es cierto para q . Entonces

(y J1 )(q)

d
X

t2LTq J2 ;:::;Jq =1 i=1

0

Por la regla de la cadena, ftJr (Jr ) =

1

q
Y

d
P
K =1

33

ftJi (Ji ) :
1

ftJr (Jr )K (y K )(1) , as que derivando

1

(y J )(q) obtenemos
1

(y J )(q+1) =
1

d
P

q
P

q
Q

t2LTq J2 ;:::;Jq =1 r=1 i=1;i6=r

q
q
d
P P
P
Q

ftJi (Ji )

t2LTq r=1 J2 ;:::;Jq+1 =1 i=1;i6=r

d
P
K =1

ftJr (Jr )K (y K )(1)

1

ftJi (Ji ) ftJr (Jr )Jq f Jq :

1

+1

+1

La clave para proseguir es observar que todo arbol de orden q + 1 se obtiene

a~nadiendo en algun lugar una nueva arista a un arbol de orden q . Esto permite
establecer la aplicacion biyectiva

LTq  Aq
(t; jr )
donde

7! LTq+1
7
!
t0 ;

t0 jAq = t;

t0 (jq+1 ) = jr :

Es decir, el arbol t0 de orden q + 1 se obtiene a partir del arbol t de orden q

a~nadiendo una arista en el nodo jr .
Veamos un ejemplo gracamente:
jq+1

7!

jr

jr

Con ayuda de esta aplicacion podemos escribir

(y J )(q+1) =
1

d
P

q
Q

f(Jti0 )

t0 2LTq+1 J2 ;:::;Jq+1 =1 i=1;i6=r

qQ
+1
d
P
P
f Ji0 1
0t 2LTq+1 J2 ;:::;Jq+1 =1 i=1 (t ) (Ji )

Jq+1
Jr
(Ji ) f(t0 ) 1 (Jr ) f(t0 ) 1 (Jq+1 )

P
t0 2LTq+1

F J (t0 )(y ):
1

Nuestro siguiente objetivo es calcular las derivadas de la solucion numerica.

Usaremos la formula (1.14), por lo que necesitamos calcular (f J (gj ))(q 1) , que
es lo que hacemos a continuacion.
34

Derivadas

Gracamente

(f J (gj ))(0) = f J
(f J (gj ))(1) =
(f J (gj ))(2) =

d
P

K =1
d
P

fKJ (gj )(gjK )(1)

K;L=1

+
(f J (gj ))(3) =

J (g )(g K )(1) (g L )(1)

fKL
j j
j

d
P
K =1

+
+
+

d
P
K;L=1
d
P
K;L=1

J (g )(g K )(2) (g L )(1)

fKL
j j
j
J (g )(g K )(1) (g L )(2)
fKL
j j
j

d
P
K;M =1
d
P
K =1

J (g )(g K )(2) (g M )(1)

fKM
j j
j

fKJ (gj )(gjK )(3)

j
j

J
fKLM
(gj )(gjK )(1) (gjL)(1) (gjM )(1)

K;L;M =1

fKJ (gj )(gjK )(2)

d
P

j
k

l
m

Ahora debera estar claro el proceso. Cada vez que diferenciamos tenemos
que:
J ; es decir, a~
(i) derivar el primer factor fK:::
nadir una nueva rama a la raz;

(ii) aumentar el numero de derivadas de cada una de las g en una unidad;

es decir, alargar la correspondiente rama.
Todos los arboles que se obtienen de esta forma son aquellos arboles \especiales" que no tienen ramicaciones salvo en la raz; y se obtienen todos.
35

Denicion 1.20. Un arbol t

tiene que #t 1 (ji ) 2 f0; 1g.

2 LTq

se dice especial si para i = 2; : : : ; q , se

Al conjunto de arboles especiales se lo denota LSq .

Tenemos el siguiente resultado, cuya prueba es inmediata empleando el
lenguaje de arboles.

Lema 1.21. (Formula de Faa di Bruno) Para q  1 se tiene

(f J (g

))(q 1)

d
X

fKJ ;:::;Kr (gj )(gjK )((u ))


(gjKr )((ur )) ;

1

K1 ;:::;Kr
u 2 LSq ;
u = [u1 ; : : : ; ur ]

donde r = r(u) es el numero de ramas que dejan la raz y (u1 ); : : : ; (ur ) es

el numero de nodos en cada una de esas ramas.
Observacion. Se tiene que q = 1 + (u1 ) +


+ (ur ).
Prueba del teorema 1.18. Debido a la semejanza de las deniciones de giJ e
ynJ+1 , probar (1.22) es equivalente a demostrar que

(giJ )(q) h=0

X
t2LTq

(t)

s
X
j =1

aij j (t)F J (t)(yn ):

(1.25)

Probaremos esta ultima formula por induccion en q .

Veamos en primer lugar que el resultado es cierto para q = 1. Para el unico
arbol  de orden 1 tenemos que j ( ) = 1 y F J ( ) = f J . As que usando (1.16)
tenemos que
P
t2LT1

(t)

s
P
j =1

aij j (t)F J (t)(yn ) =

=

s
P
j =1
s
P
j =1

y por tanto el resultado es cierto para q = 1.

36

aij j ( )F J ( )(yn )


aij f J (yn) = (giJ )(1) h=0 ;

Supongamos ahora el resultado cierto para todo numero de derivadas menor

o igual que q 1, y veamos que es cierto para q derivadas.
Usando la regla de Leibniz y la formula de Faa di Bruno tenemos que


(giJ )(q) h=0 = q

=q

s
P
j =1

aij

s
P

aij (f J (gj ))(q

j =1
P

d
P


1) 

y=yn

fKJ ;:::;Kr (gj )(gjK )((u ))


(gjKr )((ur ))

1

u 2 LSq ; K1 ;:::;Kr
u = [u1 ; : : : ; ur ]

Dado que (u1 ); : : : ; (ur ) < q , podemos aplicar la hipotesis de induccion para
obtener


(giJ )(q) h=0 =

q (t1 ) (tr )
t
2
LT
t
2
LT
r
1

(
u
)

(
u
)
r
u 2 LSq ;
1
u = [u1 ; : : : ; ur ]
d
P
P
P
P
aij ajk1 k1 (t1 )
ajkr kr (tr )
fKJ 1 ;:::;Kr F K1 (t1 )
j
K1 ;:::;Kr
k1
kr

Sea u = [u1 ; : : : ; ur ]
biyeccion

2 LSq , y tl 2 LT(ul ) , para l = 1; : : : ; r.

F Kr (tr ):
Denimos una

(u; t1 ; : : : ; tr ) 7! t 2 LTq

de forma que

(t) = q (t1 )


(tr );

P
P
P
j (t) = aij ajk k (t1 )


ajkr kr (tr );
j

F J (t) =

d
P

k1

kr

fKJ ;:::;Kr F K (t1 )


F Kr (tr ):

K1 ;:::;Kr

El arbol t se obtiene sustituyendo las ramas de u por los arboles t1 ; : : : ; tr y

tomando las etiquetas correspondientes de forma natural, es decir, en el mismo
orden.
Ejemplos.

37

p
n

m
l

7!

k
j

j
l

p
n
l

m
k

7!

k
j
p

p
m
k

q
n

j
r

l
k

7!

n
l

m
k

Usando esta biyeccion llegamos nalmente a (1.25).

1.10 Metodos explcitos. Orden obtenible

Hemos obtenido las condiciones para que un metodo de Runge-Kutta tenga
un orden dado. Ahora dirigimos nuestra atencion al problema de encontrar
soluciones a dichas condiciones de orden. Nos preguntamos que orden se puede
lograr con un metodo explcito de s etapas.
Nuestro primer resultado (debido a Butcher) da una cota superior para el
orden que se puede conseguir con un metodo explcito de s etapas. La prueba
es un ejemplo de la potencia del enfoque de arboles de Butcher.

Teorema 1.22. Un metodo de Runge-Kutta explcito de s etapas no puede

tener orden mayor que s.
Prueba. Supongamos que el metodo de s etapas tiene orden p. Consideramos

38

el arbol de orden p

jp
jp

t=

j3
j2
j1

Evidentemente (t) = p ! y por otra parte,

s
X
j1 =1

bj j (t) =
1

s
X
j1 ;:::;jp =1

bj aj j aj j
1

1 2

2 3

ajp jp :

Como el metodo es explcito, aij = 0 para j  i. As,

s
P
j1 =1

bj j (t) = 0 a menos
1

que exista una sucesion j1 ; j2 ; : : : ; jp de enteros 1; 2; : : : ; s tal que

j1 > j2 > j3 >


> jp :

Dicha sucesion solo puede existir si p  s. En caso contrario

(t)

s
X

j1 =1

bj j (t) = 0 6= 1
1

y no se cumple la condicion de orden correspondiente al arbol t.

>Se puede conseguir orden p = s? Para s = 1; 2; 3; 4 se puede ver que s,
y es relativamente facil encontrar familias de soluciones a las correspondientes
condiciones de orden.
Ejemplo. Sea s = 1. La condicion de orden 1 es en este caso b1 = 1. El unico
metodo explcito de una etapa y orden 1 es por tanto el metodo de Euler.
Ejemplo. Sea s = 2. Las condiciones de orden 1 y 2 son en este caso
b1 + b2 = 1 y b1 c1 + b2 c2 = 1 respectivamente (estamos suponiendo que se
cumple la condicion de suma por las). Si el metodo es explcito entonces
c1 = 0, y la condicion de orden 2 queda que b2 c2 = 1=2. Hay pues una familia
uniparametrica de metodos explcitos de 2 etapas y orden 2,
1
1
c2 = ;
b1 = 1
;
b2 = ;
 2 R ;  6= 0:
2
2

39

>Que pasa con s = 5? Para que un metodo sea de orden 5 se tienen que
satisfacer 17 condiciones de orden. Si s = 5 tenemos solo 15 parametros libres.
Sin embargo, el sistema es no lineal, as que podra haber solucion. Pero no es
as.

Teorema 1.23. Para p  5 no existe ningun metodo de Runge-Kutta explcito

de orden p con s = p etapas.
La prueba se puede encontrar en [HNW].
Este es un primer ejemplo de barrera de Butcher. Hay mas barreras. Por
ejemplo, se tiene el siguiente teorema.

Teorema 1.24. (i) Para p  7 no existe ningun metodo de Runge-Kutta

explcito de orden p con s = p + 1 etapas.
(ii) Para p  8 no hay ningun metodo de Runge-Kutta explcito con s = p + 2
etapas.

En la siguiente tabla recogemos el estado del arte:

Orden 1 2 3 4 5 6 7 8
9
10
s mnimo 1 2 3 4 6 7 9 11 12  s  17 13  s  17

1.11 Estimaciones de error

Supongamos que f

2 C p. Por un lado tenemos por (1.9) que

Rn = y (xn + h) y (xn) h

s
X
i=1

bi ki :

(1.26)

Por otra parte, usando la formula de Taylor se tiene que

p
y (xn + h) = y (xn ) + y 0(xn )h + y 00(xn ) h2 +


+ y (p+1) (xn + h) (hp+1)!
;
p
(p)
h
0
ki (h) = ki (0) + ki (0)h +


+ ki (i h) p! :
2

40

+1

Si el metodo (1.12) satisface las condiciones de orden hasta orden p, entonces

los terminos hasta orden p en (1.26) se anulan. Se tiene entonces que
!

s
1
1X
(
p
+1)
max ky
(xn + th)k +
jb j max kk(p)(th)k :
(p + 1)! t2[0;1]
p! i=1 i t2[0;1] i

kRnk  hp+1

Incluso suponiendo conocidas las derivadas que aparecen en la formula anterior,

esta cota es con frecuencia mucho mayor que el tama~no que realmente tiene el
error de truncacion, y no es por tanto adecuada para estimarlo. Es mucho mas
realista considerar el primer termino distinto de cero del desarrollo de Taylor
del error. El trabajo realizado hasta ahora nos da dicho termino.

Teorema 1.25. Si un metodo de Runge-Kutta satisface las condiciones de

orden p, y f 2 C p+1 , entonces el error de truncacion satisface
RnJ =

hp+1 X
(t)e(t)F J (t)(yn) + O(hp+2 );
(p + 1)! t2Tp
+1

donde

e(t) = 1 (t)

s
X
j =1

bj j (t):

Las expresiones e(t) reciben el nombre de coecientes de error.

Ejemplo. Calculemos el termino principal del error al aplicar el metodo de
tablero (1.24) al problema lineal

y 0 = y;

y (0) = ;

(1.27)

donde  es una matriz d  d con elementos IJ constantes. La verdadera

solucion es y (x) = ex  .
Se tiene que

e(t41 ) = 0; e(t42 ) = 1=3; e(t43 ) = 1; e(t44 ) = 0;

41

de manera que los unicos diferenciales elementales que hay que calcular son
los correspondientes a los arboles t42 y t43 ,

F J (t

42 )

d
X
K;L;M =1

J f K f Lf M ;
fKM
L

F J (t

43 )

d
X
K;L;M =1

K f Lf M :
fKJ fLM

En nuestro caso f es lineal; cualquier derivada segunda suya se anula, y por

tanto F J (t42 ) = 0, F J (t43 ) = 0. Concluimos que RnJ = O(h5 ), y el metodo es
al menos de orden 5 para problemas del tipo (1.27).
Veamos cual es el termino de orden 5. Todos los diferenciales que tengan
derivadas de orden superior a uno se van a anular, por ser f lineal. Habra
derivadas de orden superior cuando haya alguna bifurcacion en el arbol. As
pues, el unico arbol de los nueve de orden 5 cuyo diferencial elemental no se
anula es
o

t59 =

k
j

Teniendo en cuenta que fKJ = JK , obtenemos que

F J (t

59 )

d
X
K;L;M;O=1

fKJ fLK fML fOM f O = (5 y )J :

Por otra parte

e(t59 ) = 1 (t59 )

s
X
j =1

bj j (t59 ):

El orden del arbol lineal t59 es mayor que el numero de etapas. As pues (ver
s
P
la prueba del Teorema 1.22), bj j (t59 ) = 0 , y por tanto e(t59 ) = 1. Como
j =1

(t59 ) = 1, concluimos nalmente que

RnJ =

h5 5 J
( yn) + O(hp+2 ):
5!
42

La estimacion del error simplicada,

hp+1 X
RnJ 
(t)e(t)F J (t)(yn );
(p + 1)! t2Tp

(1.28)

+1

que se contenta con el termino principal del error, es tambien de poco interes
practico, pues requiere calcular varias derivadas parciales de orden superior. Y
precisamente la principal ventaja de los metodos de Runge-Kutta (comparados
con metodos de Taylor), es que no exigen calcular derivadas.
Sin embargo, necesitamos una estimacion practica para el error. Por un
lado para asegurarnos de que los tama~nos de paso hi son sucientemente
peque~nos como para tener la precision deseada en los resultados, y por otra
parte para asegurarnos de que los tama~nos de paso son sucientemente grandes
como para evitar trabajo computacional innecesario.

1.12 Extrapolacion de Richardson

La idea es utilizar de forma cuidadosa el comportamiento del error como funcion de h. Esta idea, debida a Richardson, se utiliza en otros contextos,
por ejemplo para estimar el error en formulas de cuadratura numerica. Expliquemosla.
Sea a(h) una aproximacion a b cuando h tiende a 0. Queremos estimar el
error b a(h) sin usar b, pues no lo conocemos. Supongamos que sabemos que
b coincide con a(h) hasta orden p,

b a(h) = Chp+1 + O(hp+2):

Si aplicamos esta formula con paso 2h, se tiene que

b a(2h) = 2p Chp + O(hp+1):

de donde, eliminando Chp+1 se tiene que
2p+1 a(h) a(2h)
+ O(hp+2);
b=
p
+1
2
1
43

y de ah la estimacion de error, conocida como estimacion de Richardson,

b a(h) =

a(h) a(2h)
+ O(hp+2):
2p+1 1

Observacion. La cantidad

2p+1 a(h)
2p+1

a(2h)
;
1

es una aproximacion a b de orden p + 1.

>Como aplicar la idea de Richardson para estimar el error de truncacion?
Supongamos dado un \valor inicial" (xn ; y (xn)) y un tama~no de paso h. Si el
metodo es de orden p, entonces

y (xn+1) yn+1 = (y (xn))hp+1 + O(hp+2);

(1.29)

donde (y (xn ))hp+1 es el termino principal del error de truncacion,

(y )hp+1 =

hp+1 X
(t)e(t)F J (t)(y ):
(p + 1)! t2Tp
+1

Damos otro paso de tama~no h. Si en un paso el error es aproximadamente

(y (xn ))hp+1 , en dos sera aproximadamente 2 (y (xn ))hp+1. En efecto, probaremos mas adelante que la aproximacion numerica yn+2 verica que

y (xn+2 ) yn+2 = 2 (y (xn))hp+1 + O(hp+2):

(1.30)

Por otra parte, si empezando en (xn ; y (xn)) damos un paso de tama~no 2h,
la aproximacion numerica resultante w satisface que

y (xn+2) w = (y (xn))2p+1 hp+1 + O(hp+2):

Restando esta igualdad de la anterior y despejando se tiene que
2 (y (xn))hp+1 =

yn+2 w
+ O(hp+2):
p
2 1
44

En resumen: hemos probado que el error de yn+2 satisface

y
w
y (xn+2) yn+2 = n+2
+ O(hp+2);
(1.31)
p
2 1
lo que da una forma simple de estimar el error. Por otra parte, podemos
obtener un valor mejorado, yn+2 , para la aproximacion de y (xn+2), mediante
y
w
yn+2 = yn+2 + n+2
p
2 1
aproximacion que resulta ser de orden p + 1.
Solo queda probar (1.30). Para ello descomponemos y (xn+2 )
suma de dos terminos,

yn+2 como

y (xn+2 ) yn+2 = (y (xn+2) y^(xn+2 )) + (^y(xn+2 ) yn+2 );

(1.32)

donde y^ es la solucion de

y^0 = f (x; y^);

y^(xn+1 ) = yn+1 :

El primer sumando de (1.32) es la diferencia de dos soluciones exactas de la

misma ecuacion con distintos valores iniciales en x = xn+1 . Tenemos que
R

y (x) = y (xn+1 ) + xxn f (s; y (s)) ds;

R
y^(x) = yn+1 + xxn f (s; y^(s)) ds:
+1

+1

Restando y tomando normas, y usando que f es Lipschitz en y se tiene que

R

ky(x) y^(x)k  ky(xn+1) yn+1k + xxn kf (s; y(s)) f (s; y^(s))k ds

R
 ky(xn+1) yn+1k + L xxn ky(s) y^(s)k ds:
+1

+1

Para  > 0 arbitrario denimos

g (x) =  + ky (xn+1) yn+1k + L

As pues

Z x
xn+1

(1.33)

ky(s) y^(s)k ds:

g0 (x) = Lky (x) y^(x)k  Lg (x):

45

Dado que g(x) > 0 para todo x  xn+1 , tenemos la desigualdad g0 =g  L,
que integrada en (xn+1 ; x) produce g(x)  g(xn+1 )eL(x xn ) . Usando una
vez mas que ky (x) y^(x)k  g(x), y pasando al lmite  ! 0+ llegamos a
+1

ky(x) y^(x)k  ky(xn+1) yn+1keL(x

xn+1 ) :

(1.34)

Observacion. Esta desigualdad prueba que el problema de valor inicial (1.1)

tiene solucion unica, y que esta dependen de manera continua del dato inicial.
Observacion. El paso de la estimacion (1.33) a la estimacion puntual (1.34)
es lo que se conoce en la literatura como Lema de Gronwall, [QV].

Procediendo como antes, pero pasando y (xn+1)

y sustituyendo x = xn+2 ,
R

yn+1 al primer miembro

ky(xn+2) y^(xn+2) (y(xn+1) yn+1)k  L xxnn ky(s) y^(s)k ds

R
 Lky(xn+1) yn+1k xxnn eL(s xn ) ds = ky(xn+1) yn+1k(eLh
+2

+1

+2

+1

+1

1)

= O(hp+2):
As, usando (1.29) tenemos que

y (xn+2 ) y^(xn+2 ) = (y (xn))hp+1 + O(hp+2):

(1.35)

Para estimar el segundo sumando observamos que

y^(xn+2 ) yn+2 = (^y(xn+1 ))hp+1 + O(hp+2):

Dado que y^(xn+1 ) = yn+1 = y (xn ) + O(h), se tiene que

y^(xn+2 ) yn+2 = (y (xn))hp+1 + O(hp+2):

Sumando (1.35) y (1.36) se obtiene nalmente (1.30).
46

(1.36)

1.13 Cambio de paso

Queremos escribir un codigo que ajuste automaticamente la longitud del paso
de forma que el error local se ajuste a una tolerancia prescrita.
Una vez elegida una longitud de paso inicial h, el programa calcula dos
pasos de tama~no h y un paso de tama~no 2h. Se estima entonces el error en
norma innito usando (1.31) en la forma
1
ERROR = p
max jy i
wi j:
2 1 i=1;:::;m n+2
Tambien son de uso frecuente otras normas, como por ejemplo la eucldea.
Una vez estimado el error se compara con la tolerancia prescrita, TOL.
Entonces, si ERROR  TOL, se aceptan los dos pasos calculados y la solucion
se avanza con yn+2 , o bien con yn+2 (en este caso se dice que se esta haciendo
extrapolacion local). Se intenta ahora un nuevo paso con una longitud de paso
h0 . Dicha longitud se intenta escoger de forma optima, es decir, la mayor que
produzca un ERROR menor que TOL.
Queremos que
Pero

2 (y (xn+2))(h0 )p+1 = TOL:

 0 p+1
 0 p+1
h
0
p
+1
p
+1
2 (y (xn+2))(h )  2 (y (xn))h
 ERROR h
:

As,

h0 = h
es la eleccion \optima".

TOL
ERROR

p+1

(1.37)

Para tener un buen codigo hay que poner un poco mas de cuidado. Se
suele multiplicar el segundo miembro de (1.37) por un factor de seguridad
FAC, normalmente FAC = 0:9, de forma que el ERROR sea aceptable la siguiente
vez con probabilidad alta. Tampoco se suelen permitir incrementos de paso
47

muy grandes, con el n de hacer el metodo mas seguro. As, por ejemplo
podemos poner


TOL
h0 = h min FACMAX; FAC
ERROR

p+1

FACMAX se suele tomar entre 1:5 y 5.

Si ERROR > TOL, los dos pasos se rechazan y se vuelven a repetir los calculos
con el nuevo tama~no de paso h0 .

1.14 Pares encajados

La estimacion de Richardson funciona bien, pero es computacionalmente cara;
si un metodo de Runge-Kutta explcito tiene s etapas, se necesitan otras s 1
evaluaciones de funcion (la primera etapa en el calculo de w ya se conoca)
para estimar el error.
Se utiliza una idea alternativa, debida a Merson, que pasamos a explicar.
Sea un Runge-Kutta que cumple las condiciones de orden p

A
bT .

Si tuvieramos otro metodo de Runge-Kutta con las mismas etapas,

c A
^bT
pero de distinto orden, p^ > p, podemos calcular
P

yn+1 = yn + h si=1 bi ki ;
P
y^n+1 = yn + h si=1 ^bi ki :
No hay que hacer nuevas evaluaciones de funcion, pues las etapas ki coinciden.
48

Como siempre, supondremos que se cumple la hipotesis de localizacion,

y (xn ) = yn. Sabemos que los errores locales satisfacen

Rn = y (xn+1) yn+1 + O(hp+2);

R^n = y (xn+1) y^n+1 + O(hp^+2 ):
Eliminando y (xn+1) y usando que R^ n = O(hp^+1 ) y que p + 2  p^ + 1, se tiene
que
Rn = y^n+1 yn+1 + O(hp+2):
Esto justica la estimacion

Rn  y^n+1

yn+1 = h

s
X
i=1

(^bi

bi )ki :

(1.38)

As, podemos tomar

ERROR = max jy^ni +1

yni +1 j;

i=1;:::;m

y se sigue como en el caso de la extrapolacion de Richardson.

Veamos dos ejemplos debidos a Fehlberg. El primero es el RKF2(3), de
tablero
0
1
1
2

yn+1
y^n+1

0
1
1
4
1
2
1
6

1
4
1
2
1
6

0
4
6:

El metodo de orden 2 que se usa para avanzar es en realidad de dos etapas,

mientras que el de orden 3 que se usa para estimar es de tres etapas. El
segundo par encajado que queremos presentar es el RKF2(3)B,
49

1
4
27
40

yn+1
y^n+1

1
4
189
800
214
891
214
891
533
2106

729
800
1
33
1
33

650
891
650
891
800
1053

1
78 .

El metodo de avance es de 3 etapas y orden 2, y el de estimacion es de 4 etapas

y orden 3.
Notacion. \Nombre p(q )" quiere decir que el orden de yn+1 (el metodo que
se utiliza para avanzar) es p, y el orden de y^n+1 (el metodo que se usa para
estimar) es q .

Comparando los metodos RKF2(3) y RKF2(3)B podramos pensar que

el primero es mucho mejor, por emplear una etapa menos para conseguir el
mismo orden. Sin embargo, el metodo RKF2(3)B tiene la propiedad FSAL
(First Same As Last),
asi = bi
i = 1; : : : ; s:
Por tanto, recordando que es un metodo explcito,

yn+1 = yn + h

s 1
X
i=1

bi ki = yn + h

y se tiene que

ks = f (t + h; yn + h

s 1
X
i=1

s 1
X
i=1

asi ki;

asi ki )

coincide con la primera etapa del paso siguiente,

k1 = f (t + h; yn+1):
Es decir, la primera etapa de cada paso coincide con la ultima del paso anterior.
Por tanto la primera etapa solo hay que calcularla en el primer paso, y desde el
punto de vista computacional es como si el metodo tuviera una etapa menos.
50

Aunque la estimacion (1.38) da el error de truncacion para el metodo del

par encajado de orden mas bajo, si h es sucientemente peque~no el error para
el metodo de orden alto sera aun menor; as que, >por que no avanzar con el
metodo de orden mas alto? Esto es lo que hacen los metodos de Dormand
y Prince DOPRI5(4) y DOPRI8(7), que son algunos de los mejores metodos
existentes hasta el momento.
Habra que tomar una precaucion. Queremos que el error de cometido al
avanzar, dado aproximadamente por (1.28), sea lo menor posible. Por tanto,
si se avanza con el metodo de orden alto, en el dise~no del algoritmo conviene
minimizar los coecientes de error de dicho metodo, en lugar de los del metodo
de orden bajo, que al n y al cabo solo se usa para controlar el paso.

1.15 Diagramas de eciencia

>Como podemos comparar la eciencia de distintos metodos? Una manera muy
conveniente es elaborar un diagrama de eciencia. En estos diagramas se representa, en escalas logartmicas, la precision conseguida, medida, por ejemplo,
con la norma innito del error en el punto nal del intervalo de integracion,
frente al coste computacional, medido en evaluaciones de funcion, o en tiempo
de CPU. En el caso de metodos de paso jo podemos obtener diferentes puntos
del diagrama cambiando el paso h. En el caso de metodos de paso variable lo
que haremos para obtener distintos puntos es cambiar la tolerancia. Veamos
algunos ejemplos.
En la Figura 1.1 representamos errornorma innito del error en b = 2
frente a fcnumero de evaluaciones de funcion, en escalas logartmicas para
el problema lineal escalar

y 0 (x) = y (x);

y (0) = 1:

Cuanta mayor precision, mas abajo estaremos en el diagrama, y cuanto

51

10

RungeKutta clsico
HIHA5(4)
DOPRI5(4)
4

10

10

error

10

10

10

12

10

14

10

16

10

10

10

Figura 1.1

10
fc

10

10

Diagrama de eciencia para un problema lineal.

mas trabajo realicemos, mas a la derecha.

El diagrama muestra que, salvo para precisiones muy exigentes, proximas
a 10 14 (donde el error de redondeo empieza a desempe~nar un papel), a mayor precision mayor trabajo. Tambien queda claro que para este problema el
metodo mas eciente es DOPRI5(4), si bien HIHA5(4) de Higham y Hall no
le va a la zaga. El metodo de Runge-Kutta clasico es obviamente peor, siendo
la diferencia mas grande cuanto mayor sea la precision requerida.
Los puntos parecen disponerse (salvo en casos extremos) a lo largo de rectas.
Veamos una explicacion cuando el metodo es de paso jo. En este caso
error = ky (xN )

yN k  Chp =

Usando que fc= Ns (o bien fc= N (s

1) si se tiene FSAL), obtenemos que

error = c fc p ;

52

C (b a)p
:
Np

10

RungeKutta clsico
HIHA5(4)
DOPRI5(4)
4

10

10

error

10

10

10

12

10

14

10

16

10

10

10

Figura 1.2

10
fc

10

10

Diagrama de eciencia para el Brusselator.

y tomando logaritmos
log10 (error) = p log10 (fc) + c:
La pendiente de la recta es p. Esto produce una forma de obtener
empricamente el orden de convergencia, incluso para metodos de paso variable. Por ejemplo, en la graca vemos que el metodo de Runge-Kutta clasico
tiene orden 4, mientras que DOPRI5(4) y HIHA5(4) tienen orden 5.
En la Figura 1.2 representamos error frente a fc, en escalas logartmicas
para el problema conocido como Brusselator, que modela ciertas reacciones
qumicas,
y10 = 1 + y12y2 4y1 ;
y1 (0) = 1:5;
y20 = 3y1 y12 y2 ;
y2 (0) = 3;
con intervalo de integracion 0  x  20.
Notese que el orden de convergencia emprico es 4 para Runge-Kutta, 6 para
53

HIHA5(4) y 7 para DOPRI5(4). Para todos los metodos, dada una precision,
necesitamos mas trabajo para resolver el Brusselator que el problema lineal.
Para metodos de paso jo hay otra forma de obtener empricamente el
orden de convergencia. Para ello integramos la ecuacion en [a; b] con paso h y
h=2. Tenemos
error(h) = ky (xN )

yh;N k  Chp ;
error(h=2) = ky (xN ) y h ;N k  Chp =2p:
2

Por tanto

error(h=2) 

error(h)

y concluimos que

p

log(error(h))

2p

log(error(h=2))
:
log 2

De hecho el lmite del segundo miembro cuando h ! 0+ es p.

1.16 Comentarios y bibliografa

Se puede encontrar una prueba del Teorema de Euler en [HNW]. El primer
captulo de dicho libro presenta la teora clasica de ecuaciones diferenciales
mostrando su evolucion historica. Contiene ademas mucha informacion sobre el
tema que no siempre se puede encontrar en otros libros. Por todo ello su lectura
resulta muy recomendable. Otras posibles referencias para la teora clasica son
los libros que el alumno habra manejado en la asignatura de Calculo IV, como
[BD] o [Si].
El metodo de Euler, predecesor de todos los demas metodos numericos para
ecuaciones diferenciales ordinarias, aparece en la mayora de los libros de texto
de introduccion al calculo numerico, por ejemplo en [A1], [A2], [JR] o [BF].
La regla del trapecio y los metodos de Taylor tambien se pueden encontrar en
cualquiera de esos textos.
54

El libro [JR] contiene una prueba de convergencia de metodos de un paso

parecida a la que damos nosotros, si bien restringida a metodos explcitos.
Contiene tambien una introduccion asequible a los metodos de Runge-Kutta.
Los metodos de Runge-Kutta aparecen por primera vez en un artculo de
Runge en 1895 [R]. Posteriormente, en 1900, Heun [Heu] introdujo nuevos
metodos. Fue Kutta en 1901 [Ku] quien formulo el esquema general de los
que hoy se conocen como metodos de Runge-Kutta. Para saber mas sobre la
historia del analisis numerico se puede consultar [Go].
Las primeras ideas de la teora de Butcher aparecen en un artculo de
Merson en 1957; su desarrollo completo es debido a John Butcher, en una
larga serie de artculos, el primero de los cuales es de 1963. En [Bu] se puede
encontrar una exposicion completa de la teora, con el encanto de que esta
escrita por su creador. El libro contiene ademas una extenssima bibliografa
sobre metodos numericos para ecuaciones diferenciales ordinarias. En [L] se
presenta una version simplicada (sin pruebas) de la teora. En [HNW] se
expone el tema de forma razonablemente clara y completa.
Hay una manera alternativa de obtener las condiciones de orden debida
a P. Albrecht. Usa algunas de las ideas que se emplean al estudiar metodos
lineales multipaso. Se puede encontrar una descripcion breve en [L].
La mejor referencia para cuestiones practicas (estimacion del error, seleccion de la longitud del paso, . . . ) es sin duda [HNW].
Existen metodos de Runge-Kutta especcos para ciertas clases de problemas, con la caracterstica de que las soluciones numericas retienen algunas
propiedades especialmente importantes de las soluciones del problema original.
Como ejemplos tenemos los metodos simplecticos para problemas hamiltonianos [SC], que preservan el hamiltoniano y los volumenes, o los metodos que
mantienen el caracter disipativo de un problema [DV].

55

Cap
tulo 2

Problemas sti

2.1 >Que es un problema sti?

Intentemos resolver el problema lineal

y 0 = Ay; y (0) = y0 ;

con A =

100
0

(2.1)

1
10

por el metodo de Euler. Se tiene que

yn+1 = (1 + hA)yn;
y por induccion que

yn = (I + hA)n y0 ;

n = 0; 1; 2; : : : :

Pasando a la base canonica de Jordan

A = V DV 1 ; donde V =

1 1
0 999
10
57

y D=

100
0

1
10

de donde se obtiene que la solucion exacta del problema viene dada por

y (x) = exA y0 = V exD V 1 y0 ; x  0; donde exD =

100x

10

En otras palabras, existen dos vectores v1 y v2 , dependientes de y0 , pero no de

x, tales que
y (x) = e 100x v1 + e x=10 v2 ;
x  0:
La funcion e 100x decae extraordinariamente mas aprisa que e x=10 . Por tanto,
incluso para x > 0 peque~no, la contribucion de v1 es despreciable a todos los
efectos, y y (x)  e x=10 v2 . >Que podemos decir de la solucion numerica fyng
obtenida con el metodo de Euler? Tenemos que

yn = V (I + hA)n V 1 y0 ;
y, dado que
(1 + hD)n =
se sigue que

yn = (1

(1 100h)n
0
0
(1 10h )n


100h)nv

n = 0; 1; : : : ;

+ 1

h
10

n

v2 ;

n = 0; 1; : : : :

Si h > 1=50, entonces j1 100hj > 1, y para n grande los iterados de Euler crecen en magnitud geometricamente, en contraste con el comportamiento
asintotico de la verdadera solucion.
Supongamos que escogemos una condicion inicial identica a un autovector
correspondiente al autovalor 1=10, por ejemplo

y0 =

999
10

(2.2)

En aritmetica exacta v1 = 0, v2 = y0 e yn = (1 10h )n y0 . La solucion numerica

converge al vector cero cuando n ! 1 para valores razonables de h, concretamente para h < 20. As que podramos esperar que todo fuera bien con
58

16

14

log kyn k

12

10

Figura 2.1

10

15

20

25

Logaritmo de la norma eucldea kyn k de los pasos de Euler, cuando se aplica

el metodo al problema (2.1) con h = 1=10 y condicion inicial igual al autovector \estable".

el metodo de Euler. No es as. Los ordenadores reales producen errores de

redondeo y antes o despues produciran una contribucion distinta de cero a un
autovector correspondiente al autovalor 100. En cuanto esto ocurre la componente \inestable" crece geometricamente como (1 100h)n , y rapidamente,
a menos que h < 1=50, eclipsa a la verdadera solucion.
En la gura 2.1 se representa log kyn k frente a n, para n = 0; : : : ; 25 con
la condicion inicial (2.2) y paso h = 1=10. El calculo se ha realizado con
Matlab. La norma de los diecinueve primeros paso decrece con la velocidad
correcta, dictada por (1 10h )n = (99=100)n. Sin embargo, a partir de ese
momento las cosas se estropean, y despues de un par de pasos la norma crece
geometricamente, como j1 100hjn = 9n . Se puede comprobar en la gura
que, en efecto, la pendiente de la curva es inicialmente log(99=100)  0:0101,
pero pasa a ser log 9  2:1972 cuando se llega al regimen inestable.
La anterior eleccion del dato inicial no es articial. Si nos enfrentamos a
un problema como (2.1) con una condicion inicial arbitraria tendremos que
59

emplear un tama~no de paso peque~no en el regimen transitorio inicial en el

que la contribucion del autovector \inestable" todava es signicativa. Una
vez pasado este corto periodo, la solucion viene esencialmente descrita por el
autovector \estable", y uno se ve tentado a incrementar h. Sin embargo hay
que resistirse a esta tentacion pues, como hemos visto antes, aunque no lo
parezca el autovector \inestable" siempre deja sentir su presencia.
Es importante comprender que este comportamiento no tiene nada que ver
con el error local del metodo numerico; hay que disminuir el paso, no por
consideraciones de precision, sino por la inestabilidad.
No todos los metodos numerico muestran tanta sensibilidad a la estabilidad.
Por ejemplo, si resolvemos (2.1) con la regla del trapecio
h
yn+1 = yn + (f (xn ; yn) + f (xn+1 ; yn+1 ))
2
obtenemos

 1

h
h
yn+1 = I
A
I + A y0 :
2
2
Dado que (I h2 A) 1 y (I + h2 A) conmutan, el orden de multiplicacion no
importa, y se tiene que

yn =

h
A
2

 1

h
I+ A
2

!n

y0 ;

n = 0; 1 : : : :

Usando como antes la forma canonica de Jordan de A y factorizando deducimos

que
!n


1 20h
1 50h n
yn =
v2 ;
n = 0; 1; : : : :
v1 +
1 + 50h
1 + 20h
Dado que


 
h
 1 50h   1 20



 1 + 50h  ;  1 + h  < 1
20
para todo h > 0, deducimos que limn!1 yn = 0. Esto recupera el comportamiento asintotico correcto para la ecuacion (2.1) independientemente del
tama~no de h.
60

En otras palabras, la regla del trapecio no requiere ninguna restriccion

en el paso para evitar la inestabilidad. Esto no quiere decir, por supuesto,
que cualquier h valga. Es necesario escoger h > 0 sucientemente peque~no
como para que el error local sea razonablemente peque~no y estemos aproximando adecuadamente la verdadera solucion. Sin embargo, no hay necesidad
de disminuir h hasta un tama~no minusculo con el n de evitar el crecimiento
incontrolado de componentes espurias de la solucion.
La ecuacion (2.1) es un ejemplo de problema sti (rgido, en ingles). En la
literatura se pueden encontrar varias deniciones de este concepto, pero quizas
sea mas interesante adoptar una denicion operativa, si bien un poco vaga.

\Denicion" 2.1. Decimos que el sistema

y 0 = f (x; y );

x  x0 ;

y (x0 ) = y0 ;

es sti si su resolucion numerica por algunos metodos exige una disminucion

signicativa del paso de integracion para evitar inestabilidades.

Es evidente que esta no es una denicion matematica correcta; pero tampoco estamos intentando probar teoremas del tipo si un sistema es sti entonces. . . . Lo importante de este concepto es que nos debe ayudar a elegir y
dise~nar metodos numericos.
En nuestro ejemplo hemos visto uno de los mecanismos (de hecho, uno de
los mas frecuentes) que hacen que un problema sea sti, el que haya modos
con escalas y tiempos de vida muy diferentes en la solucion. Es frecuente
designar al cociente de los modulos de los autovalores de mayor y menor modulo
de un sistema lineal (en el caso de un sistema general, los autovalores de la
matriz jacobiana) como coeciente de rigidez del sistema. Si dicho coeciente
es grande, el sistema suele ser sti. No obstante, haremos notar que este
analisis lineal puede fallar al aplicarlo a un sistema no lineal.
Una gran parte de los sistemas de ecuaciones diferenciales ordinarias que
aparecen en problemas reales son sti. Siempre que las ecuaciones modelen
61

varios procesos con velocidades de evolucion muy distintas, el problema sera

probablemente sti. Esto sucede por ejemplo en problemas de cinetica qumica,
combustion, prediccion meteorologica, biologa y electronica. Otra fuente muy
importante de problemas sti es el propio analisis numerico. Por ejemplo, las
ecuaciones en derivadas parciales parabolicas se aproximan con frecuencia con
sistemas de ecuaciones diferenciales ordinarias que resultan ser sti.

2.2 Dominio de estabilidad lineal y A-estabilidad

Supongamos que aplicamos un metodo numerico con paso constante h > 0 a
la ecuacion lineal escalar

y 0 = y;

x  0;

y (0) = 1;

(2.3)

donde  2 C . La solucion exacta es y (x) = ex , y por tanto limx!1 y (x) = 0

si y solo si Re  < 0.

Denicion 2.2. El dominio de estabilidad lineal (tambien llamado region

de estabilidad absoluta) de un metodo numerico es el conjunto de todos los
numeros z = h con h > 0,  2 C , tales que la solucion numerica de (2:3)
verica que nlim
y = 0.
!1 n
Es decir, el dominio de estabilidad lineal es el conjunto de los h para los
que se recupera el comportamiento asintotico de la solucion exacta cuando esta
es estable. Lo denotaremos por D.
Ejemplo. Al aplicar el metodo de Euler a (2.3) obtenemos la sucesion

yn = (1 + h)n ;

n = 0; 1; : : : ;

que es geometrica. As pues, limn!1 yn = 0 si y solo si j1+hj < 1. Concluimos

que D = fz 2 C : j1 + z j < 1g.
62

Ejemplo. Para la regla del trapecio tenemos que



1 + 12 h
yn =
1 21 h

n

n = 0; 1; : : : ;

que es una progresion geometrica. El dominio de estabilidad lineal es obviamente





 1 + 12 z 


D = z 2 C : 1 1z  < 1 :
2
Es facil vericar entonces que D = fz 2 C : Re z < 0g.
>Por que estamos subitamente interesados en una humilde ecuacion escalar?
Supongamos que queremos resolver

y 0 = Ay;

y (0) = y0 ;

con A una matriz d  d arbitraria por Euler. La solucion numerica es yn =

(1 + hA)n y0 . Si existe una base de autovectores de A, entonces A = V DV 1
con D = diag (1 ; 2 ; : : : ; d ), y se tiene que yn = V (I + hD)n V 1 y0 . Por tanto
existen vectores v1 ; v2 ; : : : ; vd 2 C d tales que

yn =

d
X
k=1

(1 + hk )n vk ;

n = 0; 1; : : : :

Supongamos que la solucion exacta es asintoticamente estable. Esto solo ocurre

si Re k < 0 para todo k = 1; 2; : : : ; d. Para que la solucion numerica de Euler
imite este comportamiento se necesita que el paso h > 0 sea tal que

j1 + hk j < 1;

k = 1; : : : ; d:

Es decir, todos los productos h1 ; : : : ; hd deben estar en el dominio de estabilidad lineal del metodo de Euler. Esto signica en la practica que el tama~no
del paso esta determinado por aquel autovalor que sea peor desde el punto de
vista de la estabilidad.
La restriccion de que A tenga una base de autovectores se puede relajar,
ya que la forma canonica de Jordan es suciente.
63

El analisis se puede extender tambien facilmente a sistemas homogeneos,

y 0 = Ay + a.

La region de estabilidad lineal D es tambien importante en sistemas no

lineales. Dado un sistema no lineal

x  x0 ;

y 0 = f (x; y );

y (x0 ) = y0 ;

se suele pedir que en el paso n-esimo

hn;1 ; hn;2 ; : : : ; hn;m 2 D;

donde n;1 ; : : : ; n;m son los autovalores de la matriz jacobiana

Jn 

@f (xn ; yn)
:
@y

Esto se basa en la suposicion de que localmente el comportamiento del sistema

esta bien modelado por la ecuacion variacional

y 0 = yn + Jn (y

yn ):

Hay que advertir que esto dista mucho de ser exacto.

Cuanto mayor sea la region de C := fz 2 C : Re z < 0g cubierta por D,
menos restricciones para h debidas a posibles inestabilidades (recuerdense los
ejemplos de Euler y la regla del trapecio). Eso motiva la siguiente denicion.

Denicion 2.3. Un metodo es A-estable si

C

 D:

Es decir, si un metodo es A-estable podemos elegir la longitud de paso h (al

menos para sistemas lineales) solo por motivos de precision, sin restricciones
debidas a la inestabilidad. Estos metodos seran buenos para problemas sti.
64

2.3 Estabilidad de metodos de Runge-Kutta

Queremos aplicar el metodo de Runge-Kutta (1.9) al problema escalar

x  0;

y 0 = y;

y (0) = 1:

(2.4)

Para ello conviene reescribir (1.9) en la forma alternativa

8
>
>
<

Y i = yn + h

>
>
: yn+1

Obtenemos que

8
>
>
<

s
P

aij f (xn + cj h; Yj );
j =1
s
P
yn + h bi f (xn + ci h; Yi):
i=1

Yi = yn + h

>
>
: yn+1

Sean Y; e 2 R s dados por

s
P

aij Yj ;
j =1
s
P
yn + h bi Yi :
i=1

Y = (Y1 ; : : : ; Ys)T ;

(2.5)

e = (1; 1; : : : ; 1)T :

Entonces podemos escribir (2.5) en la forma

(

Y = e yn + hAY;
yn+1 = yn + hbT Y:

La primera de estas ecuaciones da

Y = (I

hA) 1 eyn;

que introducido en la segunda produce

yn+1 = (1 + hbT (I
es decir,

hA) 1 e)yn ;

yn+1 = R(h)yn;
65

n = 0; 1; : : : ;

donde

R(z ) = 1 + zbT (I

zA) 1 e;

z 2 C:

(2.6)

La funcion R(z ) recibe el nombre de funcion de estabilidad del metodo, o tambien funcion de amplicacion.
Observese que yn = (R(h))n , as que de la denicion de
mediatamente que
D = fz 2 C : jR(z)j < 1g:

D se sigue in-

Por tanto un metodo sera A-estable si y solo si jR(z )j < 1 en C .

Notacion. P es el conjunto de todos los polinomios de grado menor o igual
que . P= es el conjunto de todas las funciones racionales p^=q^, con p^ 2 P y
q^ 2 P .

Un sencillo analisis permite probar el siguiente lema.

Lema 2.4. Para todo metodo de Runge-Kutta existe R 2 Ps=s tal que la solucion por dicho metodo del problema (2:4) viene dada por
yn = (R(h))n :
Es mas, si el metodo de Runge-Kutta es explcito, entonces R 2 Ps .
Prueba. Tenemos que probar que la funcion R dada por (2.6) pertenece a
Ps=s . Usamos que
(adj (I zA))T
(I zA) 1 =
:
det(I zA)
Cada componente de I zA es lineal en z . Por tanto, como cada elemento de
adj (I zA) es el determinante de una matriz (s 1)  (s 1), esta en Ps 1 .
As pues,
bT (adj (I zA))T e 2 Ps 1 :

Como det(I

zA) 2 Ps , concluimos que R 2 Ps=s .

66

Si el metodo es explcito, entonces A es estrictamente triangular inferior y

I zA es una matriz triangular inferior con unos a lo largo de la diagonal. Por
tanto det(I zA) = 1, y R es un polinomio.

Corolario 2.5. Ningun metodo de Runge-Kutta explcito puede ser A-estable.

Prueba. La funcion de estabilidad de un metodo de Runge-Kutta explcito es
un polinomio. Ademas R(0) = 1. Por otra parte, ningun polinomio, salvo la
funcion constante R(z ) = c 2 ( 1; 1) esta estrictamente acotado por uno en
C .

En el caso de metodos explcitos podemos decir aun mas sobre la funcion

de amplicacion.

Teorema 2.6. Si un metodo de Runge-Kutta explcito tiene orden de consistencia p, entonces

R (z ) = 1 + z +

z2
zp
+


+ + O(z p+1 ):
2!
p!

Prueba. La solucion exacta del (2.4) es y (x) = ex . Por ser el metodo
consistente de orden p tenemos que

y (h) y1 = O(hp+1);
as que

z2
zp
+


+ + O(z p+1);
2!
p!
donde z = h. Usando que y1 = R(z )y0 concluimos inmediatamente el resultado.
y1 = ez + O(hp+1 ) = 1 + z +

Como consecuencia tenemos que todos los metodos de Runge-Kutta explcitos con p = s tienen como funcion de estabilidad

R(z ) = 1 + z +

z2
zs
+


+ :
2!
s!

67

3i
s=4

s=3
s=2

s=1

-1

3i

Figura 2.2

Dominio de estabilidad lineal de los metodos de Runge-Kutta explcitos de

orden p = s.

En la gura 2.2 mostramos los correspondientes dominios de estabilidad lineal.

>Que podemos decir de los metodos implcitos? Veamos con un ejemplo
que pueden ser A-estables.
Ejemplo. Consideramos el metodo

0
2
3

1
4
1
4
1
4

1
4

5
12
3.
4

Su funcion de amplicacion es

R(z ) =

1 + 13 z
:
1 32 z + 16 z 2

Veamos si es A-estable. Representamos z 2 C en coordenadas polares, z =

ei , con  > 0, j  j < =2, y nos preguntamos si jR(ei )j < 1. Esto es
equivalente a
2 
2





1 + 1 ei  < 1 2 ei + 1 2 e2i  ;


3   3
6
68

y por tanto a
2
1
1 +  cos  + 2 < 1
3
9

4
1
4
 cos  + 2 cos 2 +
3
3
9

2 3
1
 cos  + 4 :
9
36

Reordenando terminos, la condicion ei 2 D pasa a ser



1
1
1
2
1
2 1 + 2 cos  < 2 (1 + cos 2) + 4 = 2 cos2  + 4 ;
9
3
36
3
36
que se satisface para todos los z 2 C , puesto que cos  < 0 para todo z de
este tipo. El metodo es por tanto A-estable.
Esto explica por que son interesantes los metodos implcitos. Son mas
caros, pero, ademas de permitir romper las barreras de Butcher, algunos de
ellos tienen regiones de estabilidad grandes, y pueden ser convenientes a la
hora de tratar problemas sti.
Sin embargo, no todos los metodos implcitos son A-estables.
Ejemplo. Consideramos el metodo

0 0 0
2
3

1
3
3
4

1
3
1.
4

Su funcion de amplicacion es

R(z ) =

1 + 23 z + 16 z 2
:
1 13 z

As pues, lim jR(z )j = 1, por lo que el metodo no es A-estable.

jzj!1

En realidad no hace falta comprobar todos los z 2 C para ver si una

funcion racional dada tiene su origen en un metodo A-estable (una R(z ) tal se
 es el contenido del siguiente lema, cuya prueba es una
dice A-aceptable). Ese
aplicacion del principio del maximo.
69

Lema 2.7. Sea R(z ) una funcion racional no constante. Entonces jR(z )j < 1
para todo z 2 C si y solo si todos los polos de R tienen parte real positiva y
jR(it)j  1 para todo t 2 R .
Prueba. Si jR(z )j < 1 para todo z 2 C , por continuidad jR(z )j  1 para
todo z 2 C . En particular R no puede tener polos en el semiplano cerrado
izquierdo, y jR(it)j  1 para todo t 2 R .

Para probar el recproco observamos que si los polos estan a la derecha del
eje imaginario, entonces la funcion racional R es analtica en el cerrado C .
Por tanto, al no ser R constante, jRj no puede tener un maximo en el interior,
y alcanza su maximo en la frontera. Parte de la frontera esta en el innito.
Sin embargo, por ser R una funcion racional existe el lmite (nito o innito)
limjzj!1 jR(z )j, y es igual a limt!1 jR(it)j. Por tanto, si jR(it)j  1, t 2 R ,
entonces jR(z )j < 1 para todo z 2 C .
Ejemplo. Consideramos el metodo de Gauss-Legendre de dos etapas (orden 4),
1
2
1
2

p3
p63
6

1
4

1
4p

3
6

1
2

1
4

p3

1
4
1.
2

Un sencillo calculo muestra que

1 + 12 z + 121 z 2
:
1 21 z + 121 z 2
p
Los polos de esta funcion, 3  i 3, estan en el semiplano abierto derecho, y
jR(it)j = 1, t 2 R , por lo que el metodo es A-estable.

R(z ) =

2.4 Comentarios y bibliografa

El captulo dedicado en [L] a los problemas sti y a la estabilidad lineal nos
parece excelente. Mas simple, pero aun mas claro, nos parece [I]. La referencia
mas completa sobre este tema es [HW].
70

Los metodos de Runge-Kutta implcitos, de los que no hemos dicho nada,

pueden ser interesantes a la hora de tratar problemas sti. Se puede encontrar
una breve introduccion en las secciones 3.3 y 3.4 de [I], donde ademas se explica
su conexion con los llamados metodos de colocacion. Si se quiere saber mas,
habra que consultar [HNW] y [HW].
Si se usan metodos implcitos (de Runge-Kutta o lineales multipaso), sera
necesario resolver sistemas no lineales de ecuaciones algebraicas. Se puede
encontrar un tratamiento exhaustivo de los metodos con que se resuelven estos
sistemas en el clasico [OR], as como en [F] y el mas moderno [K].
La teora de estabilidad no lineal es relativamente reciente. Se puede encontrar una introduccion en [L] y un estudio mas riguroso y completo en [DV].

71

Cap
tulo 3

Diferencias nitas

3.1 La formula de cinco puntos

En este captulo y el siguiente mostraremos algunos metodos para resolver
numericamente ecuaciones en derivadas parciales. Consideraremos tan solo
problemas elpticos de segundo orden. No obstante, algunas de las ideas se
pueden utilizar provechosamente para resolver otros problemas.
Nuestra primera familia de metodos usa un enfoque que ya hemos empleado con exito para resolver ecuaciones diferenciales ordinarias: sustituir
las derivadas por cocientes incrementales. Esto dara lugar a los metodos de
diferencias nitas.
Sea
 R 2 un abierto acotado del plano, tal que todo punto de la frontera
sea accesible desde el exterior por una curva continua simple. Consideramos
la ecuacion de Poisson con condiciones de frontera de tipo Dirichlet

u = f

en
;

u=g

en @
:

(3.1)

Supondremos que g es continua en

y que existen numeros C > 0 y 0 <  < 1
73

tales que jf (x) f (y )j  C jx y j para todo x; y 2

(diremos que f es Holder continua en
, y escribiremos que f 2 C 0; (
)). Se tiene el siguiente
resultado de existencia y unicidad.

Teorema 3.1. Si f 2 C 0; (
) y g 2 C (
), el problema (3:1) tiene una unica
T
solucion u 2 C 2 (
) C (
).
Una herramienta fundamental en el analisis de esta ecuacion es el principio
del maximo.
T

Teorema 3.2 (Principio del maximo). Sea u 2 C 2 (

) C (
).
(i) Si
u  0 para todo (x; y ) 2
, entonces max u  max u.

(ii) Si
u  0 para todo (x; y ) 2

, entonces min u  min u.
Supongamos que el dominio sobre el que se formula la ecuacion (3.1) es el
rectangulo

= [a1 ; a2 ]  [b1 ; b2 ]:
Se puede ajustar una red uniforme formada por los puntos de coordenadas
(xi ; yj ) = (a1 + ih; b1 + jk)
para 0  i  N y 0
determinados por

 j  M,
h=

a2

a1

donde los tama~nos de paso h y k estan

k=

b2

b1

:
N
M
Sobre los nodos de la red denimos una funcion U cuyo valor U (xi ; yj ) = Uij
es una aproximacion de la solucion u en (xi ; yj ),
Uij  u(xi ; yj ):
Nuestro objetivo es encontrar una forma de generar las aproximaciones
Uij . Para ello utilizaremos formulas de diferencias nitas que aproximen las
derivadas a partir de valores de la funcion.
74

Nota. El termino diferencias nitas es clasico, y se remonta a los tiempos

en que, siguiendo a Leibniz (1646{1716), se conceba la diferencial df de una
funcion f como diferencia \innitamente peque~na" f (x + dx) f (x) entre los
valores de f en dos puntos \innitamente" proximos x y x + dx. Entonces se
contrapona a la diferencial innitesimal la diferencia nita f (x + h) f (x).

Sea f 2 C 4 ([x h; x + h]). Queremos obtener una expresion aproximada para la derivada segunda combinando valores de f en distintos puntos.
Desarrollando por Taylor en torno a x se tiene que

f (x  h) = f (x)  hf 0 (x) +

h3
h4
h2 00
f (x)  f 000 (x) + f (4) ();
2
3!
4!

donde + y  son puntos intermedios. De esta forma

f 00 (x) =

f (x + h) 2f (x) + f (x h)
h2

h2 (4)
(f (+) + f (4) ( )):
4!

Esto permite aproximar la derivada segunda por la formula de diferencias

centrales
f (x + h) 2f (x) + f (x h)
:
f 00 (x) 
h2
Si u 2 C 4 (
), podemos aplicar la formula de diferencias centrales a u(x; y )
para aproximar uxx (x; y ) y uyy (x; y ), obteniendo

u(x; y ) =
hk u(x; y ) + O(k2 + h2 );
donde
hk es el operador Laplaciano discreto,

hk u(x; y ) =

u(x + h; y ) 2u(x; y ) + u(x h; y )

h2
u(x; y + k) 2u(x; y ) + u(x; y k)
+
;
k2

tambien conocido como formula de cinco puntos. Se representa gracamente

mediante la molecula computacional
75

k2

h2

h2

k2

h2

k2

Si u es la solucion de (3.1), en cada uno de los nodos interiores tenemos

f (xi ; yj ) 
hk u(xi ; yj ) 
hk Uij :

Deniendo fij = f (xi ; yj ) y gij = g (xi ; yj ), lo que haremos sera pedir que U
satisfaga el problema de Dirichlet discreto para la ecuacion de Poisson
(

hk Uij = fij

Uij = gij

en los nodos interiores;

en los nodos de la frontera:

(3.2)

Por tanto, para un nodo interior (xi ; yj ), 0 < i < N , 0 < j < M , tenemos

Ui+1;j + 2Uij
h2

Ui

1;j

Ui;j +1 + 2Ui;j
k2

Ui;j

= fij ;

relacion que no solo involucra a Uij , sino tambien a los valores en los nodos
adyacentes, a la izquierda y a la derecha, arriba y abajo. As, para calcular la
solucion aproximada es preciso resolver un sistema lineal de (N 1)  (M 1)
ecuaciones cuyas incognitas son todos los elementos de la matriz Uij , excepto
los que corresponden a los nodos de frontera, que ya se conocen.
El Laplaciano discreto reproduce algunas de las propiedades mas importantes del operador Laplaciano, por ejemplo el principio del maximo discreto.
76

Teorema 3.3 (Principio del maximo discreto).

(i) Si
hk Uij  0 en todos los nodos interiores, max
Uij = max Uij .
i;j
(xi ;yj )2@

(ii) Si
hk Uij  0 en todos los nodos interiores, min

U = min Uij .
i;j ij
(xi ;yj )2@

Prueba. (i) Si
hk Uij  0, entonces



2
2
1
1
+ 2 Uij  2 Ui+1;j + 2 Ui
2
h k
h
h
Se puede escribir esta desigualdad como

Uij  1 Ui+1;j + 2 Ui

1;j

1
1
Ui;j +1 + 2 Ui;j 1 :
2
k
k

1;j + 3 Ui;j +1 + 4 Ui;j 1 ;

P4
i=1 i = 1. Es decir, el

donde 0 < i < 1 para i = 1; : : : ; 4 y

valor de U en
un nodo interior es menor que una media aritmetica ponderada de los valores
en los cuatro nodos adyacentes. Por consiguiente el valor de U en el nodo
interior no puede ser mayor que en los adyacentes y por tanto no hay maximos
interiores.
(ii) Observamos que
hk ( Uij )  0, y aplicamos el apartado (i) a U .

Para calcular las aproximaciones Uij tenemos que resolver el problema (3.2). >Tiene dicho problema solucion unica?
La respuesta es armativa: si hubiera dos soluciones, al aplicar el Teorema 3.3 a la diferencia de ambas se concluye que esta es trivial.

 >Se parece la solucion numerica a la teorica cuando h; k ! 0?

Siguiendo lo que hacamos para ecuaciones diferenciales ordinarias, diremos
que el metodo converge si
lim max ju(xi ; yj )

h;k!0+ i;j

Uij j = 0:

Sea Rhk u(x; y ) el residuo o error de truncacion, es decir, lo que le falta a

la solucion teorica para satisfacer la ecuacion del metodo,

Rhk u(x; y ) =
hk u(x; y ) f (x; y ):

77

El error E = u U es solucion de
(

hk Eij = Rhk u(xi ; yj )

Eij = 0

en los nodos interiores;

en los nodos de frontera;

que es el problema (3.2) con fij = Rhk u(xi ; yj ) y gij = 0, un problema con
dato Dirichlet homogeneo.

Lema 3.4. La solucion Uij del problema (3:2) con gij = 0 satisface

jUij j  161 jf j1(a2 + b2 );

a = a2

a1 ; b = b2

b1 ;

(3.3)

donde jf j1 = max jfij j.

i;j

Prueba. Sean (c1 ; c2 ) las coordenadas del centro del rectangulo [a1 ; b1 ][a2 ; b2 ].
Consideramos la funcion auxiliar

w(x; y ) = (x c1 )2 + (y

c2 )2 ;

y denotamos wij = w(xi ; yj ). La funcion w cumple que

1
0  wij  (a2 + b2 ):
4
Por construccion, el Laplaciano discreto y el Laplaciano coinciden en la
clase de los polinomios de grado menor o igual que 2. Por consiguiente,

hk wij =
w(xi ; yj ) = 4;
y por tanto

hk (Uij + 14 jf j1wij ) = fij jf j1  0;

hk (Uij 41 jf j1wij ) = fij + jf j1  0;
Usando los principios del maximo y del mnimo discreto se obtiene

Uij + 14 jf j1  max
Uij

1
4

(xi ;yj )2@

Uij + 41 jf j1wij = 14 jf j1 max wij ;

jf j1  (ximax
Uij
;yj )2@

1
4

jf j1wij

78

1
4

(xi ;yj )2@

jf j1 (ximin
wij :
;yj )2@

De estas relaciones y de la acotacion de w (notese que jxi

jyi c2 j  b=2) se deduce que

Uij  14 jf j1 max wij 

(xi ;yj )2@

c1 j

 a=2

jf j1(a2 + b2 );
Uij  14 jf j1 min wij  161 jf j1(a2 + b2 );
(xi ;yj )2@

1
16

de donde se concluye (3.3).

Aplicando el lema al problema satisfecho por el error, obtenemos que
2

+b
jE j  a 16
jRhk uj1:
Ahora bien,

jRhk u(x; y)j = j

As pues,
y concluimos que

hk u(x; y ) +
u(x; y )j = O(h2 + k2 ):

jRhk uj1 = O(h2 + k2 );

jE j = O(h2 + k2 ):

Hemos probado de esta forma el siguiente resultado de convergencia.

Teorema 3.5. La solucion del problema de Dirichlet discreto para la ecuacion

de Poisson (3:2) converge a la solucion del correspondiente problema continuo (3:1).

 >Hay formas ecientes de resolver el sistema a que da lugar el metodo

numerico?

Lo primero que hay que hacer es numerar adecuadamente los valores Uij como incognitas del sistema. Es habitual utilizar el siguiente orden lexicograco,
(i1 ; j1 ) < (i2 ; j2 ) ,

8
>
>
< j1
>
>
:

79

< j2 ;
o
j1 = j2 e i1 < i2 :

La dimension del sistema es (N 1)(M 1)  (N 1)(M 1), y es

en general alta. No obstante, es facil ver que el sistema se puede expresar
en forma tridiagonal por bloques. Esto se puede aprovechar por un lado para
reducir el almacenamiento de datos de la matriz de coecientes, y por otro para
simplicar el calculo de la descomposicion LU de la descomposicion gaussiana
(ver [M], captulo 7). Sin embargo, para (3.2) resulta mas sencillo resolver el
sistema lineal mediante tecnicas iterativas.
El metodo de diferencias nitas que acabamos de explicar se puede aplicar
a otras ecuaciones, no necesariamente elpticas. As mismo, se pueden usar
distintas aproximaciones de diferencias nitas para las derivadas que aparezcan
en la ecuacion. De esta forma se obtienen, para un mismo problema, metodos
diversos.
>Que sucede si el dominio no es rectangular? La dicultad se presenta al
intentar incorporar las condiciones de contorno al sistema lineal. Esto se puede
conseguir haciendo modicaciones del esquema, (ver [M]); sin embargo, los
esquemas resultantes ya no son tan faciles de programar, ni son tan ecientes.
Es preferible utilizar una idea completamente diferente: los elementos nitos.

3.2 Comentarios y bibliografa

El lector que quiera saber mas sobre la teora analtica de la ecuacion de Poisson
puede consultar [P], o el libro clasico [GT]. Otra buena referencia es [H], que
tiene la ventaja de combinar aspectos analticos y numericos.
En cuanto al metodo de diferencias nitas, aunque solo lo hemos aplicado a
la ecuacion de Poisson, es util para tratar otros muchos problemas, por ejemplo
problemas de evolucion (incluyendo problemas de primer orden). Se pueden
encontrar introducciones sencillas en [I] y [M]. El segundo de estos libros tiene
otras dos ventajas para el alumno: esta en castellano y contiene bastantes
problemas resueltos. Como referencias mas completas y avanzadas tenemos
80

los libros ya clasicos de [RM] y [MG], y los un poco mas modernos [St] y [T].
El estudio detallado de los diferentes metodos de resolucion de los sistemas
lineales dispersos a que da lugar el metodo de diferencias nitas queda fuera
del alcance del curso. En la seccion 4.9 damos abundantes referencias sobre
este tema.

81

Cap
tulo 4

Elementos nitos

4.1 Introduccion
Presentaremos el metodo de los elementos nitos estudiando un problema concreto,
(

u + u = f
en
 R d ;
u=0
en @
:
Para simplicar todava mas, nos centraremos inicialmente en dimension d = 1.
Tenemos as el problema de contorno

u00 + u = f si a < x < b;

u(a) = u(b) = 0:

(4.1)

Nota. El problema con condiciones no homogeneas u(a) = , u(b) =  se

reduce a este mediante el cambio de incognita u~ = u u0 , siendo u0 cualquier
funcion regular tal que u0 (a) = , u0 (b) =  (tomar, por ejemplo, u0 afn). La
funcion u~ es solucion del problema homogeneo, eso s, con una nueva f .

Denicion 4.1. Sea f 2 C ([a; b]). Una solucion clasica del problema (4:1) es
una funcion u 2 C 2 ((a; b)) \ C ([a; b]) que verica (4:1) para todo x 2 [a; b].
83

Gracias a los resultados conocidos para ecuaciones diferenciales ordinarias

tenemos el siguiente teorema de existencia y unicidad.

Teorema 4.2. Dado f

2 C ([a; b]), existe una unica solucion clasica de (4:1).

>Como calcular la solucion numericamente? En este caso, como la geometra del dominio es sencilla, podemos utilizar un metodo de diferencias nitas. Por ejemplo, si consideramos el conjunto de puntos

xn = a + nh;

n = 0; 1; 2 : : : ; N = (b a)=h;

y aproximamos la derivada segunda usando la formula de diferencias centradas,

los valores un de las aproximaciones de la solucion en los nodos xn se obtienen
resolviendo el sistema
un+1 2un + un 1
+ un = fn si 0 < n < N;
u0 = uN = 0:
h2
Los valores un se pueden interpolar para obtener una aproximacion de u(x)
en puntos x que no pertenezcan al conjunto fxn g. Si, por ejemplo, usamos
interpolacion lineal, entonces obtenemos como aproximacion

uh(x) =
con

'n (x) =

8
>
>
>
>
>
>
>
>
>
<
>
>
>
>
>
>
>
>
>
:

N
X1
n=1

un 'n (x);

x xn 1
;
xn xn 1

x 2 [xn 1 ; xn ];

xn+1 x
;
xn+1 xn

x 2 [xn ; xn+1 ];

0;

en el resto:

(4.2)

Veamos lo que hemos conseguido desde una perspectiva distinta. Hemos

obtenido una aproximacion uh a la verdadera solucion dentro del espacio de
dimension nita


Vh = fv 2 C ([a; b]) : v (a) = v (b) = 0; v [xn

84

;xn ]

2 P1 ([xn 1 ; xn]); 1  n  N g:

 sera precisamente el punto de partida para el metodo de elementos

Este
nitos:

 Buscar una aproximacion uh a la solucion dentro de un espacio de dimen-

sion nita.

>Como hacerlo? Elegimos M funciones linealmente independientes que

satisfagan las condiciones de frontera, 'l (a) = 'l (b) = 0, l = 1; : : : ; M y
consideramos el espacio vectorial de dimension nita generado por ellas,

Vh := L(f'1 ; : : : ; 'M g):

Buscamos una aproximacion uh a u dentro de Vh ,

uh (x) =

M
X
l=1

l 'l (x);

a  x  b;

(4.3)

donde 1 ; : : : ; M son constantes reales. Notese que, dado que todos los elementos de Vh satisfacen las condiciones de frontera, uh tambien lo hace.
Elegiremos los numeros 1 ; : : : ; M de manera que uh
\lo mejor posible". >Que queremos decir con esto?

2 Vh aproxime a u

Para cualquier eleccion de 1 ; : : : ; M consideramos el defecto o residuo

dh(x) := u00h (x) + uh(x) f (x);

a < x < b:

Nota. Para que el defecto este bien denido necesitamos que uh 2 C 2 ((a; b)).
Por tanto pedimos en principio que 'l 2 C 2 ((a; b)) \ C ([a; b]), aunque mas
adelante rebajaremos esta hipotesis.

Si uh fuera la solucion de (4.1), el defecto sera identicamente cero. Por

tanto, cuanto mas proximo este dm a la funcion cero, mejor esperamos que sea
nuestro candidato a solucion.
Dado un producto escalar en C ([a; b]) (espacio que incluye al de las soluciones clasicas) se puede escribir C ([a; b]) = Vh  Vh?, lo que se traduce en
85

la descomposicion u = uh + u?h , y la formula para el defecto es simplemente

dh = (u?h )00 u?h . Como esperamos que Vh sea tan grande que contenga todas
las funciones de interes, es natural imponer, para minimizar errores que la
proyeccion de dh sobre Vh sea nula, esto es, dh 2 Vh?, o lo que es lo mismo

hdh; 'k i = 0;

k = 1; : : : ; M:

(4.4)

Estas condiciones se llaman ecuaciones de Galerkin. Tambien podemos reinterpretar lo anterior diciendo que consideramos que Vh? es el espacio de funciones
despreciables, y lo que imponemos es que el \operador de error" u?h 7! dh =
(u?h )00 u?h aplique Vh? en Vh? para que una aproximacion buena no de lugar
a un defecto grande.
Nuestro segundo principio es por tanto:

 Elegir la aproximacion de forma que el defecto sea ortogonal al espacio Vh.

En nuestro caso podemos tomar como producto escalar en C ([a; b]) el producto escalar de L2 ((a; b)),

hv; wi =

Z b
a

vw:

Por otra parte, usando la representacion (4.3) tenemos que

M
X

dh =

l=1

X
'00 + '
l l

l=1

f;

l l

que sustituido en (4.4) produce

M
X
l=1

Z b
a

( '00l 'k + 'l 'k ) =

Z b
a

f'k ;

para k = 1; : : : ; M . Es decir, tenemos el sistema lineal

M
X
l=1

akl l =

Z b
a

f'k ;
86

k = 1; : : : ; M;

(4.5)

donde

akl :=

Z b
a

( '00l 'k + 'l 'k );

k; l = 1; : : : ; M:

La resolucion del sistema nos produce los l y por tanto la aproximacion

 es el llamado metodo de Galerkin.
numerica uh . Este
Para la puesta en practica del metodo hay que evaluar M 2 integrales, bien
explcitamente, o bien por medio de formulas de cuadratura numerica, as
como resolver un sistema lineal de dimension M . Si M es un poco grande la
tarea computacional, as como las necesidades de almacenamiento en memoria,
pueden ser excesivas. Ah es donde entra el principio que distingue al metodo
de elementos nitos de otros metodos de Galerkin:

 Elegir las funciones 'k de manera que sop 'k \ sop 'l = ; para la mayor

parte de las elecciones de k; l = 1; : : : ; M .

De esta forma, habremos reducido ampliamente tanto el numero de integrales a calcular, como el numero de elementos no nulos de la matriz del sis sera dispersa, con el consiguiente ahorro de memoria y de numero
tema. Esta
de operaciones para resolver el sistema.
Una eleccion excelente en este sentido sera tomar

Vh = L(f'1 ; : : : ; 'N

g);

con 'n como en (4.2). De esta forma, sop 'k \ sop 'l = ; si jk lj > 1. Pero
hay un problema: estas funciones basicas 'n no pertenecen a C 2 ((a; b)), y de
hecho, ni siquiera a C 1 ((a; b)). Aun as, desearamos poder trabajar con este
espacio, pues la matriz de coecientes sera tridiagonal.
>Que hacer? Podemos integrar por partes las ecuaciones de Galerkin (4.5),
usando que las funciones 'l se anulan en los extremos del intervalo, para obtener
M
X
l=1

Z b
a

('0l '0k + 'l 'k ) =

87

Z b
a

f'k ;

(4.6)

para k = 1; : : : ; M . Es decir, tenemos el sistema lineal

M
X
l=1

donde

a^kl :=

a^kl l =

Z b
a

Z b
a

f'k ;

('0l '0k + 'l 'k );

k = 1; : : : ; M;
k; l = 1; : : : ; M:

>Que hemos ganado? La integracion por partes elimina una derivada, as
que ya no hay que pedir que las funciones basicas tengan dos derivadas. De
hecho, incluso pedir 'l 2 C 1 ((a; b)) es excesivo. Los coecientes a^kl tienen
sentido para funciones 'l mucho mas generales, por ejemplo para funciones
derivables a trozos como las denidas en (4.2).
La integracion por partes reduce los requisitos de regularidad sobre las 'l .
Cuanto menores sean dichos requisitos, mayor sera la clase de funciones base
que podemos considerar. Esto facilitara nuestra tarea de construir funciones
base que se anulen en casi todo el intervalo. Tenemos as un nuevo principio:

 Integrar por partes para reducir lo mas posible las exigencias de regula-

ridad sobre el espacio Vh.

4.2 Un poco de analisis funcional. Espacios de Sobolev

Llegados a este punto, nos gustara saber cual es la regularidad mnima necesaria para que las ecuaciones (4.6), o equivalentemente los coecientes a^kl ,
tengan sentido. Esto nos lleva a estudiar los llamados espacios de Sobolev.
Sea I = (a; b) un intervalo, no necesariamente acotado, y sea p 2 R con
1  p  1. Denotamos por Cc1 (I ) al conjunto de las funciones derivables en
I con soporte compacto en I .

Denicion 4.3. El espacio de Sobolev W 1;p(I ) se dene por

W 1;p(I )

= fu 2

Lp (I )

: 9g 2

Lp (I )

tal que

88

Z b
a

u'0 =

Z b
a

g'

8' 2 Cc1(I )g:

Notacion. Se suele escribir H 1 (I ) = W 1;2 (I ).

La funcion g desempe~na el papel de u0 cuando se integra por partes. As,

para cada u 2 W 1;p(I ) se escribe g = u0 , y se dice que g es la derivada debil
de u. De esta forma, W 1;p(I ) es el espacio de las funciones de Lp (I ) cuya
derivada (debil) esta en Lp (I ).
Observacion. Si u 2 C 1 (I ) \ Lp (I ) y u0 2 Lp (I ) (aqu u0 es la derivada usual
de u), entonces u 2 W 1;p(I ). Ademas la derivada usual de u coincide con la
derivada de u en el sentido de W 1;p. En particular, si I esta acotado, entonces
C 1 (I )  W 1;p(I ) para todo 1  p  1.
Ejemplo. Consideramos la funcion u(x) = jxj para x
' 2 Cc1 (I ) se tiene que
Z

u'0 =

Z 0

x'0 +

Z 1
0

x'0 =

Z 0

'+

Z 1
0

( 1; 1). Dada

';

donde la ultima igualdad proviene de integrar por partes teniendo en cuenta

que '(1) = 0. Por tanto,
Z

con

u'0 =

g';

si 1 < x < 0;
1
si 0 < x < 1:
Evidentemente g 2 Lp (I ) para todo 1  p  1. As u 2 W 1;p(I ) para todo
1  p  1, y la derivada debil de u es u0 = g . Sin embargo, la funcion
u(x) = jxj no es derivable en x = 0 en el sentido usual, y por tanto u 62 C 1 (I ).

g (x) =

De la misma forma se prueba que si el intervalo I esta acotado, entonces

toda funcion continua en I y derivable con continuidad a trozos en I pertenece
a W 1;p(I ) para todo 1  p  1.
Notaciones. El espacio W 1;p esta dotado de la norma

kukW ;p = (kukpLp + ku0kpLp )1=p :

1

89

En H 1 esta norma corresponde al producto escalar

hu; viH

= hu; v iL + hu0 ; v 0 iL :
2

Teorema 4.4. El espacio W 1;p es un espacio de Banach para 1  p  1.

Prueba. Sea fung una sucesion de Cauchy en W 1;p; entonces fung y fu0ng son
sucesiones de Cauchy en Lp . Por la completitud un ! u en Lp y u0n ! g en
Lp . Se tiene
Z b
Z b
0
un' =
u0n '
8' 2 Cc1(I );

y en el lmite

Z b

Z b

8' 2 Cc1 (I );
a
a
1
;p
0
Por tanto u 2 W , u = g y kun ukW ;p ! 0.
u'0 =

g'

Consecuentemente, H 1 (I ) es un espacio de Hilbert. Dicho espacio parece

ser el adecuado para plantear las ecuaciones (4.6), salvo porque sus elementos
no satisfacen en general las condiciones de frontera homogeneas. Las funciones
de W 1;p estan denidas en casi todo punto, y la frontera @I tiene medida cero,
as que en principio no tiene sentido pedir que las funciones de la base 'l , y
por tanto la solucion aproximada uh, valgan cero en @I . Nos tendremos que
contentar con que esa propiedad se verique en un sentido debil. Para ello
introducimos un nuevo espacio.

Denicion 4.5. Dado 1  p < 1, se designa por W01;p(I ) al cierre de Cc1 (I )

en W 1;p(I ). Se nota H01(I ) = W01;2 (I ).
El espacio W01;p esta dotado de la norma inducida por W 1;p y con esta
norma es un espacio de Banach; el espacio H01 esta dotado del producto escalar
inducido por H 1 , y con este producto escalar es un espacio de Hilbert.
Pedir que una funcion u este en W01;p(I ) es una forma debil de pedir que
u = 0 en @I , en el sentido de que podemos aproximar u tanto como queramos
(en la norma de W 1;p) por funciones de Cc1 (I ) que tienen esa propiedad.
90

En realidad no es completamente satisfactorio que las condiciones de contorno se satisfagan tan solo en un sentido debil. Afortunadamente, en dimension d = 1 las funciones de W 1;p(I ) son mejores de lo que cabra esperar, lo
que permitira dar sentido al valor que toman en un punto.

Teorema 4.6. Sea u 2 W 1;p(I ); existe una funcion u 2 C (I ) tal que u = u

c.t.p. en I ; es decir, u tiene un representante continuo.
Prueba. Basta tomar u(x)

u(y ) =

Rx 0
y u (t) dt,

8x; y 2 I .

En general se sustituira sistematicamente u por su representante continuo; y con el n de no recargar la notacion se designara tambien por u al
representante continuo. Podemos por tanto hablar del valor de una funcion u 2
W 1;p(I ) en un punto de I . Nos estaremos reriendo al valor de su representante
continuo en dicho punto. En particular, las funciones de W01;p (I ) se pueden
caracterizar de la forma siguiente.

Teorema 4.7. Sea u 2 W 1;p(I ); entonces u 2 W01;p si y solamente si u = 0

sobre @I .
As que ya tenemos el espacio que estabamos buscando: pediremos que
'l 2 H01 .
Para terminar la seccion probamos, usando la caracterizacion que acabamos
de dar, que en W01;p(I ) la cantidad ku0 kLp es una norma equivalente a la norma
de W 1;p.

Teorema 4.8 (Desigualdad de Poincare). Sea I acotado. Existe una constante C , dependiente de jI j, tal que

kukW ;p  C ku0kLp
1

Prueba. Para u 2 W01;p(I ) se tiene

ju(x)j = ju(x) u(a)j

8u 2 W01;p:

Z x 


=  u0 
a

 ku0kL ;
1

de donde se deduce el resultado gracias a la desigualdad de Holder.

91

4.3 Soluciones debiles. Formulacion de Galerkin

Conviene que revisemos lo que hemos hecho. Para tener mas exibilidad en
la eleccion del espacio de dimension nita Vh , hemos rebajado las hipotesis de
regularidad sobre las funciones 'l , pidiendo tan solo que 'l 2 H01 . Con esto
las ecuaciones de Galerkin (4.6) tienen perfecto sentido y permiten calcular
la aproximacion uh 2 Vh . Estas ecuaciones caracterizan a uh como la unica
funcion en Vh (subespacio de dimension nita de H01 ) tal que
Z b
a

u0 '0 +

Z b

uh ' =

Z b
a

f'

8' 2 Vh:

(4.7)

Observese que uh 2 H01 (I ), pero que en general uh 62 C 2 (I ). Es decir, la

aproximacion tiene menos regularidad que la verdadera solucion. Esto, que
en principio no es demasiado importante, produce sin embargo un punto debil
en nuestros argumentos. Para llegar a las ecuaciones (4.6) hemos partido de
las ecuaciones (4.4), que involucran al defecto dh; y el defecto no esta bien
denido para funciones que solo estan en H01 , pero no en C 2 . Esta situacion
parece inevitable, >acaso es posible escribir un problema en derivadas segundas
sin usar derivadas segundas? Sorprendentemente la respuesta es armativa en
muchos casos importantes, y conduce al concepto de solucion debil.
Sea u una solucion clasica de (4.1). Si multiplicamos la ecuacion por una
funcion ' 2 Cc1 (I ) e integramos por partes se tiene que
Z b
a

(u0 '0 + u') =

Z b
a

f':

(4.8)

De hecho, por ser Cc1 (I ) denso en H01 (I ), la igualdad es cierta para cada ' 2
H01 (I ). En efecto, dada ' 2 H01 (I ), existe una sucesion f'n g de funciones
de Cc1 (I ) tal que k'n 'kH (I ) ! 0 cuando n ! 1. Para cada una de las
funciones 'n se tiene que
1

Z b
a

(u0 '0n + u'n ) =

92

Z b
a

f'n :

(4.9)

Ahora bien, usando la desigualdad de Cauchy-Schwarz tenemos que

Z b

Z b


0
0
 (u0 '0 + u'n )

(
u
'
+
u'
)
n


Z a b
 a
Z b



f'n
f' = f; 'n ' L2

a
a

jh

= jhu; 'n

'iH

j  kukH k'n 'kH ;

1

i j  kf kL k'n 'kL :
2

Por tanto, pasando al lmite en (4.9) obtenemos que (4.8) es cierta para
cualquier ' 2 H01 (I ).
En realidad para que (4.8) tenga sentido basta con que f
H01 (I ), lo que motiva la siguiente denicion.

2 L2(I ) y u 2

Denicion 4.9. Sea f 2 L2 (I ). Una solucion debil de (4:1) es una funcion

u 2 H01(I ) que verica (4:8) para toda funcion ' 2 H01 (I ).
Con esta denicion las condiciones de contorno estan implcitas en el espacio
funcional.
Observacion. Sea f 2 C ([a; b]). Hemos visto que la solucion clasica es
solucion debil. Recprocamente, si una solucion debil es C 2 (I ) \ C (I ), entonces
es solucion clasica. Para verlo, integramos por partes la expresion (4.8) y
obtenemos que
Z b
a

( u00 + u

f )' = 0

8' 2 Cc1 (I ):

Supongamos que g = u00 + u f es distinto de cero para algun punto x0 2 I .

Para concretar supongamos que g (x0 ) > 0. Por continuidad, existe " > 0 tal
que g (x) > 0 en (x0 "; x0 + "). Tomando ' 2 Cc1 (I ) tal que '  0, '(x0 ) > 0,
sop '  (x0 "; x0 + "), tenemos que
Z b
a

( u00 + u f )' > 0;

una contradiccion. Concluimos que

u00 + u = f

si a < x < b:
93

Las condiciones de contorno se verican, por ser u una funcion de H01 (I ).
Sea Vh un subespacio de dimension nita de H01 (I ). Si pedimos que se satisfaga la condicion de solucion debil (4.8), pero restringiendonos a funciones uh
y ' del subespacio Vh (en lugar de considerar todo el espacio H01 ), obtenemos
(4.7), que es la llamada forma de Galerkin del metodo de elementos nitos.
Notese que para llegar a (4.7), que es equivalente a (4.6), no hemos usado para
nada el defecto dh, ni nada que nos exija una mayor regularidad para uh.
Se podra objetar que de esta forma lo que hemos obtenido es una aproximacion a una solucion debil del problema. Sabemos que la solucion clasica es
tambien solucion debil; pero, >no habra mas soluciones debiles? La respuesta
viene dada por el siguiente teorema.

Teorema 4.10. El problema (4:1) tiene una unica solucion debil.

Prueba. Sean u; v 2 H01 dos soluciones debiles. De (4.8) se tiene que
Z b
a

v )0 '0 +

(u

Z b
a

(u v )' = 0

Tomando ' = u v se tiene que

Z b
a

((u

v )0 )2 +

Z b
a

8' 2 H01:

(u v )2 = 0;

y por tanto u v = 0 c.t.p., es decir, u = v .

El metodo de elementos nitos da una aproximacion de la solucion debil
incluso cuando esta no es clasica. Pero, >sucede esto alguna vez? La respuesta
es armativa. Un ejemplo sencillo es tomar el problema (4.1) con un dato
f 2 L2 (I ), pero f 62 C ([a; b]). El problema tiene una solucion debil (ver
el teorema 4.12 mas adelante), pero no es, evidentemente, clasica. Veamos
que regularidad tiene una solucion debil que corresponda a un dato con estas
caractersticas. Sea u solucion debil, entonces u 2 H01 y ademas
Z b
a

u0 '0 =

Z b
a

(f

u)'
94

8' 2 Cc1:

Dado que f u 2 L2 , concluimos que u0 2 H 1 . La funcion u0 esta por tanto

en el espacio H 2 (I ), que se dene por recurrencia como

H 2 (I ) = fu 2 H 1 (I ) : u0 2 H 1 (I )g:
A la derivada debil de u0 , (u0 )0 se la llama derivada debil segunda de u; se
denota por u00 . En nuestro caso se tiene que u00 = u f .
Resumiendo, si f 2 L2 entonces u 2 H 2 . >Que pasa si f 2 C ([a; b])?
Tenemos entonces que u00 = u f 2 C ([a; b]), y por tanto u es solucion clasica.

4.4 Minimizacion de funcionales. Formulacion de Ritz

Un paradigma que se remonta a Euler es que el mundo fsico se explica imponiendo que ciertas cantidades (\acciones") son mnimas; por ello no es extra~no que podamos representar algunos problemas de ecuaciones ordinarias y
parciales (que tratan de representar fenomenos fsicos) a otros de minimizacon.
Sea f

2 L2 (I ). Denimos un funcional de energa J : H01 7! R mediante

J (v) = 12

Z b
a

((v 0 )2 + v 2 )

Z b
a

fv:

Teorema 4.11. (i) Si u 2 H01 es tal que J (u)  J (v ) para todo v 2 H01 (I ),
entonces u es solucion debil de (4:1).
(ii) Recprocamente, si u 2 H01 (I ) es solucion debil de (4:1), entonces J (u) 
J (v) para todo v 2 H01(I ).
Prueba. (i) Dada ' 2 H01 , denimos g () = J (u + '),  2 R . Se tiene
que g (0)  g () para todo  2 R . Si la derivada g 0(0) existe, se tendra que
g 0 (0) = 0. Ahora bien,

g () = J (u) + 

Z b
a

(u0 '0 + u'

2
f') +
2
95

Z b
a

(('0 )2 + '2 );

(4.10)

y por tanto

g 0(0) =

Z b
a

(u0 '0 + u' f'):

La condicion g 0 (0) = 0 para todo '

solucion debil de (4.1).

2 H01, es precisamente la condicion de

(ii) Si u 2 H01 (I ) es solucion debil, entonces (4.10) se convierte en

Z

2 b
J (u + ') = J (u) + 2 (('0)2 + '2 )  J (u) 8 2 R; ' 2 H01:
a
1
Para cualquier v 2 H0 (I ) tomamos ' = v u y  = 1 y obtenemos que
J (v)  J (u).

El area de la matematicas que estudia problemas de minimizacion es el

calculo de variaciones. Hemos visto que minimizar el funcional es equivalente
a encontrar una solucion debil de una ecuacion. Dicha ecuacion se conoce como
ecuacion de Euler-Lagrange del problema de minimizacion.
A partir del teorema observamos que una forma alternativa de encontrar
una solucion aproximada al problema (4.1) es minimizar J (v ) en un subespacio
Vh  H01 de dimension nita; es decir, buscamos uh 2 Vh tal que

J (uh)  J (vh)

8vh 2 Vh:
Sea f'1 ; : : : ; 'm g una base de Vh. Cualquier vh 2 Vh se puede escribir como
vh =
y el valor del funcional sera
m
1X
J (vh) = 2 l k
k;l=1

Z b
a

m
X
l=1

l 'l ;

('0l '0k + 'l 'k )

m
X
k=1

Z b
a

f'k :

Minimizar el funcional en Vh es por tanto equivalente a minimizar la forma

cuadratica

Q( 1 ; : : : ; m ) =

m
X

k;l=1

k a^kl l

m
X
k=1

Z b

96

f'k ;

a^kl =

Z b
a

('0l '0k + 'l 'k )

La condicion rQ( 1 ; : : : ; m ) = 0 es equivalente al sistema lineal (4.6). Obtenemos las ecuaciones de Galerkin partiendo de un punto de vista distinto
(aunque equivalente). Esta forma de llegar al sistema (4.6) se conoce como
formulacion de Ritz del metodo de elementos nitos.

4.5 Formulacion abstracta

Nos gustara generalizar lo que hemos hecho para el problema (4.1) a otros
problemas. No queremos repetir los mismos argumentos cada vez, as que
introduciremos una formulacion abstracta que nos permita tratar de forma
unicada una amplia clase de problemas.
Empezamos por escribir la formulacion debil del problema (4.1) de una
manera susceptible de generalizacion. Tenemos un espacio de Hilbert V = H01 ,
una forma bilineal sobre V  V dada por

a(u; v ) =

Z b
a

(u0 v 0 + uv );

y una forma lineal sobre V dada por

L(v ) =

Z b
a

fv:

Lo que hemos hecho ha sido encontrar una funcion u 2 V tal que

a(u; v ) = L(v );

8v 2 V ;

alternativamente, hemos encontrado u 2 V tal que J (u) = min J (v ), siendo

v 2V

J (v) = 12 a(v; v) L(v):

>Podemos hacer esto para otras formas a y L, denidas quiza sobre otro
espacio de Hilbert distinto? La respuesta es armativa; pero habra que pedir
que las formas a y L tengan las propiedades que nos han permitido probar
97

los teoremas en el caso del problema (4.1). Eso nos lleva a introducir algunos
conceptos nuevos.
Sea V un espacio de Hilbert con producto escalar h
;
iV y norma correspondiente k
kV . Sean a(
;
) una forma bilineal sobre V  V y L una forma
lineal sobre V tales que:

 a(
;
) es simetrica;
 a(
;
) es continua, es decir, existe una constante  > 0 tal que
ja(u; v)j  kukV kvkV

8u; v 2 V ;

 a(
;
) es coerciva, es decir, existe una constante  > 0 tal que

a(u; u)  kuk2V

8v 2 V ;

 L es continua, esto es, existe una constante  > 0 tal que

jL(v)j  kvkV

8v 2 V:

Consideramos los siguientes problemas abstractos:

Problema de minimizacion (M): encontrar u 2 V tal que

J (u) = min
J (v ) ;
v2V
donde

J (v) = 21 a(v; v) L(v):

Problema en forma debil (D): encontrar u 2 V tal que

a(u; v ) = L(v )
98

8v 2 V:

(4.11)

Teorema 4.12. Los problemas (M ) y (D) son equivalentes; esto es, u 2 V es

solucion de (M ) si y solo si es solucion de (D). Estos problemas tienen una
unica solucion, y se satisface la estimacion de estabilidad
kukV   :
(4.12)
Prueba. La existencia de solucion de (D) se obtiene usando el teorema de
Lax-Milgram, que es una variante del teorema de representacion de Riesz (para
los detalles, ver por ejemplo [B]).

La prueba de la equivalencia de (M) y (D) es analoga a la que dimos en el

Teorema 4.11 para el problema (4.1). Probemos, por ejemplo, que si u 2 V
satisface (M) entonces satisface (D). Sean v 2 V y  2 R arbitrarios. Entonces
u + v 2 V , y por tanto

J (u)  J (u + v)

8 2 R :

Deniendo g () = J (u + v ) vemos que g tiene un mnimo en  = 0. Ahora

bien,

g () = 12 a(u + v; u + v ) L(u + v )

= 12 a(u; u) + 2 a(u; v ) + 2 a(v; u) + 2 a(v; v ) L(u) L(v )
= 12 a(u; u) L(u) + a(u; v ) L(v ) + 2 a(v; v );
2

(4.13)

donde hemos usado la simetra de a(
;
). Se sigue que

0 = g 0(0) = a(u; v ) L(v );
es decir, u es solucion de (D).

Para probar la implicacion contraria observamos que si u 2 V es solucion

de (D), entonces (4.13) se convierte en
2

J (u + v) = J (u) + 2 a(v; v)  J (u)

Dado ' 2 V arbitrario, tomamos ' = v

J (u), es decir, u es solucion de (M).
99

8 2 R; v 2 V:

u y  = 1, y tenemos que J (') 

Para probar la estimacion de estabilidad tomamos v = u en (4.11) y usamos

la coercividad de a y la continuidad de L para obtener que

kuk2V

 a(u; u) = L(u)  kukV :

Dividiendo por kukV =

6 0 obtenemos (4.12).
Veamos la unicidad. Sean u1 ; u2 2 V tales que

8v 2 V; i = 1; 2:

a(ui ; v ) = L(v )
Restando tenemos que u1

u2 2 V satisface

a(u1
Tomando v = u1

u2 ; v ) = 0

8v 2 V:

u2 y usando la coercividad de a tenemos nalmente que

ku1

u2 kV

 a(u1 u2; u1 u2) = 0;

es decir, u1 = u2 .
Observacion. Se puede probar la existencia de una unica solucion u 2 V del
problema (D) sin la condicion de simetra sobre a. Sin embargo, en este caso
no habra un problema de minimizacion asociado.
Observacion. Para que el problema de minimizacion (M) este bien planteado
no es necesario que V sea un espacio de Hilbert; basta con que sea un subconjunto convexo cerrado de un espacio de Hilbert. El problema (M) tendra
entonces solucion unica, aunque no sera equivalente al problema (D) (ver el
ejemplo de condiciones de Dirichlet no homogeneas mas adelante).
Ejemplo. (Condiciones de Neumann). Consideramos el problema

u00 + u = f si a < x < b;

u0 (a) = ;

u0(b) = :

Una solucion clasica es una funcion u 2 C 2 ((a; b)) \ C 1 ([a; b]) que verica la
ecuacion y las condiciones de contorno de forma puntual. Notese que para dar
100

sentido a las condiciones de contorno hemos tenido que pedir que la funcion
sea derivable hasta la frontera.
Para escribir el problema en forma debil, multiplicamos la ecuacion por
v 2 H 1 e integramos por partes. Imponiendo las condiciones de frontera, se
tiene que
Z b
a

(u0 v 0 + uv ) =

Z b
a

fv

v (a) + v (b)

8v 2 H 1:

Una solucion debil del problema de Neumann sera una funcion u

verique (4.14).

(4.14)

2 H 1 que

Observacion. En realidad, para ser rigurosos en el razonamiento anterior

tendramos que haber multiplicado la ecuacion por una funcion v 2 Cc1 (R ) y
despues haber utilizado un argumento de densidad para probar el resultado
para v 2 H 1 (I ) general. El resultado de densidad que se necesita, cuya prueba
se puede encontrar por ejemplo en [B], es el siguiente.

Lema 4.13. Sea u 2 W 1;p(I ) con 1  p < 1. Entonces existe una sucesion
fung en Cc1(R ) tal que unjI ! u en W 1;p(I ).
En posteriores ejemplos omitiremos este tipo de argumentos de densidad.
Queremos formular el problema en forma abstracta. Para ello tomamos
R
R
V = H 1, a(u; v ) = ab (u0 v 0 + uv ) y L(v ) = ab fv v (a) + v (b). Se puede
probar que a(
;
) y L(
) estan en las hipotesis del Teorema 4.12; por tanto
existe una unica solucion debil del problema y ademas coincide con la funcion
u 2 H 1 que minimiza el funcional
Z

b
b
J (v) = 21 ((v0)2 + v2)
fv + v (a) + v (b):
a
a
Observese que no se exige explcitamente que la solucion satisfaga las condiciones de frontera. Las condiciones aparecen implcitamente en la forma lineal
L, y a traves de ella en el funcional de energa. En el caso de datos Dirichlet homogeneos, por el contrario, se imponen explcitamente los datos de contorno a

101

la solucion, exigiendo que esta pertenezca al espacio funcional H01 . Las condiciones de frontera que no necesitan ser impuestas explcitamente se llaman
condiciones de frontera naturales; las que s lo necesitan se llaman condiciones
de frontera esenciales.
Para nalizar el ejemplo, probaremos que una solucion debil que pertenezca
a C 2 (I ) \ C 1 (I ) es solucion clasica. Notese que la funcion esta por hipotesis
en el espacio correcto. Empezamos probando que se satisface la ecuacion para
todo x 2 I . Para ello, consideramos (4.14) para el caso particular de v 2 Cc1 (I ).
Integrando por partes se tiene que
Z b
a

( u00 + u f )v = 0:

Argumentando como lo hicimos en el caso de datos de Dirichlet homogeneos

se obtiene que u00 + u = f para todo x 2 I .
Una vez probado esto volvemos a (4.14), y tomamos ahora una funcion
v 2 C 1 (I ) tal que v (a) = 0, v (b) = 1. Integrando por partes se tiene que
Z b
a

( u00 + u f )v = 

v (b):

Dado que ya hemos probado que se verica la ecuacion, se tiene que  v (b) =
0, es decir, la condicion de contorno en x = b. Para probar la condicion de
contorno en x = a hay que tomar v 2 C 1 (I ) tal que v (a) = 1, v (b) = 0.
Ejemplo. (Condiciones de contorno mixtas). Consideramos el problema

u00 + u = f si a < x < b;

u(a) = 0;

u0 (b) = :

(4.15)

La condicion Dirichlet es esencial, y tiene que aparecer en el espacio funcional. As que para obtener la formulacion debil multiplicamos la ecuacion
por v 2 H := fv 2 H 1 : v (a) = 0g. Tras integrar por partes, imponiendo las
condiciones de frontera, llegamos a
Z b
a

(u0v 0 + uv ) =

Z b
a

fv + v (b)
102

8v 2 H:

(4.16)

Una solucion debil del problema (4.15) sera una funcion u

que (4.16).

2H

que veri-

Para escribir el problema en forma abstracta tomamos V = H , a(u; v ) =

R
+ uv ), L(v ) = ab fv + v (b). La condicion natural u0 (b) =  aparece
en el funcional, y no en el espacio H .
Rb 0 0
a (u v

Ejemplo. (Condicion de Robin). Consideramos el problema

u00 + u = f si a < x < b;

u0(a) ku(a) = 0;

u(b) = 0;

(4.17)

donde k  0 es un dato del problema. Necesitamos un espacio funcional que

recoja la condicion de frontera esencial; H := fv 2 H 1 : v (b) = 0g cumple ese
requisito. Una solucion debil de (4.17) sera una funcion u 2 H que verique
Z b
a

(u0v 0 + uv ) + ku(a)v (a) =

Z b
a

8v 2 H:

fv

La condicion de frontera de tipo Robin es una condicion natural, pues aparece

R
en el funcional; lo hace a traves de la forma bilineal a(u; v ) = ab (u0 v 0 + uv ) +
R
ku(a)v (a), siendo L(v ) = ab fv . Notese la diferencia con el problema de
Neumann, donde la condicion de frontera apareca en la forma lineal L.
Hasta aqu hemos tratado problemas en una dimension. Para tratar problemas en mas dimensiones tenemos que introducir los espacios de Sobolev
correspondientes.

Denicion 4.14. Sea

 R d abierto. El espacio de Sobolev H 1 (
) se dene
como
H 1 (
) = fu 2 L2 (
) : 9g1 ; : : : ; gd 2 L2 (
)
R
@' = R g ' 8' 2 C 1 (
) 8i = 1; : : : ; dg:
tales que
u @x
c

i
i
El espacio H 1 (
) esta dotado del producto escalar

hu; viH

= hu; v iL +
2

d 
X
@u
i=1

103

@v
;
@xi @xi


L2

y de la norma inducida por el

kukH

2 !1=2
d
X


@u
2 +


:
L2
@x 2
i
L
i=1

kuk

Con esta norma es un espacio de Hilbert.

Tambien necesitamos un espacio que recoja las condiciones Dirichlet (que
son esenciales).

Denicion 4.15. Sea

Cc1 (
) en H 1 (
).

 Rd

abierto. El espacio H01 (

) es la clausura de

El espacio H01 es un espacio de Hilbert con la norma que hereda de H 1 .

Para d > 1 no es cierto que las funciones de H 1 tengan un representante
continuo. >En que sentido verican entonces que u = 0 para x 2 @
? Esto lle es un asunto bastante delicado, y no lo trataremos
va a la teora de trazas. Este
aqu. Digamos simplemente que existe un operador lineal y continuo

tr : H 1 (
) 7! L2 (@
)
tal que tr(u) coincide con uj@
cuando u es una funcion continua. De esta
forma el espacio de Sobolev H01 (
) verica que

H01 (
) = fu 2 H 1 (
) : tr(u) = 0 c.t.p. en @
g:
La continuidad del operador garantiza que

ktr(u)kL (@
)  C kukH (
) :
2

Por fortuna estas complicaciones tecnicas no nos afectan en la practica, y

fundamentalmente sirven para asegurar que existe un espacio funcional en el
que construir la teora.
104

Veamos un par de ejemplos con d > 1. Para pasar a la formulacion debil

sera fundamental la formula de Green (desempe~nara el papel de la integracion
por partes)
Z
Z
Z
v
u =
rv
ru + vru
n;

donde n = n(x) es el vector normal exterior al dominio en el punto x 2 @

.
Esta formula se obtiene a partir de la igualdad
div (v ru) = rv
ru + v
u
integrando en
y aplicando el teorema de la divergencia.
Ejemplo. (Problema de Dirichlet homogeneo d-dimensional). Consideramos
el problema
(

u + u = f
en
 R d ;
u=0
en @
:
Multiplicamos por v 2 H01 , aplicamos la formula de Green y las condiciones de
contorno y obtenemos
Z

ru
rv + uv =

8v 2 H01(
);

fv

que es la forma debil del problema. Para aplicar el Teorema 4.12 tomamos
R
R
V = H01 (
), a(u; v ) =
ru
rv y L(v ) =
fv . Veamos que se cumplen las
condiciones del teorema. En efecto, a(
;
) es obviamente bilineal y simetrica.
Por otra parte,
a(u; u) = kukH (
) ;
1
0

luego a(
;
) es coerciva. La continuidad es consecuencia de la desigualdad de

Cauchy-Schwarz, pues
a(u; v ) = hu; v iH (
) :
1
0

La linealidad de L es tambien trivial. La continuidad es consecuencia de la

desigualdad de Cauchy-Schwarz en L2 (
),

jL(v)j =



fv 

 kf kL (
) kvkL (
)  kf kL (
) kvkH (
)
2

105

1
0

8v 2 H01(
):

Ejemplo. (Problema de Neumann d-dimensional). Consideramos el problema

(

en
 R d ;
en @
:

u + u = f
@u = g
@n

Multiplicamos por v 2 H 1 (recordemos que la condicion de Neumann es natural), aplicamos la formula de Green y las condiciones de contorno y obtenemos
Z

ru
rv + uv =

f+

8v 2 H 1(
);

gv

que es la forma debil del problema. Para aplicar el Teorema 4.12 tomamos
R
R
R
V = H 1 (
), a(u; v ) =
ru
rv y L(v ) =
fv + @
gv . La forma bilineal
a(
;
) cumple las condiciones del teorema (copiar del ejemplo anterior). Para
ver la continuidad de L observamos que

jL(v)j 

 Z
 
fv  + 



gv  :

Para probar la continuidad del segundo sumando aplicamos la desigualdad de

Cauchy-Schwarz en L2 (@
) y la continuidad del operador traza,
Z






gv 

 kgkL (@
) kvkL (@
)  C kgkL (@
) kvkH (
)
2

8v 2 H 1(
):

Ejemplo. (Problema de Dirichlet no homogeneo d-dimensional).

ramos el problema
(

u + u = f
u=g

en
 R d ;
en @
:

Conside(4.18)

La condicion de Dirichlet debera aparecer en el espacio funcional. As pues,

nos gustara tomar

V = fv 2 H 1 (
) : v = g en @
g:
El conjunto V no es un espacio vectorial, pero es un subconjunto convexo
cerrado del espacio de Hilbert H 1 (
), lo que sera suciente para nuestros
106

propositos. Sin embargo, si multiplicamos la ecuacion por funciones test v 2 V

R
@u ,
y aplicamos la formula de Green, aparecera un termino de la forma @
g @n
que no se puede incorporar ni a la forma bilineal ni a la lineal.
Con el n de que no aparezca el termino de borde molesto, utilizamos como
tests funciones v 2 H01 (
). Obtenemos que
Z

(ru
rv + uv ) =

fv

8v 2 H01(
):

(4.19)

Esto sugiere la siguiente denicion: una solucion debil del problema (4.18)
R
es una funcion u 2 V que verica (4.19). Tomando a(u; v ) =
ru
rv y
R
L(v ) =
fv el problema debil se puede escribir en forma abstracta: encontrar
una funcion u 2 V que verica que

a(u; v ) = L(v )

8v 2 H01(
):

>Como garantizar la existencia y unicidad de solucion debil? En H 1 (

) la
forma bilineal a(
;
) es continua y coerciva y la forma lineal L(
) es continua.
Sin embargo, no podemos aplicar el teorema 4.12, pues, por una parte la solucion u y los tests v no estan en el mismo espacio, y por otra el conjunto V no
es un espacio de Hilbert. As que adoptamos una estrategia distinta, e intentamos ver la solucion debil como solucion de un problema de minimizacion.
Para ello consideramos el funcional
J (v) = 12 a(v; v) L(v);
que esta bien denido para funciones v 2 V . Es facil ver (repitiendo los
argumentos de la prueba del teorema 4.12), que una funcion u 2 V es solucion
debil de (4.18) si y solo si

J (u) = min
J (v):
v2V

As que el problema debil es equivalente al problema de encontrar un minimizante para J en V , que es un subconjunto convexo cerrado del espacio
de Hilbert H 1 (
). Por tanto el problema de minimizacion tiene una unica
solucion. Concluimos que hay una unica solucion debil de (4.18).
107

4.6 Discretizacion. Convergencia

Sea ahora Vh un subespacio de V de dimension nita M . Podemos formular
los siguientes problemas discretos, analogos a los problemas (M) y (D):

Problema (Mh ): encontrar uh 2 Vh tal que

J (uh)  J (v)


8v 2 Vh:

Problema (Dh ): encontrar uh 2 Vh tal que

a(uh; v ) = L(v )

8v 2 Vh:

Sea f'1 ; : : : ; 'M g una base de Vh , de forma que cualquier u

representar en la forma

u=

M
X
j =1

j 'j ;

(4.20)

2 Vh se puede

j 2 R :

Usando la bilinealidad y la simetra de a, as como la linealidad de L, vemos

que ambos problemas son equivalentes a resolver el sistema lineal
M
X
i=1

a('i ; 'j ) j = L('i );

i = 1; : : : ; M:

Deniendo = ( i) 2 R M , b = (bi ) 2 R M con bi = L('i ), A = (aij ) una matriz

M  M con elementos aij = a('i ; 'j ), podemos escribir el sistema en forma
matricial como
A = b:
La matriz A recibe el nombre de matriz de rigidez y b es el vector de carga. Estos nombres provienen de las aplicaciones mas antiguas del metodo de
elementos nitos, que correspondan a calculos de estructuras mecanicas.

Lema 4.16. La matriz de rigidez A es simetrica y denida positiva.

108

Prueba. La simetra se hereda de la simetra de la forma bilineal a. Que A

sea denida positiva es consecuencia de la coercividad.

Podemos probar el siguiente resultado basico.

Teorema 4.17. Existe una unica solucion de los problemas equivalentes (Mh )
y (Dh ). Ademas se tiene la siguiente estimacion de estabilidad
kuhkV   :
Prueba. A es denida positiva, y por tanto no singular, lo que prueba la
existencia y unicidad. Para probar la estimacion de estabilidad se toma v = uh
en (4.20); la coercividad de a y la continuidad de L producen

kuhk2V

 a(uh; uh) = L(uh)  kuhkV :

Dividiendo por kuhkV =
6 0 se obtiene el resultado. Si kuhkV
trivialmente.

= 0 se cumple

A continuacion obtenemos una estimacion para el error.

Teorema 4.18 (Lema de Cea). Sea u 2 V solucion de (D) y uh

solucion de (Dh ) con Vh  V . Entonces
ku uhkV   ku wkV 8w 2 Vh:

2 Vh

Prueba. Dado que Vh  V , de (4.11) tenemos en particular que

a(u; v ) = L(v )

8v 2 Vh;

restando de aqu (4.20) llegamos a

a(u uh; v ) = 0

8v 2 Vh:

(4.21)

Para w 2 Vh arbitrario denimos v = uh w. Entonces v 2 Vh . Usando la

coercividad y la continuidad de a y la ecuacion (4.21) tenemos

ku uhk2V  a(u uh; u uh ) = a(u uh ; u uh ) + a(u uh ; v )

= a(u uh; u uh + v ) = a(u uh; u w)  ku uh kV ku wkV :
109

Dividiendo por ku uhkV =

6 0 se obtiene el resultado.
Podemos denir una nueva norma k
ka en V , llamada norma de energa,
mediante
kvk2a = a(v; v); v 2 V:
La coercividad y la continuidad de a nos dicen que esta norma es equivalente
a la norma k
kV . La norma k
ka proviene del producto escalar

hv; wia = a(v; w):

La ecuacion (4.21) se puede escribir

hu uh; via = 0

8v 2 Vh;

(4.22)

es decir, la solucion de elementos nitos uh es la proyeccion de la solucion

exacta u sobre Vh con respecto al producto escalar h
;
ia.
Por otra parte, la propiedad (4.22) nos dice que hu uh ; uh v ia = 0 para
todo v 2 Vh. Por tanto, el teorema de Pitagoras nos dice que

ku uhk2a + kuh vk2a = ku vk2a;

de donde

ku uhka  ku vka

8v 2 Vh;

es decir, uh es el elemento de Vh mas proximo a u en la norma de energa.

El lema de Cea nos da una condicion suciente para la convergencia del
metodo de elementos nitos.

Corolario 4.19. Sea fVh gh>0 una familia de subespacios de V de dimension

nita tales que para todo u 2 V se tiene que
lim inf ku uhkV = 0;

h!0 v2Vh

entonces el metodo de elementos nitos converge, es decir,

lim ku uhkV = 0:

h!0

110

(4.23)

Observacion. Para que el metodo converja basta con que la condicion (4.23)
se cumpla para la verdadera solucion.

Notese que hemos reducido la cuestion de la convergencia del metodo de

elementos nitos a una cuestion de teora de aproximacion: >Como de lejos de
Vh estan las funciones de V ?

4.7 Espacios de elementos nitos

Queremos elegir una familia de subespacios fVhgh>0 de forma que se verique
(4.23) al menos cuando u es la verdadera solucion.
Una eleccion comun es tomar espacios de funciones polinomicas a trozos
sobre subdivisiones o \triangulaciones" Th = fK g de un dominio acotado

 R d , d = 1; 2; 3, en elementos K . Para d = 1 los elementos K seran

intervalos, para d = 2, triangulos o cuadrilateros y para d = 3 tetraedros.
Si estamos estudiando un problema de contorno de segundo orden, Vh tendra que estar contenido en H 1 (
). Dado que el espacio Vh esta formado por
funciones polinomicas a trozos, se tiene que

Vh  H 1 (
) , Vh  C (
):
Para denir un espacio Vh del tipo mencionado tendremos que especicar:





4.7.1

la triangulacion Th = fK g del dominio

;
la naturaleza de las funciones v de Vh sobre cada elemento K (lineal,
cuadratica, cubica, etc);
una base de Vh .
Caso de una variable espacial

Dada una particion

a = x0 < x1 <


< xN = b;

111

del intervalo [0; L] denotaremos

hj = xj

xj 1 ;

j = 1; : : : ; N;

h = max hj :
1j N

Nuestros elementos K seran en este caso los intervalos [xi 1 ; xi ].

Para r = 1; 2; : : : , consideramos los espacios Sh;r de funciones continuas en
[a; b] que restringidas a cada intervalo [xj 1 ; xj ] son polinomios de grado menor
o igual que r, esto es

Sh;r = fv 2 C ([a; b]) : vj xj

[

1 ;xj ]

2 Pr ([xj 1; xj ]); j = 1; : : : ; N g:

La dimension de Sh;r es facil de calcular. Una funcion  2 Sh;r depende

en principio de (r + 1)N parametros, al haber N intervalos y, en cada uno
de ellos  es un polinomio de grado r, que depende de r + 1 parametros.
Para mantener la continuidad en cada nodo interno xj de la particion hay que
imponer la condicion (xj +) = (xj ). Los nodos internos son los N 1 nodos
x1 ; : : : ; xN 1 . La dimension de Sh;r es por tanto N (r + 1) (N 1) = rN + 1.
Todava tenemos que indicar que parametros vamos a utilizar para describir
las funciones de Vh. Eso es equivalente a escoger una base en Vh . Veamos como
hacerlo en el caso de elementos lineales (r = 1). En este caso consideramos los
nodos x0 ; : : : ; xN y consideramos las llamadas funciones base nodales:

'n (x) =

8
>
>
>
>
>
>
>
>
>
<
>
>
>
>
>
>
>
>
>
:

x xn 1
;
xn xn 1

x 2 [xn 1 ; xn ];

xn+1 x
;
xn+1 xn

x 2 [xn ; xn+1 ];

0;

en el resto:

'0 (x) =

8
<
:

x1 x
;
x1 x0
0;
112

n = 1; : : : ; N

x 2 [x0 ; x1 ];
en el resto:

1;

8
<

x xN 1
;
x 2 [xN 1 ; xN ];
'N (x) = xN xN 1
:
0;
en el resto:
Obviamente cualquier funcion vh 2 Sh;1 se puede escribir como
vh =

N
X
n=0

n n:

Dado que i(xj ) = ij tenemos que

n = vh (xn ):
As pues, los N + 1 parametros mediante los cuales se describe la funcion vh
son sus valores en los N + 1 nodos; de ah el nombre de base nodal.
En el caso en que r > 1 el numero de grados de libertad es mayor que el
de nodos. Lo que se hace es a~nadir nodos intermedios (denominados nodos
fraccionarios), entre los de la particion, dados por
j
xi+j=r = xi + hi ;
i = 0; : : : ; N 1;
j = 1; : : : ; r 1:
r
De esta forma se tienen tantos nodos como grados de libertad (rN + 1), y
se pueden elegir como parametros para describir cualquier funcion de Sh;r sus
valores en los nodos. Esto denira una base, la base nodal, caracterizada
porque las funciones de la base toman el valor 1 en el nodo (fraccionario o
entero) al que estan asociadas y el valor 0 en el resto de los nodos.
Estudiemos la convergencia del metodo. Nos limitaremos al caso r = 1; es
decir, Vh = Sh;1.
Observamos que inf ku v kV  ku uI kV , donde uI 2 Vh es el interpolante
v2Vh
lineal a trozos de u en los nodos x0 ; : : : ; xN de la particion. Se tiene el siguiente
resultado de aproximacion.

Lema 4.20. Si u 2 C 2 ([a; b]), entonces

1
max ju uI j  h2j max ju00j;
[xj ;xj ]
8 [xj ;xj ]
1

113

max ju0

[xj 1 ;xj ]

u0I j  hj max ju00 j;

[xj 1 ;xj ]

y por tanto

ku uI kH (I )  Ch max
ju00j:
[a;b]
1

Prueba. Consideramos la diferencia
(x) = u(x) uI (x) sobre el intervalo

[xj 1 ; xj ]. Dado que
se anula en ambos extremos del intervalo, debe existir
al menos un punto z en el que
0 (z ) = 0. Entonces para todo x,

0 (x) =

Z x
z

00 (y ) dy:

Ahora bien, como uI es lineal,
00 = u00 , y se tiene que

Z x


j
0(x)j = 



u00 (y ) dy   hj max ju00 j;

x 2 [xj 1 ; xj ];

[xj 1 ;xj ]

que es la segunda de las estimaciones.

El maximo de j
(x)j en [xj 1 ; xj ], ocurrira en un punto en el que la derivada
se anula,
0 (z ) = 0. Miramos en que mitad del intervalo esta z ; supongamos,
por ejemplo, que esta mas proximo al extremo derecho, de manera que xj z 
hj =2. Desarrollando en serie de Taylor alrededor de z ,
1

(xj ) =
(z ) + (xj z )
0 (z ) + (xj z )2
00 (w);
2
donde z < w < xj . Usando que
se anula en el extremo xj y que
00 = u00 ,

j
(z)j


1
=  (xj
2



2
00
z )
(w)

 18 h2j [xmax
j
u00 j:
j ;xj ]
1

Combinando las dos primeras estimaciones se obtiene inmediatamente la

estimacion del error de interpolacion en norma H 1

ku uI kH (I ) =
1



h4 (b a)
64

N Rx
P
j
xj
j =1

+ h2 ( b

[(u uI )2 + (u0


a) (max ju00 j)2

[a;b]

!1=2

u0I )2 ]

1=2

 Ch max
ju00j:
[a;b]

El lema garantiza la convergencia de la solucion de elementos nitos a la

verdadera solucion u si esta es tal que u 2 C 2 ([a; b]).
114

Ejemplo. Sea u la solucion (clasica) del problema (4.1) con dato f 2 C ([a; b]).
Por ser u00 = u f , y ser u 2 C ([a; b]) tenemos a posteriori que u00 2 C ([a; b]).
Por tanto tenemos convergencia de la solucion de elementos nitos a la solucion
clasica.

El resultado del lema 4.20 no es en realidad el mas adecuado para nuestros

nes. El factor h es correcto, y re eja la tasa observada de decrecimiento del
error cuando se rena la particion. La imperfeccion esta en el otro factor,
max ju00 j. Las soluciones debiles u de problemas elpticos de segundo orden no
verican en general que u 2 C 2 ([a; b]), sino que u 2 H 2 (I ). La estimacion
idonea para el error de interpolacion se da a continuacion.

Lema 4.21. Si u 2 H 2 (I ), entonces

ku uI kL (I )  Ch2 ku00kL (I ) ;
2

y por tanto

ku0 u0I kL (I )  Chku00kL (I ) ;

ku uI kH (I )  Chku00kL (I ) :
1

Se puede encontrar una prueba, basada en desarrollos en serie de Fourier,

en [SF].
Combinando este resultado con el lema de Cea tenemos el siguiente resultado de convergencia.

Corolario 4.22. Sea I = (a; b), u 2 V  H 1 (I ) solucion de (D) y uh

solucion de (Dh ) con Vh  V . Si u 2 H 2 (I ) , entonces

2 Vh

ku uhkH (I )  Chku00kL (I ) :
1

4.7.2

Caso de dos variables espaciales

Para simplicar supondremos que @

es poligonal, en cuyo caso diremos que

es un dominio poligonal (si el dominio

no es poligonal, podemos aproximarlo
115

por un dominio que s lo sea). Hacemos una triangulacion del dominio
,
t
subdividiendolo en un conjunto Th = fKj gnj =1
de triangulos Kj que satisfacen
nt
[

k=1

Kj :

siendo la interseccion Kj \ Kl , j 6= l, o bien vaca, o bien un lado comun, o

bien un vertice comun. As pues, la unica restriccion es que un vertice de un
triangulo no puede estar en mitad de una arista de otro; en otras palabras,
una conguracion como

(donde la posicion de los vertices se enfatiza con un `') no esta permitida. A

continuacion se muestra una triangulacion que se ajusta a las reglas.

Introducimos el parametro de malla

h = Kmax
diam (Kj );
2T
j

donde diam (Kj ) = diametro de Kj = longitud del mayor lado de Kj .

Para r = 1; 2; : : : , consideramos los espacios Sh;r de funciones continuas en

que restringidas a cada triangulo Kj son polinomios de grado menor o igual

116

que r, esto es

Sh;r = fv 2 C (
) : vjKj

2 Pr (Kj ); j = 1; : : : ; nt g:

Como en el caso unidimensional, queremos caracterizar las funciones de

Sh;r por sus valores en un conjunto de nodos

x 1 ; : : : ; xM :
Veamos como hacerlo en el caso de funciones lineales (r = 1). En este caso los
nodos son los vertices de la triangulacion.
Empezamos por probar que las funciones de Sh;1, y en particular las funciones base, quedan completamente determinadas por sus valores nodales. Sobre el triangulo Kl las funciones v 2 Sh;1 son de la forma

vjKl (x; y ) = al + bl x + cl y:
Supongamos que el triangulo Kl en cuestion es
xl1

xl3
xl2

de nodos xl = (xl ; yl ), xl = (xl ; yl ) y xl = (xl ; yl ). La funcion v toma

valores vl , vl , vl en los nodos xl , xl , xl , lo que da lugar al sistema lineal
1

al + bl xli + cl yli = vli ;

i = 1; 2; 3:

Como los vertices no son colineales, el sistema es no singular y se puede resolver

facilmente, produciendo al , bl y cl .
Para que la funcion v resultante este bien denida necesitamos que las
restricciones de v a triangulos Ki , Kj con una arista comun coincidan sobre
117

la arista. Dichas restricciones son funciones lineales en un segmento, y estan

unvocamente determinadas por los valores en los extremos del segmento. Como los extremos del segmento son nodos, los valores de ambas restricciones
coinciden en ellos (son valores nodales, que estan prescritos); se concluye que
las dos restricciones coinciden a lo largo de la arista.
La base conveniente para describir las funciones de Sh;1 por sus valores
nodales es la base nodal; las funciones de dicha base, '1 ; : : : ; 'M 2 Sh;1, se
denen por la propiedad
'i (xj ) = ij :
De esta forma cualquier v 2 Sh;l se expresa como

v=

M
X
j =1

v (xj )'j :

El soporte de la funcion base asociada a un nodo es la union de los triangulos

de los que dicho nodo es vertice. En efecto, para el resto de los triangulos toma
el valor cero en los nodos, y es por tanto nula. A continuacion mostramos tres
ejemplos de funciones nodales, dibujadas sobre su soporte.

.....................
....... ................
................... .........
.
...

....... ................
...
... ........................
... .........
........
.
.
.
.
.
.. .
........ ....
..
..

.............
... ............
... ............
...
..............................
. .......
..
...
..
.

Por motivos evidentes, las funciones base reciben tambien el nombre de

funciones piramide. El numero de elementos en cada soporte puede cambiar
de vertice a vertice, y la descripcion de las funciones base nodales, aunque
sencilla en teora, acaba convirtiendose en algo engorroso e inconveniente desde
el punto de vista practico.
Sin embargo, restringidas a cada triangulo las funciones base tienen una
118

expresion sencilla. Como veremos mas adelante, esto es todo lo que se necesita
en la practica.
Las unicas funciones de la base nodal que no son identicamente nulas sobre
un triangulo dado son las correspondientes a los nodos del triangulo. Estas
funciones base se expresan en terminos de las funciones base para un triangulo
 es el triangulo de
muy especial, el llamado elemento de referencia K0 . Este
vertices
" #
" #
" #
0
1
0
;
2 =
;
3 =
;
1 =
0
0
0
que representamos en la gura en el plano  ,  .


-

Para cada i = 1; 2; 3 denotamos por Ni a la funcion base del triangulo de

referencia asociada al nodo  i , y N a la aplicacion
2

N1 (;  )
7
N2 (;  ) 7
5:
N3 (;  )

N (;  ) = 6
4

Las funciones Ni son muy faciles de calcular. Se tiene que

2
6

N (;  ) = 6
4

 



7
7:
5

El triangulo K , de vertices
"

x1 =

x1
;
y1

"

x2 =

x2
;
y2

119

"

x3 =

x3
;
y3

se obtiene a partir del triangulo de referencia K0 mediante la transformacion

afn dada por

 7! x = ( ) = x1 + [x2

x1 ; x3

x1 ];

donde, obviamente, [x2 x1 ; x3 x1 ] es la matriz cuyas columnas son los vectores x2 x1 y x3 x1 , aplicacion que expresada en funcion de sus componentes
se escribe como
" #
"
# "
# "
#" #

x
x1
x2 x1 x3 x1

=
7! x =
=
+
:

y
y1
y2 y1 y3 y1

Gracamente,


La aplicacion

x3

x2
x1

-
2

-x

'1 (x; y )
6
7
7
'(x; y ) = 6
'
(
x;
y
)
4 2
5
'3 (x; y )
de las funciones base correspondientes a nodos del triangulo K se expresa
entonces como
' = N 1;
cuando nos restringimos a dicho triangulo.
Recordando que Ni ( j ) = ij , se tiene que se puede expresar tambien en
terminos de las funciones Ni , mediante
( ) = x1 N1 ( ) + x2 N2 ( ) + x3 N3 ( ):

K,

Sea ahora r = 2; las funciones v 2 Sh;2 son cuadraticas sobre cada triangulo

vjT (x; y ) = a + bx + cy + dx2 + exy + fy 2:

120

Hay por tanto seis parametros; para que una funcion quede completamente
determinada sobre el triangulo habra que especicar su valor en seis puntos.
Ademas de los vertices habra que considerar otros tres nodos. Se suelen tomar
los puntos medios:
x3

x5

x6

x1

x2
x4

La numeracion estandar de los nodos es la de la gura: primero los tres

vertices en sentido positivo, y luego, en el mismo sentido, los tres puntos
medios, empezando por el del primer lado (el de los vertices primero y segundo).
En este caso habra seis funciones base nodales en el elemento de referencia,
Ni , i = 1; : : : ; 6. La aplicacion
2
6

'(x; y ) = 6
4

'1 (x; y )
..
.
'6 (x; y )

3
7
7
5

de las funciones base correspondientes a nodos del triangulo K se expresa

tambien en este caso como ' = N 1 , cuando nos restringimos a dicho
triangulo, donde
2
3
N1 (;  )
6
7
..
7:
N (;  ) = 6
.
4
5
N6 (;  )
Al ir tomando r mas grande hara falta tomar cada vez mas nodos. Por
ejemplo, si r = 3 hay diez parametros en cada triangulo, y habra que especicar
los valores de la funcion en diez nodos; se toman los de la gura:
121

. .......... 7
.
.
.
.
.
.
.
8 ....
... ........ . . . . ........
. .10. . ...... ...
..... . . . ....... ......5.
... ...
.4

6
2

La numeracion estandar es la de la gura. Observese que cada arista queda

dividida en tres partes iguales, y que el baricentro del triangulo tambien es un
nodo.
Para nalizar esta seccion estudiamos la convergencia en el caso r = 1.
Para K 2 Th denimos

hK = diametro de K;
K = diametro de la circunferencia inscrita en K;
h = Kmax
h :
2T K
h

En realidad no vamos a usar una triangulacion, sino una familia de triangulaciones fThg, caracterizadas por el parametro h. Supondremos que existe
una constante positiva  independiente de la triangulacion Th 2 fTh g, esto es,
independiente de h, tal que
K
  8K 2 Th:
hK
Esta condicion signica que no se permite que los angulos de los triangulos se
hagan arbitrariamente peque~nos; la constante  es una medida del angulo mas
peque~no en cualquier K 2 Th para cualquier Th 2 fTh g.
Sean xi , i = 1; : : : ; M , los nodos de Th . Dada una funcion u 2 C (
),
denimos el interpolante lineal a trozos uI como la unica funcion de Sh;1 tal
que
u(xi ) = uI (xi )
8i = 1; : : : ; M:
Sea

H 2 (
) = fu 2 H 1 (
) :

@u
@xi
122

2 H 1(
); i = 1; : : : ; dg:

Se tiene la siguiente estimacion para el error de interpolacion.

Lema 4.23. Si u 2 H 2 (
), entonces

ku uI kL (
)  Ch2
2

ku0 u0I kL (
)  Ch

2

!1=2
d Z 
2 u 2
X
@


;
 @x @x 
i
j

i;j =1
!1=2
d Z 
2 u 2
X
@


;
 @x @x 
i
j

i;j =1

y por tanto


ku uI kH (
)  h C1
1

h2

C
+ 22


1=2

!1=2
d Z 
2 u 2
X
@


:
 @x @x 
i
j

i;j =1

Se puede encontrar una prueba de este lema (con hipotesis adicionales sobre
u) en [J].
Observacion. Los experimentos numericos muestran que la estimacion da
la dependencia de  correcta. Esto tiene consecuencias practicas importantes.
Cuanto menor sea  , peor sera la aproximacion. Habra que intentar, por tanto,
que la triangulacion no contenga triangulos muy \delgados".

Gracias al lema se tiene convergencia de las soluciones del metodo de elementos nitos a la verdadera solucion, siempre y cuando esta tenga suciente
regularidad; concretamente se necesita que u 2 H 2 (
). Esta regularidad se
tiene para muchos problemas elpticos si el dominio
es convexo.

4.8 Calculos con elementos nitos. Programacion

En esta seccion concretaremos como se usan todos los conceptos anteriores a
la hora de hacer calculos con el ordenador. Nos limitaremos al caso de dos
dimensiones y funciones lineales sobre triangulos. Ilustraremos las distintas
ideas que vayamos introduciendo por medio de un problema modelo. Entre
123

otras cosas indicaremos instrucciones (o conjuntos de instrucciones) de MatLab

que realicen las diversas tareas que vayamos explicando. No intentaremos hacer
el mejor programa posible, sino uno que permita ver claramente cuales son las
ideas del metodo de elementos nitos. Concluiremos la seccion explicando
sucintamente que hay que hacer para tratar otros problemas.
Sea
el dominio de la gura,
D

La frontera de
esta descompuesta en dos partes disjuntas, D y N .
La primera de ellas aparece marcada en la gura por medio de un trazo mas
grueso.
Nuestro problema modelo es el siguiente:
Dados f 2 L2 (
) y g 2 L2 (
la solucion debil del problema
8
>
>
<
>
>
:

N ),

encontrar una aproximacion numerica de

u + u = f
u = 0;
@u
@n = g

en
;
en D ;
en N ;

(4.24)

donde   0 es un parametro real. Si  = 0 recuperamos la ecuacion de

Poisson, que es el caso que trataremos en los experimentos numericos.
Segun hemos visto, una solucion debil de (4.24) es una funcion u de

V = fv 2 H 1 (
) : v = 0 c.t.p. en
tal que

a(u; v ) = L(v )
124

8v 2 V;

donde

a(u; v ) =

(ru
rv + uv );

L(v ) =

fv +

Z
N

gv:

Para resolver el problema de forma aproximada lo primero que tenemos

que hacer es triangular el dominio. Veamos una posible triangulacion.
12

13
14

12
13

15

14
5

11

11
9

10

1
15

16

17

10

4
3

5
4

En esta gura hemos numerado los nodos y los triangulos; los ndices de
los triangulos estan recuadrados.
Hay un total de M = 15 nodos,

x 1 ; : : : ; xM ;
y nt = 17 triangulos,

K1 ;


; Knt :

La base nodal consta as de M funciones, 'i , i = 1; : : : ; M , que verican que

'i (xj ) = ij :
La aproximacion de elementos nitos es una funcion uh 2 Vh  V ,

Vh = fv 2 C (
) : v (x) = 0 si x 2

D;

125

vjKj

2 P1 (Kj ); j = 1; : : : ; nt g;

tal que

8v 2 Vh:

a(uh; v ) = L(v )

Cualquier v 2 Vh se puede representar en la forma

v=

M
X
j =1

v (xj )'j :

Ahora bien, si v 2 Vh, entonces v (xj ) = 0 para xj 2

v=

X
j 62ID

v (xj )'j ;

D,

ID = fj : xj 2

y por tanto
D

g:

As que podemos tomar como base de Vh las funciones 'i , i = 1; : : : ; M con
i 62 ID .
Denotando uj = uh(xj ) tenemos que la solucion de elementos nitos viene
dada por la solucion del sistema lineal
X
j 62ID

a('i ; 'j )uj = L('i );

i 62 ID :

(4.25)

Las incognitas son los valores nodales uj con j 62 ID .

Deberamos por tanto construir la matriz de rigidez de elementos

i; j 62 ID :

a('i ; 'j );

Sin embargo, como veremos mas adelante, la matriz que se puede construir de
forma sistematica y sencilla es la matriz \completa",

A = (aij );

i; j = 1; : : : ; M;

aij = a('i ; 'j ):

As pues, lo razonable sera construir la matriz completa A y extraer de ella

la matriz de coecientes del sistema (4.25) eliminando las las y columnas
de A correspondientes a ndices de ID . Ahora bien, para el ordenador esta
extraccion es extraordinariamente costosa, pues obliga a trasvasar datos de un
126

sector de memoria a otro. Por ello, lo que se hace en la practica es completar

las ecuaciones (4.25) a~nadiendo las ecuaciones

j 2 ID :

uj = 0;

(4.26)

Denotemos por D0 a la matriz diagonal M  M cuyas componentes diagonales dii son 1 si i 62 ID y 0 si i 2 ID . Entonces para un vector v 2 R M , el
producto D0 v es simplemente hacer 0 las componentes de v correspondientes
a ndices de ID .
Sea uh = (u1 ; : : : ; uM )T . El sistema de ecuaciones (4.25){(4.26) se escribe
como
A0 uh = D0 b;
(4.27)
donde A0 = D0 AD0 + (I D0 ), e I es la matriz identidad de R M , esto es,
A0 es la matriz que se obtiene al hacer cero las las y las columnas de A
correspondientes a ndices de ID y a continuacion dando valor 1 a los elementos
de la diagonal principal que correspondan a columnas cuyos ndices esten en
ID .
Teniendo en cuenta todo esto, descompondremos la programacion en las
siguientes tareas:






4.8.1

Suministrar los datos al programa, en particular la triangulacion.

Calcular la matriz de rigidez.
Calcular el vector de carga.
Construir el sistema reducido (4.27) a partir del sistema completo, para
tomar en cuenta las condiciones de Dirichlet.
Resolver el sistema resultante.
Datos de la triangulaci
on

>Como se almacena la triangulacion en el ordenador? La idea es muy simple.

Se utilizan dos matrices. La primera de ellas, llamemosla T, es una matriz de
127

elementos enteros que tiene una columna por cada triangulo, y en la columna
j -esima tres componentes que guardan los ndices de los vertices del triangulo
j -esimo.
En la triangulacion de nuestro ejemplo vemos que los ndices de los vertices
del primer triangulo son 2, 5 y 6, luego las tres componentes de la primera
columna de T seran
2 3
2
6 7
6 5 7:
4 5
6
Para el segundo triangulo los ndices de sus vertices son 2, 6 y 7, luego las tres
componentes de las dos primeras columnas de T seran
2
6
6
4

2 2
5 6
6 7

7
7:
5

Siguiendo con este proceso llegamos a que la matriz T de la triangulacion de

nuestro ejemplo es
2
2 2 2 4 4
6
4 5 6 7 7 8

4 4 3 3 3 1 1 1 1 1 1 3
9 10 10 11 12 12 13 14 15 5 2 2 7
5:
6 7 4 8 9 10 3 11 12 1 13 14 15 5 2 3 4

Hay varios aspectos que se~nalar. En cada triangulo hemos colocado los
ndices ordenados en sentido positivo (antihorario); esta ordenacion no tiene
por que coincidir con el orden creciente de los nodos. Por ejemplo, la tercera
columna de T es
2 3
2
6 7
6 7 7;
4 5
4
en vez de [2; 4; 7]T . Observese tambien que para los triangulos que tienen un
lado en la frontera, este es siempre el que corresponde al lado cuyos extremos
son el segundo y el tercer vertice. Por otra parte, no hay triangulos con dos
128

lados de frontera; as el primer vertice de cada triangulo siempre es interior a

.
No todas estas convenciones son habituales. Sin embargo, si se mantienen
se pueden utilizar provechosamente para simplicar la programacion.
La segunda matriz necesaria para describir la triangulacion es una matriz
que denotamos por Z, que tiene una columna por cada nodo, y que en la
columna j -esima guarda las coordenadas x; y del nodo j -esimo. En el ejemplo
de la gura la matriz Z sera
"

Z=

x1 x2
y1 y2

x15 :



y15

Con estas dos matrices queda perfectamente descrita la triangulacion, al

menos desde el punto de vista de los calculos.
Tambien tenemos que almacenar de alguna forma los datos de frontera, en
nuestro caso los lados en los que se tienen datos de tipo Dirichlet (los lados
que pertenecen a D ) y los lados en los que se tienen datos de tipo Neumann
(los que pertenecen a N ). Esta informacion se guarda en dos vectores, que
llamaremos ID e IN.
El vector ID (resp. IN) contiene los ndices de los triangulos que tienen un
lado en D (resp. N ).
En nuestro ejemplo las matrices seran
ID =

5 6 8 9 ;

IN =

1 2 4 11 12 13 14 :

Conociendo el triangulo al que pertenece una arista, se conocen sus nodos

de frontera gracias al convenio que hemos seguido: si el triangulo es el j-esimo,
los nodos de frontera seran T(j,2) y T(j,3).
Una vez construidas las matrices se le pasan al programa como argumentos.
MatLab. La primera lnea de nuestro programa de MatLab podra ser

129

function u=mef(T,Z,ID,IN,alpha,fsm,dn)

Los argumentos alpha, fsm y dn se utilizan para pasar al programa el

valor del parametro  y los nombres de la funcion f del segundo miembro de
la ecuacion y de la funcion g del dato de Neumann.
El programa, que se almacenara en un chero con extension .m, devolvera
como resultado la matriz columna u, que contendra los valores nodales de la
solucion de elementos nitos.
4.8.2

Construcci
on de la matriz de rigidez

Veamos ahora como construir la matriz de rigidez A, de elementos

aij = a('i ; 'j );

i = 1; : : : ; M:

Como vimos, las funciones de la base nodal, 'i , admiten una expresion
sencilla cuando nos restringimos a cada triangulo. Esto nos lleva a nuestra
primera idea basica.

Calcular la matriz de rigidez como suma de contribuciones de cada

triangulo.
Para ello expresamos la forma bilineal a como suma de formas bilineales
restringidas a los triangulos:

a=

nt
X
j =1

aKj ;

donde aKj es la forma bilineal a restringida al triangulo Kj . En nuestro ejemplo,

Z
Z

a(u; v ) =

de manera que

(ru
rv + uv ) =

aKj (u; v ) =

Z
Kj

nt
X

j =1 Kj

(ru
rv + uv );

(ru
rv + uv );

130

En terminos de la matriz de rigidez la descomposicion se traduce en

A=

nt
X
j =1

AKj

donde las matrices AKj son las matrices de las formas bilineales aKj en la base
'1 ; : : : ; 'M , y son por tanto matrices M  M . Calcularemos las matrices AKj
para luego ensamblarlas para formar la matriz A.
Sobre cada triangulo Kj (de nodos xj , xj , xj ) solo hay 3 funciones base
no nulas, las correspondientes a los nodos del triangulo ('j , 'j , 'j ). Por
tanto, la matriz AKj tiene a lo sumo 32 elementos no nulos, que son los que se
encuentran en la interseccion de las las y columnas cuyos ndices corresponden
a los nodos del triangulo.
1

Los 32 elementos de la matriz AKj que pueden ser no nulos conforman una
matriz Ej dada por
2

Ej =

6
6
4

aKj ('j ; 'j ) aKj ('j ; 'j ) aKj ('j ; 'j )

aKj ('j ; 'j ) aKj ('j ; 'j ) aKj ('j ; 'j )
aKj ('j ; 'j ) aKj ('j ; 'j ) aKj ('j ; 'j )
1

7
7:
5

Observese que Ej es la matriz de la forma bilineal aKj restringida al espacio

generado por las funciones base que son no nulas en el triangulo Kj , y expresada
en esa base.
Las matriz Ej (de dimension 3  3) es la matriz elemental del triangulo
(elemento) Kj ; la correspondiente matriz AKj (de dimension M  M ) es la
contribucion de la matriz elemental Ej a la matriz A.
En nuestro ejemplo los unicos elementos que pueden ser no nulos de la
matriz AK son los que ocupan la posicion (i; j ) para i; j = 2; 5; 6, y la matriz
E1 correspondiente viene dada por
1

2
6

E1 = 6
4

aK ('2 ; '2 ) aK ('2 ; '5 ) aK ('2 ; '6 )

aK ('5 ; '2 ) aK ('5 ; '5 ) aK ('5 ; '6 )
aK ('6 ; '2 ) aK ('6 ; '5 ) aK ('6 ; '6 )
1

131

7
7:
5

Observacion. Los elementos de la matriz Ej no tienen por que coincidir con los
correspondientes elementos de la matriz A; si un lado [xi ; xj ] corresponde a dos
triangulos, el elemento de A correspondiente, a('i ; 'j ), recibira contribuciones
de las matrices Ej de ambos. As, en nuestro ejemplo,

a25 = a('2 ; '5 ) =

N
X
j =1

aKj ('2 ; '5 ) = aK ('2 ; '5 ) + aK ('2 ; '5 ):

1

15

Analogamente, si un nodo xi pertenece a varios triangulos, el elemento correspondiente a('i ; 'i ) tiene contribuciones de todos ellos.
Observacion. La matriz Ej no se calcula necesariamente con la misma ordenacion con que se suma a A. En nuestro ejemplo, en el caso del tercer triangulo
los programas de elementos nitos calcularan la matriz
2
6
6
4

aK ('2 ; '2 ) aK ('2 ; '7 ) aK ('2 ; '4 )

aK ('7 ; '2 ) aK ('7 ; '7 ) aK ('7 ; '4 )
aK ('4 ; '2 ) aK ('4 ; '7 ) aK ('4 ; '4 )
3

3
7
7;
5

pero, a la hora de sumar estos elementos a los de A se reordenan para obtener la

matriz E3 . Afortunadamente, usando MatLab no nos tendremos que preocupar
por estos detalles.
Veamos como se hace el calculo de Ej en nuestro ejemplo. Para no recargar la notacion, abandonamos la j . Los elementos de la matriz elemental E
correspondiente al triangulo K son

aK ('i ; 'k ) =

Z
K

('i;x'k;x + 'i;y 'k;y + 'i 'k );

i; k = 1; 2; 3:

En principio podramos calcular cada uno de estos nueve elementos por separado. Sin embargo:

 El calculo de la matriz E no se efectua componente a componente, sino

que se calcula toda la matriz a la vez.

132

Con el n de expresar E en una forma adecuada para este proposito denimos

2
3
'1
6
7
7;
'=6
'
2
4
5
'3
donde 'i , i = 1; 2; 3, es la funcion de la base nodal asociada al vertice xi ,
i = 1; 2; 3, del triangulo K . Esto nos permite escribir

E=

('x 'Tx + 'y 'Ty + ''T ):

Como hemos visto, las funciones base restringidas al triangulo K se pueden

expresar de forma muy sencilla en terminos de las funciones de base en el
elemento de referencia, mediante

' = N 1:
Esto motiva la siguiente idea.

 Referir las integrales, mediante un cambio de variables, al elemento de

referencia.

Aplicando la regla de la cadena

'x = N x + N x ;

'y = N y + N y ;

por lo que los productos 'x 'Tx , 'y 'Ty se expresan como

'x 'Tx = x2 N NT + xx (N NT + N NT ) + x2 N NT ;
'y 'Ty = y2 N NT + y y (N NT + N NT ) + y2 N NT ;
que sumados dan
'x 'Tx + 'y 'Ty = (x2 + y2 )N NT +(xx + y y )(N NT + N NT )+(x2 + y2 )N NT :
Para calcular x; y ; x ; y observamos que son precisamente las componentes de D 1 . Recordando que
"

( ) =

x
y

"

x1
y1

"

x
+ 2
y2
133

x1 x3
y1 y3

x1
y1

#"


;


tenemos que

"

D =

x x
y y

"

x2
y2

x1 x3
y1 y3

x1
;
y1

y de aqu que
"

D 1 =

x y
x y

1
=
J

"

y3 y1
(y2 y1 )

(x3 x1 )
;
x2 x1

donde J es el jacobiano de la aplicacion ,

J = (x2

x1 )(y3

y1 ) (x3

x1 )(y2

y1 ):

Vemos que tanto J como x, y , x y y son constantes.

Deniendo
R

S1 = K N NT ;
R
S3 = K N NT ;
0
0

tenemos que

S2 = K (N NT + N NT );

R
S4 = K NN T ;
0
0

E = aS1 + bS2 + cS3 + dS4 ;

donde

y1 )2 + (x3 x1 )2
;
jJ j
(y3 y1 )(y2 y1 ) + (x3
b = jJ j(xx + y y ) =
jJ j
(y2 y1 )2 + (x2 x1 )2
2
2
c = jJ j(x + y ) =
;
jJ j
d = jJ j:
a = jJ j(x2 + y2 ) =

(y3

x1 )(x2

x1 )

Observacion. Si se han numerado los nodos con la orientacion positiva (antihoraria), la aplicacion mantiene la orientacion, y por tanto J > 0.

Es decir, para cada triangulo K , la matriz elemental de aK (u; v ) es combinacion lineal de las mismas cuatro matrices Si , i = 1; : : : ; 4, con coecientes
a; b; c; d que dependen exclusivamente de las coordenadas x; y de los vertices
134

del triangulo K . Por tanto, para cada triangulo en una triangulacion, solo hay
que calcular los coecientes a, b, c y d, y con ellos combinar las matrices Si .
Estas matrices reciben el nombre de matrices basicas.
Las matrices basicas son
2
6

S1 = 21 6
4
2
6

S3 = 12 6
4

1
1
0

1 0
1 0
0 0

1 0
0 0
1 0

1
0
1

7
7;
5

S2 =

6
16
24
2

7
7;
5

S4 =

6
1 6
24 4

2
1
1

1
0
1

2 1 1
1 2 1
1 1 2

1
1
0

3
7
7;
5

3
7
7:
5

Observacion. Si en el funcional aparecen terminos cruzados, ux uy , es necesario

T , donde
disponer de las matrices S12 y S12
2

S12 =

Z
K0

N NT

1 0
16
6
= 4 1 0
2
0 0

1
1
0

3
7
7:
5

MatLab. Al principio del programa, despues del encabezamiento, habra que

introducir las matrices basicas. Por ejemplo, para S1 escribimos
S1=[1,-1,0;-1,1,0;0,0,0]/2;

De la misma forma se introducen las demas matrices, que denotaremos

dentro del programa por S2, S3 y S4.
Tambien habra que incluir un bucle que ensamble la matriz de rigidez. La
secuencia de instrucciones puede ser parecida a esta:
% reservamos espacio para la matriz de rigidez
M=size(Z,2); % n
umero de nodos
A=sparse(M,M); % creo la matriz de rigidez (dispersa)

135

% ensamblado de la matriz de rigidez

for j=1:size(T,2); % recorremos los tri
angulos
I=T(1:3,j); % obtenemos los v
ertices del tri
angulo
x=Z(1,I); % coord. x de los v
ertices
y=Z(2,I); % coord. y de los v
ertices
% calculamos el jacobiano
J= (x(2)-x(1))*(y(3)-y(1))-(x(3)-x(1))*(y(2)-y(1));
% calculamos los coeficientes a, b, c, d.

Para b utilizo

% el nombre bb para no confundirlo con el vector de carga

a=((y(3)-y(1))*(y(3)-y(1))+(x(3)-x(1))*(x(3)-x(1)))/J;
bb=-((y(3)-y(1))*(y(2)-y(1))+(x(3)-x(1))*(x(2)-x(1)))/J;
c=((y(2)-y(1))*(y(2)-y(1))+(x(2)-x(1))*(x(2)-x(1)))/J;
d=J;
Se=a*S1+bb*S2+c*S3+alpha*d*S4; % matriz elemental
% se suma la matriz elemental a las componentes
% de A que correspondan
A(I,I)=A(I,I)+Se;
end

4.8.3

Construcci
on del vector de carga

A continuacion veremos como construir el vector de carga

2

b=

L('1 )
..
.
L('M )

6
6
4

7
7:
5

fv +
2
3
'1
6 . 7
6 .. 7 ;
4
5

En nuestro problema modelo L(v ) =

=

'M

136

gv , por lo que deniendo

se tiene que

b = f h + gh ;

fh =

con

f ;

gh =

Z
N

g :

Calcularemos por separado la contribucion al vector de carga de cada uno de

los sumandos.
Contribucion de f . Las ideas son similares a las que se aplicaron al construir
la matriz de rigidez. El calculo se hace como suma de contribuciones de cada
triangulo:
N Z
X

fh =
Cada uno de los vectores

j =1 Kj

Z
Kj

f :

f

tiene todas sus componentes no nulas, salvo quizas las 3 componentes que
corresponden a las 3 funciones base no nulas sobre el elemento Kj . De esta
forma, el procedimiento para calcular f h consiste en para cada triangulo Kj ,
cuyos nodos tienen ndices j1 , j2 , j3 , calcular el vector de 3 componentes
2

fKj =

Z
Kj

6
f6
4

'j
'j
'j

1
2

7
7;
5

y sumar estos valores a las componentes j1 , j2 , j3 de f h (en la practica se

suman a las componentes de b).
Sea pues el triangulo K de vertices x1 , x2 , x3 . Queremos expresar la integral
sobre K como una integral sobre el triangulo de referencia K0 . Recordemos
que ' = N 1 , y que por tanto

' = N;
por lo que

fK =

Z
K

f' = jJ j

Z
K0

137

f ( ( ))N ( ) dd;

donde J es el determinante jacobiano de la transformacion , que es positivo

bajo nuestra hipotesis acerca de la orientacion de los nodos.
Normalmente estas integrales no pueden expresarse en terminos de funciones elementales, por lo que es necesario recurrir a formulas de cuadratura
numerica. En general usamos una formula de s nodos

c1 ; : : : ; cs ;
y sus correspondientes pesos, que agrupamos en el vector

w = [w1 ; : : : ; ws ]T :
Con esta formula,
Z

2
6

K0

fg  f (cj )wj g (cj ) = [f (c1 ); : : : ; f (cs )] 6

4

w1

32

...

ws

76
76
54

g (c1 )
..
.
g (cs )

3
7
7:
5

Denotando entonces por N la matriz de dimension 3  s cuyas columnas son

N (cj ), esto es,
2
6

N = [N (c1 ); : : : ; N (cs )] = 6
4

N1 (c1 ) : : : N1 (cs)
N2 (c1 ) : : : N2 (cs)
N3 (c1 ) : : : N3 (cs)

7
7;
5

y por fq al vector cuyas componentes son los valores de f en los transformados

(cj ) de los nodos de cuadratura, esto es
2
6

fq = 6
4

f ( (c1 ))
..
.
f ( (cs ))

3
7
7
5

podemos pues aproximar

fK = J

K0

N ( )f ( ( )) dd  JN diag(w)fq ;
138

(4.28)

donde diag(w) es la matriz diagonal cuyos elementos diagonales son las componentes del vector w.
Los transformados de los nodos de cuadratura se calculan usando que
( ) = x1 N1 ( ) + x2 N2 ( ) + x3 N3 ( );
de donde

2
6

(ci ) = [x1 x2 x3 ] 6
4
es decir,

N1 (ci )
N2 (ci )
N3 (ci )

3
7
7;
5

[ (c1 ) : : : (cs )] = [x1 x2 x3 ]N:

MatLab. Tenemos un procedimiento sistematico:

1. Se calcula la matriz N .
2. Para cada triangulo Kj , j=1,. . . ,N,
(a) se transforman los nodos de cuadratura;
(b) se evalua f en los transformados de los nodos de cuadratura para
formar el vector fq ;
(c) se aproxima fKj mediante (4.28)
Usaremos la formula de cuadratura gaussiana denida por los siguientes
coecientes:

i
1
2
3
4
5
6
7

i
1/3
p
(6 + 15)=21
p
(9 2 15)=21
p
(6 + 15)=21
p
(6
15)=21
p
(9 + 2 15)=21
p
(6
15)=21

i
1/3
p
(6 + 15)=2
p
(6 + 15)=2
p
(9 2 15)=21
p
(6
15)=21
p
(6
p15)=21
(9 + 2 15)=21
139

wi
0.1125
p
(155 + 15)=2400
p
(155 + 15)=2400
p
(155 + 15)=2400
p
(155
15)=2400
p
(155
p15)=2400
(155
15)=2400

Esta formula, que se usa con frecuencia, es exacta en el triangulo de referencia K0 para polinomios de grado menor o igual que 5.
En el programa principal tendremos que incluir un conjunto de instrucciones parecido al siguiente:
M=size(Z,2); % n
umero de nodos
b=zeros(M,1); % reservo sitio para el vector de carga
% pesos de la cuadratura
w=[0.1125,...
(155+sqrt(15))/2400,(155+sqrt(15))/2400,...
(155+sqrt(15))/2400,(155-sqrt(15))/2400,...
(155-sqrt(15))/2400,(155-sqrt(15))/2400];
% matriz Nc:

valores de N en los nodos de la cuadratura

Nc=[1/3,-(96+23*sqrt(15))/42,...
-(102+17*sqrt(15))/42,(6+sqrt(15))/21,...
(9+2*sqrt(15))/21,(6-sqrt(15))/21,(6-sqrt(15))/21;
1/3,(6+sqrt(15))/21,(9-2*sqrt(15))/21,(6+sqrt(15))/21,...
(6-sqrt(15))/21,(9+2*sqrt(15))/21,(6-sqrt(15))/21;
1/3,(6+sqrt(15))/2,(6+sqrt(15))/2,(9-2*sqrt(15))/21,...
(6-sqrt(15))/21,(6-sqrt(15))/21,(9+2*sqrt(15))/21];
% ensamblado del vector de carga
for j=1:size(T,2); % recorremos los tri
angulos
I=T(1:3,j); % obtenemos los nodos del tri
angulo
x=Z(1,I); % coord. x de los v
ertices
y=Z(2,I); % coord. y de los v
ertices
% calculamos el jacobiano
J= (x(2)-x(1))*(y(3)-y(1))-(x(3)-x(1))*(y(2)-y(1));
% imagenes por Psi de los nodos de cuadratura
Psicx=x*Nc; % coord. x de las im
agenes
Psicy=y*Nc; % coord. y de las im
agenes

140

% fq:

valor de f (denotada por fsm dentro del programa)

% en las transformadas de los nodos.

f se debe definir

% externamente y debe devolver un vector columna

fq=feval(fsm,Psicx,Psicy);
% sumamos la contribuci
on del tri
angulo a las
% componentes de b que correspondan
b(I)=b(I)+J*Nc*diag(w)*fq;
end

Como se ha indicado, la funcion f del lado derecho de la ecuacion se debe

denir externamente. Esto se hara por medio de un chero .m independiente.
Si, por ejemplo, el lado derecho es, f (x; y ) = xy entonces el chero podra ser
function func=fsm(x,y); % los argumentos son matrices
func=(x.*y)'; % la multiplicaci
on es componente a componente
% la salida debe ser un vector columna

Al llamar al programa principal habra que introducir como penultimo argumento el nombre del chero, que no tiene por que ser fsm, entre comillas.
Contribucion de los datos Neumann. A continuacion estudiamos la contribucion al vector de carga del termino procedente de los datos de Neumann, gh .
Escribiendo N como union de segmentos, N = l1 [ l2


[ lns , podemos
expresar g h como suma de contribuciones de los distintos segmentos,

gh =
Cada uno de los vectores

Z
N

g =
Z
lj

ns Z
X
j =1 lj

g :

g

tiene todas sus componentes nulas excepto quiza aquellas que correspondan a
los nodos xj , xj que estan sobre el segmento lj . As pues, para construir el
1

141

vector gh calculamos el vector de 2 componentes

"

'
glj = g j
'j
lj

1
2

para cada segmento lj , y sumamos los valores obtenidos a las componentes j1 ,

j2 de gh (en la practica se suman a las componentes j1 , j2 de b).
Abandonamos la j para no recargar la notacion. Sea pues un segmento l
de N , de extremos x1 , x2 , con funciones de base asociadas '1 , '2 . Denimos
"

'=

'1
:
'2

Parametrizando l mediante la aplicacion : [0; 1] 7! R 2 ,

( ) = x1 +  (x2
se tiene que

gl =

Z
l

g' ds = kx2

x1 k

Z 1
0

x1 );
g ( ( ))L( ) d;

(4.29)

donde L = [L1 ; L2 ]T se dene como L = ' . Notese que L1 y L2 son funciones

lineales de  , y que L1 (0) = L2 (1) = 1, L1 (1) = L2 (0) = 0. Por tanto,

L1 ( ) = 1 ;

L2 ( ) = :

Normalmente las integrales (4.29) no pueden expresarse en terminos de

funciones elementales, por lo que de nuevo es necesario recurrir a formulas de
cuadratura. Sean c1 ; : : : ; cs , los nodos de una de dichas formulas de cuadratura
en [0; 1], y w1 ; : : : ; ws los correspondientes pesos, que agrupamos en el vector

w = [w1 ; : : : ; ws ]T :
Denotando por L la matriz cuyas columnas son los valores de L en cada uno
de los nodos de cuadratura cj , esto es,

L = [L(c1 );


; L(cs )];

142

y por gq el vector

g ( (c1 ))
6
7
..
7;
gq = 6
.
4
5
g ( (cs))
tenemos entonces que el vector gl se aproxima como
gl = kx2

x1 k

Z 1
0

g ( ( ))L( ) d  kx2

x1 kL diag(w)gq :

(4.30)

Para calcular las transformadas (ci ) de los nodos de cuadratura observamos que ( ) = x1 L1 ( ) + x2 L2 ( ), de donde
[ (c1 ) : : : (cs )] = [x1 x2 ]L:
Una situacion bastante frecuente es aquella en que g es constante. En ese
caso

gl =

Z
l

g' ds = g kx2

x1 k

Z 1
0

L d = g kx2

x1 k

"

1
2
1
2

(4.31)

MatLab. Teniendo la matriz IN (con los lados correspondientes a datos

Neumann) y la matriz Z (con las coordenadas de cada nodo), el procedimiento
sistematico para sumar el vector g h al vector b es como sigue:

1. Se calcula la matriz L.
2. Para cada segmento lj , j = 1; : : : ; ns ,
(a) se transforman los nodos de cuadratura cj ;
(b) se evalua g en los transformados para formar el vector gq ;
(c) se aproxima glj mediante (4.30) (o (4.31) si g es constante).
Como regla de cuadratura usamos la de Simpson, con nodos y pesos

cj 0 1/2 1
wj 1/6 4/6 1/6
143

Las instrucciones de MatLab podran ser:

wborde=[1/6 4/6 1/6]; % pesos de la cuadratura
% les llamamos wborde para no confundirlos con los pesos w
% utilizados para calcular la contribucion de f.
Lc=[1,1/2,0;0,1/2,1]; % valores de L en nodos de cuadratura
% ensamblado del vector de carga
for j=1:length(IN); % recorremos los lados
I=T(2:3,IN(j)); % obtenemos los nodos del lado
x=Z(1,I); % coord. x de los v
ertices
y=Z(2,I); % coord. y de los v
ertices
% imagenes por gamma de los nodos de cuadratura
gammacx=x*Lc; % coordenada x de la transformada
gammacy=y*Lc; % coordenada y de la transformada
% gq:

valor de g (denotada por dn dentro del programa)

% en las transformadas de los nodos.

g se debe definir

% externamente, y debe devolver un vector columna

gq=feval(dn,gammacx,gammacy);
% longitud del intervalo
lg= sqrt((x(2)-x(1))*(x(2)-x(1))+(y(2)-y(1))*(y(2)-y(1)));
% sumamos la contribuci
on del segmento
% a las componentes de b que correspondan
b(I)=b(I)+lg*Lc*diag(wborde)*gq;
end

La funcion g tambien debera ser denida externamente por medio de un

chero .m independiente. Al utilizar el programa principal habra que pasarle
el nombre de dicho chero (entre comillas) como ultimo argumento.

144

4.8.4

Condiciones Dirichlet y resoluci

on del sistema lineal

Una vez construidos la matriz de rigidez A y el vector de carga b, tenemos que

formar a partir de ellos el sistema (4.27), que tiene en cuenta las condiciones
de Dirichlet homogeneas. Para ello:

Se hacen cero las las y columnas de A, y las las de b correspondientes

a ndices de ID.
Se coloca en los elementos diagonales de A correspondientes a ndices de
ID el valor 1.

MatLab.
% vector con los 
ndices de nodos con datos Dirichlet
I=reshape(T(2:3,ID),2*length(ID),1)
% hacemos ceros las filas y columnas de A de 
ndices de ID
A(I,:)=zeros(size(A(I,:))); % las filas
A(:,I)=zeros(size(A(:,I))); % las columnas
% creamos un vector que tenga unos en 
ndices de ID, y ceros
% en las dem
as posiciones.
v=zeros(M,1);
v(I)=1;
% unos en la diagonal en las posiciones de los 
ndices de ID
A=A+spdiags(v,0,M,M); % matriz diagonal dispersa
% ceros en las componentes de b de 
ndices de ID
b(I)=0;

La resolucion del sistema lineal

Au = b
es en general una tarea compleja, difcil de programar, pues hay que aprovechar
el hecho de que la matriz A es dispersa, lo que obliga a usar algoritmos especiales para ese tipo de matrices. Sin embargo, nosotros no nos tendremos que
preocupar por eso, pues el MatLab realizara la tarea por nosotros.
145

MatLab. El sistema se resuelve por medio de la instruccion

u=Anb

En [I] se describen de forma sencilla algunos metodos para resolver sistemas

lineales dispersos.
4.8.5

Resultados numericos

Uniendo adecuadamente los sectores de codigo que hemos ido suministrando,

se construye sin mucha dicultad un programa de MatLab que permite resolver nuestro problema modelo (y otros muchos). Veamos algunos resultados
numericos obtenidos con dicho programa.
Empezamos por asegurarnos de que el programa funciona correctamente.
Para ello resolvemos un problema cuya solucion conocemos. Concretamente
consideramos la funcion u(x; y ) = x(x2 x y 2 + y ), que es la solucion de
8
>
>
<
>
>
:

con
= (0; 1)  (0; 1),
 = 1, f = 2 4x y

g (x; y ) =

(x; y ) 2
;
(x; y ) 2 D ;
(x; y ) 2 N ;

u + u = f;
u(x; y ) = 0;
@u
@n = g;
D

= f(x; y ) 2 R 2 : x = 0; 0  y  1g,

y2 y
x

= @
n

D,

si x = 0; 0 < y < 1;
si y = 0 o y = 1; 0 < x < 1:

En la gura 4.1 mostramos las soluciones numericas correspondientes a dos

triangulaciones. La segunda triangulacion es un renamiento de la primera,
obtenida descomponiendo cada triangulo en cuatro, mediante la union de los
puntos medios de cada lado:

146

0.15

0.2

0.1

0.15
0.1

0.05

0.05
0

0.05
0.1

0.05

0.15
0.2
1

0.1
0.8

1
0.6

0.8

0.15

0.6

0.4
0.4

0.2

0.2
0

0.2

0.25
0.15
0.2
0.15
0.1
0.1
0.05
0.05

0
0.05

0.1
0.15

0.05

0.2
1
0.8

1
0.6

0.8
0.4

0.2

0.15

0.2
0

Figura 4.1

0.1

0.6

0.4

Ejemplo en un cuadrado. Soluciones numericas.

0.25

0.2
0.25
0.15

0.2
0.15

0.1

0.1
0.05

0.05

0
0.05

0
0.1
0.15
1

0.05
0.8

1
0.6

0.8
0.4

0.2

0.2
0

Figura 4.2

0.1

0.6

0.4

Ejemplo en un cuadrado. Solucion teorica.

147

800

700
900
800

600

700
600

500

500
400

400

300
200

300

100
0
60

200
50

80
40

100

60
30

40
20

20
10

Figura 4.3

Problema modelo. Supercie solucion.

55

800

50
700
45
600
40
500
35
400
30
300
25
200

20

100

15

10
10

20

30

Figura 4.4

40

50

60

70

80

Problema modelo. Representacion plana.

148

En la gura 4.2 se muestra la verdadera solucion. Comparando las dos

guras observamos que la solucion numerica se parece, a medida que vamos
renando la triangulacion, cada vez mas a la teorica, y que el parecido entre
ambas es muy alto.
A continuacion resolvemos nuestro problema modelo con  = 0, f = 4
p
y g = 0:6 j(100 x)(100 y )j. En este caso no conocemos la solucion
explcitamente, y tenemos que recurrir al calculo numerico si queremos tener
una idea cuantitativa de como es. En la gura 4.3 se representa la supercie
solucion. Este problema se puede ver como un modelo para una distribucion
estacionaria de temperaturas en un medio en el que hay una fuente de calor
(representada por f ), un ujo de calor en la frontera (representado por g ), y
una zona de la frontera, D , que se mantiene a cero grados. Los problemas
de este tipo admiten una representacion plana, mediante un codigo de colores,
muy conveniente. Damos un ejemplo en la gura 4.4.
4.8.6

Condiciones Dirichlet no homogeneas

Supongamos ahora que en lugar de condiciones Dirichlet homogeneas sobre la

parte D de la frontera tenemos condiciones Dirichlet mas generales,

u(x) = (x);

x2

D:

Lo que se hace es miminimizar el funcional en el convexo cerrado

Vh = fv 2 C (
) : vjKj

2 P1 (Kj ); v(xj ) = (xj ); j 2 ID g:

Denotando uj = uh(xj ) tenemos que la solucion de elementos nitos viene

dada por la solucion del sistema lineal
M
X
j =1

a('i ; 'j )uj = L('i );

149

i 62 ID :

Los valores nodales uj con j 2 ID no son incognitas, as que los pasamos al
segundo miembro, quedando que
X
j 62ID

aij uj = bi

X
j 2ID

aij uj ;

i 62 ID :

As pues, hay dos diferencias con el caso homogeneo: (i) hay que modicar el
vector de carga, y (ii) los valores en los nodos xj , j 2 ID no son cero, sino que
vienen dados por uj = (xj ). Por lo demas el procedimiento a seguir es como
en el caso homogeneo.
MatLab. Una vez construidos la matriz de rigidez completa A, de elementos
aij = a('i ; 'j ), i; j = 1; : : : ; M y el vector de carga b, de elementos bi = L('i ),
i = 1; : : : ; M , habra que hacer lo siguiente:







Se construye un vector u (no hay inconveniente en usar el mismo nombre

que para el vector solucion) cuyas componentes son uj = 0 si j 62ID,
uj = (xj ) si j 2ID.
Se resta a b el vector Au.
Se hacen cero las las y columnas de A correspondientes a ndices de ID.
Se coloca en los elementos diagonales de A correspondientes a ndices de
ID el valor 1.
Se colocan los valores de u de ndices de ID en las posiciones correspondientes del vector b.

El conjunto de instrucciones de MatLab podra ser:

% vector con los 
ndices de nodos con datos Dirichlet
I=reshape(T(2:3,ID),2*length(ID),1)
% construimos el vector u de datos Dirichlet
u=zeros(M,1);
u(I)=rho(Z(1,I),Z(2,I)); % rho se define externamente
% se modifica el vector de carga
b=b-A*u;

150

% hacemos ceros las filas y columnas de A de 

ndices de ID
A(I,:)=zeros(size(A(I,:))); % las filas
A(:,I)=zeros(size(A(:,I))); % las columnas
% creamos un vector que tenga unos en 
ndices de ID, y ceros
% en las dem
as posiciones.
v=zeros(M,1);
v(I)=1;
% unos en la diagonal en las posiciones de los 
ndices de ID
A=A+spdiags(v,0,M,M); % matriz diagonal dispersa
% ceros en las componentes de b de 
ndices de ID
b(I)=u(I);

El sistema se resuelve como en el caso homogeneo.

4.8.7

Datos de Robin

Supongamos que estamos estudiando un problema con condiciones de Robin

en una parte R de la frontera, por ejemplo,
8
>
>
<
>
>
:

en
 R 2 ;
en R ;
en @
n R :

u + u = f
@u
@n + u = g
u=0

El espacio de Hilbert adecuado para tratar este problema es

V = fv 2 H 1 (
) : v = 0 c.t.p. en @
n

g:

Multiplicando la ecuacion por una funcion de V e integrando por partes se

tiene, tras imponer las condiciones de contorno, que
Z

(ru
rv + uv ) +

uv =

En la forma bilineal aparece una integral sobre

r(u; v ) =

151

uv;

fv +
R

gv:

que contribuira a la matriz de rigidez con una matriz R, de dimension M  M

con elementos
Z
rij =
'i 'j :
R

Escribiendo

como union de segmentos,

rij =

nr
X
k=1

blk ('i ; 'j );

= l1 [ l2 [


[ lnr , tenemos que

donde blk (u; v ) =

Z
lk

uv:

Para cada segmento lk todos los valores blk ('i ; 'j ) son cero excepto quizas
aquellos que correspondan a las funciones base 'k , 'k asociadas a los nodos
xk , xk del segmento lk . As, para calcular la contribucion del segmento lk a
la matriz R solo es necesario calcular la matriz
1

Rk =

Z "
lk

'k
'k

['k 'k ] ds;

1

y a~nadir sus valores a la submatriz de R que se obtiene al intersecar las las

y columnas de ndices k1 , k2 . Parametrizando l mediante : [0; 1] 7! R 2 ,
( ) = x1 +  (x2 x1 ), la matriz Rk se expresa como

Rk = kx2

x1 k

Z 1
0

LLT d;

donde L = ' . Por tanto, el procedimiento sistematico para sumar R a los

otros terminos de la matriz de rigidez, y obtener as A, consiste en calcular
primero
#
"
Z 1
2
1
1
;
G=
LLT d =
6 1 2
0
y despues, para cada lado lk , multiplicar esta matriz por la longitud del lado,
y sumar sus elementos a las correspondientes componentes de A.
MatLab. Habra que crear una matriz IR con los ndices de los triangulos
tienen lados en R . Las instrucciones para recoger las contribuciones de los
segmentos con condiciones de Robin podran ser las siguientes:

152

G=[2 1
1 2]/6;
for j=1:length(IR); % recorremos los tri
angulos
I=T(2:3,IR(j)); % obtenemos los v
ertices del triangulo
x=Z(1,I); % coord. x de los v
ertices
x=Z(2,I); % coord. y de los v
ertices
% calculamos la longitud del lado
lg= sqrt((x(2)-x(1))*(x(2)-x(1))+(y(2)-y(1))*(y(2)-y(1)));
Se=lg*G; % matriz elemental
% se suma la matriz elemental a las componentes
% de A que correspondan
A(I,I)=A(I,I)+Se;
end

4.8.8

Funciones polinomicas a trozos de orden mas alto

Hemos descrito como se programan los elementos nitos en el caso en que

se utilicen elementos triangulares y funciones continuas, lineales sobre cada
triangulo. En el caso de funciones continuas, cuadraticas o cubicas sobre cada
triangulo, las cosas se hacen esencialmente de la misma forma, aunque los
calculos son bastante mas complejos.
La contribucion E del triangulo K a la matriz de rigidez se calcula como
combinacion de matrices basicas,

E = aS1 + bS2 + cS3 + dS4 ;

igual que se hacia en el caso de funciones lineales, tomando los coecientes
a; b; c; d los mismos valores que en aquel caso. Las matrices basicas tambien se
denen igual,
R

S1 = K N NT ;
R
S3 = K N NT ;
0
0

S2 = K (N NT + N NT )

R
S4 = K NN T :
0
0

153

La diferencia estriba en que cambia el numero de nodos en cada triangulo y

en consecuencia el numero de componentes de N , que pasa a ser 6 en el caso
de funciones cuadraticas y 10 en el caso de funciones cubicas. De esta forma
las matrices elementales son matrices 6  6 en el caso de funciones cuadraticas
y 10  10 en el caso de funciones cubicas. De la misma forma, la matriz
N = [N (c1 ); : : : ; N (cs )] que aparece al construir el vector de carga tendra mas
las, no habiendo por lo demas cambios especialmente signicativos.
Las adaptaciones necesarias para tratar las condiciones de frontera son
tambien mas o menos obvias.

4.9 Comentarios y bibliografa

El analisis funcional que utilizamos en el captulo (espacios de Sobolev y teora
variacional de ecuaciones elpticas) esta muy bien explicado en [B].
Probablemente el texto de introduccion a los elementos nitos mas facil de
leer y con una crtica excelente sea [J]. Vino a sustituir al clasico \StrangFix" [SF], quizas el primer libro que explicaba los elementos nitos \para
matematicos", y que sigue resultando magnco como libro de consulta .
Otro dos libros muy accesibles son [BCO] y [CO1], que son los dos primeros
volumenes de una \enciclopedia" de seis peque~nos tomos dedicada a los elementos nitos. Los restantes volumenes, [CO2], [CO3], [OC1], [OC2], mas
especializados, tambien son recomendables. El ultimo volumen [CO3] trata
sobre elementos nitos aplicados a la mecanica de uidos. Este tema se trata
de una forma sencilla en [Pi].
Algo mas antiguo que [J], aunque con aplicaciones mas cuidadas y mas
completo es [MW]. El libro de Ciarlet [C] trata en profundidad el analisis de
la convergencia del metodo de elementos nitos. El libro de Zienkiewicz [Z],
ademas de ser un clasico, esta muy orientado a las aplicaciones, y es extraordi154

nariamente completo en la discusion de todo tipo de elementos para el calculo

de estructuras. Tiene ademas el atractivo del punto de vista originario de los
ingenieros. Hay una version abreviada y simplicada, [ZM]. Otros dos textos con buena crtica son [AB] y [BS]. Para captar la amplitud del tema se
recomienda [QV].
Sobre la programacion del metodo hay unas notas excelentes (no publicadas) del profesor Bosco Garca Archilla (U. Sevilla). El libro de Schwartz [Sc],
que es muy completo en este aspecto, es indispensable si se quiere de verdad
aprender a programar el metodo de elementos nitos.
Por falta de tiempo y por ser un tema excesivamente tecnico, no se trata en
el curso la resolucion de los sistemas lineales dispersos que aparecen al aplicar
el metodo de elementos nitos. Se puede encontrar una introduccion clara y
breve, tanto a los metodos directos como a los iterativos para la resolucion
de este tipo de sistemas en [I]. Como referencias mas avanzadas tenemos
[DER] y [GL] para metodos directos, y [Ax], [HY], [N] o el clasico [V] para
metodos iterativos. Dentro de los metodos iterativos se pueden utilizar tecnicas
multimalla. De nuevo [I] proporciona una breve introduccion. Hay varios
buenos libros introductorios sobre la materia, por ejemplo [BR], [H2] o [W].
Como referencia general para temas de algebra lineal numerica recomendamos el excelente libro [GV], y a un nivel mas basico el estupendo [TB].

155

Cap
tulo 5

Metodos lineales multipaso

5.1 Introduccion
Segun hemos visto en el captulo anterior, el metodo de Euler se puede obtener
partiendo de la aproximacion de la derivada y 0(x) por el cociente incremental

y (x + h) y (x)
:
h
Esta aproximacion es de orden 1. Una idea para obtener un metodo mejor es
usar una aproximacion de la derivada mas precisa, por ejemplo, la formula de
diferencias centradas
y (x + h) y (x h)
y 0(x) 
;
2h
que es una aproximacion de orden 2. Aplicandola en x = xn+1 se tiene que
y (x ) y (xn)
f (xn+1 ; y (xn+1 )) = y 0(xn+1 )  n+2
;
2h
lo que da lugar a la iteracion
y 0 (x) 

yn+2 = yn + 2hf (xn+1 ; yn+1 );

157

(5.1)

que se conoce como metodo leap-frog. Como vemos, para calcular yn+2 se
utilizan dos valores anteriores, yn e yn+1 , y no uno solo, como en los metodos
del captulo anterior. Diremos que es un metodo de dos pasos.
Este metodo se puede obtener tambien a partir de una formula de cuadratura. Para ello observamos que

y (xn+2 ) y (xn) =

Z xn+2
xn

y 0(x) dx;

(5.2)

y usamos la regla del punto medio para aproximar la integral, obteniendo que

y (xn+2) y (xn)  2hy 0(xn+1 );

lo que da lugar al metodo (5.1), que tambien se conoce como regla del punto
medio.
Si aproximamos la integral (5.2) usando la regla de Simpson obtenemos que

h
y (xn+2) y (xn)  (y 0 (xn ) + 4y 0(xn+1 ) + y 0(xn+2 )):
3
Esta aproximacion conduce al metodo numerico,
h
yn+2 yn = (f (xn ; yn) + 4f (xn+1 ; yn+1) + f (xn+2 ; yn+2));
3
conocido como regla de Simpson, que tambien es un metodo de dos pasos.
Tanto la regla del punto medio como la regla de Simpson son metodos lineales de dos pasos. En general, un metodo lineal de k-pasos obtiene las aproximaciones yn a la solucion del problema de valor inicial (1.1) como solucion de
la recurrencia
k
X
j =0

j yn+j = h

k
X
j =0

j f (xn+j ; yn+j );

n = 0; : : : ; N

k:

Las constantes j y j son los coecientes del metodo y satisfacen

k = 1;

j0j + j0j 6= 0:
158

(5.3)

El metodo se dice lineal porque yn+k es una combinacion lineal de yn+j y

f (xn+j ; yn+j ); y se dice de k pasos porque se necesitan k valores anteriores
yn ; : : : ; yn+k 1 para calcular yn+k .
La condicion k = 1 elimina la arbitrariedad que surge del hecho de que
podemos multiplicar ambos lados de (5.3) sin cambiar el metodo. La condicion
j0j+j0j 6= 0 elimina la posibilidad de que 0 y 0 sean ambas cero, excluyendo
as metodos como por ejemplo

yn+2

yn+1 = hf (xn+1 ; yn+1 )

que es esencialmente de 1 paso, y no de 2, y que se puede escribir como

yn+1

yn = hf (xn ; yn):

Notacion. Para abreviar se suele escribir fn = f (xn ; yn ).

Una forma alternativa de especicar un metodo lineal multipaso es dando

su primer y segundo polinomio caractersticos,  y  , que se denen mediante

( ) =

k
X
j =0

j  j ;

 ( ) =

donde  2 C es una variable muda.

k
X
j =0

j  j ;

Si k = 0, el metodo se dice explcito. En este caso se puede despejar

explcitamente cada yn+k en terminos de los k valores anteriores yn ; : : : ; yn+k 1.
Si se han almacenado los valores fn ; fn+1 ; : : : ; fn+k 2, entonces una unica evaluacion de funcion, fn+k 1, producira yn+k .
Si k 6= 0, el metodo es implcito. Para obtener yn+k en cada paso tendremos que resolver la ecuacion no lineal

yn+k = hk f (xn+k ; yn+k ) + g;

donde g es una funcion conocida de los valores ya calculados,

g=

k 1
X
j =0

(hj f (xn+j ; yn+j ) j yn+j ):

159

(5.4)

Si el lado derecho de la igualdad (5.4) tiene una constante de Lipschitz M con

respecto a yn+k tal que M < 1, entonces este sistema de ecuaciones tiene una
unica solucion yn+k , que se puede aproximar tanto como se desee por medio
de la iteracion
[ ]
yn[++1]
k = hk f (xn+k ; yn+k ) + g;

yn[0]+k arbitrario:

(5.5)

Por el Teorema de la Aplicacion Contractiva se cumple que

lim y [ ]
 !1 n+k

= yn+k :

Notese que siempre que f sea Lipschitz podemos forzar que se cumpla la condicion M < 1 tomando un valor de h sucientemente peque~no. Si este fuera
inaceptable se debera optar por otro metodo iterativo, como el de Newton.
Los metodos lineales multipaso presentan una dicultad computacional.
Para aplicarlos es necesario disponer de los k valores de arranque y0 ; : : : ; yk 1,
y sin embargo el problema de valor inicial (1.1) solo proporciona y0 . Si k 
2 habra que obtener los k 1 valores de arranque restantes por algun otro
procedimiento, como puede ser por desarrollo de Taylor, utilizando un metodo
de Runge-Kutta, etc. Estos valores de arranque no coinciden por tanto en
general con los valores de la solucion teorica; tendran un cierto error, y no hay
por que esperar que un metodo converja si no son sucientemente precisos.
Como ahora hay k valores de arranque, tenemos que modicar la denicion
de convergencia que se introdujo en el captulo anterior.

Denicion 5.1. Un metodo lineal de k pasos se dice convergente (respectivamente convergente de orden p) si para todo problema (1:1) con solucion
sucientemente regular se tiene que
lim max ky (xn) ynk = 0;

h!0+ knN

si lim max

h!0+ 0nk 1

ky(xn) ynk = 0;

(respectivamente
max ky (xn ) yn k = O(hp);

knN

si max

0nk 1

160

ky(xn) ynk = O(hp)

cuando h ! 0+ ).
Observacion. La condicion lim max ky (xn) yn k = 0 sobre los valores de
h!0 0nk 1
arranque es equivalente a pedir que lim yn = y (x0 ) para n = 0; : : : ; k 1.
+

h!0+

>Que condiciones deben satisfacer los coecientes de un metodo lineal multipaso (5.3) para que este converja? Si converge, >con que orden lo hace?
Intentaremos responder a estas preguntas en las siguientes secciones.

5.2 Consistencia
Queremos proceder como lo hacamos en el caso de los metodos de un paso.
La situacion ideal sera que la solucion y del problema de valor inicial (1.1)
satisciese las ecuaciones del metodo (5.3), pero en general esto no sera as, y
habra un residuo o error de truncacion

Rn =
=

k
P

(j y (xn+j )

hj f (xn+j ; y (xn+j ))

(j y (xn+j )

hj y 0(xn+j ));

j =0
k
P
j =0

n = 0; : : : ; N

k;

que, como en el caso de metodos de un paso, es lo que le falta a la solucion

teorica para satisfacer la ecuacion del metodo.

Denicion 5.2. Un metodo lineal multipaso se dice consistente de orden p si

para todo problema de valor inicial (1:1) con solucion y sucientemente regular
se tiene que
max kRn k = O(hp+1)
cuando h ! 0+.

0nN k

Determinar las condiciones de orden es mucho mas sencillo para un metodo

lineal multipaso que para un metodo de Runge-Kutta. En efecto, si tomamos
161

desarrollos de Taylor para y (xn+j ) = y (xn + jh) e y 0 (xn+j ) = y 0(xn + jh)

alrededor del punto x = xn ,

y (p+1) (xn )(jh)p+1

+ O(hp+2);
(p + 1)!
(p+1) (x )(jh)p
y
n
y 0 (xn+j ) = y 0 (xn ) + y 00 (xn )jh +


+
+ O(hp+1);
p!
y (xn+j ) = y (xn ) + y 0(xn )jh +


+

obtenemos que

Rn = C0 y (xn ) + C1 y 0 (xn )h +


+ Cp+1y (p+1) (xn )hp+1 + O(hp+2);

donde

C0 =
C1 =
..
.

k
X
j =0
k
X
j =0

1
Cq =
q!

j ;
j j
k
X
j =0

k
X
j =0

j j q

j ;

k
X
j =0

j j q

q > 1:

As pues, un metodo lineal multipaso sera consistente de orden p si y solo si

Cj = 0, j = 0; : : : ; p. El valor Cp+1 recibe el nombre de constante de error.
Ejemplo. El metodo

yn+2 + yn+1

h
2yn = (fn+2 + 8fn+1 + 3fn)
4

es un metodo lineal de dos pasos con

0 = 2; 1 = 1; 2 = 1;

0 = 3=4; 1 = 8=4; 2 = 1=4:

Un sencillo calculo muestra que C0 = C1 = C2 = C3 = 0 y que C4 = 1=4!. As

pues, el metodo es de orden 3 y la constante de error es C4 = 1=4!.
No es necesario efectuar el desarrollo de Taylor en x = xn . Se puede (y se
debe) efectuar en otro punto si hay posibilidad de aprovechar simetras.
162

Ejemplo. Consideramos la regla de Simpson,

h
yn+2 yn = (f (xn ; yn) + 4f (xn+1 ; yn+1) + f (xn+2 ; yn+2)):
3
Si desarrollamos alrededor del punto medio x = xn+1 queda un desarrollo en
potencias impares de h,

Rn = y (xn+2 ) y (xn ) h3 (y 0(xn ) + 4y 0(xn+1 ) + y 0(xn+2 )) =

= (2y 0(xn+1 )h + 13 y 000 (xn+1 )h3 + 5!2 y (v) (xn+1 )h5 +


)
2 (v)
h (6y 0 (x
000
2
4
n+1 ) + y (xn+1 )h + 4! y (xn+1 )h +


):
3
Concluimos que Rn = O(h5 ), y el metodo es por tanto consistente de orden 4.
Ejemplo. Consideramos la regla del trapecio
h
yn+1 yn = (fn+1 + fn ):
2
Si desarrollamos alrededor del punto medio x = xn+ := xn + h2 , de nuevo
queda un desarrollo en potencias impares de h,
1
2

Rn = y (xn+1) y (xn ) h2 (y 0(xn ) + y 0(xn+1 )) =

= (y 0(xn+ )h + 241 y 000(xn+ )h3 +


)
1 000
h (2y 0 (x
2
n+ ) + 4 y (xn+ )h +


):
2
1
2

1
2

1
2

1
2

Concluimos que Rn = O(h3 ), y el metodo es por tanto consistente de orden 2.

Lo que tienen en comun los dos ejemplos anteriores es que ambos metodos
son simetricos, en el sentido de que

j = k j ;

j = k j :

Es facil comprobar, desarrollando alrededor del punto medio, que los metodos
simetricos tienen orden par. As pues, si un metodo es simetrico solo sera
necesario considerar las condiciones de orden impar, pues las de orden par se
cumplen necesariamente.
El siguiente teorema, que caracteriza el orden en terminos de los polinomios
caractersticos puede resultar util.
163

Teorema 5.3. El metodo lineal multipaso (5:3) es consistente de orden p  1

si y solo si existe una constante c 6= 0 tal que cuando  ! 1
1)p+1 + O(j

( )  ( ) log  = c(

Prueba. Sea  = ew ; por tanto 
en serie de Taylor tenemos que

(ew ) w (ew ) =
=
=
Por tanto,

k
P

j ejw

j =0
k
X

j =0
1
X
q=0

j
1
q!

! 1 corresponde a w ! 0.

k
P

j ejw

j =0
1 q q!
X
jw

q=0
k
X
j =0

1jp+2):

q!

j j q wq

k
X
j =0
1
X

(q
q=1

1 q q
X
jw

q=0

1)!

Desarrollando

q!

k
X
j =0

j j q

wq :

(ew ) w (ew ) = cwp+1 + O(wp+2)

para algun c 6= 0 si y solo si el metodo es de orden p. Volviendo a la variable

 y usando que
log  = (

1) + O((

para  ! 1;

1)2 )

se concluye el resultado.
Ejemplo. Consideramos el llamado metodo de Adams-Bashforth de dos pasos,
h
yn+2 yn+1 = (3fn+1 fn ):
2
Al aplicar el teorema convendra expresar las cosas en terminos de la variable
 =  1. Tenemos que

( )  ( ) log  = ( 2  ) 32  12 log 

= ( + 2 ) 1 + 32  ( 12 2 + 13 3 +


)
= 125 3 + O(4);
el metodo es consistente de orden 2.
164

Como en el caso de los metodos de Runge-Kutta, la consistencia es necesaria

para la convergencia.

Teorema 5.4. Si un metodo es convergente, entonces es consistente de orden

al menos 1.
Prueba. Si es convergente, en particular lo es para el problema

y 0(x) = 0;

y (0) = 1;

cuya solucion es y (x) = 1. El metodo lineal multipaso (5.3) aplicado a este

problema produce la ecuacion en diferencias
k
X
j =0

j yn+j = 0;

n = 0; : : : ; N

k:

(5.6)

Tomamos valores de arranque yn = 1, n = 0; : : : ; k 1. Por ser el metodo

convergente tenemos que lim yn = 1, j = k; : : : ; N k. As pues, pasando al
h!0
lmite h ! 0+ en (5.6) obtenemos que
+

k
X
j =0

j = 0;

o lo que es lo mismo, C0 = 0.
Para ver que C1 tambien es cero consideramos el problema

y 0(x) = 1;

y (0) = 0;

cuya solucion es y (x) = x. El metodo, aplicado a este problema, se reduce a

k
X
j =0

j yn+j = h

k
X
j =0

j ;

que es una ecuacion en diferencias no homogenea. Es facil ver que una solucion
de esta ecuacion es la sucesion yn = nhM , n = 0; 1; : : : , donde
k
P

M=

j =0
k
P
j =0

165

jj

(Veremos en breve que para un metodo convergente

k
P
j =0

jj =
6 0, de manera

que M esta bien denido.) Si tomamos como valores de arranque yi (h) = ihM ,
i = 0; : : : ; k 1, entonces el metodo producira la sucesion yn = nhM = xn M .
Por tanto,
y (xn) yn = xn xn M = xn (1 M );
y para que haya convergencia se tiene que tener que M = 1, es decir
k
X
j =0

j =

k
X
j =0

jj ;

o lo que es lo mismo, C1 = 0.
Observacion. En terminos del primer y segundo polinomio caracterstico la
condicion de consistencia de orden al menos 1 se escribe

0 (1) =  (1):

(1) = 0;

5.3 0-estabilidad
La consistencia de orden 1 es necesaria para la convergencia. >Es tambien
suciente? Veamos que no con un ejemplo.
Consideramos el metodo

yn+2 + yn+1

2yn = h(5fn+1

2fn );

que es consistente de orden 1. Se lo aplicamos al problema mas simple que

uno pueda imaginar,
y 0 = 0;
y (0) = 0;
cuya solucion es y (x) = 0. El metodo en este caso queda

yn+2 + yn+1
166

2yn = 0:

La solucion general de esta ecuacion en diferencias se puede hallar facilmente.

Se buscan soluciones de la forma yn = n , y se tiene que  debe vericar que
2 +  2 = 0. As  debe ser 1 o 2, lo que da lugar a las soluciones yn = 1,
yn = ( 2)n . La solucion general del problema sera una combinacion lineal de
estas dos soluciones:

yn = c1 + c2 ( 2)n :
Consideramos la solucion

yn = h( 2)n :
Observese que

y0 = h ! 0 = y (x0); y1 = 2h ! 0 = y (x0 )
Sin embargo,

cuando h ! 0+ :

ky(xN ) yN k = h2N = b N a 2N ! 1:

As pues, al contrario de lo que sucede para metodos de un paso, en el caso

de metodos multipaso la consistencia de un metodo no garantiza su convergencia.
La divergencia se debe a la inestabilidad del metodo numerico, es decir,
de la ecuacion en diferencias. La solucion numerica con datos de arranque
y0 = y1 = 0 coincide con la teorica; sin embargo, una peque~na perturbacion
en los datos hace que la solucion numerica explote.
Para que un metodo numerico sea de alguna utilidad debera ser estable,
imitando la estabilidad del problema (1.1) al cual intenta aproximar (recordemos que, bajo las condiciones del teorema de existencia y unicidad, las soluciones de (1.1) dependen de forma continua de los datos).

Denicion 5.5. Un metodo lineal multipaso se dice 0-estable si para cada

problema (1:1) con f Lipschitz, y para cada dos sucesiones fungNn=0 y fvn gNn=0
167

satisfaciendo para 0  n  N
k
P
j =0
k
P
j =0

j un+j

j vn+j

se tiene que

max kun

vn k = O

knN

k
P
j =0
k
P
j =0

j f (xn+j ; un+j ) = hn ;
j f (xn+j ; vn+j ) = h n ;


max

0j k 1

kuj vj k + 0max
kj j k :
j N k

La 0-estabilidad es la version discreta de la dependencia continua de los

datos.

Teorema 5.6. Si un metodo lineal multipaso es consistente de orden p y 0estable, entonces es convergente de orden p.
Prueba. Por un lado tenemos que
k
X
j =0

j yn+j = h

k
X
j =0

j f (xn+j ; yn+j );

y por otro que

k
X
j =0

j y (xn+j ) = h

k
X
j =0

j f (xn+j ; y (xn+j )) + hn ;

donde n = Rn =h. Pero el metodo es 0-estable, as que

max ky (xn)

knN

yn k = O

max

0j k 1

ky(xj ) yj k + 0max
knk
nN k

cuando h ! 0+ . Como el metodo es consistente de orden p, se tiene que

knk = O(hp). Concluimos que
max ky (xn ) yn k = O(hp)

knN

si max

0j k 1

ky(xj ) yj k = O(hp);

es decir, el metodo es convergente de orden p.

Nos preguntamos ahora que condiciones debe satisfacer un metodo lineal
multipaso para ser 0-estable.
168

Denicion 5.7. Se dice que un metodo lineal multipaso satisface la condicion

de la raz si todas las races del primer polinomio caracterstico tienen modulo
menor o igual que 1, y aquellas que tienen modulo 1 son simples.
Observacion. Si un metodo es consistente, por el Teorema 5.3 una de las
 es la raz principal, a la que se etiqueta como
races de ( ) debe ser 1. Esta
1 . Las restantes races, i , i = 2; : : : ; k, son las llamadas races espurias; surgen
porque decidimos representar un sistema diferencial de primer orden por un
sistema en diferencias de orden k. Obviamente, los metodos lineales multipaso
con k = 1 no tienen races espurias, y satisfacen por tanto la condicion de la
raz.

Teorema 5.8. Un metodo lineal multipaso es 0-estable si y solo si satisface la

condicion de la raz.
Prueba. Consideramos el problema (1.1) con f = 0 y tomamos dos soluciones
un y vn de la ecuacion en diferencias sin perturbar
k
X
j =0

j un+j = 0:

(5.7)

Si el metodo es 0-estable, existe una constante C tal que

max kun

kj N

vn k  C max

0j k 1

kuj vj k:

Veamos que esto es imposible si no se satisface la condicion de la raz.

(i) Supongamos que existe  tal que j j > 1 y ( ) = 0. Las sucesiones
un = 0 y vn = h n son solucion de la ecuacion en diferencias (5.7). Tenemos
que
kuN vN k = hj jN = b N a j jN ! 1;
mientras que,
b a k 1
max kuj vj k = max hj jj 
j j :
0j k 1
0j k 1
N
169

As pues, el metodo no es 0-estable.

(ii) Supongamos que existe  tal que j j = 1, ( ) = 0 = 0 ( ). Las
p
sucesiones un = 0 y vn = hn n son solucion de la ecuacion en diferencias (5.7)
(notese que  es raz doble). Entonces

p
p p
kuN vN k = hN j jN = b a N ! 1;

mientras que
max

0j k 1

kuj vj k = 0max
j k 1

p
b a
hj j j  p (k
j

1);

y el metodo no es 0-estable.
Veamos ahora que si el metodo cumple la condicion de la raz, entonces es
0-estable. Damos la prueba (que se puede omitir en una primera lectura) en
el caso escalar. El caso vectorial es analogo.
La idea es escribir el metodo como metodo de un paso. Para ello denotamos

Yn = (yn+k 1; : : : ; yn)T ;
P
F (Yn ) = ( kj=0 j fn+j ; 0; : : : ; 0)T ;

n = 0; : : : ; N
n = 0; : : : ; N

k + 1;
k + 1:

Notese que estamos diciendo que

=

k
X
j =0

j f (xn+j ; yn+j )

es funcion solo de yn ; : : : ; yn+k 1, y no de yn+k (los valores xn se consideran

jos). Para verlo observamos que

=

k 1
X
j =0

j f (xn+j ; yn+j ) + k f xn+k ; h

k 1
X
j =0

j yn+j ;

ecuacion implcita que dene a  como funcion de yn ; : : : ; yn+k 1.

Con esta notacion podemos escribir el metodo como

Yn+1 = AYn + hF (Yn);

170

n = 0; : : : ; N

k;

donde A es la matriz compa~nera del primer polinomio caracterstico ( ) =

0 + 1  +


+ k 1  k 1 +  k ,
2

k
1
0
..
.
0

6
6
6
=6
6
6
6
4

k
0
1
...
:::

:::
:::
...
...
0

0
0

...

..
.
0

3
7
7
7
7
7:
7
7
5

Iterando el proceso tenemos que

Yn =

An Y

0+h

n 1
X
l=0

An

1 l F (Y ):
l

Para las un y vn de la denicion de estabilidad, llamando

n ; 0; : : : ; 0)T ;

Rn = (n

n = 0; : : : ; N

k;

y restando tendremos que

Un

Vn =

An (U

V0 ) + h

n 1
X
l=0

An 1 l (F (Ul )

F (Vl )) + h

n 1
X
l=0

An

1 lR

l:

Si las potencias de A estan todas ellas acotadas, esto es, si

kAnk  M;

n = 0; 1; : : : ;

tomando normas y utilizando que f es Lipschitz con constante L tenemos que

kUn Vnk  M kU0 V0k + hLMB

n 1
X
l=0

kUl Vl k + (b a)M 0max

kl lk;
lN k

siendo B es una cota de los modulos de los j , desigualdad que se puede escribir
como

n  A + hLMB
171

n 1
X
l=0

l ;

deniendo n y A de la manera obvia.

Sea

= A + hLMB

nP1
l=0

n+1

l . Se tiene que
n

= hLMBn  hLMB n ;

de donde concluimos que

n 

 (1 + hLMB )n

 e(n

1)hLMB (A +

kU0 V0 k);

y de aqu la 0-estabilidad.
Para que las potencias de A esten todas acotadas es necesario y suciente
que las races del polinomio  sean de modulo menor que la unidad, y aquellas
que sean de modulo 1 sean simples, es decir, que se cumpla el criterio de la
raz. Para ver esto no hay mas que observar que los autovalores de A son las
races de .
Finalmente nos preguntamos si es cierto el recproco del Teorema 5.6.

Teorema 5.9 (Teorema de equivalencia de Dahlquist). Un metodo lineal

multipaso es convergente si y solo si es consistente de orden al menos 1 y 0estable.
Prueba. Por los teoremas 5.4, 5.6 y 5.8, basta con probar que

convergencia ) criterio de la raz:

Supongamos que no se satisface el criterio de la raz. Consideramos el problema

y 0 = 0;

y (0) = 0;

cuya solucion es y (x) = 0.

(i) Supongamos que existe  tal que j j > 1 y ( ) = 0. La sucesion yn =
n = 0; 1; : : : , es solucion de la ecuacion en diferencias y tiene condiciones

h n ,

172

de arranque yn = h n, n = 0; : : : ; k
0; : : : ; k 1. Sin embargo,

1; por tanto yn

ky(xN ) yN k = hj jN = b N a j jN ! 1

0 = y (0), n =

cuando h ! 0+ ;

y el metodo no es convergente.
(ii) Supongamos que existe  tal que j j = 1, ( ) = 0 = 0 ( ). La sucesion
p
yn = hn n , n = 0; 1; : : : , es solucion de la ecuacion en diferencias y tiene
p
condiciones de arranque yn = hn n , n = 0; : : : ; k 1; por tanto yn ! 0 =
y (0), n = 0; : : : ; k 1. Sin embargo,

p
p p
ky(xN ) yN k = hN j jN = b a N ! 1

cuando h ! 0+ ;

y el metodo no es convergente.
>Cual es el mejor orden que podemos conseguir para un metodo, para un
numero de pasos k dado? Puesto que un metodo lineal de k pasos viene determinado por 2k + 1 coecientes, y las condiciones de orden de consistencia son
relaciones lineales de estos coecientes, existen metodos lineales de k pasos con
orden de consistencia p = 2k; reciben el nombre de maximales. Sin embargo
no sirven para nada, puesto que no son convergentes.

Teorema 5.10 (Primera barrera de Dahlquist (1956)). El orden de convergencia p de un metodo lineal de k pasos 0-estable satisface:

 p  k + 2 si k es par;
 p  k + 1 si k es impar;
 p  k si k =k  0 (en particular si es explcito).
La prueba se puede encontrar en [HNW].
Veremos en el proximo captulo que los metodos con p = k + 2, llamados
optimales, tampoco funcionan bien en la practica. As pues, el mejor orden que
se puede alcanzar con un metodo de k pasos que sirva para algo es p = k + 1.
173

5.4 Comentarios y bibliografa

La prueba del Teorema de Equivalencia y el establecimiento de la primera
barrera por Germund Dahlquist, en 1956 y 1959 respectivamente, fueron hitos
importantes en la historia del analisis numerico. No solo son resultados de
gran importancia intrnseca, sino que sirvieron para dar al analisis numerico
el rango de disciplina matematica de verdad, e impulsaron el rigor matematico
en el pensamiento numerico. La monografa clasica de Henrici [He], dedicada
a explicar los resultados de Dahlquist, sigue siendo un modelo de exposicion
clara y bonita.
Los libros [I] y [JR] tratan el tema de forma mucho mas sencilla (tambien
menos completa). Un poco mas avanzado resulta [L], donde se estudia el
Teorema de equivalencia (sin dar todas las pruebas) en una forma mas general,
aplicable a metodos multipaso no necesariamente lineales. Ese mismo punto
de vista se adopta en el tambien clasico [IK]. De nuevo es en [HNW] donde
podemos encontrar mas informacion. Ademas de pruebas completas (aunque
a veces diferentes) de todo el material que presentamos, contiene informacion
adicional sobre, por ejemplo, los metodos de Nordsieck (una clase de metodos
lineales multipaso especialmente adecuada para cambiar la longitud del paso)
o aspectos practicos de programacion.
En [L] se puede encontrar todo un captulo dedicado a los pares predictorcorrector, de los que aqu apenas hemos dicho lo que son.
Los metodos lineales multipaso de hoy en da no solo tienen paso variable,
sino tambien orden variable. Son los llamados algoritmos VSVO (Variable
Step Variable Order). En [L] se explica la estructura de este tipo de algoritmos. El lector interesado en conocer detalles sobre algoritmos concretos podra
encontrarlos por ejemplo en [G] y [SG].

174

Bibliografa
[A1]

Atkinson, K.E. (1989), An Introduction to Numerical Analysis (2nd

ed.), John Wiley & Sons, New York.
[A2] Atkinson, K.E. (1993), Elementary Numerical Analysis (2nd ed.),
John Wiley & Sons, New York.

[ABD] Aubanell, A.; Bensey, A.; Deshalms, A. (1993) Utiles
basicos del
calculo numerico, Labor, Barcelona.
[Ax] Axelsson, O. (1994), Iterative Solution Methods, Cambridge University
Press, Cambridge.
[AB] Axelsson, O.; Barker, V.A. (1984), Finite Element Solution of Boundary Value Problems: Theory and Computation, Academic Press, Orlando, FL.
[BCO] Becker, E.B.; Carey, G.F.; Oden, J.T. (1981), Finite elements. An
Introduction. Volume I, Prentice Hall, New Jersey.
[B]
Brezis, H. (1984), Analisis Funcional, Alianza Editorial, Madrid.
[BR] Briggs, W.L. (1987), A Multigrid Tutorial, SIAM, Philadelphia.
[BD] Boyce, W.E.; DiPrima, R.C. (1986), Elementary Dierential Equations and Boundary Value Problems (4th ed.), Wiley, New York.
[BS] Brenner, S.C.; Scott, L.R. (1994), The Mathematical Theory of Finite
Element Methods, Springer-Verlag, New-York.
[BF] Burden, R.L.; Faires, J.D. (1998), Analisis Numerico (6a ed.), International Thompson Editores, Madrid.
[Bu] Butcher, J.C. (1987), The Numerical Analysis of Ordinary Dierential
Equations, John Wiley, New York.
175

[CO1] Carey, G.F.; Oden, J.T. (1983), Finite elements: A Second Course.
Volume II, Prentice Hall, New Jersey.
[CO2] Carey, G.F.; Oden, J.T. (1984), Finite elements: Computational Aspects. Volume III, Prentice Hall, New Jersey.
[CO3] Carey, G.F.; Oden, J.T. (1986), Finite elements: Fluid Mechanics.
Volume VI, Prentice Hall, New Jersey.
[C]
Ciarlet, P.G. (1978), The Finite Element Method for Elliptic Problems,
North-Holland, Amsterdam.
[CB] Conte, S.D.; de Boor, C. (1981), Elementary Numerical Analysis. An
algorithmic approach (3rd ed.), McGraw-Hill, Singapore.
[CM] Crouzeix, M.; Mignot, A.L. (1986), Exercises d'analyse numerique des
equations dierentielles, Masson, Paris.
[DB] Dahlquist, G.; Bjork, 
A. (1974), Numerical Methods, Prentice Hall,
Englewood Clis, NJ.
[Da] Datta, B.N. (1995), Numerical Linear Algebra and Applications,
Brooks/Cole, Pacic Grove CA, USA.
[DV] Dekker, K.; Verwer, J.G. (1984), Stability of Runge-Kutta Methods for
Sti Nonlinear Dierential Equations, North-Holland, Amsterdam.
[DER] Du, I.S.; Erisman, A.M.; Reid, J.K. (1986), Direct Methods for Sparse
Matrices, Oxford University Press, Oxford.
[F]
Fletcher, R. (1987), Practical Methods of Optimization (2nd ed.), Wiley, London.
[G]
Gear, C.W. (1971), Numerical Initial Value Problems in Ordinary Differential Equations, Prentice-Hall, Englewood Clis, NJ.
[GL] George, A.; Liu, J.W.-H. (1981), Computer Solution of Large Sparse
Positive Denite Systems, Prentice-Hall, Englewood Clis, NJ.
[GT] Gilbarg, D.; Trudinger, N.S. (1983), Elliptic Partial Dierential Equations of the Second Order (2a ed.), Springer-Verlag.
[Go] Goldstine, H.H. (1977), A History of Numerical Analysis, SpringerVerlag, Berlin.
176

[GV]

Golub, G.H.; Van Loan, C.F. (1996) Matrix Computations (3rd ed.),
John Hopkins U. Press, Baltimore.
[H]
Hackbusch, W. (1992), Elliptic Dierential Equations: Theory and
Numerical Treatment, Springer-Verlag, Berlin.
[H2] Hackbusch, W. (1985), Multigrid Methods and Applications, SpringerVerlag, Berlin.
[HY] Hageman, L.A.; Young, D.M. (1981), Applied Iterative Methods, Academic Press, New York.
[HNW] Hairer, E.; Nrsett, S.P.; Wanner, G. (1991), Solving Ordinary Differential Equations I: Nonsti Problems (2nd ed.), Springer-Verlag,
Berlin.
[HW] Hairer, E.; Wanner, G. (1991), Solving Ordinary Dierential Equations II: Sti Problems and Dierential-algebraic Equations, SpringerVerlag, Berlin.
[HH] Hammerlin, G.; Homan, K.-H. (1991), Numerical Mathematics,
Springer-Verlag, Berlin.
[He] Henrici, P. (1962), Discrete Variable Methods in Ordinary Dierential
Equations, Wiley, New York.
[Heu] Heun, K. (1900), Neue Methode zur approximativen Integration der
Dierentialgleichungen einer unabhangigen Veranderlichen, Z. Math.
Phys. 45, 23{38.
[I]
Iserles, A. (1996), A First Course in the Numerical Analysis of Differential Equations, Cambridge University Press, Cambridge.
[IK]
Isaacson, E.; Keller, H.B. (1966), Analysis of Numerical Methods, John
Wiley & Sons, New York.
[J]
Johnson, C. (1987), Numerical Solution of Partial Dierential Equations by the Finite Element Method, Cambridge University Press,
Cambridge.
[JR] Johnson, L.W.; Riess, R.D. (1982), Numerical Analysis, Addison Wesley, London.
177

[K]

Kelley, C.T. (1995), Iterative Methods for Linear an Nonlinear Equations, SIAM, Philadelphia.
[KC] Kincaid, D.; Cheney, W. (1996), Numerical Analysis, Brooks/Cole,
Pacic Grove CA, USA.
[Ku] Kutta, W. (1901), Beitrag zur naherungsweisen Integration totaler Differentialgleichungen, Z. Math. Phys. 46, 435{453.
[L]
Lambert, J.D. (1991), Numerical Methods for Ordinary Dierential
Systems, Wiley, London.
[MF] Mathews, J.H.; Fink, K.D. (2000), Metodos Numericos con Matlab,
Prentice Hall, Madrid.
[MG] Mitchell, A.R.; Griths, D.F. (1980), The Finite Dierence Method
in Partial Dierential Equations, Wiley, London.
[MW] Mitchell, A.R.; Wait, R. (1977), The Finite Element Method in Partial
Dierential Equations, Wiley, London.
[M]
Moreno Gonzalez, C. (1999), Calculo Numerico II, Universidad Nacional de Educacion a Distancia, Madrid.
[N]
Nevanlinna, O. (1993), Convergence of Iterations for Linear Equations, Birkhauser, Basel.
[OC1] Oden, J.T.; Carey, G.F. (1983), Finite elements: Mathematical Aspects. Volume IV, Prentice Hall, New Jersey.
[OC2] Oden, J.T.; Carey, G.F. (1984), Finite elements: Special Problems in
Solid Mechanics. Volume V, Prentice Hall, New Jersey.
[OR] Ortega, J.M.; Rheinholdt, W.C. (1970), Iterative Solution of Nonlinear
Equations in Several Variables, Academic Press, New York.
[P]
Peral, I. (1995), Primer Curso de Ecuaciones en Derivadas Parciales,
Addison-Wesley/U. Autonoma de Madrid, Madrid.
[Pi]
Pironneau, O. (1989), Finite Element Methods for Fluids, John Wiley
& Sons, Chichester.
[Po]
Pozrikidis, C. (1988), Numerical Computation in Science and Engineering, Oxford University Press, Oxford.
178

[QV]
[RT]
[RM]
[R]
[S]

[SC]
[SG]
[Sc]
[Sc2]
[Si]
[Ste]
[SB]
[Str]
[SF]

Quarteroni, A.; Valli, A. (1997), Numerical Approximation of Partial

Dierential Equations (2nd ed.), Springer-Verlag, Berlin.
Rabier, P.; J.-M. Thomas (1993), Exercises d'analyse numerique des
equations aux derivees partielles, Masson, Paris.
Richtmyer, R.D.; Morton, K.W. (1967), Dierence Methods for InitialValue Problems, Interscience, New York.
Runge, C. (1895), Ueber die numerische Au osung von Dierentialgleichungen, Math. Ann. 46, 167{178.
Sanz-Serna, J.M. (1998), Diez lecciones de calculo numerico, Secretariado de Publicaciones e Intercambio Cientco, Universidad de Valladolid.
Sanz-Serna, J.M.; Calvo, M.P. (1994), Numerical Hamiltonian Problems, Chapman-Hall, London.
Shampine, L.F.; Gordon, M.K. (1975), Computer Solution of Ordinary
Dierential Equations, W.H. Freeman, San Francisco.
Schwartz, H.R. (1998), Finite Element Methods, Academic Press, London.
Schwartz, H.R. (1989), Numerical Analysis. A Comprehensive Introduction, John Wiley & Sons, Chichester.
Simmons, G.F. (1993), Ecuaciones Diferenciales con Aplicaciones y
Notas Historicas (2a ed.), McGraw-Hill, Madrid.
Stewart, G. (1996) Afternotes on Numerical Analysis, SIAM, Philadelphia.
Stoer, J.; Bulirsch, R. (1993), Introduction to Numerical Analysis (2nd
ed.), Springer-Verlag, Berlin.
Strang, G. (1982), Algebra lineal y sus aplicaciones, Fondo Educativo
Interamericano, Mexico.
Strang, G., Fix, G.J. (1973), An Analysis of the Finite Element
Method, Prentice-Hall, Englewood Clis, NJ.
179

[St]
[T]
[TB]
[V]
[W]
[Z]
[ZM]

Strikwerda, J.C. (1989), Finite Dierence Schemes and Partial Differential Equations, Brooks/Cole Publishing, Pacic Grove, CA.
Thomas, J.W. (1995), Numerical Partial Dierential Equations, Springer-Verlag, Berlin.
Trefethen, L.N., Bau, D. III (1997), Numerical Linear Algebra, SIAM,
Philadelphia.
Varga, R.S. (1962), Matrix Iterative Analysis, Prentice-Hall, Englewood Clis, NJ.
Wesseling, P. (1992), An Introduction to Multigrid Methods, Wiley,
Chichester.
Zienkiewicz, O.C. (1977), The nite Element Method (3rd ed.),
McGraw-Hill, London.
Zienkiewicz, O.C.; Morgan, K. (1983), Finite Elements and Approximation, John Wiley & Sons, New York.

180