Tema 7
Resoluci
on Num
erica de
Problemas de Contorno.
Indice
1. Introduccion
2. Existencia y unicidad de solucion
3. Metodo de Tiro
4. Metodo de Diferencias Finitas
5. Metodos de Galerkin y Ritz
Introducci
on
Un problema de contorno para una ecuacion diferencial es aquel en el que a la
variable dependiente o incognita se le exigen condiciones en dos o mas puntos
distintos. Nos centraremos en los problemas de contorno para ecuaciones
diferenciales lineales de segundo orden, es decir, problemas del tipo
y 00 = f (x, y, y 0 ), x [a, b],
(1)
siendo [a, b] R y f : [a, b] R2 R, con condiciones del tipo
y(a) = ,
y(b) = ,
(2)
a0 y(a) a1 y 0 (a) = ,
b0 y(b) + b1 y 0 (b) = ,
(3)
o, en general, del tipo
con |ao | + |a1 | 6= 0 y |b0 | + |b1 | 6= 0.
1
Del mismo modo que ocurre para problemas de valores iniciales, solo
una minora de ecuaciones lineales pueden resolverse mediante metodos que
aparecen en las obras dedicadas a su estudio. Por otro lado, no siempre es
posible explicitar una solucion que verifique ciertas condiciones del tipo (2)
o, en general, del tipo (3). Es por ello que se hace necesaria, en la practica
la aplicacion de metodos numericos de resolucion aproximada.
Entre los mas conocidos y utilizados se encuentran el metodo de tiro, el
de diferencias finitas, los metodos de Ritz y Galerkin, o los metodos tipo
spline entre los que destaca el metodo de los elementos finitos.
Existencia y unicidad de soluci
on
Consideremos el problema de contorno (1)(2). Nos proponemos estudiar,
en principio, bajo que condiciones podemos asegurar la existencia y unicidad
de solucion para este problema. Especial atencion merecera el caso en el que
la ecuacion diferencial (1) sea lineal.
Teorema 1 Sea D = [a, b] R2 y supongamos que f : D R verifica las
siguientes condiciones:
1. f ,
2.
f f
y
son continuas en D,
u v
f
> 0 para todo (x, u, v) D,
u
3. existe una constante M > 0 tal que
f
M, (x, u, v) D.
u
Entonces el problema (1)(2) tiene solucion u
nica.
En particular, como consecuencia de este teorema, se obtiene el siguiente
resultado para el caso en que la ecuacion diferencial (1) sea lineal, es decir,
si f (x, y, y 0 ) = p (x) y 0 + q (x) y + r (x), para x [a, b].
Corolario 2 Consideremos el problema de contorno
y 00 = p (x) y 0 + q (x) y + r (x) , x [a, b]
y (a) = , y (b) = .
(4)
Si se verifica que
1. p, q y r son continuas en [a, b] ,
2. q (x) > 0 para todo x [a, b],
entonces el problema (4) tiene solucion u
nica.
M
etodo de tiro
Vamos es estudiar un metodo de aproximacion de la solucion de un problema
de contorno lineal de la forma (4), o mas general, de la forma (1)(2), que
en las condiciones del corolario 2 o del teorema 1, respectivamente, existe y
es u
nica, mediante la resolucion de problemas de valores iniciales asociados
a los mismos.
Nos ocupamos en principio del problema lineal. Por tanto, consideremos
el problema de contorno (4) y supongamos que se verifican las condiciones
del corolario 2. Asociados a este problema, consideremos los problemas de
valores iniciales
y 00 = p (x) y 0 + q (x) y + r (x) , x [a, b]
y (a) = , y 0 (a) = 0.
y 00 = p (x) y 0 + q (x) y, x [a, b]
y (a) = 0, y 0 (a) = 1.
(5)
(6)
Bajo las condiciones del corolario 2, los problemas (5) y (6) admiten
soluciones u
nicas, definidas en [a, b].
Teorema 3 Supongamos que se verifican las hipotesis del corolario 2. Sean
y1 (x) e y2 (x) las u
nicas soluciones de los problemas (5) y (6), respectivamente. Entonces
y1 (b)
y (x) = y1 (x) +
y2 (x)
(7)
y2 (b)
3
es la u
nica solucion del problema (4).
El metodo de tiro para el problema de contorno (4) consiste, por tanto, en
hallar una aproximacion del problema de valores (5), otra del (6), y construir
la solucion y dada por (7), tal y como se indica en el siguiente algoritmo:
ALGORITMO DEL METODO
DE TIRO LINEAL
Entrar p, q, r, a, b, ,
Hallar y1 como una aproximacion de la solucion de (5)
mediante un metodo numerico para problemas de valores iniciales
Hallar y2 como una aproximacion de la solucion de (6)
mediante un metodo numerico para problemas de valores iniciales
y1 (b)
y = y1 +
y2
y2 (b)
Escribir y.
En los pasos segundo y tercero de este algoritmo se necesita incluir el algoritmo de alguno de los metodos de resolucion aproximada de problemas de
valores iniciales estudiados en el captulo anterior, tanto para el problema
(5) como para el (6). El error del metodo de tiro para problemas lineales
vendra determinado por el error del metodo de aproximacion empleado para
estos problemas de valores iniciales, ya que la construccion de la solucion del
primero es lineal con respecto a las soluciones de los otros dos.
Consideremos ahora el caso general de un problema de contorno de segundo orden, no necesariamente lineal. Es decir, se trata del problema de
contorno
y 00 = f (x, y, y 0 ) , x [a, b]
y (a) = , y (b) = .
(8)
En este caso, la solucion del problema no puede expresarse como una
combinacion lineal de las soluciones de dos problemas de valores iniciales,
debido a la no linealidad de f . Para hallarla hara falta realizar un proceso
iterativo.
4
Consideremos el problema de valores iniciales
y 00 = f (x, y, y 0 ) , x [a, b]
y (a) = , y 0 (a) = m.
(9)
para un cierto valor de m dado, que indicara la pendiente de salida de la
solucion del problema (de ah el nombre de metodo de tiro).
Sea y (x, m) la solucion de este problema. Coincidira con la solucion de
(8) si
y(b, m) = .
(10)
Por tanto, se trata de resolver (en la variable m) esta u
ltima ecuacion,
para lo cual se puede aplicar cualquier metodo numerico de los estudiados
en el captulo 2. Pero el empleo de un metodo u otro necesitara de un mayor
o menor conocimiento de la solucion de (9) y, quiza, de alguna derivada, lo
cual puede hacer complicada su aplicacion.
Uno de los metodos mas sencillos que se puede utilizar para resolver (10)
es el metodo de la secante, para lo que se necesita dar dos valores iniciales,
m0 y m1 , y hallar las siguientes aproximaciones mediante la iteracion
mk = mk1
(y (b, mk1 ) ) (mk1 mk2 )
, k 2.
y (b, mk1 ) y (b, mk2 )
(11)
A continuacion, conocido el valor de mk en la iteracion kesima, se vuelve
a resolver el problema (9) para m = mk y, mediante (11) para k = k + 1,
obtener un nuevo valor de la pendiente del tiro m = mk+1 , reiterando el
proceso hasta que se verifique un cierto criterio de parada.
Por tanto, el metodo de tiro, usando el metodo de la secante, consiste en
hallar una solucion aproximada de (8) mediante el siguiente algoritmo (en el
que la variable niter indica el n
umero de iteraciones que deseamos realizar):
ALGORITMO DEL METODO
DE TIRO NO LINEAL (SECANTE)
Entrar f , a, b, , , m0 , m1 , niter
Hallar y (x, m0 ) una aproximacion de la solucion
de (9) para m = m0
Hallar y (x, m1 ) como una aproximacion de la solucion
de (9) para m = m1
Para k = 2, . . . , niter
k1 ))(mk1 mk2 )
mk = mk1 (y(b,m
y(b,mk1 )y(b,mk2 )
Hallar y (x, mk ) como una aproximacion de la solucion
de (9) para m = mk
Escribir y (x, mniter ).
M
etodo de diferencias finitas
Los metodos de tiro, tanto para problemas lineales como no lineales, presentan problemas de inestabilidad. Para evitar esto se puede hacer uso de
un metodo de diferencias finitas como los que describiremos a continuacion,
que, aunque son menos precisos que el metodo de tiro, presentan mayor estabilidad.
Estos metodos se basan en la sustitucion de las derivadas que aparecen
en el problema de contorno por formulas de derivacion numerica adecuadas,
con el orden de error deseado.
Para ello, sea N > 0 y consideremos en el intervalo [a, b] una particion
uniforme en N subintervalos, es decir, la particion de nodos
xi = a +
ba
i, i = 0, . . . , N.
N
Sean y0 , . . ., yn valores reales que suponemos que aproximan a y (x0 ), . . .,
y (xN ), siendo y (x) la solucion de un problema de contorno dado. Entonces,
se sabe que, para i = 1, 2, . . . , N 1,
y (xi+1 ) y (xi1 )
+ O h2 ,
2h
2
y (xi+1 ) 2y (xi ) + y (xi1 )
y 00 (xi ) =
+
O
h ,
h2
y 0 (xi ) =
y, por tanto,
yi+1 yi1
,
2h
yi+1 2yi + yi1
y 00 (xi ) '
.
h2
y 0 (xi ) '
En el caso lineal, es decir, si consideramos el problema 4, efectuando las
aproximaciones antes indicadas, para cada uno de los nodos interiores de la
particion fijada, se obtienen las igualdades
y0 = ,
yi+1 2yi + yi1
yi+1 yi1
= p (xi )
+ q (xi ) yi + r (xi ) , i = 1, . . . , N 1,
2
h
2h
yN = .
que se pueden escribir en forma matricial como
Ay = b,
donde y = (y1 , y2 , . . . , yN 1 )T ,
2 + q1 h2 1 + h2 p1
1 h p2 2 + q2 h2 1 + h p2
2
2
..
..
..
.
.
.
A=
...
...
...
(12)
...
...
1 h2 pN 1
b=
r1 h2 + 1 + h2 p1
r2 h2
..
.
..
.
rN 2 h2
rN 1 h2 + 1 h2 pN 1
7
1 + h2 pN 2
2 + qN 1 h2
(13)
(14)
habiendo utilizado la notacion pi = p (xi ), qi = q (xi ) y ri = r (xi ), i =
1, 2, . . . N 1.
Por tanto, el algoritmo es el siguiente:
ALGORITMO DEL METODO
DE DIFERENCIAS FINITAS LINEAL
Entrar p (x), q (x), r (x), a, b, ,
A = (0)(N 1)(N 1) , b = (0)(N 1)1
Para i = 1, . . . , N 1,
aii = 2 + h2 q (xi )
2
bi = h
r (xi )
h
b1 = b1 + 1 + p (x1 )
2
bN 1 = bN 1 + 1 h2 p (xN 1 )
Para i = 1, . . . , N 2
h
ai+1,i = 1 p (xi )
2
h
ai,i+1 = 1 + p (xi+1 )
2
Resolver Ay = b
Escribir y.
Observese que la salida del algoritmo es el vector de incognitas (y1 , . . . , yN 1 ).
Por otra parte, la resolucion del sistema Ay = b se puede hacer mediante el
metodo de descomposicion LU para el caso de matrices tridiagonales.
Los siguientes resultados, cuya demostracion omitimos, establecen condiciones suficientes tanto para la unisolvencia de este sistema lineal, como para
la determinacion de un orden de convergencia cuadratico.
Teorema 4 Supongamos que p, q y r son continuas en [a, b] y que p (x) > 0,
2
para todo x [a, b]. Sea L = max |p (x)|. Si h <
entonces el sistema
x[a,b]
L
tridiagonal (12)-(14) tiene una u
nica solucion.
Teorema 5 Bajo las hipotesis del teorema anterior, sean y (x) la u
nica
solucion de (4) e y1 , . . . , yN 1 la u
nica solucion de (12)-(14). Si y C 4 ([a, b])
entonces
yi = y (xi ) + O h2 , i = 1, . . . , N 1.
8
5
5.1
M
etodos de Galerkin y Ritz
Introducci
on
Los metodos de Galerkin y Ritz son metodos numericos de resolucion aproximada de problemas de contorno que consisten, en esencia, en discretizar el
problema reemplazando el espacio de funciones de dimension infinita (normalmente un espacio de funciones con cierta regularidad) donde se plantea
el problema original por otros espacios de dimension finita, que estaran generados por funciones que verifican las condiciones de contorno.
Principalmente, los espacios de discretizacion son espacios de funciones
polinomicas (qmetodos), o funciones trigonometricas (metodos espectrales)
o espacios de funciones splines (lineales, cuadraticos o c
ubicos, fundamentalmente).
Se considera en el espacio de discretizacion una base de funciones que
sean ortogonales dos a dos respecto de cierto producto escalar definido en el
espacio de funciones considerado. La base de estos metodos, como veremos,
es una buena eleccion de una sucesion de funciones, llamadas funciones coordenadas, que constituya un sistema independientes de funciones del espacio
considerado.
En lo que sigue, continuamos con problemas de contorno lineales de segundo orden con condiciones de Dirichlet, Newmann o Mixtas.
5.2
Homogeneizaci
on de condiciones de contorno
Consideremos el Problema de Contorno:
Ly = f (x), y = y(x), x [a, b]
y(a) = ya , y(b) = yb ,
(15)
donde L es un operador diferencial lineal de segundo orden, es decir,
L(y)(x) = a0 (x)y 00 (x) + a1 (x)y 0 (x) + a2 (x)y(x),
con a0 , a1 , a2 C[a, b]. Mediante un cambio de funcion incognita transformamos el problema (15) en otro equivalente con condiciones homogeneas.
9
En concreto, sea
yb ya
z(x) = y(x) ya + (x a)
ba
Entonces, z(x) es solucion del problema
Lz = f(x), z = z(x), x [a, b]
z(a) = 0, z(b) = 0,
yb ya
yb ya
siendo f (x) = f (x)
a1 (x) ya + (x a)
a2 (x).
ba
ba
Para condiciones de tipo Neumann y 0 (a) = y a , y 0 (b) = y b , el cambio
1
2 yb ya
z(x) = y(x) y a (x a) + (x a)
2
ba
transforma en problema en uno homogeneo, ya que z 0 (a) = z 0 (b) = 0.
Para las condiciones mixtas y(a) = ya , y 0 (b) = y b , el cambio que transforma las condiciones en homogeneas puede ser
1
2 yb
z(x) = y(x) ya + (x a)
.
2
ba
5.3
Sistema de funciones coordenadas
Sea D = {y C 1 [a, b] : y(a) = 0, y(b) = 0} (para el caso de condiciones
de Dirichlet) y definamos en el espacio C 1 [a, b], y por tanto en cualquier
subespacio suyo, la distancia
d1 (y1 , y2 ) = max (|y1 (x) y2 (x)| + |y10 (x) y20 (x)|) ,
x[a,b]
para cualesquiera y1 , y2 C 1 [a, b].
Una familia de funciones {n }n1 D es un sistema de funciones coordenadas de D si
Para todo n 0, {1 , . . . , n } son linealmente independientes;
10
Para toda y D y todo > 0, existe n N y 1 , . . . , n R tales que
n
X
si =
i i entonces
i=1
d1 (y, ) < ,
es decir, cualquier funcion de D se puede aproximar tanto como se desee
por una combinacion lineal finita de funciones del sistema de funciones
coordenadas.
Se puede comprobar que las siguientes familias numerables de funciones
constituyen sistemas de funciones coordenadas en los casos que se indican:
Para condiciones de Dirichlet homogeneas, y(a) = y(b) = 0,
sistema de funciones coordenadas de tipo polinomico:
n (x) = (x a)(b x)xn1 , n 1;
sistema de funciones coordenadas de tipo trigonometrico
xa
n (x) = sen n
, n 1.
ba
Para condiciones de Neumann homogeneas, y 0 (a) = y 0 (b) = 0,
sistema de funciones coordenadas de tipo polinomico:
1 (x) = 1,
1
1
2 (x) = x3 + (a + b)x2 abx,
3
2
2
n (x) = (x a) (b x)2 xn3 , n 3;
sistema de funciones coordenadas de tipo trigonometrico
xa
n (x) = cos (n 1)
, n 1.
ba
Analogamente, se pueden calcular sistemas de funciones coordenadas para
otros tipos de condiciones, por ejemplos las mixtas.
11
5.4
M
etodo de Galerkin
Consiste en discretizar el problema de contorno dado en el espacio de dimension finita generado por las n primeras funciones de un sistema de funciones coordenadas dado adecuado a las condiciones del contorno de problema.
Por tanto, consideremos el problema de contorno dado por la ecuacion
diferencial
Ly = f (x), x [a, b],
con ciertas condiciones de contorno (de Dirichlet, Neumann o Mixtas) homogeneas.
Sea {n }n1 un sistema de funciones coordenadas de C 1 [a, b] verificando
las condiciones dadas y Xn el espacio generado por las n primeras funciones
del sistema considerado.
Si y D es solucion de Ly = f (x), x [a, b], entonces, para cualquier
D se verifica que
Z b
Z b
Ly(x)(x)dx =
f (x)(x)dx.
a
El Metodo de Galerkin consiste en discretizar este problema en Xn , es
decir, hallar yn Xn tal que
Z b
Z b
Lyn (x)(x)dx =
f (x)(x)dx,
a
para toda Xn .
Teniendo en cuenta que yn Xn entonces yn =
n
X
i i y, por linealidad,
i=1
el problema se reduce a calcular 1 , . . . , n R tales que
!
Z b
Z b X
n
f (x)j (x)dx, j = 1, . . . , n,
L
i i (x)j (x)dx =
a
i=1
es decir, se trata de un sistema de n ecuaciones lineales con n incognitas
1 , . . . , n que, en forma matricial, se puede escribir como
A = b,
12
L(i )(x)j (x)dx
siendo A =
a
, = (1 , . . . , n ) y b =
1i,jn
f (x)i (x)dx
a
A la solucion discreta obtenida yn se le llama aproximacion nesima, por
el metodo de Galerkin, del problema de contorno dado, a partir del sistema
de funciones coordenadas considerado.
La matriz A no tiene por que ser, en principio simetrica.
5.5
Transformaci
on a forma autoadjunta
Si un operador diferencial lineal de segundo orden es autoadjunto entonces
se puede escribir de la forma
Ly(x) = (p(x)y 0 (x))0 + q(x)y(x).
En este caso, como veremos en la siguiente seccion, el metodo de Galerkin
va a producir un sistema de ecuaciones lineales cuya matriz de coeficientes
es simetrica lo cual facilita los calculos.
Por ello, es interesante trabajar con operadores que sean autoadjuntos. A
continuacion describimos un metodo para transformar el problema diferencial
lineal de segundo orden
a0 (x)y 00 (x) + a1 (x)y 0 (x) + a2 (x)y(x) = f (x), x [a, b],
en uno problema autoadjunto.
Se puede comprobar que, llamando
Z x
a1 (x)
a2 (x)
f (x)
p(x) = exp
dx , q(x) = p(x)
, r(x) = p(x)
,
a0 (x)
a0 (x)
a a0 (x)
entonces el problema anterior se puede escribir en la forma
(p(x)y 0 (x))0 + q(x)y(x) = r(x), x [a, b].
5.6
M
etodo de Ritz
Consideremos el problema de contorno lineal de segundo orden, dado en
forma autoadjunta y con condiciones homogenas (por ejemplo de Dirichlet,
analogamente se razona para los otros tipos):
13
.
1in
(p(x)y 0 (x))0 + q(x)y(x) = r(x), x [a, b],
y(a) = y(b) = 0.
Sea {n }n1 un sistema de funciones coordenadas verificando las condiciones de contorno que aparecen en el problema dado.
Aplicando el metodo de Galerkin, la nesima aproximacion de la solucion
n
X
de este problema es la funcion yn (x) =
i i (x), siendo = (1 , . . . , n =t
i=1
la solucion del sistema de ecuaciones lineales A = b, con A = (aij )1i,jn y
b = (bi )1in , siendo, para cada i, j = 1, . . . , n,
aij =
Rb
((p(x)0i (x))0 + q(x)i (x)) j (x)dx =
Rb
Rb
a (p(x)0i (x))0 j (x)dx + a q(x)i (x)j (x)dx
a
y, aplicando partes a la primera integral, obtenemos
Z b
aij =
p(x)0i (x)0j (x) + q(x)i (x)j (x) dx.
a
El metodo as obtenido se llama metodo de Ritz. En conclusion, la n
esima aproximacion de la solucion de problema de contorno dado, mediante
n
X
el metodo de Galerkin es yn (x) =
i i (x), siendo = (1 , . . . , n ) la
i=1
solucion del sistema
A = b
A=
a
p(x)0i (x)0j (x)
Z
+ q(x)i (x)j (x) dx
r(x)i (x)dx
a
,
1i,jn
b=
.
1in
Observese que, frente al metodo de Galerkin, en este caso la matriz de
coeficientes es simetrica.
14