DIFERENCIAS FINITAS Y ECUACIONES DIFERENCIALES

PARCIALES
L. H
ector Ju
arez V.
Departamento de Matem
aticas, UAM-I, M
exico, D. F.
Guanajuato, Gto, Junio 2006

Introducci
on

Este minicurso constituye una introduccion a la solucion numerica de ecuaciones diferenciales

parciales del tipo parabolico. Se discuten los aspectos basicos fundamentales como son:
derivacion de los esquemas de diferencias finitas, estudio del error, de la convergencia y de la
estabilidad de los metodos. Ademas se incluyen programas en el ambiente MATLAB para
ilustrar algunos resultados numericos con los diferentes esquemas.

Ejemplo Cl
asico Sencillo

Comezamos con el estudio del modelo mas sencillo de una ecuacion diferencial parcial del
tipo parabolico, a saber la ecuacion que modela el flujo o propagacion de calor en una barra
homogenea de longitud finita sin fuentes de calor. La ecuacion es la siguiente:
k 2u
u
=
, en 0 < x < L, t > 0
t
Cp x2
u(0, t) = u(L, t) = 0, t > 0,
u(x, 0) = u0 (x) , 0 x L ,

(2.1)
(2.2)
(2.3)

donde u = u(x, t) es la temperatura de la barra en la posicion x al tiempo t, k es el coeficiente

de conductividad termica, es la densidad de la barra, Cp es el calor especfico, y L es la
longitud de la barra. Suponemos que los parametros k, y Cp son constantes.
Adimensionalizaci
on. Con el objeto de simplificar la ecuacion diferencial anterior se inx 0
u
0
troducen las siguientes variables adimensionales: x = , u = , donde U = max u0 (x)
L
U
u
u0
con 0 x L. Con lo anterior se obtienen las siguientes relaciones
= U
y
t
t
0
2u
U 2u
=
. Sustituyendo estas u
ltimas expresiones en la ecuacion (2.1) se obtiene
x2
L2 x0 2

Diferencias Finitas y E.D.P.

2
0

L2 Cp u0
2u
0
la ecuacion
=
. Si se toma t = t k/L2 Cp , entonces las ecuaciones (2.1)
02
k
t
x
(2.3) se tranforman, despues de quitar la comilla de las variables x y u, en la ecuacion mas
sencilla:
ut (x, t) = uxx (x, t) en 0 < x < 1, t > 0

(2.4)

u(0, t) = u(1, t) = 0, t > 0,

(2.5)

u(x, 0) = u0 (x), 0 x 1,

(2.6)

en donde ut y uxx representan derivadas de u respecto a t y x, respectivamente. En las

u(x,t)

1
0
Figura 1. Distribucion de temperatura en una barra uniforme

secciones siguientes consideraremos la ecuacion adimensional (2.4)(2.6), a menos que se

especifique lo contrario. La Figura 1 ilustra el problema modelo.

M
etodo de separaci
on de variables

La ecuacion anterior tiene una solucion en forma de series de Fourier que se puede calcular
por medio del metodo de separacion de variables. El metodo de separcion de variables es un
metodo con limitaciones ya que no es aplicable a casos un poco mas generales. Sin embargo,
proporciona soluciones u
tiles para propositos de comparacion, ademas de que proporciona
un medio para el analisis de la estabilidad de los metodos de diferencias finitas. En el metodo
de separcion de variables se buscan soluciones de la forma u(x, t) = X(x) T (t). Sustituyendo
00
en la ecuacion (2.4), obtenemos X T = X T , en donde el punto indica derivacion respecto al
00
X
T
=
tiempo y las comillas derivacion respecto a la posicion. De lo anterior se tiene que
T
X
debe ser igual a una constante k 2 , debido a que el lado izquierdo depende solo de t y el lado
derecho depende solo de x. Nombramos k 2 a esta constante para enfatizar que debe ser
2

negativa. Luego entonces, se tiene X(x) = sen kx, y T (t) = ek t . Por lo tanto, la solucion
debe de ser de la forma
2

u(x, t) = ek t sen kx.

(3.1)

Diferencias Finitas y E.D.P.

Ahora es claro que el valor negativo k 2 se toma porque no es posible obtener soluciones
con crecimiento exponencial en el problema fsico. En ausencia de fuentes de calor, y debido
a las condiciones de frontera u = 0 en los extremos de la barra, el fenomeno de difusion
presente en el problema provoca que la temperatura disminuya a partir del valor inicial.
Condiciones de frontera. Las condiciones de frontera sirven para deteminar los valores
que puede tomar la constante k. La primera condicion de frontera u(0, t) = 0 es satisfecha
automaticamente por (3.1). La segunda condicion de frontera
2

u(1, t) = ek t sen k = 0,
es satisfecha por (3.1) si k = m para m = 1, 2 . . ..
Modos de Fourier. De acuerdo a lo anterior, el problema (2.4)(2.5) tiene una infinidad
de soluciones denominadas modos de Fourier, y estas son de la forma
2 2 t

um (x, t) = em

sen mx .

(3.2)

Por lo tanto, la solucion general de la ecuacion diferencial (2.4), y que satisface las condiciones
de frontera (2.5), es una superposicion de estos modos de Fourier. Es decir
u(x, t) =

2 2 t

am em

sen mx

(3.3)

m=1

en donde las constantes am se encuentran de las condiciones iniciales (2.6) y de las condiciones
de ortogonalidad de las funciones propias sen mx.
Ortogonalidad. Las funciones sen mx, m = 1, 2, . . . son ortogonales bajo el producto
R1
R1
escalar de funciones hf, gi = 0 f g dx. Ademas 0 sen2 mx dx = 1/2. La condicion inicial
(2.6) debe ser satisfecha por (3.3):
u(x, 0) =

am sen mx = u0 (x).

m=1

Multiplicando por sen nx esta u

ltima ecuacion e integrando se obtiene
Z 1
an = 2
u0 (x) sen(nx) dx, n = 1, 2, . . .

(3.4)

Por lo tanto la solucion del problema (2.4)(2.6) en terminos de series de Fourier es

Z 1

X
(m)2 t
u(x, t) =
am e
sen(mx),
am = 2
u0 (x) sen(mx) dx.
m=1

(3.5)

Esta solucion puede considerarse como una solucion analtica del problema, pero tiene el
inconveniente que debemos calcular una infinidad de integrales para encontrar los coeficientes

Diferencias Finitas y E.D.P.

am y despues sumar un n
umero infinito de terminos, lo cual no es posible en la practica. Para
encontrar una solucion concreta se deben evaluar solo un n
umero finito de terminos de la
serie, lo cual la convierte en una solucion aproximada desde el punto de vista practico. Una
de las principales limitacion del metodo es que no es generalizable a ecuaciones diferenciales
parciales ligeramente mas complicadas.

Esquema explcito de diferencias finitas

Los metodos de aproximacion constituyen una alternativa cada vez mas viable para encontrar
soluciones de ecuaciones diferenciales. De entre los metodos mas populares se encuentran
los metodos de diferencias finitas debido a su simplicidad. El esquema mas sencillo para
el problema parabolico unidimensional es el esquema explcito de diferencias finitas. A
continuacion presentamos como se encuentra el esquema:
1. Construcci
on de la malla Se divide el intervalo (0, 1) en J subintervalos y el intervalo
de tiempo (0, tf ), donde tf es el tiempo final, en N subintervalos
x = 1/J
xj = j x, j = 0, 1, ..., J.
tf
t =
tn = n t, n = 0, 1, ..., N.
N
El conjunto de puntos (xj , tn ) as construidos constituyen una malla o latice en donde

t n+1
tn

1111
0000

1
0
1
0

11
00

11111111
00000000
1
0
1
0
1
0
0
1
1
0
1
0

11
00

xj1 x j x j+1

Figura 2. Malla y molecula computacional para el esquema explcito

queremos encontrar los valores aproximados Ujn de la solucion exacta u(xj , tn )
2. Aproximaci
on de las derivadas por medio de diferencias finitas. Para encontrar
una aproximacion de las derivadas que aparecen en la ecuacion diferencial, en el punto

Diferencias Finitas y E.D.P.

(xj , tn ), utilizamos las siguientes expansiones en serie de Taylor (suponiendo que la solucion
es suficientemente suave):
(t)2
(t)3
utt (x, t) +
uttt (x, t) + ...
2!
3!
(x)2
(x)3
u(x + x, t) = u(x, t) + x ux (x, t) +
uxx (x, t) +
uxxx (x, t) + ...
2!
3!
(x)2
(x)3
u(x x, t) = u(x, t) x ux (x, t) +
uxx (x, t)
uxxx (x, t) + ...
2!
3!
u(x, t + t) = u(x, t) + t ut (x, t) +

(4.1)
(4.2)
(4.3)

Aproximaci
on de ut en el punto (xj , tn ). Se obtiene evaluando (4.1) en dicho punto:
(x)2
utt (xj , ), (tn , tn+1 ) ,
2!
y posteriormente despejando ut (xj , tn ), para obtener
u(xj , tn+1 ) = u(xj , tn ) + t ut (xj , tn ) +

ut (xj , tn ) =

u(xj , tn+1 ) u(xj , tn ) t

utt (xj , ).
t
2

(4.4)

Ujn+1 Ujn
con error O(t).
t

(4.5)

Por lo tanto, podemos escribir:

ut (xj , tn )

Aproximaci
on de uxx en el punto (xj , tn ). Se obtiene evaluando (4.2)(4.3) en dicho
punto y sumando para obtener:
u(xj+1 , tn ) + u(xj1 , tn ) = 2 ut (xj , tn ) + (x)2 uxx (xj , tn ) +

(x)4
uxxxx (, tn )
12

con (xj1 , xj+1 ). Despejando uxx (xj , tn ), se obtiene:

uxx (xj , tn )

u(xj+1 , tn ) 2 ut (xj , tn ) + u(xj1 , tn ) (x)2

uxxxx (, tn ).
(x)2
12

(4.6)

Por lo tanto, podemos escribir:

uxx (xj , tn )

n
n
Uj+1
2 Ujn + Uj1
con error O((x)2 ).
(x)2

(4.7)

3. Aproximaci
on de la ecuaci
on diferencial
De acuerdo a (4.5) y (4.7) el problema (2.4)(2.6) se puede sustituir por el problema
aproximado
n
n
Ujn+1 Ujn
2Ujn + Uj1
Uj+1
=
,
t
(x)2
U0n = UJn = 0, 1 n N,

Uj0 = u0 (xj ), 0 j J.

(Ecuacion en diferencias)
(Condiciones de frontera)
(Condiciones iniciales)

(4.8)
(4.9)
(4.10)

Diferencias Finitas y E.D.P.

La ecuacion en diferencias (4.8) puede escribirse como

n
n
Ujn+1 = Ujn + (Uj+1
2Ujn + Uj1
), donde =

t
,
(x)2

(4.11)

para j = 1, ..., J 1, y n = 0, ..., N 1. Claramente cada valor en el nivel tn+1 puede

calcularse independientemente de los demas en el mismo nivel, y solo evaluando valores en el
nivel tn . De ah el nombre de esquema explcito. Las condiciones de frontera (4.9) e iniciales
(4.10) son valores conocidos, y por tanto es posible calcular la solucion aproximada en todos
los puntos interiores para valores sucesivos de n a partir de la formula (4.11) (ver Figura 2).
Ejemplo 1 Para ilustrar el comportamineto del esquema se escoge la distribucion inicial de
temperaturas dada por
(
2x
si 0 x 1/2,
u0 (x) =
2 2x
si 1/2 x 1 .
Se utiliza el siguiente parametro de discretizacion: x = 0.05, el cual corresponde a J = 20.
Ademas se consideran dos casos t = 0.0012 y t = 0.0013, es decir = 0.48 y = 0.52.
Para los calculos se utilizo el programa expllicito1.m escrito para el ambiente MATLAB y
que se incluye al final de este trabajo. Los resulltados se muestran en la Figura 3.
Despues de realizar los experimentos numericos, se observa que con t = 0.0012 se obtienen resultados numericos aceptables. Por otro lado, con t = 0.0013 la solucion numerica
muestra oscilaciones que crecen rapidamente conforme t aumenta. Esto muestra que no
es posible escoger x y t de manera arbitraria. Es sorpredente la diferencia entre las
soluciones a pesar de la diminuta diferencia entre los dos valores de t. Este fenomeno
es un resultado tpico de estabilidad o inestabilidad que aparece frecuentemente cuando
se utilizan esquemas explcitos para resolver numericamente ecuaciones diferenciales. Mas
adelante estudiaremos detenidamente el fenomeno de estabilidad y veremos su dependencia
respecto de la razon de mallas de los parametors de discretizacion por medio de .

Error de truncamiento

El error de truncamiento es un concepto u

til para estudiar el error de los metodos de aproximacion as como la convergencia de los mismos.
El error de truncamiento T (x, t) para el esquema (4.8) se obtiene cuando la aproximacion

Diferencias Finitas y E.D.P.

7
n=0

0.9

0.9

0.8

0.8

0.7

0.7

0.6

0.6

0.5

0.5

0.4

0.4

0.3

0.3

0.2

0.2

0.1
0

0.1
0

0.2

0.4

0.6

0.8

n = 10

0.9

0.8

0.8

0.7

0.7

0.6

0.6

0.5

0.5

0.4

0.4

0.3

0.3

0.2

0.2

0.4

0.6

0.8

0.2

0.4

0.6

0.8

0.2

0.4

0.6

0.8

0.2

0.4

0.6

0.8

0.1
0

0.2

0.4

0.6

0.8

n = 30

1
0.9

1
0.9

0.8

0.8

0.7

0.7

0.6

0.6

0.5

0.5

0.4

0.4

0.3

0.3

0.2

0.2

0.1
0

0.2

0.9

0.1

0.1
0

0.2

0.4

0.6

0.8

n = 50

1
0.9

1
0.9

0.8

0.8

0.7

0.7

0.6

0.6

0.5

0.5

0.4

0.4

0.3

0.3

0.2

0.2

0.1
0

0.1
0

0.2

0.4

0.6

t = 0.0012

0.8

t = 0.0013

Figura 3. Solucion del problema por medio del metodo explcito.

Ujn se reemplaza por la solucion exacta de la ecuacion diferencial:

T (x, t) :=

u(x, t + t) u(x, t)
u(x + x, t) 2u(x, t) + u(x x, t)

t
(x)2

(5.1)

Sustituyendo u(x, t+t), u(x+x, t), y u(xx, t) por sus expansiones en series de Taylor
(4.1), (4.2), y (4.3), respectivamente, y simplificando posteriormente, se obtiene
1
1
utt (x, t) t
uxxxx (x, t) (x)2 + . . .
2
12
1
1
uxxxx (x, t) (x)2 + T OA,
= utt (x, t) t
2
12

T (x, t) = ut uxx +

(5.2)
(5.3)

Diferencias Finitas y E.D.P.

donde T OA indica terminos de orden alto. Frecuentemente es conveniente truncar la serie

infinita introduciendo un termino de residuo:
1
1
T (x, t) = utt (x, ) t
uxxxx (, t) (x)2
2
12
donde (x x, x + x), (t, t + t).

(5.4)

Consistencia del esquema. Para el esquema explcito vemos que

T (x, t) 0

cuando t, x 0,

(x, t) ,

independientemente de la relacion entre t y x. Gracias a esta propiedad se dice que el

esquema es consistente con la ecuacion diferencial.
Orden del esquema. Tambien podemos observar que
1
1
t
Mxxxx
Mtt t +
Mxxxx (x)2 =
(Mtt +
),
2
12
2
6
donde Mtt es una cota superior para |utt |, y Mxxxx es una cota superior para |uxxxx |. Por
|T (x, t)|

lo tanto, para valores fijos de observamos que |T (x, t)| se comporta asintoticamente como
O(t). Con la excepcion de ciertos valores especiales de , esta sera la mayor potencia de t
para la cual se satisface esta propiedad. Por lo tanto, se dice que el esquema tiene precision
de primer orden en t, y se escribe T (x, t) O(t).
Subiendo el orden del esquema. Si u es suficientemente suave, entonces ut = uxx implica
utt = uxxxx . En este caso la relacion (5.4) se puede escribir como
1
1
T (x, t) = (1 ) uxxxx t + O(t2 )
2
6
Por lo tanto, si tomamos = 1/6 se obtiene que el esquema es de segundo orden en t. Es
decir, en este caso se tiene
T (x, t) v O(t2 )
t
1
= es muy severa, pues el paso del tiempo debe tomarse
2
(x)
6
muy peque
no. Por ejemplo
Sin embargo, la condicion

x = 0.1 = t = 1.67 103 ,

x = 0.01 = t = 1.67 105 .

Convergencia del m
etodo explcito

Sin entrar en detalles tecnicos, podemos decir que el esquema explcito de diferencias finitas
es convergente si para alg
un valor fijo de la solucion aproximada Ujn tiende al valor

Diferencias Finitas y E.D.P.

exacto de la solucion cuando ambos t y x tienden a cero. Probaremos que el esquema es

convergente si 1/2.
Definimos
en+1
:= Ujn+1 u(xj , tn+1 ) .
j

(6.1)

Evaluando (5.1) en (xj , tn ) y despejando u(xj , tn+1 ), obtenemos

u(xj , tn+1 ) = u(xj , tn ) + {u(xj+1 , tn )2u(xj , tn )+u(xj1 , tn )}+t T (xj , tn ).

(6.2)

Sustituyendo (6.2) y (4.11) en (6.1), y simplificando, se obtiene

en+1
= (1 2) enj + enj+1 + enj1 t Tjn .
j

(6.3)

1
Si , entonces los primeros tres coeficientes del lado derecho de esta igualdad son no
2
negativos y suman 1. Por lo tanto
|en+1
| = (1 2) |enj | + |enj+1 | + |enj1 | + t |Tjn | .
j

(6.4)

E n := max{ |enj |, j = 0, 1, ..., J } ,

(6.5)

Definiendo
y
T =

t
Mxxxx
(Mtt +
)
2
6

cota superior de |T (x, t)| ,

se obtiene
E n+1 E n + t T .
Si ademas suponemos E 0 = 0 (valores iniciales exactos), entonces
E n N t T tf T 0

si t 0.

(6.6)

Para lograr lo anterior se ha supuesto que Mtt y Mxxxx se satisfacen uniformemente sobre
toda la region. Esto solo es posible si suponemos que los datos iniciales y de frontera son
consistentes, y si los datos inciales son suficientemente suaves. De otra manera debemos se
debe depender del efecto suavizador del operador de difusion para asegurar cotas de esta
forma para todo tiempo t > 0
Ejemplo2. En el siguiente ejemplo donde u(0, t) = u(1, t) = 0 y u(x, 0) = 0.9, es un caso
donde las condiciones inciales no son consistentes con las condiciones de frontera, pues es las
esquinas x = 0 y x = 1 hay una discontinuidad. La Figura 4 muestra como el operador de
disusion casi inmediatamente suviza la solucion

Diferencias Finitas y E.D.P.

10

0.9

0.9

0.8

0.8

0.7

0.7

0.6

0.6

n=0

0.5
0.4

0.4

0.3

0.3

0.2

0.2

0.1
0

0.1
0

0.2

0.4

0.6

0.8

0.9

0.9

0.8

0.8

0.7

0.7

0.6

0.6

0.5

0.2

0.4

0.6

0.8

0.8

0.5

n=5

0.4

n = 15

0.4

0.3

0.3

0.2

0.2

0.1
0

n=1

0.5

0.1
0

0.2

0.4

0.6

0.8

0.2

0.4

0.6

Figura 4. Efecto suavizador del operador de difusion

La propiedad de convergencia de la solucion aproximada a la solucion exacta cuando los

parametros de discretizacion t y x tienden a cero es fundamental para cualquier esquema
de aproximacion. Muestra que puede obtenerse la precision deseada con el uso de una malla
lo sificientemente fina. Sin embargo, no es conveniente hacer un refinamiento excesivo pues,
aparte que el volumen de calculos aumenta dramaticamente, tambien se acumula el error
de redondeo. Es bien sabido que a partir de cierto refinamiento el error de redondeo puede
dominar al error de truncamiento.

An
alisis de Fourier del error

Hemos expresado la solucion exacta de la ecuacion diferencial como una serie de Fourier,
la cual es una suma infinita de modos como la expresada en la ecuacion (3.5). En forma
compleja esta serie se puede escribir como
u(x, t) =

Z
am e

(m)2 t imx

am = 2

u0 (x)eimx dx.

En este caso los modos de Fourier son um (x, t) = e(m) t eimx = ek t eikx .

(7.1)

Diferencias Finitas y E.D.P.

11

N
umero de onda. Al n
umero k := m se le denomina el n
umero de onda del modo. Si el
n
umero de onda k es grande se dice que el modo es de alta frecuencia. En un punto de la
malla xj = jt, tn = nt, se tiene
um (xj , tn ) = (ek

2 t

)n eik(jx) .

(7.2)

La amplitud de cada modo es el primer termino del lado derecho, y este tiende a cero con el
tiempo:
(ek

2 t

)n 0 cuando n .

Veremos que un modo similar de Fourier es la solucion exacta de la ecuacion en diferencias.

n
n
) tiene un
2Ujn + Uj1
Caso Discreto. La ecuacion en diferencias Ujn+1 = Ujn + (Uj+1

modo de fourier discreto de la forma

Ujn (k) := n eik(jx) .

(7.3)

Factor de amplificaci
on del modo: Al termino = (k) se le demonima el factor de
amplificacion del modo de Fourier discreto con n
umero de onda k, y sustituye al termino
exponencial ek

2 t

Sustituyendo el modo de Fourier discreto (7.3) en la ecuacion en diferencias (4.11), se

obtiene
n+1 eik(jx) = n eik(jx) + (n eik(j+1)x 2n eik(jx) + n eik(j1)x ) .
Simplificando esta expresion, se obtiene
(k) = 1 + ( eikx 2 + eikx )
= 1 4 sen2 (kx/2),

t
.
(x)2

Por lo tanto, la solucion de la ecuacion en diferencias o solucion aproximada del problema

se puede escribir como la serie (al poner k = m)
X
Am (m)n eim(jx)
Ujn =

(7.4)

mZ

con (m) = 1 4 sen2 (mx/2).

Comparaci
on entre la soluci
on exacta y la soluci
on discreta. Para efectuar una
comparacion entre las soluciones exacta y discreta solo comparamos el factor de ampliicacion
2

(k) con el termino exponencial ek t , debido a que esta es la u

nica diferencia entre las
expresiones para los modos de Fourier dicretos y los modos de Fourier continuos. Haciendo

Diferencias Finitas y E.D.P.

12

una expansion en series de Taylor de ambos, encontramos que

1 4
k (t)2 . . .
2
1 4
(k) = 1 k 2 t +
k t (x)2 . . .
12
y observamos que ambas son iguales hasta el segundo termino, que es un termino de primer
ek

2 t

= 1 k 2 t +

orden en t. Por lo tanto concluimos que la solucion aproximada es de primer orden

t
t. Ademas observamos que la aproximacion es de segundo orden si =
=
(x)2
Observese que hemos recobrado las mismas conclusiones que con el analisis del error

en
1
.
6
de

truncamiento y convergencia. Sin embargo, con el analisis de Fourier podemos decir un poco
mas. Por ejemplo, los terminos de baja frecuencia (k peque
no) son bi
en aproximados
pues

1
1
) + . . . C() k 4 (t)2 .
12 2
Esto no sucede con los terminos de alta frecuencia como veremos en seguida. Por un lado,
2
para valores grandes de k en la solucion exacta ek n t 0 rapidamente cuando n .
1
Pero en la solucion numerica |(k)|n cuando n si > , debido a que en este
2
caso |(k)| > 1. En particular para m = J, k x = , y
|(k) ek

2 t

| = k 4 (t)2 (

1
|(k)| = |1 4 sen2 | = 4 1 > 1 si > .
2
2

Estabilidad

Para el modelo que hemos considerado decimos que el metodo es estable si existe K > 0,
independiente de k, tal que
|(k)n | K,

n,

k.

(8.1)

Condici
on de Von Neumann. En la expresion anterior, si K > 1 entonces la funcion K s
es convexa en el intervalo [0, 1], de tal forma que queda por debajo del segmento de recta
dado por 1 + (K 1)s como se muestra en la Figura 5.
1
Lugo entonces para s = con n = 1, 2 . . . se tiene
n
1
t
|(k)| K 1/n 1 + (K 1) = 1 + (K 1)
.
n
tf
0

Si llamamos K a la constante

k1
obtenemos la celebre condici
on de Von Neumann
tf

Diferencias Finitas y E.D.P.

13

K
1+(K1)s
s

0<s<1

Figura 5. Graficas de K y 1 + (K 1)s, con 0 < s < 1

para estabilidad:
0

|(k)| 1 + K t

k.

(8.2)

Esta condicion es necesaria y suficiente para la convergencia de un esquema de diferencias

consistente que aproxime a una ecuacion diferencial. Como la anterior condicion debe satisfacerse para t > 0 arbitrario, entonces para que el esquema sea estable basta con que el
factor de amplificacion de todo modo de Fourier discreto sea menor que uno, es decir
|(k)| 1,

k.

(8.3)

En particular, para nuestro modelo, el metodo explcito de diferencias es estable si

1
|(k)| = |1 4 sen2 kx| 1
2
21
1 1 4 sen kx 1
2
21
0 4 sen kx 2.
2
Por lo tanto, el m
etodo explicito es
Estable si 1/2 ,
Inestable si > 1/2 .
Numero finito de modos discretos de Fourier. Hemos usado una serie infinita para
representar Ujn en terminos de los modos de Fourier como se muestra en (7.4). Sin embargo,
sobre una malla discreta solamente puede haber un n
umero finito de modos distintos. De
hecho, dos modos con n
umero de onda k1 y k2 son iguales si (k1 k2 )x es un multiplo de
2. Por lo tanto, los modos distintos asociados a la malla discreta son aquellos asociados a

Diferencias Finitas y E.D.P.

14

los numero de onda k = m con m = (J 1), (J 2), . . . , 1, 0, 1, . . . , J. El modo mas

alto sobre la malla es el asociado a k = J. De acuerdo a la ecuacion (7.3) este modo es:
Ujn (J) = (J)n eiJ(jx) = (1 4)n cosj = (1 4)n (1)j
De hecho, este es el modo m
as inestable y es el responsable de las oscilaciones cuando
> 1/2, como se observa en el Ejemplo 1 de la Seccion 4.

M
etodo implcito

La restriccion de estabilidad para el metodo explcito 1/2 es muy severa, especialmente

cuando queremos resultados en tiempos grandes. Ademas, si se quisiera mejorar la precision
de los resultados disminuyendo x = 1/J, la cantidad de trabajo aumentara proporcionalmente a J 2 , lo cual representa un gran costo computacional. Uno de los remedios mas
utilizados para remontar estos inconvenientes es utilizar un esquema de tipo implcito, es
decir, utilizar equemas de diferencias finitas de tiempo haca atras. En nuestro modelo para
el flujo de calor en la barra, esto se logra evalundo el lado derecho de la ecuacion (4.8) en
n + 1 en lugar de n. El esquema de diferencias implcito es, por lo tanto
n+1
n+1
Ujn+1 Ujn
Uj+1
2Ujn+1 + Uj1
=
t
(x)2

La molecula computacional se muestra en la siguiente Figura 6.

t n+1
tn

1111
0000

11111111
00000000
1
0
0
1
0
1
0
1
0
1
0
1

11
00

1
0
0
1

11
00

xj1 x j x j+1

Figura 6. Molecula computacional del esquema implcito

(9.1)

Diferencias Finitas y E.D.P.

15

La ecuacion en diferencias puede escribirse en la forma siguiente

n+1
n+1
Uj1
+ (1 + 2) Ujn+1 Uj+1
= Ujn ,

j = 1, . . . , J 1 ,

(9.2)

complementada por las condiciones de frontera (4.9) y las condiciones iniciales (4.10).
Estabilidad. El esquema es incondicionelamente estable pues al sustituir U (k)nj =
n ei k j x en la ecuacion en diferencias (9.2), simplificando y despejando , en forma analoga
a como se hizo en el caso explcito, se obtiene el factor de amplificacion
(k) =

1
,
1 + 4 sen2 ( 12 k x)

(9.3)

el cual es menor que 1 en valor absoluto para toda k y para todo valor de . Por lo tanto
t y x se escogen tomando solo criterios de precision.
Costo adicional. Sin embargo, el esquema es un poco mas complicado debido a que para
obtener la solucion en cada nivel tiempo tn+1 es necesario resolver un sistema de J 1
ecuaciones con J 1 incognitas. Es decir, el costo adicional, es resolver en cada nivel de
tiempo un sistema de ecuaciones de tama
no (J 1) (J 1) con matriz constante

2 1
1 2 1

A = I + . . . . . . . . . . . . . . . . . . . . . .

1 2 1
1 2
donde I indica la matiz identidad de tama
no (J 1) (J 1).
Ejemplo 3. Para ilustrar el comportamiento del esquema implcito, consideremos el problema resuelto en el Ejemplo 1 (Seccion 4). Se utilizan el mismo parametro de discretizacion
x = 0.05 (J = 20). Se consideran los casos t = 0.0013 ( = 0.52) y t = 0.0125 ( = 5).
Para los calculos se utilizo el programa implicito1.m escrito para el ambiente MATLAB y
que se incluye al final de este trabajo. Los resultados se muestran en la Figura 7.
El lector puede hacer experimentos con algun valor fijo de , y varios valores de J, por
ejempo J = 200, 400, 1000, etc., para verificar que la memoria se agota rapidamente, y en
consecuencia la eficiencia del metodo disminuye dramaticamente.
Algoritmo de Thomas. Realmente no hemos hecho gran cosa, excepto que hemos utilizado

Diferencias Finitas y E.D.P.

16

0.9

0.9

0.8

0.8

0.7

0.7

0.6

0.6

0.5

0.5

0.4

0.4

0.3

0.3

n=0

0.2

n=0

0.2

0.1

0.1
0

0.2

0.4

0.6

0.8

0.9

0.9

0.8

0.8

0.7

0.7

0.6

0.6

0.5

0.5

0.4

0.4

0.3

0.3

n = 10

0.2

0.2

0.4

0.6

0.8

1
0.9

0.8

0.8

0.7

0.7

0.6

0.6

0.5

0.5

0.4

0.4

0.3

0.8

0.2

0.4

0.2

0.4

0.2

0.4

0.6

0.8

0.3

0.2

0.2

n = 30

0.1
0

0.2

0.4

0.6

n=3

0.1
0.8

0.9

0.9

0.8

0.8

0.7

0.7

0.6

0.6

0.5

0.5

0.4

0.4

0.3

0.3

0.2

0.2

n = 50

0.1
0

0.6

0.1

0.9

0.4

n=1

0.2

0.1
0

0.2

0.2

0.4

0.6

t = 0.0013 ( = 0.52 )

0.8

0.8

n=5

0.1
0.8

0.6

0.6

t = 0.0125 ( = 5)

Figura 6. Resultados numericos del problema modelo con el metodo implcito

un esquema que es estable para todo valor de los parametros de discretizacion. Sin embargo,
el costo computacional es demasido alto a
un, debido a la utilizacion de memoria (por la
necesidad de almacenar la matriz) como por el volumen de operaciones que hay que realizar,
al tener que resolver un sistema de ecuaciones en cada iteracion de tiempo.
Afortunadamente la matriz A del sistema es una matriz rala (sparse), es decir, la
inmensa mayora de su coeficientes es cero. De hecho, la matriz A tiene una estructura muy
regular, debido a que es una matiz tridiagonal y ademas diagonalmente dominante.
Esta estructura se puede explotar inteligentemente para ahorrar memoria y resolver el sistema
de manera muy eficiente. El metodo mas eficiente para el presente caso es el denominado
algoritmo de Thomas. Este metodo consiste en reducir el sistema tridiagonal a un sistema

Diferencias Finitas y E.D.P.

17

bidiagonal superior de la forma

Uj ej Uj+1 = fj

j = 1, 2, . . . , J 1,

por medio de eliminacion del termino Uj1 en cada ecuacion de la forma

Uj1 + (1 + 2) Uj Uj+1 = dnj ,

j = 1, 2, . . . , J 1 .

en donde identificamos dnj con Ujn , el lado derecho de (9.2) (obsevese que hemos suprimido
los superndice n y n + 1 en los terminos con U para simplificar la notacion).
Para lograr la reduccon anterior se procede como sigue: en la primera ecuacion del
sistema tridiagonal j = 1 y Uj1 = U0 = 0, por los que esta ecuacion se transforma en

dn1
, f1n =
.
1 + 2
1 + 2
Despejando U1 de esta ecuacion, sustituyendo en la segunda ecuacion del sistema tridiagonal,
U1 e1 U2 = f1n ,

con e1 =

y reordenando, se obtiene
U2 e2 U3 = f2n ,

con e2 =

,
1 + 2 e1

f2n =

d2 + f1n
.
1 + 2 e1

Observese que e1 y e2 no dependen de n, pero f1n y f2n s. Continuando inductivamente de

esta forma, obtenemos el sistema bidiagonal
Uj ej Uj+1 = fjn ,

con ej =

,
1 + 2 ej1

fjn =

n
dj + fj1
,
1 + 2 ej1

(9.4)

para j = 1, . . . , J 1, con e0 = f0 = 0. Este sistema bidiagonal resultante se puede resolver

entonces de manera muy facil por medio de sustituci
on hacia atr
as:
UJ = 0
Para j = J 1, . . . , 1, hacer
Uj = fjn + ej Uj+1
finalizar
La anterior metodologa es parte de los algoritmos que se conocen come como resolvedores
rapidos o fast solvers por su nombre en el idioma Ingles. La tecnica de solucion es un caso
particular de una tecnica denominada de reduccon cclica. Debido a que la matriz original
es diagonalmente dominante cada coeficiente ej en (9.4) es menor a uno en valor absoluto,
lo cual garaniza que el uso de las relaciones de recurrencia como Uj = fj + ej Uj+1 = fj ,
para encontrar la solucion, es n
umericamente estable. El trabajo computacional del metodo
implcito es apenas tres veces aproximadamente el del metodo explcito cuando se utiliza

Diferencias Finitas y E.D.P.

18

el algoritmo de Thomas. La verdadera importancia del metodo implicito radica en que los
pasos de tiempo pueden tomarse significativamente mayores a los permitidos por el metodo
explcito, con lo cual al final hay un ahorro real en el costo computacional.
Ejemplo 4. Consideremos de nuevo el problema resuelto en el Ejemplo 1 (Seccion 4). Se
pueden utilizar los parametros de discretizacion: x = 0.05 (J = 20) y = 1. Los resultados
son muy semejantes a los ya obtenidos anteriormente. El progama implicito1t.m, que
contiene el algorimto de Thomas para resolver los sistemas tridiagonales, se incluye al final
de las notas.
El lector puede repetir el experimeto del Ejemplo 3 con varios valores de J, y verificar
que con el algoritmo de Thomas hay un ahorro real en memoria y un aumento considerable
en la eficiencia para encontrar la solucion numerica.

10

Metodo (Promedios pesados)

Los metodos explcito e implcito difieren solamente en que el primero aproxima uxx por tres
punto en el nivel de tiempo tn , y el segundo por tres puntos en el nivel de tiempo tn+1 . Una
generalizacion natural de estos metodos sera una aproximacion que use los seis puntos. Esto
puede lograrse mediante un promedio pesado de las dos ecuaciones en diferencias:
n+1
n+1
n
n
Ujn+1 Ujn
Uj+1
2Ujn+1 + Uj1
Uj+1
2Ujn + Uj1
=
+
(1

)
,
t
(x)2
(x)2

(10.1)

con 0 1.
Con = 0

Obtenemos el metodo explcito.

Con = 1

Obtenemos el metodo implcito.

En la siguiente figura ilustramos la molecula computacional del metodo

La ecuacion en diferencias constituye un sistema tridiagonal de J 1 ecuaciones con J 1
incognitas en cada paso del tiempo n = 1, 2, . . .
n+1
n+1
n
n
),
2Ujn + Uj+1
Uj1
+ (1 + 2)Ujn+1 Uj+1
= Ujn + (1 )(Uj1

(10.2)

con j = 1, . . . , J 1, y complementada por las condiciones de frontera (4.9) y las condiciones

iniciales (4.10). Los coficientes de este sitema de ecuaciones satisfacen las condiciones para
poder aplicar el algoritmo de Thomas, y al hacer eliminacion obtenemos el sistema bi

Diferencias Finitas y E.D.P.

19

1111
0000

t n+1

1
0
0
1

1
0
0
1

11
00

tn

11111111
00000000
1
0
0
1
0
1
0
1
0
1
0
1

11
00

1
0
0
1

xj1 x j x j+1

Figura 7. Molecula computacional del esquema

diagonal
Uj ej Uj+1 = fj ,

ej =

,
1 + (2 ej1 )

fj =

dnj + fj1
,
1 + (2 ej1 )

(10.3)

j = 1, . . . , J 1, con e0 = f0 = 0 y dnj es el lado derecho de (10.2) para cada j. Este sistema

se resuelve eficientemente por medio de sustitucion regresiva.
Al final de las notas se incluye el programa teta1t.m con este esquema.
Estabilidad. Sustituyendo el modo discreto Ujn = n eik(jx) en la ecuacion (10.2), y simplificando obtenemos

ikx

ikx
21
1 = { + (1 )} e
2+e
= { + (1 )} 4 sen kx .
2

Despejando de la ecuacion anterior, obtenemos

1 4(1 ) sen2 ( 21 kx)
(k) =
,
1 + 4 sen2 ( 21 kx)

(10.4)

la cual es siempre menor a uno debido a que > 0 y 0 1. Por lo tanto puede haber
inestabilidad solo si < 1, es decir si
1
4(1 2) sen2 ( kx) > 2.
2
El modo mas inestable, con n
umero de onda k = J, maximiza el termino en el lado
izquierdo, pues en este caso 21 kx = /2. As que la condicion de inestabilidad es 4(1

Diferencias Finitas y E.D.P.

20

2) > 2. Por lo tanto

1
= inestabilidad
2
1
(1 2) = estabilidad
2
De lo anterior se tienen dos casos:
(1 2) >

Para 0 < < 1/2 , el metodo es estable si (1 2) 1/2. (Estabilidad condicional).

Para 1/2 < 1 , el metodo es estable para toda > 0 .

(Estabilidad incondicional).

Para cualquier > 0 se debe resolver un sistema tridiagonal y, a primera vista no se ve

1
ninguna ventaja de usar 0 < < , pues en estos casos el esquema es condicionalmente
2
1
estable, a menos que los resultados sean mas precisos que en los casos 1. Para esto
2
es necesario estudiar el error de truncamiento del metodo.
Error de truncamiento en el m
etodo . Para calcular el error de truncamiento del
esquema , que involucra seis puntos, se escoge el punto (xj , tn+ 1 ) como punto de expansion
2

en serie de Taylor. Este es el mejor punto, por simetra, y porque es el que cancela el mayor
de n
umero de terminos en la expansion. Con estas consideraciones el error de truncamiento
es:
n+ 12

Tj

n+1
n
n
n
un+1
unj
(un+1
+ un+1
j
j1 2uj
j+1 ) + (1 (uj1 2uj + uj+1 )

.
t
(x)2

Para expander alrededor del punto (xj , tn+ 1 ) se hace

2

un+1
j

unj

n+ 21
1
1 1
1 1
2
3
= u + tut + ( t) utt + ( t) uttt + . . .
2
2 2
6 2
j

n+ 21
1
1 1
1
1
u tut + ( t)2 utt ( t)3 uttt + . . .
2
2 2
6 2
j
n+ 12

1
1
(t)t uttttt . . .
.
= tut + (t)3 uttt +
24
1920
j

(10.5)

Diferencias Finitas y E.D.P.

21

Por otro lado, se tiene

un+1
j1

2un+1
j

un+1
j+1

n+1

1
1
4
6
2
(x) uxxxxxx + . . .
= (x) uxx + (x) uxxxx +
12
360
j

n+ 21
1
1
= (x)2 uxx + (x)4 uxxxx +
(x)6 uxxxxxx + . . .
12
360
j

n+ 21
1
1
1
2
4
6
+ t (x) uxxt + (x) uxxxxt +
(x) uxxxxxxt + . . .
2
12
360
j
n+ 12

1 1
1
1
2
2
4
6
+ ( t) (x) uxxtt + (x) uxxxxtt +
(x) uxxxxxxtt + . . .
2 2
12
360
j
+ ...

Analogamente

unj1

2unj

unj+1

n+ 21
1
1
4
6
= (x) uxx + (x) uxxxx +
(x) uxxxxxx + . . .
12
360
j

n+ 12
1
1
1
2
4
6
t (x) uxxt + (x) uxxxxt +
(x) uxxxxxxt + . . .
2
12
360
j

n+ 12
1 1
1
1
2
2
4
6
+ ( t) (x) uxxtt + (x) uxxxxtt +
(x) uxxxxxxtt + . . .
2 2
12
360
j
2

...
Entonces
n+1
n
n
n
(un+1
+ un+1
j1 2uj
j+1 ) + (1 )(uj1 2uj + uj+1 )

n+ 21
1
1
2
4
6
= (x) uxx + (x) uxxxx +
(x) uxxxxxx + . . .
12
360
j

n+ 12
1
1
1
2
4
6
+ ( )t (x) uxxt + (x) uxxxxt +
(x) uxxxxxxt + . . .
2
12
360
j

n+ 12
1
1
1
2
2
2
4
6
(x) uxxxxxxtt + . . .
+ (t) (x) (x) uxxtt + (x) uxxxxtt +
8
12
360
j

+ ...

Diferencias Finitas y E.D.P.

22

Por lo tanto, el error de truncamiento es

n+ 12
n+ 12
1
1
1
1
n+ 12
2
Tj
=
uxxt t
uxxx (x)
+
uttt uxxtt
(t)2
2
12
24
8
j
j

n+ 12
1 1
2
+
uxxxxt t
uxxxxxx (x)2
(x)2 + T OA .
12 2
6!
j

(10.6)

Claramente el esquema es consistente para toda y para toda . Ademas el error de

truncamiento muestra que normalmente se obtiene orden de convergencia lineal en t, pero
1
que en el caso simetrico especial = , el metodo es de segundo orden tanto en t como en
2
x
Esquema de CrankNicolson. El caso especial =
CrankNicolson

1
da origen al esquema que de
2

n+1
n+1
n
n
Ujn+1 Ujn
(Uj+1
2Ujn+1 + Uj1
) + (Uj+1
2Ujn + Uj1
)
=
,
t
2(x)2

(10.7)

el cual es incondicionalmente estable y de segundo orden tanto en x como en t. El error

de truncamineto de este metodo es

n+ 21
1
n+ 21
2
2
Tj
=
+ T OA .
uxxxx (x) + uttt (t)
12
j
Por lo tanto, podemos explotar las propiedades adicionales de estabilidad para tomar pasos
en el tiempo mayores, como por ejemplo t = x, y seguir obteniendo precision de segundo
orden.
Convergencia del metodo . Se puede hacer un analisis de convergencia para el metodo
semejante al hecho en la Seccion 6 para el metodo explcito.
Como en aquel caso definimos
en+1
:= Ujn+1 un+1
.
j
j

(10.8)

Despejando Ujn+1 de (10.1) y un+1

de (10.5) y sustituyendo en (10.8), se obtiene
j

n+1
n+1
n
n
2Ujn + Uj1
2Ujn+1 + Uj1
) + (1 )(Uj+1
en+1
= Ujn + (Uj+1
j

n+ 12
n+1
n
n
n
n
+ un+1
(un+1
j1 ) + (1 )(uj+1 2uj + uj1 ) uj tTj
j+1 2uj

n+ 12
n+1
n+1
n
n
n
= enj + (en+1

2e
+
e
)
+
(1

)(e

2e
+
e
)

tT
.
j+1
j
j1
j+1
j
j1
j

Diferencias Finitas y E.D.P.

23

Reordenando, se obtiene
n+ 12

n+1
n
n
n
(1 + 2)en+1
= (en+1
j
j1 + ej+1 ) + (1 )(ej1 + ej+1 ) + [1 2(1 )]ej tTj
n+ 21

Definiendo E n como en (6.5) y T como la cota superior de Tj

. (10.9)

sobre todos los puntos

(xj , tn+ 1 ), se obtiene

2

(1 + 2)E n+1 2E n+1 + E n + tT ,

debido a que todos los coeficientes en (10.9) excepto el u

ltimo son no negativos. Simplificando
la expresion anterior se obtiene
E n+1 E n + tT ,
y suponiendo que el error incial es cero, E 0 = 0, se obtiene el resultado deseado
E n ntT < tf T = 0

cuando t 0 ,

1
sobre toda trayectoria de refinamiento que satisfaga (1 ) .
2

11

Principio del m
aximo

Si se considera otras propiedades debe tener un esquema de aproximacion de la ecuacion

diferencial ut = uxx aparte de la convergencia, estabilidad, y orden razonable de presicion,
lo siguiente sera pedir que satisfaca un principio del m
aximo. Se sabe tanto desde el
punto de vista matematico y como fsico que la soluci
on u(x, t) esta acotada inferiormente y
superiormente por los extremos de los valores iniciales y los valores de frontera. Es deseable
que el esquema discretizacion tenga esta propiedad. Los valores extremos sobre la malla
computacional se muestran en la siguiente Figura 8.
1
Resultado: El metodo con 0 1 y (1 )
produce la sucesion {Ujn } con
2
Umin Ujn Umax , donde
Umin = min { U0m , UJm , 0 m n; Uj0 0 j J },
Umax = max { U0m , UJm , 0 m n; Uj0 0 j J }.
Demostraci
on. Despejando Ujn+1 de (10.1), se obtiene

n+1
n+1
n
n
Ujn+1 = Ujn + (Uj+1
2Ujn+1 + Uj1
) + (1 ) (Uj+1
2Ujn + Uj1
) .
Reagrupando, se obtiene
n+1
n+1
n
n
(1 + 2)Ujn+1 = (Uj1
+ Uj+1
) + (1 )(Uj1
+ Uj+1
) + [1 2(1 )]Ujn

Diferencias Finitas y E.D.P.

24

1111
0000

t n+1

11
00

tn

11
00

xj1 x j x j+1

Figura 8. Valores extremos sobre la malla computacional.

Observese que todos los coeficientes en la relacion anterior son no negativos ya que 0,
> 0, y 2(1 ) 1. Supongase ahora que U toma su maximo en un punto interno, y que
n+1
n+1
n
n
este es Ujn+1 , y sea U el mayor de los cinco valores adyacentes Uj1
, Uj+1
, Uj1
, Ujn , Uj+1

en la ecuacion anterior. Entonces Ujn+1 U debido a que los coeficientes en el lado derecho
suman 1 + 2. Por otro lado U Ujn+1 , pues Ujn+1 es el valor maximo por hipotesis. Se
concluye que U = Ujn+1 . El argumento puede repetirse en cada uno de los puntos U as
obtenidos, obteniendo una sucecion de puntos hasta que se toma un punto frontera. Un
argumento similar muestra que el mnimo tambien se toma en un punto frontera.
1
La condicion del resultado anterior, (1 ) , para garantizar el principio del
2
1
maximo es mucho mas restrictiva que la obtenida por analisis de Fourier, (1 2) ,
2
1
para garantizar estabilidad. Por ejemplo, el esquema de CrankNicolson ( = ) es estable
2
para toda , pero satisface el principio del maximo solo si 1.
El analisis del principio del maximo puede usarse como medio alternativo para obtener
condiciones de estabilidad. Este tipo de analis tiene la ventaja sobre el analisis de Fourier de
que es facilmente aplicable a problemas con coeficientes variables, entre otros. Sin embargo,
con el principio del maximo solo pueden derivarse condiciones suficientes para estabilidad.
Ejemplo 5. Con el objeto de ilustrar el efecto de las condiciones para que se satisfaga un
principio del maximo, consiremos el problema donde la distribucion incial de temperaturas
en la barra tiene la forma de un pico, es decir u0 (x) es 1 en el punto medio de una malla
con J = 20 y 0 en cualquier otro punto de la malla. Si tomamos = 1/2, enonces para

Diferencias Finitas y E.D.P.

25

0.9

0.9

0.8

0.8

0.7

0.7

0.6

0.6

0.5

0.5

0.4

0.4

n=0

0.3

0.2

0.1
0

0.1
0

0.2

0.4

0.6

0.8

n=2

0.2

0.4

0.6

0.8

0.4

0.6

0.8

0.4

0.6

0.8

0.4

0.6

0.8

n=1

0.2

0.4

0.6

0.8

0.5

0.2

0.5

0.5

n=4

n=2

0.2

0.4

0.6

0.8

0.5

0.2

0.5

0.5

n=6

n=3
0

0.5

0.5

0.5

0.5

0.5

n=0

0.3

0.2

0.5
0

0.2

0.4

0.6

0.8

t = 0.0025 ( = 1 )

0.2

t = 0.005 ( = 2 )

Figura 9. Resultados con = 1/2, J = 10 y = 1 y = 2.

t = 0.0025 ( = 1) se satisface la condicion para que el esquema satisfaga el principio del
maximo, mientras que para t = 0.005 ( = 2) no se satisface la condicion. Los resultados
se muestran el la Figura 9. Como puede observarse, en el segundo caso, en el primer nivel
del tiempo la solucion es negativa en el centro, despues tiende a corregirse conforme avanza
el tiempo. Esto no sucede en el primer caso, donde si se tiene la condicion necesaria 1
para que el esquema satisfaga el principio del maximo. Se utilizo el programa teta1t.m.

Diferencias Finitas y E.D.P.

12

26

Problemas m
as generales

12.1

Coeficientes variables

En el modelo de la ecuacion de flujo de calor en la barra, se ha supuesto que las propiedades

del material conductor son constantes con el tiempo e independientes de x. En realidad,
en muchos problemas estas propiedades no son constantes y dependen cuando menos de x.
Consideremos el problema mas general
ut = a(x, t) uxx ,

con a(x, t) > 0,

x, t.

(12.1)

k
, donde k denota el coeficiente de conductividad termica, c el calor especfico, y
c
la densidad del material.
con a =

An
alogo del m
etodo . En el caso de coeficientes variable el esquema de promedios
pesados o es
n+1
n+1
n
n
Ujn+1 Ujn
Uj+1 2Ujn+1 + Uj1
Uj+1
2Ujn + Uj1

=a
+ (1 )
,
t
(x)2
(x)2

(12.2)

donde a puede escogerse de diferentes maneras. Algunas de las elecciones adecuadas son

1111
0000

t n+1

1
0
0
1

1
0
0
1

11
00

tn

11111111
00000000
1
0
0
1
0
1
0
1
0
1
0
1

11
00

1
0
0
1

xj1 x j x j+1

Figura 10. El coeficiente a se evalua en el punto central.

n+ 21

a = aj

). Con la primera eleccion , a(x, t) se evalua en el punto

, y a = 12 (anj + an+1
j

intermedio de la molecula computacional (ver Figura 10), y se puede probar que el error de

Diferencias Finitas y E.D.P.

27

truncamiento no se altera excepto por la inclusion del extrafactor a:

n+ 12
n+ 12
1
a
1
a
n+ 12
2
Tj
=
uxxt t
uxxx (x)
+
uttt uxxtt
(t)2
2
12
24
8
j
j

n+ 12
1 1
2a
+
uxxxxt t
uxxxxxx (x)2
(x)2 + T OA .
12 2
6!
j
Se observa inmediatamente que el esquema es consistente, y que el orden del mismo es O(t)
1
en el caso general, y O((t)2 + (x)2 ) en el caso que = (CrankNicolson).
2
La ecuacion en diferencias (12.2) se puede escribir en la forma
n+1
n+1
n
n
+Uj1
) . (12.3)
(1+2 a )Ujn+1 = a (Uj1
+Uj+1
)+[12 a (1) ] Ujn +a (1) (Uj+1

De esta ecuacion se observa que los coeficientes del lado derecho suman lo mismo que el coeficiente del lado izquierdo, y que se satisface un principio del maximo si todos los coeficientes
son no negativos. Esto u
ltimo se logra si 1 2 a (1 ) 0. Por lo tanto, la convergencia
del m
etodo se puede realizar por medio del principio del maximo, pero se debe adoptar la
condici
on de estabilidad
(1 ) a(x, t) 1/2 x, t .

En ocasiones es conveniente usar la segunda alternativa a =

(12.4)

1 n
(a
2 j

+ an+1
), la cual
j

produce un error de truncamiento del mismo orden. Por el contrario, el siguiente esquema
n+1
n
U n+1 2Ujn+1 + Uj1
Ujn+1 Ujn
U n 2Ujn + Uj1
n+1 j+1
n j+1
= aj
+ (1 ) aj
.
t
(x)2
(x)2

produce un error de truncamiento O(t) independientemente del valor de . Es decir, no es

posible aumentar el orden para alg
un valor especial de , motivo por el cual este esquema
no es recomendable.

12.2

Ecuaci
on parab
olica en forma autoadjunta

En las aplicaciones, muy frecuentemente la ecuacion de difusion aparece en la forma

ut = (a ux )x ,

a = a(x, t) > 0 ,

(12.5)

Diferencias Finitas y E.D.P.

28

denominada forma autoadjunta. Es posible escribir esta ecuacion como ut = a uxx + ax ux ,

y tratar de construir un esquema de aproximacion a partir de esta expresion. Sin embargo,
es mejor construir una aproximacion de la ecuacion respetando su forma original. Para
construir una aproximacion, hacemos la siguiente definicion
x Uj := Uj+ 1 Uj 1 .
2

(12.6)

donde x es el operador de diferencias centrales con respecto a la variable de posicion x.

Entonces, el esquema explcito para el problema en forma autoadjunta es
Ujn+1 Ujn
x [anj x (Ujn )]
=
,
t
(x)2
el cual, como puede observarse, respeta la simetria del problema. Desarrollando los operadores de diferencias centrales x obtenemos el esquema explcito
n
n
anj+ 1 (Uj+1
Ujn ) + anj 1 (Ujn Uj1
)
Ujn+1 Ujn
2
2
=
.
t
(x)2

(12.7)

La molecula computacional de este esquema se muestra en la Figura 11, donde se han

dibujado simbolos para indicar los puntos donde se evalua la funcion a(x, t).

1111
0000

11
00

t n+1
tn

11111111
00000000
1
0
0
1
0
1
0
1
0
1
0
1

xj1 x j x j+1

11
00

Figura 11. Molecula computacional para el problema en forma adjunta.

El error de truncamiento de este metodo viene dado por

n
n
a uxxxx ax uxxx axx uxx axxx ux
1
n
utt t
+
+
+
(x)2 + . . .
Tj =
2
12
6
8
24
j
j
donde, como se esperaba, el error es de primer orden en t.

(12.8)

Diferencias Finitas y E.D.P.

29

Para estudiar la estabilidad del metodo reescribimos (12.7) en la forma

h
i
n
n
Ujn+1 = 1 (anj+ 1 + anj 1 ) Ujn + anj+ 1 Uj+1
+ anj 1 Uj1
.
2

Esta ecuacion satisface un principio del maximo si los coeficientes son acotados y no negativos. Como hemos supuesto que a(x, t) > 0 para toda x y t, entonces el principio del
maximo se satisface si.
1
donde A = max{a(x, t)} .
A ,
2
Esta condicion es una condicion necesaria para la estabilidad del metodo explcito, lo cual a
su vez asegura la convergencia del metodo, gracias a la consistencia del mismo.
Un esquema analogo al metodo es
Ujn+1 Ujn
x [an+1
x (Ujn+1 )]
x [anj x (Ujn )]
j
=
+ (1 )
.
t
(x)2
(x)2
Desarrollando los operadores de diferencias centrales y reordenando el resultado, se obtiene
h
i
i
h
n+1
n+1
n+1
n
n
1 + (aj 1 + aj+ 1 ) Uj = 1 (1 )(aj 1 + aj+ 1 ) Ujn + (an+1
U n+1 + an+1
U n+1 )
j 1 j1
j+ 1 j+1
2

+ (1 )

(anj 1
2

n
Uj1

n
+ anj+ 1 Uj+1
).
2

Este esquema satisface un principio del maximo si

(1 )A

12.3

1
,
2

donde A = max{a(x, t)} .

Problemas no lineales

Considerese la ecuacion diferencial

ut = a(u) uxx ,

a(u) > 0 ,

en donde el coeficiente a es una funcion de la incognita u. El coeficiente a tambien podra

depender de x y t directamente, pero este caso no es significativamente mas difcil. El
esquema mas sencillo para este problema es el siguiente esquema explcito
n
n
Ujn+1 Ujn
2Ujn + Uj+1
Uj1
n
= a(Uj )
.
t
(x)2

(12.9)

Diferencias Finitas y E.D.P.

30

El error de truncamiento de este metodo esta dado por

n
1
1
n
n
2
Tj =
utt t
a(Uj ) uxxxx (x) + . . . O(t) ,
2
12
j
de donde se concluye que el esquema es consistente y de primer orden en t. Despejando
Ujn+1 de (12.9), se obtiene
Ujn+1

n
n
n
n
n
= 1 2 a(Uj ) Uj + a(Uj ) Uj+1 + Uj1 .

(12.10)

Los coeficientes en el laso derecho suman uno y son no negativos si 1 2 a(Ujn ) 0 para
toda j y n. Por lo tanto de satisface un principio del m
aximo si
1
donde A = max{a(Ujn )}.
(12.11)
2
Sin embargo, el analisis del error no es triival en este caso, debido a que aparece un problema
si

que no surge en el caso de ecuaciones lineales: el error se propaga de manera nolineal

cuando n crece. Mostraremos brevemente las dificultades que aparecen al intentar hacer el
analisis del error. De la ecuacion (12.10) y de la definicion del error de truncamiento se
obtiene
Ujn+1 = Ujn + a(Ujn ) x2 Ujn
un+1
= unj + a(unj ) x2 unj + t Tjn .
j
y no pueden restarse estas ecuacione en forma directa para obtener el error enj = Ujn+1 un+1
j
debido a que a(Ujn ) 6= a(unj ). Pero se puede linealizar a(unj ) por medio de su serie de Taylor
del primer orden
a(unj ) = a(Ujn ) + (un+1
Ujn+1 )
j
a n
=

( ) ,
u j
en donde jn se encuentra entre unj y UJn . Por lo tanto
a(Ujn )

a n
( )
u j

(12.12)

enj

a
(x2 unj ) enj t Tjn
u
(12.13)
n

a n n
n
n
n
n
n
( ) e t Tj .
= 1 2 a(Uj ) ej + a(uj ) (ej1 + ej+1 ) +
u j j
a
Para poder garantizar convergencia se debe controlar el crecimiento de
. Suponiendo que
u
a
existe K > 0 tal que | | K, y que |unj+1 2unj + unj1 | Mxx (x)2 , donde Mxx es una
u
en+1
= enj + a(Ujn ) x2 enj +
j

Diferencias Finitas y E.D.P.

31

cota superior de |uxx |, entonces

E n+1 [1 + K Mxx t] E n + t T
[1 + K Mxx t]n+1 E 0
+ {1 + [1 + K Mxx t] + [1 + K Mxx t]2 + . . . + [1 + K Mxx t]n }t T
(12.14)
[1 + K Mxx t]n+1 1

t T
K Mxx t
eK Mxx tf

T 0 cuando t, x 0 .
K Mxx
Se observa que la estimacion del error es mucho peor que en el caso lineal, a menos que K
y Mxx sean peque
nos.
Esquemas implcitos. Muchos autores prefieren utilzar esquemas implcitos o tambien
esquemas semiimplcitos para aproximar problemas no lineales. El esquema implcito mas
simple es
n+1
Ujn+1 = Ujn + a(Ujn+1 ) (Ujn+1 2 Ujn+1 + Uj1
).

(12.15)

Se podra hacer una aproximacion semiimplcita tomando a(Ujn ) en lugar de a(Ujn+1 ) en la

ecuacion anterior. El sistema resultante (12.15) es un sistema algebraico no lineal, y debe
utilizarse alg
un metodo iterativo o alg
un tipo de linealizacion para resolverlo. Dentro de los
metodos iterativos mas comunes se encuentran los metodos de tipo Picard (o punto fijo), y el
metodo de Newton. El valor inicial para comenzar las iteraciones comunmente es la solucion
en el paso de tiempo anterior. Por ejemplo, un metodo tipo Picard para resolver (12.15) es
el siguiente
Dado Ujn+1, 0 = Ujn
Para k = 1, 2, . . . hasta convergencia
n+1, k
Calcular Ujn+1, k = Ujn + a(Ujn+1, k1 ) (Ujn+1, k 2 Ujn+1, k + Uj1
)

Finalizar
kU n+1, k U n+1, k1 k
, donde es un
kU n+1, k k
o menor. Observese que en cada iteracion del tiempo

Se paran las iteracione con alg

un criterio como
numero peque
no, por ejemplo 106

se deben resolver kmin sistema de ecuaciones lineales, donde kmin es el numero mnimo de
itreaciones para obtener la tolerancia en el criterio de paro.
Algunos metodos alternativos estan basados en algoritmos tipo predictorcorrector. Por

Diferencias Finitas y E.D.P.

32

ejemplo
n+1/2

Predictor:
Corrector:

n+1/2

x2 Uj
Ujn
= a(Ujn )
t/2
(x)2
2
n+1
Ujn+1 Ujn
+ Ujn )
n+1/2 x (Uj
= a(Uj
)
t
2 (x)2

Uj

(semiimplcito)
(CrankNicolson)

Se pueden utilizar otras variantes. Observese que en este caso tambien hay que resolver dos
sistemas algebraico linelaes en cada iteracion del tiempo.

13

Problemas bidimensionales

Consideremos la caja rectangular D = 0, X] [0, Y ] y denotemos por D a su frontera. Sea

a > 0 un coeficiente de difusion, que consideraremos constante por el momento. El problema parabolico bidimensional (problema de difusion), con condiciones de frontera Dirichlet,
consiste en calcular u = u(x, y, t) tal que
ut = M u = a (uxx + uyy ),

en D, t > 0,

(13.1)

u(x, y, 0) = u0 (x, y)

en D,

(13.2)

u(x, y, t) = g(x, y, t)

sobre D, t 0.

(13.3)

Construimos una malla uniforme sobre el dominio rectangular D de la siguiente manera:

Se divide [0, X] en I subintervalos y [0, Y ] en J subintervalos Sean x = X/I y y = Y /J
los pasos de discretizacios. Los puntos de la malla sobre D son
xi = i x,

i = 0, 1, . . . , I,

yj = j y,

j = 0, 1, . . . , J.

Dividimos el intervalo de tiempo [0, tf ] en N subintervalos de tal manera que los niveles de
tiempo discreto son tn = n t, n = 0, 1, . . . , N.
n
la aproximacion a la solucion exacta u(xi , yj , tn ).
Esquema explcito. Denotamaos por Ui,j

El exquema explcito bidimensional es

n

n+1
n
n
n
n
n
n
Ui,j
Ui,j
+ Ui+1,j
+ Ui,j+1
Ui1,j 2Ui,j
2Ui,j
Ui,j1
=a
+
.
t
(x)2
(y)2

(13.4)

Diferencias Finitas y E.D.P.

33

1111
0000

yj+1

1
0
0
1

11
00

yj

11111111
00000000
1
0
0
1
0
1
0
1
0
1
0
1

11
00

yj1

1
0
0
1

xi1 x i x i+1

Figura 12. Molecula computacional para el esquema explcito bidimensional

Con el objeto de simplificar la exposicion introducimos la siguiente notacion:
U := Ui,j
x2 U := Ui1,j 2Ui,j + Ui+1,j
y2 U := Ui,j1 2Ui,j + Ui,j+1
Con la notacion anterior el esquema explcito se puede escribir como
2 n

y2 U n
x U
U n+1 U n
=a
+
,
t
(x)2 (y)2

(13.5)

o bien
U n+1 = U n + x x2 U n + y y2 U n ,

(13.6)

donde x = a t/(x)2 y y = a t/(y)2 .

Casi todo el analisis en una dimension puede extenderse facilmente al caso de dos dimensiones.
Error de truncamiento.

a
1
2
2
T (x, t) = utt t
uxxxx (x) + uyyyy (y) + . . . O(t)
2
12
Por tanto, el esquema es consistente y de orden t.
Convergencia. Bajo la condicion x + y 1/2, se satisface

a2 t
1
n
Mtt t +
E
Mxxxx + Myyyy tf 0 ,
2
12

(13.7)

Diferencias Finitas y E.D.P.

34

cuando t 0.
Estabilidad. La estabilidad tambien puede analizarse mediante series de Fourier, suponiendo
que a es constante, e ignorando las condiciones de frontera. En este caso los modos de
Fourier son de la forma
U n = n ei(k1 x+k2 y) ,

(13.8)

donde k1 y k2 son los n

umeros de onda en las direcciones x y y, respectivamente. Sustituyendo (13.8) en (13.6), se obtiene que el factor de amplificaci
on es

2 1
2 1
= (k1 , k2 ) = 1 4 x sen ( k1 x) + y sen ( k2 y) .
(13.9)
2
2
El peor modo discreto es aquel en el que k1 x = k2 y = , es decir k1 = I y k2 = J. De
aqu que

1
(13.10)
|| 1 si x + y .
2
Esta condici
on de estabilidad es m
as severa que en el caso unidimensional. Por ejemplo,
a t
1
si x = y = h, entonces = 2 . Ademas, cuando a sea variable, en cualquier pico
h
4
de ella sera necesario disminuir el paso del tiempo para mantener la estabilidad. Entonces,
el esquema es algo impractico en 2D, y mas a
un su generalizacion a 3D. Por lo tanto, debe
buscarse equemas implcitos o semiimplcitos para relajar la condicion de estabilidad.

14

El m
etodo

Una posible alternativa para superar las limitaciones del esquema explcito bidimensional es
el metodo de promedios pesados (metodo ):
y2 ( U n+1 + (1 ) U n )
U n+1 U n
2 ( U n+1 + (1 ) U n )
=a x
+
a
.
t
(x)2
(y)2

(14.1)

Para el caso unidimensional hemos visto que el esquema es incodicionalmente estable para
1/2, y es de segundo orden tanto en t como en x si = 1/2. Entonces vale la pena
considerar el metodo
CrankNicolson. Haciendo = 1/2 y reordenando (14.1), obtenemos la ecuacion en
diferencias

1
1
1
1
(14.2)
x x2 y y2 ) U n+1 = (1 + x x2 + y y2 ) U n .
2
2
2
2
Para verificar estabilidad, sustituimos el modo (13.8) en esta ecuacion, y despues de simpli(1

Diferencias Finitas y E.D.P.

35

ficar y despejar, encontramos que el factor de amplificacion en este caso es

(k1 , k2 ) =

1 2 x sen2 (k1 x/2) 2 y sen2 (k2 y/2)

.
1 + 2 x sen2 (k1 x/2) + 2 y sen2 (k2 y/2)

(14.3)

Inmediatamente se sigue que |(k1 , k2 )| 1 para todo k1 , k2 , y para todo x, y, t. Por

lo tanto, el esquema bidimensional de CrankNicolson es incondicionalmente estable.
En una dimension la gran ventaja de este esquema fue relajar la condicion de estabilidad
a un extra costo relativamente peque
no, debido a que se cuenta con un resolvedor rapido
y eficiente para sistemas tridiagonales, que es el algoritmo de Thomas. Sin embargo, en el
caso de dos dimensiones, a pesar que el algoritmo sigue siendo incondicionalmente estable, el
trabajo extra respecto al caso explcito es considerable. En el problema bidimensional, para
cada n = 1, 2, . . ., se debe resolver un sistema lineal con (I 1) (J 1) ecuaciones con el
n+1
mismo n
umero de incognitas: Ui,j
. Para ilusttrar, supongamos que escogemos x = y,

es decir x = y , entonces la matriz del sistema es:

A I
4 1
I A I

1 4 1

con A = . . . . . . . . . . . . . . . . . . . . . .
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
A = I +

I A I
1 4 1
I A
1 4
Como puede observarse la matriz del sistema A viene dada por bloques, y cada ecuacion
involucra cinco incognitas, y no hay forma de permutar renglones y columnas para que los
coeficientes nocero formen una banda delgada, de tal forma que podamos hacer eliminacion
y reduccon como en el caso unidimensional. En conclusion, no es posible resolver el sistema
en forma eficiente en el caso bidimensional, y el problema es todava mas severo en el caso
tridimensional.

15

M
etodos implcitos de direcciones alternantes (ADI)

Como los metodos implcitos son muy eficientes en una dimension, es natural buscar metodos
bidimensionales que sean implcitos en una direccion pero no en ambas. El siguiente esquema,
es un ejemplo de este tipo de esquemas

Diferencias Finitas y E.D.P.

36

PeacemanRachford (1955)
y2 U n
U n+1/2 U n
2 U n+1/2
=a x
+
a
t/2
(x)2
(y)2
y2 U n+1
2 U n+1/2
U n+1 U n+1/2
=a x
+
a
t/2
(x)2
(y)2

n t

(n+ 12)t

(
(

IMPLICITO EN x
EXPLICITO EN y
EXPLICITO EN x
IMPLICITO EN y

(15.1)
(15.2)

(n+1) t

Este esquema puede escribirse en la forma

1
1
x x2 ) U n+1/2 = (1 + y y2 ) U n ,
(15.3)
2
2
1
1
(1 y y2 ) U n+1 = (1 + x x2 ) U n+1/2 ,
(15.4)
2
2
el cual constituye un sistema de dos problemas implcitos unidimensionales. La solucion de
cada uno de estos sistemas involucra solo conjuntos de ecuaciones tridiagonales. El trabajo
(1

total, en un paso del tiempo, necesita de la solucion de J 1 sistemas tridiagonales de orden

I 1, seguido de la solucion de I 1 sistemas tridiagonales de orden J 1. Este proceso
de solucion de conjuntos de sistemas tridiagonales, es mucho mas rapido que la solucion del
problema completo de (I 1)(J 1) ecuaciones en el metodo de CrankNicolson. El trabajo
requerido es aproximadamente tres veces el del metodo explcito. La siguiente figura ilustra
el orden en el que se pueden resolver estos sistemas tridiagonales. A este tipo de esquemas se
les conoce comunmente como esquemas implcitos de direcciones alternantes o ADI
por sus siglas en el idioma Ingles (alternating direction implicit).
Para hacer el analisis del esquema, eliminamos U n+1/2 en las ecuaciones (15.3), (15.4).
Se obtiene la siguiente ecuaci
on en diferencias
1
1
1
1
(15.5)
(1 x x2 )(1 y y2 ) U n+1 = (1 + x x2 )(1 + y y2 ) U n .
2
2
2
2
Sustituyendo el modo de Fourier (13.8) en (15.5), simplificando, y despejando, se obtiene el
factor de amplificaci
on
(1 2 x sen2 (k1 x/2)) (1 2 y sen2 (k2 y/2))
(k1 , k2 ) =
.
(1 + 2 x sen2 (k1 x/2)) (1 + 2 y sen2 (k2 y/2))
de donde se sigue inmediatamente la estabilidad incondicional del metodo ADI.

(15.6)

Diferencias Finitas y E.D.P.

37

1111
0000

11
00

11
00

yj+1
yj

1
0
0
1

1
0
0
1

1
0
0
1

1
0
0
1

J1
sistemas
tridiagonales
de orden
I1

1
0
0
1

yj1

xi1 x i x i+1

I1 sistemas tridiagonales de orden J1

Figura 13. Ilustracion del metodo ADI para problemas bidimensionales

Principio del m
aximo: En la ecuacion (15.3), al desarrollar, se obiene
y
x n+1/2
n+1/2
n+1/2
n
(1 + x ) Ui, j
= (1 y ) Ui,nj + (Ui,nj1 + Ui+1,
(Ui+1, j + Ui1, j ) .
j)
2
2
n+1/2

De esta expresion vemos que si y 1, entonces Ui,j

n

vacinos de U y U

n+1/2

es combinacion lineal de valores

con coeficientes no negativos que suman 1. Se obtiene una relacion

analoga para la ecuacion (15.4). Por lo tanto si

max{x , y } 1 ,
entonces se asegura que el esquema ADI satisface el principio del maximo.
Truncamiento. Observese que el esquema ADI (15.5) puede escribirse como
1
1
x y 2 2 n+1
1
1
1
(1 x x2 y y2 +
x y ) U
= (1 + x x2 + y y2 + x y x2 ty 2 ) U n . (15.7)
2
2
4
2
2
4
As que el error de truncamiento es
n+1
n+1/2
u
un ax2 (un+1 + un ) ay2 (un+1 + un ) t a2 x2 y2 (un+1 un )
n+1/2

+
.
Ti,j
=
t
2 (x)2
2 (y)2
4(x)2 (y)2
i,j
(15.8)

Diferencias Finitas y E.D.P.

38

Se puede verificar que

n+1/2
un+1 un
(t)2 uttt
= ut +
+ ...
t
24
i,j

n+1/2
2 n+1
n
2
x (u
+u )
(x) uxxxx (t)2 uxxtt
+
+ ...
= uxx +
2 (x)2
12
12
i,j

n+1/2
2 n+1
n
2
2
y (u
+u )
(y) uyyyy (t) uyytt
= uyy +
+
+ ...
.
2
2 (y)
12
12
i,j

(15.9)
(15.10)
(15.11)

Por otro lado, gracias a (15.9) se tiene

n+1/2
t x2 y2 (un+1 un )
(t)2 uttt
(t)2
2 2
=
ut +
+ ...
(x)2 (y)2
(x)2 (y)2 x y
24
i,j
2 v
y dado que para cualquier funcion v = v(x, y, t) suficientemente suave se satisface x 2 =
(x)

n+1/2
2
(x)
entonces
vxx +
vxxxx + . . .
12
i,j

t x2 y2 (un+1 un )
2
2
2
2 n+1/2
=
(t)
u
+
O(t)
+
O(x)
+
O(y)
.
xxyyt
i,j
(x)2 (y)2

(15.12)

Por lo tanto sustituyendo las expresiones (15.9)(15.12) en (15.8) se obtiene que el error de
truncamiento es

n+1/2
1 1
n+1/2
2
2
2
uttt (t) uxxxx (x) uyyyy (y)
Ti,j
=
12 2
i,j

n+1/2
1 1
1
2
2
2
(15.13)
+ T OA

uxxtt (t) + uyytt (t) uxxyyt (t)

4 2
2
i,j

2
2
2
= O (t) + (x) + (y) .
Se concluye entonces, que el esquema implicito de direcciones alternantes es de segundo
orden respecto a t, x y y.
Ejemplo 6. Se considera el problema bidimensional con a = 1 y condiciones inciales dadas
por la funcion peaks de MATLAB. La Figura 14 ilustra esta funcion cuando n = 0.
Las condiciones de frontera son u = 0 sobre D. Se toman I = J = 48, y t = 0.0025.
Observese que las condiciones iniciales toma valores positivos y negativos. En la misma figura
se muestra la solucion para varios niveles de tiemp n. Se utilizo el programa adicn.m, que
consiste del metodo ADI con = 1/2 (CrankNicolson).

Diferencias Finitas y E.D.P.

39

n=0

n=3

10

10
1

5
1

0.8

0.5

0.8

0.5

0.6

0.6

0.4
Y

0.4

0.2
0

0.2
0

n = 25

n=6
5

5
1

5
1
1
0.8

0.5

0.6

1
0.8

0.5

0.6

0.4
Y

0.4

0.2
0

0.2
0

Figura 14. Resultados numericos para un problema bidimensional obtenidos con el metodo
ADI

Ejemplo 7. En este u
ltimo ejemplo se toma como condicion inicial una funcion con dos
torres en forma simetrica (Ver Figura 15 con n = 0), y condiciones de frontera cero. Se
toma I = J = 40 y dos casos: t = 0.0005, y t = 0.0025. El primer caso proporciona
x = y = 0.8, y satisface la condicion para el principio del maximo para el esquema ADI.
El segundo caso, por el contrario, proporciona x = y = 4, y no se satisface el principio
del maximo. Notese que, a pesar de que el metodo es estable en ambos casos, la solucion se
degrada en el caso que no se satisface la condicion para el principio del maximo.

16

Programas

Nota. En los comentarios de estos programas no se incluyen acentos.

Diferencias Finitas y E.D.P.

40

n=0

n=0

0.8

0.8

0.6

0.6

0.4

0.4

0.2

0.2

0.2
1

0.2
1
1

0.8

0.5

0.8

0.5

0.6

0.6

0.4
Y

0.4

0.2
0

0.2
0

n=5
1

0.8

0.8

0.6

0.6

0.4

0.4

0.2

0.2

0.2
1

0.2
1

n=1

1
0.8

0.5

0.8

0.5

0.6

0.6
0.4

0.4
Y

0.2
0

n = 30

0.2
0

n=2

0.8

0.8

0.6

0.6

0.4

0.4

0.2

0.2

0.2
1

0.2
1
1

1
0.8

0.5

0.8

0.5

0.6

0.6
0.4

0.4
Y

0.2
0

n = 100

0.2
0

n=5

0.8

0.8

0.6

0.6

0.4

0.4

0.2

0.2

0.2
1

0.2
1
1

1
0.8

0.5

0.6

0.8

0.5

0.6
0.4

0.4
Y

0.2
0

0.2
0

Figura 15. Ilustracion del principio del maximo en el metodo ADI. Lado izquierdo:
t = 0.0005 (x = y = 0.8). Lado derecho t = 0.0025 (x = y = 4)

Diferencias Finitas y E.D.P.

function M = explicito1(J,dt,N)
% METODO DE EULER EXPLICITO PARA EL MODELO PARABOLICO UNIDIMENSIONAL
%
%
%
%
%
%
%

ARGUMENTOS DE ENTRADA
J = numero de subdivisiones de [0,1]
dt = paso del tiempo
N = numero de pasos del tiempo
ARGUMENTO DE SALIDA
M = almacena los distintas instantaneas graficas para hacer una
animacion de la solucion

dx = 1/J;
x = 0:dx:1;

% Paso de discretizacion del intervalo [0, 1]

% Conjunto de valores x_j en la malla

% Condiciones iniciales: u = 2x
para x entre 0 y 0.5
%
u = 2 2x para x entre 0.5 y 1
j = 1:J/2; u(j) = 2*x(j);
j = J/2+1:J+1; u(j) = 2 2*x(j);
plot(x,u,r) % Grafica de las condiciones iniciales
axis([0 1 0 100])
pause
nu = dt/dx^2;
for n = 1:N
% Iteraciones en el nivel del tiempo
Unew(1) = 0; % Condicion de frontera
for j = 2:J % Esquema de diferencias.
Unew(j) = u(j) + nu*( u(j1)2*u(j)+u(j+1) );
end
Unew(J+1) = 0; % Condicion de frontera
u = Unew; % Actualizacion de la solucin
plot(x,u)
axis([0 1 0 100])
M(n) = getframe
end

41

Diferencias Finitas y E.D.P.

function A = implicito1(J,dt,N)
% METODO DE EULER IMPLICITO PARA EL MODELO PARABOLICO UNIDIMENSIONAL
% CON CONDICIONES DE FRONTERA DIRICHLET HOMOGENENAS
%
%
%
%
%
%
%
%

ARGUMENTOS DE ENTRADA
J = numero de subdivisiones de [0,1] (J debe ser par)
dt = paso del tiempo
N = numero de pasos del tiempo
ARGUMENTO DE SALIDA
A = la matriz del sistema lineal en la ecuacion en diferencias
Condiciones de frontera homogeneas en x = 0 y x = 1

dx = 1/J; % Paso de la discretizacion en el intervalo [0, 1]

x = 0:dx:1; % Conjunto de valores x_j en la malla
% Condiciones iniciales u = 2x
para x entre 0 y 0.5
%
u = 22x para x entre 0.5 y 1
u = zeros(J+1,1);
j = 1:J/2+1; u(j) = 2*x(j);
j = J/2+2:J+1; u(j) = 2 2*x(j);
plot(x,u,r) % Grafica de las condiciones iniciales
axis([0 1 0 1])
pause
% CONSTRUCCION DE LA MATRIZ DEL SISTEMA
nu = dt/dx^2;
A = (1+2*nu)*eye(J1);
A(1,2) = nu;
for j = 2:J2;
A(j,j1) = nu;
A(j,j+1) = nu;
end
A(J1,J2) = nu;
% Iteraciones en el tiempo
unew = zeros(J+1,1);
for n = 1:N,
unew(1) = 0;
% Condicion de frontera en x = 0
unew(2:J) = A\u(2:J); % SOLUCION DEL SISTEMA LINEAL
unew(J+1) = 0; % Condicion de frontera en x = 1
u = unew; % Actualizacion para avanzar en el tiempo
plot(x,u)
axis([0 1 0 1])
M(n) = getframe;
end

42

Diferencias Finitas y E.D.P.

function implicito1t(J,dt,N)
% METODO DE EULER IMPLICITO PARA EL MODELO PARABOLICO UNIDIMENSIONAL
% CON CONDICIONES DE FRONTERA DIRICHLET HOMOGENENAS
% SE UTILIZA EL METODO DE THOMAS
%
%
%
%
%
%

ARGUMENTOS DE
J = numero de
dt = paso del
N = numero de
ARGUMENTOS DE
Nunguno

ENTRADA
subdivisiones de [0,1] (J debe ser par)
tiempo
pasos del tiempo
SALIDA

dx = 1/J; % Paso de la discretizacion en el intervalo [0, 1]

x = 0:dx:1; % Conjunto de valores x_j en la malla
% Condiciones iniciales u = 2x
para x entre 0 y 0.5
%
u = 22x para x entre 0.5 y 1
j = 1:J/2+1; u(j) = 2*x(j);
j = J/2+2:J+1; u(j) = 2 2*x(j);
plot(x,u,r) % Grafica de las condiciones iniciales
axis([0 1 0 1])
pause % Pausa
nu = dt/dx^2;
% Coeficientes e(j) del sistema normalizado
e(1) = 0;
for j = 2:J
e(j) = nu/(1+nu*(2e(j1)));
end
% Iteraciones en el tiempo
unew = zeros(J+1,1);
for n = 1:N
f(1) = 0;
for j = 2:J % Lados derechos del sistema normalizado
f(j) = (u(j)+nu*f(j1))/(1+nu*(2e(j1)));
end
unew(J+1) = 0; % Condicion de frontera homogenea en x = 1
% Solucion por sustitucion regresiva del sistema normalizado
for j = J:1:2
unew(j) = f(j)+e(j)*unew(j+1);
end
unew(1) = 0;
% CONDICION DE FRONTERA HOMOGENEA en x = 0
u = unew; % Actualizacion para avanzar en el tiempo
plot(x,u)
axis([0 1 0 1])
M(n) = getframe;
end

43

Diferencias Finitas y E.D.P.

function teta1t(J, dt, N, theta)
% METODO TETA PARA EL MODELO PARABOLICO UNIDIMENSIONAL
% CON CONDICIONES DE FRONTERA DIRICHLET HOMOGENENAS
% SE UTILIZA EL AGORITMO DE THOMAS
%
%
%
%
%
%
%

ARGUMENTOS DE
J = numero de
dt = paso del
N = numero de
theta = valor
ARGUMENTOS DE
Nunguno

ENTRADA
subdivisiones de [0,1] (J debe ser par)
tiempo
pasos del tiempo
de teta entre 0 y 1 (CrankNicolson si theta = 0.5)
SALIDA

dx = 1/J;
x = 0:dx:1;
% Condiciones iniciales u = 2x
para x entre 0 y 0.5
%
u = 22x para x entre 0.5 y 1
j = 1:J/2+1; u(j) = 2*x(j);
j = J/2+2:J+1; u(j) = 2 2*x(j);
plot(x,u,r)
pause
nu = dt/dx^2;
nut = theta*nu;
% Coeficientes e(j) del sistema normalizado
e(1) = 0;
for j = 2:J
e(j) = nut/(1+nut*(2e(j1)));
end
% Iteraciones en el tiempo
unew = zeros(J+1,1);
for n = 1:N
% Lados derechos del sistema normalizado
f(1) = 0;
for j = 2:J
d = u(j) + (1theta)*nu*(u(j1)2*u(j)+u(j+1));
f(j) = (d + nut*f(j1))/(1+nut*(2e(j1)));
end
unew(J+1) = 0; % Condicion de frontera homogenea en x = 1
% Solucion por sustitucion regresiva del sistema normalizado
for j = J:1:2
unew(j) = f(j)+e(j)*unew(j+1);
end
unew(1) = 0;
% CONDICION DE FRONTERA HOMOGENEA en x = 0
u = unew; % Actualizacion para avanzar en el tiempo
plot(x,u)
axis([0 1 0.2 1])
M(n) = getframe;
%pause
end
hold off

44

Diferencias Finitas y E.D.P.

function adicn(I,J,dt,N)
% Metodo implicito de direcciones alternantes en 2D
% Condiciones de frontera Dirichlet homogeneas
dx = 1/I;
dy = 1/J;
nux = dt/dx^2;
nuy = dt/dy^2;
% Construccion de la malla bidimensional
x = (0:dx:1);
y = (0:dy:1);
% Condiciones iniciales
% Condicion inicial suave con valores postivos y negativos.
% En este caso debera escogerse necesariamente I = J = 48
u = peaks;
% Grafica de la condicion inicial
surf(x,y,u)
xlabel(X)
ylabel(Y)
pause
% Calculo de los coeficientes e(i) en el algoritmo de Thomas.
% ex(i) son los coef. para los sist. tridiagonales en la direccion x
% ey(j) son los coef. para los sist. tridiagonales en la direccion y
nuxm = 0.5*nux;
ex(1) = 0;
for i = 2:I
ex(i) = nuxm/(1+nuxm*(2ex(i1)));
end
nuym = 0.5*nuy;
ey(1) = 0;
for j = 2:J
ey(j) = nuym/(1+nuym*(2ey(j1)));
end
% Arreglo auxiliar para almacenar calculos intermedios
unew = zeros(I+1,J+1);
% ITERACIONES EN EL TIEMPO
for n = 1:N
% Solucion de los SISTEMAS TRIDIAGONALES EN LA DIRECCION X
for j = 2:J % Se resuelve un sistema tridiagonal para cada j
fx(1) = 0;
for i = 2:I % Algoritmo de Thomas
d = (1nuy)*u(i,j) + 0.5*nuy*(u(i,j1) + u(i,j+1));
fx(i) = (d + nuxm*fx(i1))/(1+nuxm*(2ex(i1)));
end
unew(I+1,j) = 0; % Condicion de frontera homogenea en x = 1
for i = I:1:2
unew(i,j) = fx(i)+ex(i)*unew(i+1,j);
end
unew(1,j) = 0;
% Condicion de frontera homogenea en x = 0
end
u = unew; % Actualizacion para ir a la otra direccion
% Solucion de los SISTEMAS TRIDIAGONALES EN LA DIRECCION Y

45

Diferencias Finitas y E.D.P.

for i = 2:I % Se resuelve un sistema tridiagonal para cada i
fy(1) = 0;
for j = 2:J % Algoritmo de Thomas
d = (1nux)*u(i,j) + 0.5*nux*( u(i1,j)+u(i+1,j) );
fy(j) = (d + nuym*fy(j1))/(1+nuym*(2ey(j1)));
end
unew(i,J+1) = 0; % Condicion de frontera homogenea en y = 1
for j = J:1:2
unew(i,j) = fy(j)+ey(j)*unew(i,j+1);
end
unew(i,1) = 0;
% Condicion de frontera homogenea en y = 0
end
u = unew; % Actualizacion para avanzar en el tiempo
% Grafica de la solucion
% contour3(x,y,u,50)
surf(x,y,u);
%shading flat
axis([0 1 0 1 5 5]);
%axis([0 1 0 1 0 1]);
xlabel(X)
ylabel(Y)
M(n) = getframe;
end

46