Luis Lara MDF
Luis Lara MDF
Luis Lara MDF
Per
u, 2008
Indice general
1. Introducci
on
1.4.
Metodos de Discretizacion . . . . . . . . . . . . . . . . . . . . . . . 10
1.4.1.
Serie de Taylor . . . . . . . . . . . . . . . . . . . . . . . . . 11
1.4.2.
Metodo Variacional . . . . . . . . . . . . . . . . . . . . . . 11
1.4.3.
Pesos Residuales . . . . . . . . . . . . . . . . . . . . . . . . 11
1.4.4.
Volumen de Control . . . . . . . . . . . . . . . . . . . . . . 12
1.4.5.
2. M
etodo de Diferencia Finitas
. . . . . . . . . . 56
ii
123
. . . . . . . . . . . . . . . . . . . . . . . . 125
. . . . . . . . . . . . . . . . . . . . . . . 153
iii
. . . . . . . . . . . . . . . . . . . . . . 155
181
iv
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 246
5. La Ecuaci
on del Transporte
252
281
296
Bibliografa
320
vi
Indice de figuras
2.1. Discretizacion del dominio . . . . . . . . . . . . . . . . . . . . . . .
. . . . . . . . . . . . . . . . . . . . . . . . . 128
. . . . . . . . . . . . . . . . . . . . . . . . . 130
. . . . . . . . . . . 135
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 149
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 150
3.10. Lax-Friedrichs
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 150
vii
. . . . . . . . . . . . . . . . . . . . . 155
. . . . . . . . . . . . . . . . . . . . . 179
viii
. . . . . . . . . . . . . . . . . . . . . . . . . . . . 181
. . . . . . . . . . . . . . . . . . . . . . . . . 182
. . . . . . . . . . . . . . . . . . . 182
. . . . . . . . . . . . . . . . . . . . . . . . 183
. . . . . . . . . . . . . . . . . . . . . . . . . . . . 184
. . . . . . . . . . . . . . . . . . . . . . . . . . . 184
. . . . . . . . . . . . . . . . . . . . . . . . . 212
. . . . . . . . . . . . . . . . . . . . . . . 214
. . . . . . . . . . . . . . 215
. . . . . . . . . . . . . . . . . . . . . . . . . . . . 250
. . . . . . . . . . . . . . . . . . . . . . . . 262
. . . . . . . . . . . . . . . . . . . . . 262
. . . . . . . . . . . . . . . . . . . . . . . . 262
. . . . . . . . . . . . . . . . . . . . 262
. . . . . . . . . . . . . . . . . . . . . . . . 272
. . . . . . . . . . . . . . . . . . . . . . . . 319
. . . . . . . . . . . . . . . . . . . . . . . . . . 319
Captulo 1
Introducci
on
En los u
ltimos a
nos, se ha notado un creciente interes de matematicos y especialistas en computacion cientfica, por tratar problemas que involucran modelos
con parametros de un tama
no arbitrario. Muchos de estos modelos interpretan
problemas que ocurren en la realidad como en flujos de fluidos, flujo en medios porosos, transferencia de calor y masa, electromagnetismo, contaminaci`on del medio
ambiente. Las ecuaciones diferenciales parciales (EDP), [6][7][13], que provienen
de este tipo de problemas se caracterizan por tener mas variables dependientes y
mas de una variables independientes, cuando una de las variables independientes
es la variable temporal las ecuaciones generalmente tienen una parte convectiva y
una parte difusiva, as como tambien interpretan los estados transitorios y estados
estacionarios, esto a su vez da origen a trabajar con ecuaciones de tipo hiperbolico,
ecuaciones de tipo parabolico o una convinacion de ambas ecuaciones. La conveccion y difusion ocurren simultaneamente en muchas situaciones de la fsica.
Por la variedad de tipo de ecuacion, no todo metodo numerico de resolucion
(1.1)
donde t representa el tiempo, x la variable espacial, a y son la velocidad convectiva y el coeficiente de difusividad respectivamente, y u es alguna cantidad escalar (temperatura por ejemplo). Decidimos trabajar con este modelo puesto que
a traves de ella, se puede interpretar una amplia y variada gama de problemas, y
puede ser vista como combinacion de la ecuacion de conveccion
u
u
+a
=0
t
x
(1.2)
(1.3)
La ecuacion del transporte por conveccion-difusion reune dos procesos diferentes: el difusivo, descrito matematicamente por la ecuacion parabolica (1.3); y el
convectivo descrito por la ecuacion hiperbolica (1.2). Teniendo en cuenta estos dos
procesos, la expresion que los rige de forma conjunta sera la ecuacion diferencial
(1.1) llamada ecuacion del transporte, la cual se puede ubicar en el area de la
Dinamica de Fluidos Computacional (DFC).
Por la facilidad en su aplicabilidad el Metodo de Diferencias Finitas ha sido
usando en muchos problemas de la DFC, obteniendo gran exito en algunos casos,
pero tambien numerosas dificultades en otros casos [9]. Se requiere de un analisis de
estabilidad a fin de determinar si un esquema de diferencias finitas (EDF), simple
o complicado resulta u
til pero este proceso a su vez llega a ser lento y voluminoso
que intentarlo por la experimentacion numerica llega a ser mas realista y facil de
implementar.
Como sabemos las EDP estan presentes en casi todas las ramas de la ciencia, asi
que cada tipo de EDP se adecuara en mayor o menor grado a un Esquema de Diferencias Finitas diferente que es una estrategia especial o particular de aproximacion
numerica. Pero para que un esquema de diferencias sea aceptable tenemos que garantizar que la solucion aproximada sea bastante cercana a la solucion exacta, esto
es convergencia. Pero como sabemos para llegar a establecer la convergencia de las
soluciones aproximadas a la solucion exacta es demasiado difcil o materialmente
imposible, por lo cual son necesarios dos conceptos importantes, consistencia y
estabilidad, [9],[16],[22],[23], que nos ayudan a probar la convergencia de manera
indirecta utilizando el Teorema de Lax [22], [26].
n
umero de Courant Cr a traves del n
umero de Peclet desarrollado para el caso de
la conveccion. Otro ejemplo que podemos citar es el esquema Upwind que es cont
dicionalmente estable con Cr = a x
y en combinacion con el esquema ATCE para
se consideran definidas en un n
umero finito de puntos o nodos dentro del dominio
a considerar, cuyas dimensiones estan dadas por el n
umero de variables independientes. Los nodos forman una red computacional y estan identificados mediante
un conjunto de ndices, por ejemplo en una dimension se puede denotar por (j, n),
relacionados con la coordenada espacial x y temporal t de los nodos, por medio de
los intervalos de discretizacion en cada dimension, por ejemplo (x, t). El conjunto de nodos los cuales son necesarios para escribir la ecuacion de diferencias y
los coeficientes en las ecuaciones diferenciales del modelo, relacionados a un nodo,
forman la llamada molecula computacional de un esquema de diferencias finitas.
Los terminos implcitos y explcito estan relacionados a la dinamica del esquema de diferencias finitas y depende de la forma en que se calcule el avance en el
tiempo. Los esquemas explcitos producen soluciones explcitas para las variables
dependientes en cada nivel de tiempo, explcitamente en funcion de valores conocidos en el nivel de tiempo anterior. Por otra parte, en los esquemas implcitos el
calculo de la solucion en un nodo en un nivel, involucra valores de otros nodos en
el mismo nivel conduciendo a la solucion de un sistema de ecuaciones simultaneas
en cada nivel de tiempo.
Al aplicar un esquema seleccionado, suele obtenerse un sistema de ecuaciones
algebraicas que puede ser lineales o no lineales, donde los coeficientes dependen
de las incognitas. En este u
ltimo caso su solucion puede obtenerse usando un
metodo iterativo como el de Newton-Raphson, en el cual el sistema de ecuaciones
se transforma en un sistema lineal para cada iteracion. El sistema suele tener una
matriz de coeficientes con una estructura dispersa y puede resolverse haciendo uso
de algoritmos especiales.
Para ubicar el Metodo de Diferencias Finitas dentro de los Metodos numericos
presentamos una descripcion general epistemologica y logica de como abordar un
problema con un metodo numerico.
1.1.
La Naturaleza de los M
etodos Num
ericos
La solucion numerica de una ecuacion diferencial consiste de un conjunto de
n
umeros a partir de los cuales la distribucion de la variable dependiente puede ser
construida. En este sentido, un metodo numerico es semejante a un experimento
de laboratorio, en el cual el conjunto de instrumentos leidos nos permiten establecer la distribucion de la cantidad medible en el dominio de investigacion. Tanto
el analista numerico como el experimentador de laboratorio deben quedarse con
u
nicamente un n
umero finito de valores numericos como resultado. Por tanto un
metodo numerico es un proceso de discretizacion de un problema a fin de lograr
los objetivos planteados.
Finalmente podemos decir que un metodo numerico trata como sus incognitas
basicas los valores de las variables dependientes en un n
umero finito de localizaciones llamados puntos o nodos los cuales se distribuyen formando una malla en el
dominio de calculo. El metodo incluye la tarea de proveer un conjunto de ecuaciones algebraicas para estas incognitas y prescribir un algoritmo para resolver estas
ecuaciones.
1.2.
El Concepto de Discretizaci
on
Teniendo en cuenta que un problema de ecuaciones diferenciales parciales
1.3.
Precisi
on de las Aproximaciones
Jhon von Neumann y H.H. Goldstein han identificado cuatro principales
fuentes de errores que son casi inevitables cuando estamos describiendo sistemas
fsicos.
1. Errores de modelaci
on: Los modelos matematicos de por s son generalmente dise
nados con varias simplificaciones.
2. Errores de medici
on: La mayora de las descripciones de fenomenos naturales requieren de datos medidos, estos datos de por s tienen un error
1.4.
M
etodos de Discretizaci
on
Para una ecuacion diferencial dada, los metodos de discretizacion requeridos
1.4.1.
Serie de Taylor
1.4.2.
M
etodo Variacional
1.4.3.
Pesos Residuales
se hacerca a cero. Por tanto si se elige una sucesion de funciones peso, podremos
generar todas las ecuaciones que se requieren para encontrar los parametros.
1.4.4.
Volumen de Control
1.4.5.
Este metodo, es el mas moderno, como es natural, esta basado en una teora
mucho mas compleja, como es la Teora del Potencial, que consiste en formular las
soluciones de ecuaciones dieferenciales como operadores de integrales de contorno,
muchas veces singulares, de modo que la discretizacion ya no es en todo el dominio,
sino u
nicamente en el contorno y la ecuacion se discretiza siguiendo la tecnica de
Galerkin y de los elementos finitos, pero con una funcion de prueba, muy especial,
como son el caso de las funciones de Green de los operadores diferenciales, que
formalmente, no son otra cosa mas que los operadores inversos de los operadores
diferenciales y son expresados como integrales sobre el contorno del dominio.
En esta primera parte de los metodos numericos para las ecuaciones diferenciales parciales esta dedicado al estudio de los esquemas de diferencias finitas. Nuestra
intencion es hacer nuestro ingreso al campo de la Dinamica de Fluidos Computacional(DFC) de una manera accesible, desde luego que no dejamos de mencionar
otros metodos que son importantes y que seran tratados en su debido tiempo, tales
como los elementos finitos, metodos espectrales elementos de contorno y vol
umenes
finitos. Utilizamos las Ecuciones Diferenciales Parciales(EDP), ya que esta clase de
ecuaciones son las que generalmente modelan los fenomenos naturales, mecanicos
y otros tipos de procesos que ocurren en las ciencias e Ingeniera.
Suponemos que el lector tiene un conocimiento basico de las ecuaciones diferenciales parciales asi como algunos conocimientos del analisis, si se desea probar
algunos de nuestros teoremas presentados. A fin de ser mas eficientes en nuestro
aprendizaje es necesario citar algunos libros de consulta como, Weinberger [27],
Farlow[8] y Tijonov [24]. Despues de esta rapida presentacion comezamos a presentar los conceptos basicos de los esquemas de diferencias finitas (EDF). Tambien
presentamos los conceptos importantes de convergencia, consistencia y estabilidad
los cuales son relacionados por el Teorema de Lax.
Captulo 2
M
etodo de Diferencia Finitas
El Metodo de Diferencias Finitas (MDF) es un metodo de caracter general que
permite la resolucion aproximada de ecuaciones diferenciales definidas en dominios
finitos. Probablemente es el primer metodo numerico utilizado en la resolucion de
problemas en dinamica de fluidos y tranferencia de calor, as como en problemas electromagneticos; existe documentacion en la que se menciona que Gauss
utilizo este metodo. Su uso se generalizo con la aparicion de los primeros ordenadores, y la bibligrafa sobre el mismo es abundante en los a
nos 60, especialmente en
relacion con el analisis de guas de ondas. Este metodo obtiene una solucion aproximada de las ecuaciones diferenciales definidas en un domnio o region de trabajo.
Sobre dicho dominio estaran definidas las condiciones de contorno o frontera y las
condiciones iniciales que marcaran el punto de partida en la solucion de problemas
concretos.
El primer paso para la aplicacion del metodo es definir el dominio donde ha de
calcularse el valor de la funcion incognita a resolverse. Dicho dominio, que en este
2.1.
M
etodo de Diferencias Finitas
Unidimensional
En esta seccion presentamos en forma general, y basica el concepto de
2.1.1.
Discretizaci
on del Dominio
}|
xn1
xn
xn+1
2.1.2.
Discretizaci
on de la Variable
Definici
on 2.1.4 Una funcion discreta v, es aquella que esta definida en una
malla , tal que a cada punto xn se le asocia un n
umero vn .
Un ejemplo facil se obtiene considerando una funcion continua u definida en
, a la cual se le puede asociar una funcion discreta v definida sobre una malla
de tal que vj = u(xj ).
A las funciones discretas tambien se les puede asociar una norma y por tanto
ubicarlas en un espacio adecuado.
Como R puede ser un intervalo finito, infinito o R. La funcion discreta puede contener un n
umero finito o infinito de valores de la forma v =
(v1 , v2 , . . . , vN ), v = (v1 , v2 , . . .), v = (vN , vN +1 , . . . , v1 , v0 , v1 , . . . , vN , vN +1 )
o v = (. . . , v1 , v0 , v1 , . . .). En cualquier caso podemos considerar a una funcion
discreta definida para j Z, ademas usamos la notacion
h = x = max {xj }
j
Espacios y Normas
Las diferentes normas finito dimensional que usaremos incluyen la norma euclidiana dependiendo de la malla
v
u N
uX
kvk2 = t
|vj |2
(2.1)
j=1
(2.2)
(2.3)
1jN
(2.4)
1jN
l2 = {v = (. . . , v1 , v0 , v1 , . . .)
|vj |2 < }
(2.5)
j=
con norma
v
uX
u
kvk2 = t
|vj |2
(2.6)
j=
(2.7)
j=
sup
<j<+
|uj | < }
(2.8)
con norma
kvk =
sup
<j<+
|vj |
(2.9)
y norma discreta
kvk,x =
sup
<j<+
|vj | x
(2.10)
Y finalmente, cuando tomamos transformaciones, encontramos el espacio infinito dimensional lineal sobre los complejos, las funciones cuadrado integrables de
Lebesgue, L2 (IR),
Z
|v(x)|2 dx < }
L2 (IR) = {v : IR IC :
(2.11)
IR
con norma
Z
kvk22
2.1.3.
|v(x)|2 dx
(2.12)
IR
Discretizaci
on de los Operadores Diferenciales
y x tal
que
f (x) = Pn (x) + Rn (x)
(2.13)
donde
Pn (x) = f (x0 ) + f 0 (x0 )(x x0 ) +
=
n
X
f (k) (x0 )
k=0
k!
f (n) (x0 )
f 00 (x0 )
(x x0 )2 + . . . +
(x x0 )n
2!
n!
(x x0 )k ,
f (n+1) ((x))
(x x0 )(n+1)
(n + 1)!
h2
hn
+ . . . + f (n) (x0 )
2
n!
(2.14)
y el residuo por
Rn (x0 + h) =
1 2
f (y0 )(x x0 ) + f (y0 )(t t0 ) +
f (y0 )(x x0 )2
x
t
2! x2
1 2
2
+
f (y0 )(x x0 )(t t0 ) +
f (y0 )(t t0 )2 + . . .
xt
2! t2
1
1
N
+
f (y0 )(x x0 )N m (t t0 )m + RN +1
(N m)! m! xN m tm
j
N N
X
X
1 j+n
f (y0 )(x x0 )j (t t0 )n + RN +1
(2.15)
=
j tn
j!n!
x
j=0 n=0
f (x, t) = f (y0 ) +
donde,
RN (x, t)
es el residuo multidimensional o error de truncamiento.
f (2) ()h2
,
2!
x0 < < x 0 + h
(2)
f (x0 + h) f (x0 )
f ()h
f (x0 ) =
+
h
2!
0
(2.16)
f (x0 + h) f (x0 )
.
h
(2.17)
f (2) ()h
= O(h)
2!
donde O(h), se lee, orden de h, y significa que la funcion del lado derecho se
comporta en forma menor o igual que la funcion h.
En general decimos que una funcion g(h) es de orden O(hp ), p > 0, si y solo si
existe una constante M 0 tal que
|g(h)| M hp
Tambien podemos usar la formula de Taylor para aproximar la funcion en el
punto x = x0 h
f (x0 h) = f (x0 ) f 0 (x0 )h +
de donde
f (2) (1 )h2
,
2!
x0 h < 1 < x 0
(2)
f (x0 ) f (x0 h)
f (1 )h
f (x0 ) =
+
h
2!
0
(2.18)
f (x0 ) f (x0 h)
h
(2.19)
f (2) (1 )h
= O(h)
2!
Se
nalaremos que en ambos casos se han involucrado dos puntos. En el caso de
la aproximacion (2.17), se utilizaron los puntos x0 y x0 + h. En la aproximacion
(2.19), se utilzaron x0 h y x0
Tambien podemos obtener otra aproximacion involucrando tres puntos, tales
como x0 , x0 + h, x0 h, para ello sumamos (2.16) y (2.18) y despejando f 0 (x0 )
obtenemos
f 0 (x0 ) =
f (x0 + h) f (x0 h) h 00
00
+
f (1 ) f ()
2h
2
(2.20)
f (x0 + h) f (x0 h)
2h
h 00
00
f (1 ) f ()
2
(2.21)
Esto significa que la pendiente de la secante a la grafica de f que pasa por los
puntos (x0 h, f (x0 h)) y (x0 + h, f (x0 + h)) es mucho mejor aproximacion a la
pendiente de la tangente que pasa por (x0 , f (x0 )), que la pendiente de las secantes
que pasan por (x0 h, f (x0 h)), (x0 , f (x0 )) y por (x0 , f (x0 )), (x0 + h, f (x0 + h)).
Esto quiere decir que el error de truncamiento en (2.21) debe ser el de mayor
orden, lo cual se logra suponiendo que f es tres veces diferenciable.
Aplicando el Teorema del Valor Medio en (2.21) tenemos;
i
h h 00
(x0 , h) =
f ()(1 ) ,
2
= O(h2 ),
1 < <
pues 1
es de orden O(h2 )
El calculo exacto del orden de truncamiento (x0 , h), se obtiene usando la formula
de Taylor, cuyo orden es el que indica el exponente de h, en el caso que n = 2
tenemos
(x0 , h) =
h2 (3)
f (1 )(2 1 ) , 2 < 1 < 1
3!
= O(h2 )
el cual es un error de truncamiento de orden dos.
Presentamos la aproximacion por diferencias finitas para la primera derivada. Para ello revisemos las aproximaciones presentadas anteriormente en forma
detallada. Tomemos la Formula de Taylor, para n=2.
1. Sea f una funcion de clase C (2) que cumple las condiciones del teorema de
Taylor
f (x) = f (x0 )(x x0 )0 + f 0 (x0 )(x x0 )1 +
f (2) ((x) )
(x x0 )2
2
(2.22)
x0 < x0 +h < x0 + h
(2.23)
donde
h
f (2) ((x0 +h) ) O(h)
2
es decir
h
f (2) ((x0 +h) ) O(h)
2
entonces la ecuacion (2.23) se puede escribir en la forma
f 0 (x0 ) =
f (x0 + h) f (x0 )
+ O(h)
h
f (x0 + h) f (x0 )
h
(2.24)
La aproximacion (2.24) es llamada aproximacion por diferencias finitas progresiva o hacia adelante para f 0 donde O(h) es llamado error de truncamiento de orden p = 1, este orden depende de la potencia de h y como
aparece es O(hp ) e involucra dos puntos.
2. Usando la formula (2.22) con x = x0 h , h > 0, tenemos
f (x0 h) = f (x0 ) f 0 (x0 )h +
f 0 (x0 ) =
f (2) ((x0 h) ) 2
h
2
f (x0 ) f (x0 h)
h
+ f (2) ((x0 h)
h
2!
donde
h
+ f (2) ((x0 h) O(h)
2
(2.25)
es decir
h
+ f (2) ((x0 +h) ) O(h)
2
entonces tenemos que la expresion (2.25) se puede escribir en la forma
f 0 (x0 ) =
f (x0 ) f (x0 h)
+ O(h)
h
f (x0 ) f (x0 h)
h
(2.26)
f 00 (x0 )
(x x0 )2
2!
f (3) ((x) )
+
(x x0 )3
3!
(2.27)
h2 00
h3
f (x0 ) + f (3) ((x0 +h) )
2!
3!
(2.28)
h2 00
h3
f (x0 ) f (3) ((x0 h) )
2!
3!
(2.29)
h3 (3)
f ((x0 +h) ) + f (3) ((x0 h) )
3!
donde
h2 (3)
f ((x0 h) ) O(h2 )
3!
h2 (3)
f ((x0 h) ) O(h2 )
3!
es decir
f (x0 + h) f (x0 h)
+ O(h2 )
2h
f (x0 + h) f (x0 h)
2h
(2.31)
La aproximacion de diferencias (2.31) se llama aproximacion por diferencia central para f 0 donde O(h2 ) es llamado error de truncamiento
de orden p = 2, esta aproximacion incluye tres puntos.
Usando la formula de Taylor se puede hallar otras aproximaciones, para la
primera derivada, las mismas que pueden involucrar mas de tres puntos, y ser de
mayor orden, como por ejemplo
f 0 (x0 )
(2.32)
la cual es de tercer orden. Dejamos como ejercicio al lector, verificar estos ordenes.
Ahora definamos las aproximaciones para la segunda derivada. Sea f una funcion de clase C (3) , que satisface las condiciones del teorema de Taylor para n = 3,
alrededor de un punto x0
f (x) = f (x0 )(x x0 )0 + f 0 (x0 )(x x0 )1 +
f 00 (x0 )
(x x0 )2
2!
(2.33)
h3
h4
h2 00
f (x0 ) + f 000 (x0 ) + f (4) ((x0 +h) )
2!
3!
4!
(2.34)
h2 00
h3
h4
f (x0 ) f 000 (x0 ) + f (4) ((x0 h) )
2!
3!
4!
(2.35)
h4 (4)
f ((x0 +h) ) + f (4) ((x0 h) )
4!
observando que
h2 (4)
f ((x0 h) ) O(h2 )
3!
(2.37)
h2 (4)
f ((x0 h) ) O(h2 )
3!
(2.38)
es decir
i1
i+1
(2.39)
(2.40)
2.2.
Discretizaci
on del Problema de EDP
Un problema de EDP
(P )
en Rn , t > 0
en R
u(x, 0) = u0
en t = 0
infinito
+ infinito
2.2.1.
Discretizaci
on del Problema (P)
t 6
tn+1
(j, n + 1)
w
tn
(j 1, n)
(j, n)
(j + 1, n)
xj1
xj
xj+1
Sean x = h y t = k n
umeros positivos; entonces la malla esta formada
por los puntos (xj , tn ) = (jx, nt) para enteros arbitrarios j y n. Para
una funcion discreta v definida sobre la malla, denotamos sus valores por vjn ,
en los puntos malla (xj , tn ). Tambien usaremos la notacion unj para u(xj , tn )
donde u es una funcion continua en el plano xt. El conjunto de puntos (xj , tn )
para un valor de n fijo es llamado nivel de malla n. La idea basica para
el esquema de diferencias finitas es remplazar derivadas por sus diferencias
finitas.
u(jh, (n + 1)k) u(jh, nk)
u
(jh, nk)
t
k
u(jh, (n + 1)k) u(jh, (n 1)k)
u
(jh, nk)
t
2k
Similares formulas aproximan las derivadas con respecto a x.
2. Discretizaci
on de la variable
t u
(jx, nt)
1! t
t2 2 u
(jx, nt) + . . .
2! t2
(2.41)
(2.42)
donde se asume que las derivadas de mayor orden de u en (jx, nt) son
acotadas. Entonces una aproximacion de ut (jx, nt) puede ser dada por
vjn+1 vjn
t
(2.43)
uxx (x, t) = lm
(2.44)
(2.45)
Asi las expresiones dadas en (2.43) y (2.45) pueden ser usadas para aproximar
la ecuacion diferencial
u
2u
= 2
t
x
(2.46)
n
n
vjn+1 vjn
vj+1
2vjn + vj1
=
t
x2
(2.47)
t n
n
(vj+1 2vjn + vj1
)
2
x
(2.48)
Finalmente, es facil ver que las condiciones iniciales y de frontera son aproximados por
vj0 = f (jx), j = 0, . . . , M
(2.49)
(2.50)
n+1
= b((n + 1)t), n = 0, . . .
vM
(2.51)
2.3.
S j = j1 ,
(S+ I)j
j+1 j
=
x
x
(I S )j
j j1
=
x
x
j+1 n1
(S+ S )j
=
2x
2x
02 j = D D+ j =
y se cumple
D D+ j = D+ D j
Ahora cuando la funcion discreta es de dos variables independientes (j, n) se
utiliza superndices con el nombre de la variable continua que se va ha discretizar,
por ejemplo, para aproximar la primera derivada con respecto a la variable t se
utiliza
t n
D+
vj
(S+t I)vjn
=
t
o
x n
D+
vj =
2.4.
(S+x I)vjn
x
Convergencia
El concepto de convergencia es un concepto matematico familiar, conocido
en la convergencia de sucesiones de n
umeros, sin embargo, aqu se refiere al hecho
de que son sucesiones de soluciones discretas obtenidas como soluciones de los
problemas discretos que aproximan a la solucion exacta de un problema continuo.
Las sucesiones de soluciones de los problemas discretos, que son soluciones aproximadas se obtienen cuando x, t tiende a cero, [22], [23]. Por tanto es suficiente
utilizar los puntos de la malla, por lo que sera necesario hablar de convergencia
puntual y/o convergencia uniforme.
Una importante pregunta concerniente a las soluciones discretas obtenidas computacionalmente es, Que garanta puede ser dada para que la solucion computacional sea proxima a la solucion exacta de la ecuacion diferencial parcial y bajo
que circunstancias la solucion computacional coincide con la solucion exacta ?. La
segunda parte de esta pregunta puede ser respondida requiriendo que la solucion
aproximada (computacional) converja a la solucion exacta cuando x, t tiende a
cero. Sin embargo, convergencia es muy difcil de establecer directamente por tanto
una ruta indirecta es utilizar dos conceptos muy relacionados, y que juntos implican convergencia mediante el Teorema de Lax; estos conceptos son consistencia y
estabilidad.
Ecuacin
Diferencial
Discretizacin
Sistema
de
Consistencia
Ecuaciones
Mtodos
Analticos
Estabilidad
.
Solucin
Exacta
Solucin
Convergencia
Aproximada
Esta ruta indirecta requiere que el operador de diferencias discreto debe ser
consistente con el operador diferencial parcial, esto implica la inversa del proceso
de discretizacion.
A traves del desarrollo de la serie de Taylor se recupera las ecuaciones diferenciales gobernantes cuando x, t tienen a cero. Ademas, el algoritmo usado para
resolver el sistema de ecuaciones algebraicas para obtener la solucion aproximada
debe ser estable.
El Teorema de Lax es de gran importancia, puesto que es relativamente facil
demostrar la estabilidad de un algoritmo y la consistencia de la discretizacion con
la ecuacion diferencial parcial original, mientras que por lo general es muy difcil
de mostrar la convergencia de las soluciones.
En realidad la convergencia es necesaria medirla en alguna norma, veamos la
siguiente.
Definici
on 2.4.1 La sucesi
on de soluciones discretas {vjn } converge en la norma
k.kx si y solo si
lm (
max
x0 j=0,1,...,T /t
kej kx ) = 0
(2.52)
(2.53)
2
=
, 0 < x < 1, t > 0,
t
x2
(2.54)
(x, 0) = f (x), 0 x 1,
(2.55)
(2.56)
(2.57)
(2.58)
(2.59)
(2.60)
De ello se obtiene
en+1
= renj1 + (1 2r)enj + renj+1 + n+1
nj
j
j
+ r(2nj nj1 nj+1 )
(2.61)
h2 2
)j,n + ( 2 )(xj +1 h,tn )
x
2! x
(2.62)
h2 2
h( )j,n + ( 2 )(xj 2 h,tn )
x
2! x
(2.63)
nj+1 = (xj + h, tn ) = nj + h(
nj1
= (xj h, tn ) =
nj
y
n+1
= (xj , tn + k) = nj + k(
j
)(x ,t + k)
t j n 3
(2.64)
(2.65)
2
)(xj ,tn +3 k) ( 2 )(xj +4 h,tn ) |
t
x
h
k2
(2.66)
se tiene que todos los coeficientes de e son positivos o cero, por consiguiente se
obtiene que
|en+1
| r|enj1 | + (1 2r)|enj | + r|enj+1 | + kMh,k
j
rEn + (1 2r)En + rEn + kMh,k
= En + kMh,k ,
(2.67)
(2.68)
(2.69)
para todo n.
Luego recursivamente tenemos la cadena de desigualdades
En+1 En + kMh,k (En1 + kMh,k ) + kMh,k = En1 + 2kMh,k
(2.70)
y por consiguiente
En E0 + nkMh,k = tn Mh,k
(2.71)
donde tn = nk y puesto que los valores iniciales de v y son las mismas se cumple
E0 = 0.
)(x,t) |
t
x2
(2.72)
Ahora puesto que es solucion de (2.54), se tiene que el valor limite de Mh,k es
cero. Por lo tanto, cuando h, k 0;
|nj nj | = |enj | max |enj | = En
(j,n)
es decir
|nj nj | En
(2.73)
para todo n. Por tanto se cumple que vjn (x, t) cuando (jh, nk) (x, t).
Se ha demostrado que v converge en la norma del maximo a cuando h 0,
r 1/2 y t es finito. Cuando r > 1/2 puede demostrarse que el esquema no genera
una sucesion de soluciones convergente.
Para esquemas de multipasos la definicion asume que alg
un procedimiento inicial es usado para computar los primeros niveles de tiempo, necesarios para emplear los esquemas de multipasos. Para el caso que los datos son especificados
sobre estos primeros niveles de tiempo, las definiciones se modifican para requerir
vjm , 0 m J converga a v0 (xj ).
(2.74)
2.4.1.
Un esquema de diferencias finitas tal como aquellos que discutiremos son usados
porque sus soluciones aproximan a las soluciones de ciertas ecuaciones diferenciales
parciales. Pero que realmente es necesario para que la solucion de las ecuaciones
en diferencias puedan aproximar a la solucion de la ecuacion diferencial parcial
con cualquier precicion deseada. Asi es necesatio alguna clase de convergencia de
la solucion de la ecuacion en diferencias a la solucion de la ecuacion diferencial
parcial.
Empezamos considerando problemas de valor inicial. Asi consideremos la ecuacion diferencial parcial
Lu = F
(2.75)
donde las funciones u y F estan definidas en toda la recta real y ademas se tiene
la condicion inicial u(x, 0) = f (x), x R. Asumimos que hemos obtenido una solucion aproximada de (2.75) por el esquema de diferencias finitas que denotaremos
por Lnj donde, n corresponde al paso del tiempo y j al punto espacial de la malla.
Esta solucion aproximada la denotaremos por vjn la cual esta definida sobre toda
la malla
= {(jx, nt), j = , . . . , +, n = 0, . . .}.
y por u la solucion exacta de nuestro problema de valor inicial.
Empezamos dando la siguiente definicion
Definici
on 2.4.2 Un esquema de diferencias finitas Lnj vjn = Gnj que aproxima a
la ecuaci
on diferencial parcial Lu = F es un esquema puntualmente convergente
(2.76)
(2.77)
(2.78)
(2.79)
Soluci
on:
Debemos entender que estamos considerando un problema de valor inicial sobre
todo R, el ndice j sobre vjn sera el rango sobre todo los enteros, < j < +.
Denotamos la solucion exacta del problema de valor inicial (2.78)-(2.79) por u
tal que
enj = vjn u(jh, nk)
(2.80)
(2.81)
Entonces sustrayendo la ecuacion (2.81) de la ecuacion (2.76), vemos que enj satisface
en+1
= (1 2r)enj + r(enj+1 + enj1 ) + O(k 2 ) + O(kh2 )
j
(2.82)
Si 0 < r < 1/2, los coeficientes del lado derecho de la ecuacion (2.82) son
no-negativos y
|enj | (1 2r)|enj | + r|enj+1 | + r|enj1 | + A(k 2 + kh2 )
(2.83)
donde A es una constante asociada con el termino big O que depende de las
hipotesis asumidas sobre las derivadas de orden superior. Por lo tanto si tomamos
la norma del supremo en el nivel n
En =
sup
<j<+
{|enj |}
(2.84)
llegamos a
E n+1 E n + A(k 2 + kh2 )
(2.85)
Hay que notar que en el supremo sobre el lado derecho se incluyen los terminos
que contienen k y h. En este caso asumimos que la constante A es una cota de
la segunda derivada con respecto al tiempo y la cuarta derivada con respecto al
espacio en toda la recta real. Asi hemos asumido que estas derivadas de la solucion
son uniformente acotadas en toda la recta.
Aplicamos repetidamente (2.85) y llegamos a
E n+1 E n + A(k 2 + kh2 ) . . . E 0 + (n + 1)A(k 2 + kh2 )
(2.86)
(2.87)
v
u(jh, (n + 1)k) (n + 1)kA(k 2 + kh2 ) 0
j
(2.88)
cuando k, h 0.
Asi vemos que para cualquier x y t, cuando k y h se aproximan a 0 en semejante
forma que (jh, (n + 1)k) (x, t), vjn se aproxima a u(x, t).
Deberiamos notar que puesto que (n + 1) en la ecuacion (2.88) puede causar
problemas en nuestra convergencia ( no es bueno que tengamos terminos que van
al infinito en una expresion que deseamos que tienda a cero), esto fue un hecho
que en el u
ltimo paso se tenga (n + 1)k t.
Ahora usamos la norma del m
aximo sobre el espacio de sucesiones acotadas, l ,
kuk =
max
<j<+
|uj |
(2.89)
n
Si un = (. . . , un1 , un0 , un1 , . . .)T y vn = (. . . , v1
, v0n , v1n , . . .)T , entonces tenemos que
probar que para t tal que (n + 1)k se aproxima a t, v n+1 se aproxima a u(., t) en
la norma del maximo (2.89).
Puesto, que en general, la convergencia puntual es difcil de probar, cambiaremos la definicion de convergencia en terminos de la norma del maximo.
Por conveniencia, denotamos el vector de valores de la solucion de la ecuacion
de diferencias vjn en el nivel n por vn y el vector de valores en nk de la solucion
de la ecuacion diferencial parcial evaluada en los puntos malla, u(jh, nk), por un ,
generalizamos la definicion de convergencia en la siguiente forma.
Definici
on 2.4.3 Un esquema de diferencias finitas Lnj vjn = Gnj que aproxima a
la ecuaci
on diferencial parcial Lu = F es un esquema convergente en la norma k.k
(2.90)
cuando h 0 y k 0
Dede especificarse que en la Definicion (2.93) no fue especificada la norma
porque en diferentes situaciones, diferentes normas seran utilizadas.
Se probo que un esquema explcito (2.76) fue convergente de acuerdo a la
Definicion (2.4.2) en la norma del supremo. En otras ocaciones veremos que la
norma natural elegida sera una variacion de la norma L2 y la norma L2,h . Asi,
siempre que trabajemos con normas estas deben ser especificadas.
Debe quedar claro que la Definicion (2.4.3) difiere de la Definicion (2.4.2).
Usando la Definicion (2.4.2), la razon de convergencia de vjn a u(x, t) (razon en
terminos de hx 0 y kt 0) puede variar considerablemente para diferentes
valores de x. Puede ocurrir, seg
un la Definicion (2.4.2), que un esquema converge
en algunos valores de (x, t) y no en otros. Sin embargo, via la Definicion (2.4.3)
(es decir k un vn k es peque
no), entonces se obtiene que vjn esta proxima a unj
para todo j.
A veces deseamos discutir convergencia en terminos la rapidez en la cual la
solucion del esquema de diferencias converge a la solucion de la ecuacion diferecial
parcial. Para este proposito utilizamos una definicion de convergencia de orden
(p, q) como sigue
Definici
on 2.4.4 Un esquema de diferencias finitas Lnj vjn = Gnj que aproxima a
la ecuaci
on diferencial parcial Lu = F es un esquema convergente de orden (p, q)
(2.91)
cuando h 0 y k 0
Cuando usamos la notacion big O, debemos recordar siempre que existe una
constante involucrada, es decir, la ecuacion (2.91) es realmente una notacion corta
para existe una constante C tal que
kun+1 vn+1 k C(O(hp ) + O(k q ))
(2.92)
2.4.2.
no o M es
y encontramos que vjn , j = 1, . . . , M 1, n = 1, . . .. Cuando h es peque
grande, el vector desconocido sera grande pero finito. Esto hace que no podamos
(2.93)
cuando i y t 0
En una forma similar, se puede definir convergencia de orden (p, q).
Aunque para satisfacer la Definicion (2.4.5) debemos considerar todas las particiones que converjan a cero, el modelo para usar la Definicion (2.4.5) es usar el
dominio [0, 1] y h = 1/M . Entonces el espacio obtenido sera finito dimensional,
(2.94)
n+1
v0n+1 = vM
=0 n0
(2.95)
vj0 = f (jx), j = 0, . . . , M
(2.96)
(2.97)
(2.98)
(2.99)
Soluci
on: Probamos convergencia del esquema (2.94)-(2.96) en la misma forma
que se prueba para un problema de valor inicial, excepto en esta caso, nuestro
vector tiene longitud finita y el tama
no ira modificandose. Tambien, porque hacemos esto en forma similar, eventualmente tenemos que hacer varias suposiones
para hacer el trabajo de prueba. Empezamos escribiendo xm como el incremento
sup
1jMm 1
|vj |
Haciendo
enj = vjn u(jxm , nt)
Como vimos anteriormente enj satisface
en+1
= (1 2r)enj + r(enj1 + enj+1 )
j
+O(t2 ) + O(t(xm )2 ), j = 1, . . . , Mm 1
(2.100)
Si 0 < r 1/2, los coeficientes del lado derecho de la ecuacion (2.100) son no
negativos, donde en0 y enMm son ambos ceros, y
|en+1
| (1 2r)|enj | + r|enj+1 | + r|enj1 | + A(O(t2 ) + O(t(xm )2 ))
j
ken kMm 1, + A(O(t2 ) + O(t(xm )2 ))
donde
en = (en1 , . . . , eMm 1 )T
Si tomamos el supremo sobre el lado izquierdo de la desigualdad anterior y aplicamos repetidamente esta desigualdad, entonces tenemos
ken+1 kMm 1, (n + 1)tA(t + (xm )2 )
Por lo tanto, cuando t 0, (n + 1)t t, y m ,
kun+1 vn+1 kMm 1, 0,
puesto que xm 0, se tiene que el esquema es convergente.
(2.101)
2.5.
Consistencia
2.5.1.
Definici
on 2.5.1 (Consistencia Puntual) El esquema de diferencias finitas
Lnj vjn = Gnj es puntualmente consistente con la ecuaci
on diferencial parcial
Lu = F en el punto (x, t) si para cualquier funcion suave = (x, t), se tiene que
(Lu F )|nj [Lnj (jh, nk) Gnj ] 0
(2.102)
(2.103)
donde
n
vn = (. . . , v1
, v0n , v1n , . . .)T ,
(2.104)
(2.105)
cuando h, k 0.
Una simple variacion de la Definicion (2.5.2) es dado cuando se quiere precisar
de consistencia.
Definici
on 2.5.3 El esquema de diferencias finitas (2.103) se dice que es consistente de orden de precisi
on (p, q) con la ecuaci
on diferencial parcial Lu = F
si
k n k = O(hp ) + O(k q )
(2.106)
(2.107)
con la ecuaci
on diferencial parcial
ut = uxx , < x < , t > 0
(2.108)
Soluci
on: Si denotamos la solucion de la ecuacion diferencial parcial (2.108) por
u, e insertando u en la expresion (2.107) se obtiene
(un+1
unj ) r(unj+1 2unj + unj1 ) = O(k 2 ) + O(kh2 )
j
(2.109)
donde r = k/h2 . Para aplicar las Definiciones (2.5.2)-(2.5.3) debemos ser cuidadosos de lo que esta contenido en O(k) + O(h2 ). Si asumimos que la segunda
(2.110)
+
+
=
x2
t
(uxxxx (x1 , nt) + uxxxx (x2 , t))
2
24
(2.111)
k
h2
+ sup |uxxxx |
Ak + Bh2
2
12
i
(2.112)
j=
(2.113)
j=
2.5.2.
La consistencia puntual sera la misma que aparece en la Definicion (2.5.1) excepto que ahora tenemos que analizar cualquier condicion de frontera que contiene
una aproximacion.
Para la norma considermos una sucesion de particiones del intervalo [0, 1]
definiendo una sucesion de incrementos espaciales {xj } y una sucesion de espacios
apropiados, {Xj }, con normas k.kj Entonces las Definiciones (2.5.2) y (2.5.3) se
vj0 = f (jx)
(2.114)
(2.115)
n+1
vM
= 0
(2.116)
(2.117)
j = 1, ...M 1
(2.118)
(2.119)
u(x, 0) = f (x),
(2.120)
u(1, t) = 0, t > 0
(2.121)
ux (0, t) = 0, t > 0
(2.122)
x [0, 1]
(2.123)
Soluci
on: La condicion de frontera (2.122) se puede aproximar por
n
v1n v1
=0
2x
(2.124)
(2.125)
con componentes
vjn+1 = [(1 2r)I + r(+ + )] vjn
Si u es la solucion exacta del problema de valor inicial-frontera continuo, entonces
podemos reemplazarlo en (2.125) y obtenemos
un+1 = Qun + t n
(2.126)
Antes de empezar nuestro analisis de consistencia, debemos decidir sobre nuestra sucesion de espacios. Es claro que los espacios que usaremos seran espacios de
dimension M que constan de vectores (v0 , . . . , uM 1 )T . Ya conocemos que jn es
O(t) + O(x)2 , j = 1, . . . , M 1.
Para examinar la consistencia de la ecuacion primero examinemos la consistencia de la discretizacion en la frontera (2.117), sea u la solucion del problema de
t2
+ ...
2
[(1 2r)un0
x2
x3
+ (uxxx )n0
+ . . .}]
2
6
t2
= [(ut )n0 (uxx )n0 ]t 2rx(ux )n0 + (utt )n0
2
tx
(uxxx )n0
+ ...
3
+ 2r{un0 + (ux )n0 x + (uxx )n0
Usando el hecho que (ux )n0 = 0 y (ut uxx )n0 = 0, vemos que
0n =
(2.127)
Ejemplo 6
Soluci
on: En efecto, el esquema es
n
vjn+1 vjn
vjn
vj+1
+a
= 0,
k
h
+a
t
x
y sea una funcion suave = (x, t). El operador de diferencias finitas esta dado
por
Lnj
n+1
nj
nj+1 nj
j
=
+a
k
h
(2.128)
donde nj = (jh, nk). Usando las series de Taylor de la funcion (x, t) en el nodo
(jh, nk),
1 2
k tt + O(k 3 )
2!
1
= nj + hx + h2 xx + O(h3 )
2!
n+1
= nj + kt +
j
(2.129)
nj+1
(2.130)
donde las derivadas son evaluadas en los nodos (xj , tn ). Reemplazando (2.129) y
(2.130) en (2.128) obtenemos
1
1
Lnj = t + ktt + O(k 2 ) + a{x + hxx + O(h2 )}
2
2
1
a
= t + ktt + ax + hxx + O(k 2 ) + O(h2 )
2
2
Entonces
1
a
Lnj Lnj = t + ax t ktt ax xx
2
2
O(k 2 ) O(h2 )
1
a
= kO(h2 )itt hxx O(k 2 )
2
2
Por lo tanto
Lnj Lnj 0; h, k 0
y as concluimos que el esquema es puntualmente consistente.
(2.131)
con la condici
on inicial
u0 (x) =
1 , 1 x 0
0 , en otro lugar
(2.132)
(2.133)
Soluci
on: El esquema puede reescribirse en la forma
k
vjn+1 = vjn (vj+1 vjn )
h
n
= (1 + Cr )vjn Cr vj+1
t 6
u 6= 0
w
u 6= 0
v=0
u=0
u0j =
1 , 1 j h 0
0 , en otro lugar
vjn es cero.
Note que concluimos que el esquema es no convergente sin especificar el tipo de
convergencia, pero ciertamente, una sucesion de funciones que son todas ceros esto es, vjn para j > 0 - no puede converger, bajo cualquier definicion razonable de
convergencia, a la funcion diferente de cero u.
En seguida demostraremos que el esquema de Lax-Friedrichs es consistente
bajo ciertas condiciones.
Ejemplo 8
ecuci
on diferencial ut + aux = 0.
Soluci
on: En efecto, sea el esquema de Lax-Friedrichs
n
n
n
n
vjn+1 21 (vj+1
+ vj1
)
vj+1
vj1
+a
=0
k
2h
n+1
12 (nj+1 + nj1 )
nj+1 nj1
j
=
+a
k
2h
2.6.
Estabilidad
Note que si un esquema es convergente, entonces vjn converge a u(x, t), ciertamente vjn es acotada en alg
un sentido. En este contexto, estos conceptos que
estamos utilizando son aplicables a los llamados Problemas bien puestos.
La siguiente definicion de estabilidad es para problemas de valor inicial homogeneos.
Definici
on 2.6.1 (Estabilidad de PVI Homog
eneos) Un esquema de diferencias finitas Lnj vjn = 0 para una ecuaci
on diferencial de primer orden, Lu = 0
es estable si existe un entero J y n
umeros positivos h0 y k0 , tal que para cualquier
tiempo positivo T , existe una constante CT tal que
h
|vjn |2
CT h
j=
X
X
|vjn |2
(2.134)
j=0 j=
(2.135)
Definici
on 2.6.2 (Problemas bien Puestos) El problema de valor inicial para
una ecuaci
on diferencial parcial de primer orden Lu = 0 es bien puesto si para
cualquier tiempo T 0, existe una constante CT tal que para cualquier solucion
u(x, t) satisface
|u(x, t)| dx CT
|u(x, 0)|2 dx
0tT
(2.136)
2.6.1.
An
alisis de Estabilidad Matricial
2
=
, 0<x<1
t
x2
(2.137)
con
(0, t) = a, (1, t) = b, 0 t T
y
(x, 0) = g(x), 0 < x < 1
Sea el esquema de diferencias finitas de Euler definido en un punto interior general
n+1
= nj + r(nj1 2nj + nj+1 )
j
para j = 1, ..., N y N h = 1, n = 0, ..., J con Jk = T .
En efecto, sobre las fronteras se tiene
= n1 + r(a 2n1 + n2 )
n+1
1
y
n
n
n
n+1
N 1 = N 1 + r(N 2 2N 1 + b)
donde N h = 1.
(2.138)
n+1
n+1
n+1
2
..
= .
n+1
N 2
n+1
N 1
n+1
n
r
1 2r
1 ra
r
n 0
1
2r
r
.. ..
..
.. ..
=
. + .
.
.
.
0
r 1 2r
r
N 2
r
1 2r
nN 1
rb
(2.139)
(2.140)
(2.141)
(2.142)
(2.143)
El esquema de diferencias finitas sera estable cuando los restos de en seran acotados cuando n se incrementa indefinidamente. En otras palabras e0 se propaga de
acuerdo a
en = Aen1 = . . . = An e0
(2.144)
(2.145)
Lax define el esquema de diferencias como estable si existe M > 0 (independiente de n, h y k) tal que kAn k M, n = 1, . . . , T . Esta condicion ciertamente
limita la amplificacion de cualquier perturbacion inicial y por lo tanto cualquier
error de redondeo puesto que
ken k M ke0 k
Desde que
kAn k = kAAn1 k kAk kAn1 k . . . kAkn
(2.146)
Nota
La norma uno de la matriz A es la suma maxima de los modulos de las
columnas de A, y se denota por kAk1 . La norma infinita de la matriz A es
la suma maxima de los modulos de las filas de A, y se denota por kAk . La
norma dos de la matriz A es la raz cuadrada del radio espectral de AH A donde
AH = (A)T (transpuesta de la conjugada de A), y se denota por kAk2 .
(A) 1 desde que (A) kAk
(el recproco no es cierto).
(2.147)
Ejemplo 10 Existen dos formas para analizar la estabilidad por medio de matrices: encontrar los autovalores de A explcitamente o usar normas para obtener
cotas sobre los autovalores.
Cuando 0 < r 1/2 entonces
kAk = r + (1 2r) + r = 1
(2.148)
(2.149)
Asi concluimos que el esquema es estable para 0 < r 1/2 e inestable para
r > 1/2.
Otra alternativa es usar el hecho de que A es real y simetrica, entonces se tiene
kAk2 = (A) = max |j | donde j son los autovalores de A. Si A = I + rS donde
I = ...
2 1
1 2 1
.. .. ..
y S=
.
.
.
1
1 2 1
1 2
(2.150)
j
), j = 1, . . . , N 1
2N
(2.151)
j
), j = 1, . . . , N 1
2N
(2.152)
j
)| 1
2N
(2.153)
o
1 1 4r sen2 (
j
)1
2N
(2.154)
1
1)
2 sen2 ( (N2N
)
(2.155)
y como
h 0, N y sen2 (
(N 1)
)1
2N
(2.156)
2.6.2.
(2.157)
(2.158)
Definici
on 2.6.4 El esquema de diferencias finitas (2.157) es estable con respecto
a la norma k.k si existen constantes positivas x0 y t0 , y constantes no negativas
K y tal que
kvn+1 k K et kv0 k,
(2.159)
kv k Ke
J
X
kvi k.
(2.160)
i=0
(2.161)
(2.162)
(2.163)
kQn+1 v0 k
K et
0
||v ||
(2.164)
(2.165)
(2.166)
(2.167)
kvn k
Si tomamos el supremo sobre ambos lados de la inecuacion, con respecto a j,
obtenemos
kvn+1 k kvn k
(2.168)
Note que para la estabilidad del esquema (2.167) se ha requerido que r 1/2.
En este caso decimos que el esquema es condicionalmente estable.
2.6.3.
Al igual como se ha trabajado anteriormente asumimos que tenemos una sucesion de particiones de nuestro intervalo [0, 1] descrito por la sucesion de incrementos, {xj }, y una sucesion de espacios {Xj } con normas k.kj . Entonces decimos
(2.169)
(2.170)
u(0, t) = u(1, t) = 0, t 0
(2.171)
(2.172)
n+1
v0n+1 = vM
= 0, n = 0, . . .
(2.173)
vjn = f (jx), j = 0, . . . , M
(2.174)
2.7.
El Teorema de Lax
2.7.1.
(2.175)
de orden precisi
on (p, q) en la norma k.k para un problema de valor inicial lineal bien puesto y estable con respecto a la norma k.k, entonces el esquema es
convergente de orden de precisi
on (p, q) con respecto a la norma k.k.
Prueba: Sea u = u(x, t) la solucion exacta del problema de valor inicial. Entonces
puesto el esquema es de orden de precision (p, q), tenemos
un+1 = Qun + t Gn + t n
(2.176)
(2.177)
n
X
Qj nj
j=0
n+1
= t
n
X
Qj nj
(2.178)
j=0
(2.179)
n
X
kQj k k nj k
j=0
t K
n
X
ej t k nj k
j=0
t K e
(n+1)t
n
X
k nj k
j=0
(n + 1) t K e
(n+1)t
(2.180)
se tiene
(n + 1) t K et C1 (t) (O(xp ) + O(tq ) t K et C1 (t)(0) = 0
(2.181)
Note que el termino Gn que aparece en las expresiones (2.175) y (2.176) son aproximaciones del termino fuente. Asi, cuando sustraemos para formar una ecuacion
para w, ellos se sustraen. Cualquier error del truncamiento debido a como Gn se
aproxima el termino de la fuente esta contenido en el termino
2.7.2.
Como es constumbre, despues de tratar los problemas de valor inicial, volvemos a tratar los problemas de valor inicial-frontera. Empezamos asumiendo que
tenemos las mismas hipotesis que teniamos para las definiciones de convergencia, consistencia y estabilidad para problemas de valor inicial-frontera. Es facil
ver entonces que para una version del Teorema de Lax para un problema de valor
inicial-frontera se agrega la subscripcion j a cada norma k.kj . Por lo tanto tenemos
el siguiente teorema
2.7.3.
(2.182)
(2.183)
(2.184)
(2.185)
para todo w IR
Prueba: Ver [22] pag. 171
Observaci
on
Suponemos que q = 0 en (2.185) y que se utiliza la condicion de estabilidad
|g()| 1, esto es, suponemos que
|etq() | 1
|g(h)| 1
(2.186)
|g(h)n | CT
(2.187)
ku0 Sv k =
||>
h
|b
u0 ()|2 d,
(2.188)
(2.189)
en h y k, es decir 0 cuando h, k 0.
Se consideran en L2 la norma de u(., tn ) Sv n
ku(., tn ) Sv n k.
(2.190)
|u(x, tn ) Sv (x)| dx =
|eq()tn g(h)|2 |b
u0 ()|2 d
+
||>
h
|eq()tn u
b0 ()|2 d
(2.191)
Se considera que el lado derecho de (2.191) como una integral sobre IR, con la
particularidad que el integrando esta dado en dos partes discretas. Es decir, el
integrando es la funcion
h () =
|eq()tn g(h)n |2 , ||
|eq()tn .b
u0 ()|2
, || >
(2.192)
(2.193)
Se concluye que el lado derecho de (2.191) converge a cero para cada valor de
cuando h, k 0. As se tiene que el conjunto de funciones h que estan en L1 (IR)
tienden a cero para cada punto cuando h, k 0.
Antes de concluir que las normas de las funciones convergen a cero se necesita
lo siguiente
|eq()tn g(h)n |2 |b
u0 ()|2 < (2CT )2 |b
u0 ()|2
(2.194)
Esto demuestra que las funciones h son uniformemente acotadas por una funcion
en L1 (IR), es decir, 4CT2 |b
u0 ()|2 . Por el teorema de la convergencia dominante de
Lebesgue se concluye que
Z
h () d =
c |2 d
|b
u(, tn ) Sv
(2.195)
(2.196)
(2.197)
(2.198)
w
bM () =
M , | M | < M
en otro lugar
P
M =1
wM (x),
se demostrara que u0 esta en L2 (IR). Como los intervalos IM son disjuntos se tiene
que
Z
|u0 (x)| dx =
=
=
|b
u0 ()|2 d
XZ
|w
bM ()|2 d
M =1
X
2
2
M
M
M =1
X
2
= 2
M =1
2
,
3
CT 1
T
M
8
(2.199)
(2.200)
As,
Z
n
kSv u(., tn )k
|g(h)n eq()nk |2 |b
u0 ()|2 d
|M |M
1
2
[(1 + M kM )n CT ]2 M
2 M
2
1 + 12 M kM )n CT 2
= 2[
]
M
Se tiene por (2.199),
1
CT 1 2 T 2
kSv n u(., tn )k2 2( kM
)
90
2
M
32
(2.201)
2.8.
An
alisis de Estabilidad
En esta seccion presentamos y desarrollamos la importante herramienta
del Analisis de Fourier, con el cual analizamos los esquemas de diferencias finitas
y sus soluciones. El analisis de Fourier es una herramienta usada para estudiar
importantes propiedades de los esquemas de diferencias finitas y sus soluciones.
Usamos el analisis de Fourier desde el principio para estudiar los esquemas de
diferencias finitas y las ecuaciones diferenciales parciales.
2.8.1.
An
alisis de Estabilidad para Problemas de Valor
Inicial
(2.202)
(2.203)
(2.204)
eiwx u
b(w, t) dw
(2.205)
1 X i j
vb() =
e
vj
2 j=
(2.206)
para [, ].
L2 ([, ]) = {f : [, ] IR :
|f ()|2 d < }
ei j u
b() d
(2.207)
(2.208)
[
, ].
h h
La transformada discreta inversa de Fourier esta dada por
1
vm (x) =
2
2.8.2.
Z+ h
eimh vb() d
(2.209)
La Transformada de Fourier
Integral de Fourier.
Transformada de Fourier de una funci
on: Consideremos la aplicacion integral
Z
T (f (x)) =
(2.210)
llamada Transformaci
on Integral donde K(x, t) es el n
ucleo.
La Integral de Fourier: Sea f una funcion periodica de periodo T = 2L definida
en [L, L], entonces (f se asume uniformemente convergente)
f (x)
a0 X
n
n
+
(an cos
x + bn sen x)
2
L
L
n=1
(2.211)
Z L
Z L
1X
n
1
n
f (t)dt +
f (t) cos
f (x) =
t. cos
x dt
2L L
L n=1
L
L
L
Z L
n
n
f (t)sen t.sen x dt]
+
L
L
L
Z L
Z
1X L
n
1
f (t)dt +
f (t) cos( (t x)) dt
=
2L L
L n=1 L
L
Ahora como
a0
1
=
2
2L
1
f (t)dt
2L
L
yL
1
2L
n
L
y =
|f (t)|dt
L
|f (t)|dt 0
L
1X
f (x) = lm
L L
n=1
Si hacemos =
Asi que
f (t) cos
L
n
L
(t x) dt
(2.212)
1X
f (x) = lm
f (t) cos (n (t x)) dtd
L
n=1
(2.213)
Si hacemos
cos (t x) =
ei(tx) + ei(tx)
2
(2.214)
Z Z
Z Z
1
1
i(tx)
f (t)e
dtd +
f (t)ei(tx) dtd
2 0
2
Z 0 Z
Z Z
1
1
f (t)ei(tx) dtd +
f (t)ei(tx) dtd
2
2 0
Z Z
1
f (t)ei(tx) dtd
2
Z Z
1
f (t)eit eix dtd
2
Z
Z
1
1
f (t)eix dt eix d
2 2
|
{z
}
Si
1
F () =
2
entonces
1
f (x) =
2
f (t)eit dt
(2.215)
F ()eix dx
(2.216)
1 , |x| < a
f (x) =
0 , |x| > a
entonces la Transformada de Fourier esta dada por
F () =
f (u)e
iu
+a
du =
Para = 0 se tiene
F (0) =
Asi que
F () =
sen a
,
6= 0
+a
du =
2.8.3.
(1)eiu du = 2
du = 2a
a
2 sen a , 6= 0
2a
, =0
Analisis de Fourier
La herramienta que usaremos mas extensamente en nuestro estudio de estabilidad y problemas bien puestos es el analisis de Fourier.
Usaremos el analisis de Fourier sobre los reales IR o sobre una malla de enteros,
Z, o h Z, el cual es definido por h Z = {hm : m Z}.
La transformada inversa de Fourier muestra que u puede ser recuperada de u
b.
La transformada inversa de Fourier es una formula que expresa la funcion u(x)
como una superposicion de ondas, dadas por eiwx , con diferentes amplitudes u
b(w).
Si por ejemplo, la variable x representa el tiempo, entonces, la variable real w
representa la frecuencia y la funcion transformada u
b, esta definida en el espacio
de frecuencias. Por otro lado si x representa una variable espacial, entonces, la
variable w representa el n
umero de onda, en este caso la funcion transformada u
b
esta definida en el espacio de n
umero de onda.
En la practica se define la transformada de Fourier para una funcion cuadrado
integrable
|u(x)|2 dx <
(2.217)
1 , |xj | 1
vj =
1/2 , |xj | = 1
0 , |x | 1
j
donde h = x = M 1 para alg
un entero M .
Tomemos la transformada de Fourier a vj
1 X i j h
e
vj h
vb() =
2 j=
1
1
1
1
1
= ei j h ( ) h + ei j h ( ) h +
2
2
2
2
2
M
1
X
ei j h
j=(M 1)
h
h sen(M 12 )h
= cos +
sen 12 h
2
2
h sen(M h ) cos( 12 h ) cos(M h ) sen( 12 h )
h
= cos +
sen 12 h
2
2
h
h
h
1
= cos + sen cot( h ) cos
2
2
2
2
1
h
= sen cot( h )
2
2
h2
=
2
/4
1
sen2 () cot( h) d
2
/4
= kuk2h
1
= h[( )2 +
2
M
+1
X
1
(1 + ( )2 )]
2
j=(M +1)
1
1
= h[ + M 1 (M + 1) + 1] = h[ + 2M 1]
2
2
1
= h[2M ]
2
1
=2 h
2
Ejemplo 15 Sea vj =
1
,
j2
1 X 1 ij
1 X 2 eij + eij
vb() =
e
=
=
2
2 j= j 2
2 j=1 j 2
2X 1
cos(j)
j=1 j 2
Proposici
on 2 La suseci
on vn es estable en la norma l2,x si y solo si la sucesi
on
bn es estable en L2 ([, ])
v
2.8.4.
Factor de Amplificaci
on
vj+1 ei j h
j=
vj ei(j1)h
j=
ih
= e
vj ei j h
j=
ih
= e
v.
(2.218)
1 X ij
F(v) =
e vj
2 j=
(2.219)
S+ v = {uj+1 }
(2.220)
S v = {uj1 }
(2.221)
F(S v) = ei F(v)
(2.222)
2.8.5.
(2.223)
u(1, t) = 0, t > 0
(2.224)
(2.225)
(2.226)
n+1
vM
= 0
(2.227)
vj0 = f (jx), j = 0, . . . , M.
(2.228)
Soluci
on : Observe que por la Proposicion (3), si consideramos la ecuacion de
diferencias (2.226) como un esquema de diferencias para un problema de valor
inicial, la estabilidad de este esquema sera necesario para la estabilidad del esquema
(2.226)-(2.228).
(2.229)
(2.230)
2.8.6.
(2.231)
donde 0 k M , n = 0, . . . y i2 = 1.
Dado un esquema diferencias finitas para la ecuacion de difusion es posible
llegar a determinar el factor de amlificacion
g = 1 4rsen2
jx
, j = 1, . . . , M
2
(2.232)
2.8.7.
Condici
on de Estabilidad de von Neumann
(2.233)
k
vn k2h
|
v ()| d =
(1 + Ck)2n
|g|2n |
v 0 ()|2 d
|
v 0 ()|2 d = (1 + Ck)2n k
v 0 k2h
x > 0, tenemos
, 0 < k k0 , 0 < h h0
|g|
s
1 + Ck
para todo
[1 , 2 ]
0,
/ [1 , 2 ]
0
=
v () =
q
h
, [ , ]
2 1
Observamos que k
v 0 k = 1.
0,
q
/ [ h1 , h2 ]
h
,
2 1
/ [ h1 , h2 ]
k
vn k2h
Z
=
Z
=
h
2
h
h1
|
v n ()|2 d
|g|2n |
v 0 ()|2 d
|g|2n |
v 0 ()|2 d
Z
(1 + Ck)
2n
2
h
h1
|v 0 ()|2 d = (1 + Ck)2n
Por tanto para cada CT = 21 e2CT existe T > 0 tal que 1 + CT CT y se cumple
kv n k2h CT kv 0 k2h
que es la no estabilidad del esquema.
(2.234)
Proposici
on 5 Es esquema de diferencias finitas
vn+1 = Qvn
(2.235)
(2.236)
Proposici
on 6 El esquema de diferencias
vn+1 = Qvn
(2.237)
(2.238)
(2.239)
entonces
|g()|n+1 e(n+1)Ct .
(2.240)
(2.241)
(2.242)
Por lo tanto |g()|n+1 no puede ser acotado. Esto contradice la Proposicion (6).
(2.243)
|vjn+1 |2
j=
n
|2
|vjn + vj+1
j=
j=
n
n
||2 |vjn |2 + 2|||||vj+1
||vjn | + ||2 |vj+1
|2
n
n
||2 |vjn |2 + ||||(|vj+1
|2 + |vjn |2 ) + ||2 |vj+1
|2
j=
(|vjn |2
j=
n
|vj+1
|2 )
=2
|vjn |2
j=
tenemos
=
=
j=
X
j=
j=
= (|| + ||)2
|vjn |2
j=
X
j=
|vjn+1 |2
(|| + ||)
X
j=
|vjn |2
(2.244)
|vjn+1 |2
(|| + ||)
j=
|vjn |2
j=
X
2
(|| + ||)
|vjn1 |2
j=
(|| + ||)2n
|vj0 |2
j=
Entonces
|vjn+1 |2
(|| + ||)
j=
2n
|vj0 |2
j=
j=
X
X
|vj0 |2
j=0 j=
Si hacemos
CT = (|| + ||)2n
|| bf v
n+1
||v
n+1
||h CT (h
||h
|vj0 )|2 )
j=
CT ||v0 ||h
Asi que el esquema es estable si |CT | 1, esto es, ||+|| 1. En el caso particular
del esquema de forward-time forward-space, tenemos que = 1 + Cr y = Cr ,
asi que
|| + || = |1 + Cr | + | Cr | 1
Entonces
|1 + Cr | + |Cr | 1
Por lo tanto, para a < 0 se tiene 0 Cr 1 o 1 Cr 0.
Mostraremos que la condicion Cr 1 es necesaria para la estabilidad de muchos
de los esquemas explcitos para la ecuacion de conveccion (??).
Ejemplo 18 Aqui incluimos una alternativa mucho mas facil de analizar la estabilidad del esquema de diferencias finitas
1
n
)
vjn+1 = vjn Cr (vj+1 vj1 ) + r(vj+1 2vjn + vj1
2
(2.245)
(2.246)
(2.247)
tal que
Si tomamos una constante r0 tal que r0 1/2, entonces (2.247) se puede escribir
r0
= (1 4r0 sen2 )2 + a2 tsen2
2
C = a2
r0
(2.248)
De manera similar se puede probar que no existe esquemas consistentes explcitos para ecuaciones hiperbolicas que sean estables para todos los valores de r (con
constante cuando h, k 0).
En este sentido presentamos el siguiente teorema; probado primeramente por
Courant,Friedrichs y Lewy.
Teorema 2.8.3 No existe esquemas de diferencias finitas consistentes, incondicionalmente estables explcitos para sistemas hiperb
olicos de ecuaciones diferen-
ciales parciales.
Ejemplo 19 Mostraremos que el esquema implcito
n+1
vjn+1 vjn
vjn+1 vj1
+a
=0
k
h
(2.249)
(2.250)
Elevando al cuadrado
n+1 2
(1 + Cr )|vjn+1 |2 = |vjn + Cr vj1
|
n+1
n+1 2
|vjn |2 + 2Cr |vjn ||vj1
| + Cr2 |vj1
|
n+1 2
n+1 2
|vjn |2 + Cr |(vjn |2 + |vj1
| ) + Cr2 |vj1
|
n+1 2
= (1 + Cr )|vjn |2 + (Cr + Cr2 )|vj1
|
X
j=
|vjn+1 |2 (1 + Cr )
|vjn |2 +
j=
X
n
(Cr + Cr2 )
|vj1
|2
j=
X
= (1 + Cr )
|vjn |2 +
j=
Cr (1 + Cr )
= ((1 + Cr ) + a(1 + Cr ))
j=
X
j=
= (1 + 2Cr + Cr2 )
X
j=
|vjn |2
|vjn |2
|vjn |2
Entonces
(1 + Cr )
|vjn+1 |2
(1 + Cr )
j=
|vjn |2
j=
de donde se obtiene
|vjn+1 |2
j=
|vjn |2
j=
j=
|vjn1 |2
|vjn2 |2
j=
...
|vj0 |2
j=
Entonces
|vjn+1 |2
j=
|vj0 |2
j=
kvn+1 kh
0
X
kvj kh
j=0
o lo que es lo mismo
h
X
m=
|vjn |2
CT h
X
X
|vjn |2
j=0 m=
para J = 0. Por lo tanto el esquema es estable para todo valor de h/k cuando
a > 0.
Nota
Siempre es posible elegir el n
umero de Courant Cr arbitrariamente grande
para que el esquema (2.249) se estable como se demuestra teoricamente, pero es
Nota
1. Las soluciones de los esquemas de diferencias finitas muestran que la inestabilidad esta relacionado con oscilaciones de alta frecuencia.
2. La precision sera menos buena cuando hay discontinuidades en las condiciones iniciales o en sus derivadas. Puede demostrarse analticamente que
cuando los valores en la frontera son constantes el efecto de las discontinuidades en los valores iniciales y sus derivadas iniciales en la solucion de una
ecuacion parabolica disminuye cuando t se incrementa.
3. Inestabilidad es el crecimiento rapido de altas frecuencia en la solucion del
esquema de diferencias finitas.
4. La inestabilidad es esencialmente un fenomeno de caracter local, esto es,
las oscilaciones acurren en los puntos donde la derivada de la solucion es
discontinua. Es claro que, las oscilaciones causadas por la inestabilidad se
propaga a otras regiones, que pueden hacer que la perturbacion finalmente
parezca global en la extencion.
5. El termino de error es definido por
Error = (
tama
nos de paso. Naturalmente, cuando el valor del paso del tiempo se incrementa,
el total de puntos malla en el dominio computacional decrece y, como resultado, el
tiempo computacional disminuye. Sin embargo, estas ventajas son acompa
nadas
por un aumneto de error.
Considerar que para todo el analisis de estabilidad se impone limitaciones sobre
algunos esquemas de diferencias. En segundo lugar, la precision de la solucion y
el tiempo computacional requerido para generar la solucion juega un importante
rol en la eleccion del tama
no de paso adecuado. Una vez que se ha obtenido la
solucion, esta debera ser comparada con otras soluciones, analticas o numericas,
y para datos experimentales si es posible. Cuando se gana experiencia en estos
metodos y en su comportamiento, facilmente se puede elegir el tama
no correcto
de paso y la tecnica numerica apropiada.
Seleccionar una tecnica numerica depende del problema planteado y la naturaleza de las ecuaciones gobernantes del problema y ademas las condiciones iniciales
y de frontera que son impuestas. Cada uno de los metodos descritos tienen sus
propias ventajas y desventajas. As cuando un problema es planteado, las ventajas
y desventajas del metodo numerico disponible debe ser cuidadosamente pensado
antes de seleccionar un algoritmo en particular.
Si hacemos un resumen sobre aplicaciones y limitaciones del metodo de estabilidad de von Neumann
1. El analisis de estabilidad de von Neumann puede ser aplicado solo a ecuaciones lineales.
2. La influencia de la condicion de frontera sobre la estabilidad de la solucion
no es incluida.
3. Para una EDP escalar la cual es aproximada por un esquema de diferencias
finitas, el requerimiento matematico es impuesto sobre el factor de amplificacion como sigue:
a. Si g es real, entonces |g| 1
b. si g es compleja, entonces |g|2 1, donde |g|2 = g g
4. Para una EDP escala la cual es aproximada por un EDF de tres niveles,
el factor de amplificacion es una matriz. En este caso el requerimiento es
impuesto sobre los autovalores de la matriz de amplificacion G como sigue:
a. Si es un autovalor real, entonces || 1
b. Si es un autovalor complejo, entonces ||2 1
5. El metodo de diferencias finitas puede ser faclmente extendido a problemas
bidimensionales.
6. El procedimiento puede ser usado para analizar la estabilidad de un sistema
de EDPs. El requerimiento es impuesto sobre los autovalores mas grandes de
la matriz de amplificacion.
7. Los valores de estabilidad para problemas uni dimensionales no estacionarios
puede ser establecida como sigue:
a. Para formulaciones explcitas: N
umero de Courant Cr 1, N
umero de
Difusion r 0.5 y N
umero de celda de Reynolds Rec
2
Cr
2.9.
An
alisis de Error
En el metodo de diferencias finitas tenemos esquemas de primir orden de precision. Un ejemplo de un esquema de primer orden de precision para la ecuacion
modelo
u
u
+a
=0
t
x
si utilizamos diferencias progresivas en el tiempo y diferencias regresivas en el
espacio obtenemos
un+1
unj
unj unj1
j
+a
=0
t
x
Note que, en el proceso de aproximacion, la serie de Taylor fue truncada en la
segunda derivada, es decir,
un+1
unj
u
t 2 u (t)2 3 u
j
=
+
+
+ ...
t
t
2! t2
3! t3
| O(t)
Para una formulacion de segundo orden de precision, el termino del error de truncamiento esta en la tercera derivada; por ejemplo,
un1
un+1
u
2(t)2 3 u
j
j
=
+ ...
t
2t
3!
t3
| O(t)2
2.10.
Ecuaci
on Modificada
Para determinar el termino de error dominante de una ecuacion de diferencias, el desarrollo de la serie de Taylor fue reemplazada por una ecucacion en
diferencias y, despues de alguna manipulacion algebraica, la as llamada ecuacion
modificada es obtenida. Para ilustrar el procedimiento consideremos dos ejemplos.
Los ejemplos seleccionados representan un orden de precision de primer y segundo.
Estos ejemplos mostraran el termino de error dominante de estos esquemas y su
relacion con los errores de disipacion y dispersion.
Ejemplo 20 Un esquema de diferencias de primer orden para la ecuaci
on
modelo
ut + aux = 0
es dado por el esquema Upwind
n
vjn+1 vjn
vjn vj1
+a
=0
t
x
n
Los terminos vjn+1 y vj1
son expandidos en su serie de Taylor como sigue
v (t)2 2 v (t)3 3 v
+
+
+ O(t)4
t
2! t2
3! t3
(2.251)
v (x)2 2 v (x)3 3 v
+
+ O(x)4
t
2! t2
3! t3
(2.252)
+
a
t
x
2 t2
2 x2
2 3
2 3
(t) v
(x) v
a
+ O(t)3 + O(x)3
3
3
6 t
6 x
(2.253)
Analizamos la ecuacion (2.253) y comparamos esta con la ecuacion diferencial parcial original, las derivadas de orden superior en el tiempo deben ser reemplazadas
por las derivadas espaciales. Esta sustitucion requiere la determinacion de
un orden de error [O(t)2 , O(x)2 ] y
3v
t3
2v
t2
con
2v
,
t2
2v
2 v t 3 v
x 3 v
= a 2
+
a
xt
x
2 xt2
2 x3
2
4
2 4
(t) v
(x) v
3
3
a
+
O(t)
,
O(x)
6 xt3
6 x4
(2.254)
+
a
t2
tx
2 t3
2 tx2
6 t4
x 4 v
3
3
a
+
O(t)
,
O(x)
6 tx3
(2.255)
3
2v
2v
t 3 v
2v
2 t v
+ 2 = a2 2 + a
a
xt t
x
2 xt2
2 x3
4
4
t v
2v
t 3 v
x 3 v
2 t v
+ a
+a
a
+a
6 xt3
6 x4
tx
2 t3
2 tx2
2 4
x
(t) v
a
+ O(t)3 , O(x)3
4
6 t
6
Simplificando
2
t
3v
3v
2v
2 v
= a
+
a
3
t2
x2
2
xt2
t
3
3
x
v
2 v
2
2
+
a
a
+
O(t)
,
O(x)
2
tx2
x3
Esta ecuacion requiere que determinemos
3
3v
, v
t3 tx2
3v
,
t2 x
(2.256)
v
t
3v
t 4 v
x 4 v
+
a
t2 x
2 t4
2 x2 t2
2 5
2
(t) v
(x) 5 v
3
3
a
+
O(t)
,
O(x)
(2.257)
6 t5
6 t2 x3
= a
3
3v
t
4v
4v
2 v
= a
+
a 2 2
xt2
x3
2
x t
xt3
x
4v
2 v
+
a
+ O(t)2 , O(x)2
a
3
4
2
tx
x
(2.258)
3
4
t
4v
3v
4v
3 v
2 v
= a
+
4
a 2 2 a
t3
x3
2
t x
tx3
t
4
4
4
t
v
v
3 v
2
2
+
+
a
+
O(t)
,
O(x)
a 2 2 a2
2
t x
tx3
x4
Puesto que estamos interesados en la relacion de presicion de primer orden, podemos escribir
3
3v
3
=
a
+ [O(t), O(x)]
t3
x3
(2.259)
(2.260)
3
3v
2 v
=
a
+ (O(t), O(x))
t2 x
x3
(2.261)
2
3
2v
t
2 u
2 v
= a
+a
a
+ [O(t), O(x)]
t2
x2
2
x3
3
t
3 v
a
+ [O(t), O(x)]
2
x3
x
3v
+ a
a 3 + [O(t), O(x)]
2
x
3
x v
+ O(t)2 , O(x)2
a2
3
2 x
Por lo tanto,
2
2v
3v
2 v
3
2
=
a
+
(a
t
a
x)
+ O[(t)2 , t, x, (x)2 ]
t2
x2
x3
(2.262)
y
3
3v
3 v
= a
+ O[t, x]
t3
x3
(2.263)
(t) v
(x) v
3
2
2
3
+ a3
a
+
O
(t)
,
(x)(t)
,
(x)
(t),
(x)
6 x3
2 x3
(2.264)
(2.265)
v
v (t)2 3 v
(x)2 3 v
4
4
= a
a
+
O
(t)
,
(x)
t
x
6 t3
6 x3
(2.266)
(2.267)
3
v
v
(x)2 2 t 2
v
= a
+a
a(
) 1
+ O[(t)4 , (x)3 (t), (x)4 ] (2.268)
t
x
6
x
x3
Esta es la ecuaci
on modificada de la ecuacion modelo
u
u
= a
t
x
(2.269)
Claramente se indica el t`ermino del error dominate. Note que su termino de error
incluye derivadas de tercer orden.
2.11.
Viscosidad Artificial
El termino de error en el esquema (2.264) de primer orden de precision es claramente indicado con el segundo termino del lado derecho, empezando el termino
dominante del error. El coeficiente que lo denotamos por c
h
c = a (1 Cr )
2
(2.270)
+
a
+ O((t)2 , (x)2 )
t
x
2 t2
2 x2
(2.271)
(2.272)
u
= a
a
x
x
2
u
= a2 2
x
Por lo tanto
2
2u
2 u
=
a
t2
x2
(2.273)
u
u
2 u t
x 2 u
= a
a2 2
+a
+ O[t2 , x2 ]
t
x
x 2
2 x2
(2.274)
Por lo tanto,
u
u
x
t 2 u
= a
+ (a
a2 ) 2 + O[t2 , x2 ]
t
x
2
2 x
u
x
t 2 u
= a
+a
(1 a
)
+ O[t2 , x2 ]
x
2
x x2
o, en terminos del n
umero de Courant
u
x
2u
u
= a
+a
(1 Cr ) 2 + O[t2 , x2 ]
t
x
2
x
(2.275)
x
(1 Cr )
2
(2.276)
Captulo 3
Ecuaci
ones Diferenciales
Hiperb
olicas
3.1.
(3.1)
(3.2)
3.1.1.
Clasificaci
on de las ecuaciones de primer orden
aij
uij X uj
+
bij
= ci
x
y
j=1
i = 1, ..., n
(3.3)
(3.4)
Pn () = det(A B)
El sistema (3.5) es
1. Elptico si Pn () no tiene races reales.
(3.5)
2. Hiperb
olico si Pn () tiene raices reales y distintas, o si Pn () tiene n raices
en donde por lo menos una se encuentra repetida y el problema generalizado
de valor propio (At B t )V = 0 proporciona n vectores propios linealmente
independientes.
3. Parab
olico si Pn () tiene n raices reales en donde por lo menos una se
encuentra repetida y el problema generalizado anterior de valor propio proporciona menos de n vectores propios linealmente independientes.
En el caso que Pn () tiene raices reales y complejas, no se puede llevar a cabo una
clasificacion exhaustiva del sistema (3.5).
Ejemplo 22 Los sistemas (3.1) y (3.2) son hiperb
olicos.
3.2.
M
etodo de las Caractersticas
En esta parte nos referimos a las curvas caractersticas de una EDP hi-
3.2.1.
Problema Lineal
(P V I)
u(x, 0) = u (x),
0
x IR
x IR, t > 0
(3.6)
(3.7)
(3.8)
Entonces comparando (3.6) con (3.7) podemos seguir los siguientes pasos para
hallar la solucion del PVI:
Paso 1: Solucionar las dos ecuaciones diferenciales ordinarias
dt
= a(x, t), t(0) = 0
ds
dx
= b(x, t), x(0) = x0
ds
(3.9)
u(0) = u ( ), IR
0
(3.10)
u(x, 0) = sen(x), x IR
t(0) = 0
x(0) = x0
resolviendo tenemos
t(s) = s,
x(s) = as + x0 ,
s>0
s>0
x0 IR
s>0
u(0) = sen(x0 )
la cual tiene por solucion general u(s) = Ce2s . Verificamos que se cumplan las
condiciones iniciales se tiene que
u(s, ) = e2s sen(x0 )
Regresando a las variables originales tenemos la solucion del PVI
u(x, t) = e2t sen(x at)
u(x, 0) = u (x), x IR
0
Resolviendo las ecuaciones diferenciales ordinarias
dt
= 1, t(0) = 0
ds
dx
= a, x(0) = x0
ds
tenemos
t(s) = s,
x(s) = as + x0 ,
s>0
s>0
IR
a >0
s>0
t>0
la cual representa que la onda inicial. Si a > 0 es desplazada a la derecha con una
velocidad constante a.
Para examinar los problemas que incluyen conveccion, podemos tomar la ecuacion que modela una Ley de conservacion lineal
ut + aux = 0
(3.11)
, )
t x
(3.12)
(a, 1)
1
(x,t)=(xo+at,t)
vector
direccion
x=at+xo
(Xo, 0)
Ahora, fijamos un punto sobre el eje x, por ejemplo (x0 , 0). La recta que pasa a
traves de este punto, paralela a x = at es dado por x = x0 + at. Desde que nuestra
solucion es constante a lo largo de esta recta tenemos
u(x, t) = u(x0 + at, t) = u(x0 , 0)
(3.13)
(3.14)
(3.15)
(3.16)
u(x, 0)
= u0 (x),
< x < .
(3.17)
= t,
t = ,
o
(3.18)
= x at,
x = + a.
luego definiendo u(, ) = u(x, t), donde (, ) y (x, t) estan relacionadas por
el cambio de variable definida previamente, entonces la ecuacion en (3.17), se
transforma en
u
t
x
=
ut +
ux = ut + aux = bu + f ( + a, )
es decir
u
= b
u + f ( + a, )
u(x, t) = u0 (x at)e
(3.19)
(3.20)
Notar que u(x, t) depende solamente de los valores de (x1 , t1 ) tal que x at =
x1 at1 , es decir, la solucion depende u
nicamente de los valores de u y f sobre
las caractersticas que pasan por (x, t), para 0 t1 t. Este metodo de solucion
puede extenderse facilmente a ecuaciones no lineales de la forma
ut + aux = f (x, t, u).
3.3.
Problema no lineal
Sea el PVI
(P V I)
ut + g(u)ux = 0,
x IR, t > 0
u(x, 0) = u (x),
0
x IR
(3.21)
en donde por analoga con el problema lineal podemos suponer que una partcula
de agua se empieza a mover desde x0 con una velocidad de g(u). Entonces, despues
de t segundos la posicion de la partcula x sera:
x = x0 + g(u)t
Como la concentracion dada por la funcion u(x, t) permanece constante a lo largo
de la curva caracterstica, podemos escribir
g(u) = g[u(x0 , 0)]
por lo tanto la curva caracterstica que comienza en (x0 , 0) tiene por ecuacion
x = x0 + g[u(x0 , 0)]t
Entonces la solucion del PVI no lineal esta dada por la solucion implcita de la
funcion
u(x, t) = u0 (x) = u0 (x g[u(x0 , 0)]t)
entonces
u(x, t) = u0 (x g[u]t)
(3.22)
1 , x>0
0 , x<0
Las curvas caractersticas que comienzan en el punto (x0 , 0) estan dadas por
x = x0 + g[u(x0 , 0)],
g(u) = 3u
luego
para x0 < 0
x = x0 + g[0]t
x = x0
para x0 > 0
x = x0 + g[1]t
x = x0 + 3t
por tanto la grafica de las curvas caractersticas estan dadas en la figura.
t
t1
xo
xo
X1
3.4.
conservacion de una cierta cantidad dada por ejemplo por la funcion u, la cual
representa la energa, masa, momento, temperatura, etc.
Sea IR. La forma general de un sistema de leyes de conservacion en varias
ut +
fi (u) = 0
xi
i=1
(3.23)
donde
m
u : IRn IR+
0 IR ;
fi : IRm IRm
u es llamado el vector de las variables de conservacion. Las funciones fi diferenciables son llamadas funciones de flujo, las cuales representan un mecanismo de
transporte de la cantidad u.
m
Definici
on 3.4.1 Sean las funciones u : IRn IR+
0 IR ;
fi : IRm
n
X
[fi (u)]xi = 0
(3.24)
i=1
(3.25)
f : IRm IRm
blema de Valor Inicial. Este sistema tiene solucion suave o clasica si la solucion u
es por lo menos de clase C 1 y si ademas la matriz Jacobiana de f , funcion de flujo,
tiene solo valores propios reales y es diagonalizable para cada valor de u, entonces
el sistema (3.24) es Hiperbolico.
ut + 2ux + vx = 0
v + u + 2v = 0
t
x
x
se puede escribir como
2 1 u
u
= 0.
+
1 2
v
v
t
3.4.1.
Condiciones de Contorno
(P V IC)
ut + aux = 0,
0 < x < 1, t > 0
u(x, 0)
= u0 (x),
< x < .
u(0, t)
= g(x, t), t , 0 < x < 1.
(3.26)
a>0
a<0
condiciones
condicion
de
de
frontera
frontera
u(x,0)=f(x)
1
u ( x,0 ) = f ( x )
u(x, t) =
u0 (x at)
, si x at > 0
g(t a1 x) ,
x at < 0
3.5.
Derivaci
on de Leyes de Conservaci
on
x2
x1
x2
m=
(x, t)dx
(3.27)
x1
(3.28)
luego, la razon de cambio de masa en [x1 , x2 ] viene dado por la diferencia de flujos
en x1 y x2 , esto es,
dm
= (x1 , t)v(x1 , t) (x2 , t)v(x2 , t)
dt
(3.29)
x2
(3.30)
x1
difcil resolver.
Otra forma integral de (3.30) es expresar la masa en [x1 , x2 ] en el tiempo t2 ,
(t2 > t1 ) en terminos de la masa en el tiempo t1 , y el flujo total en las fronteras
para este perodo. Par obtener esta nueva forma integral, basta integrar (3.30) en
[t1 , t2 ] dando lugar a la expresion
Z
x2
x2
(x, t2 )dx =
(x, t1 )dx +
x1
x1
t2
t2
(3.31)
Ahora para obtener una forma diferencial de la ley de conservacion, asumiremos
que las funciones (x, t) y v(x, t) son diferenciables. Entonces, la variacion de la
densidad en el intervalo de tiempo [t1 , t2 ] esta dado por
Z
t2
(x, t)dt
t
(x, t2 ) (x, t1 ) =
t1
x2
x2
x2
Entonces
Z
x2
x1
t2
t1
(x, t)dtdx =
t
t2
t1
x2
x1
t2
t1
x2
x1
(x, t)dxdt =
t
t2
t1
x2
x1
factorizando
t2
x2
[
t1
x1
(x, t) +
(x, t)v(x, t)]dxdt = 0
t
x
Esta u
ltima ecuacion es valida para cualquier seccion [x1 , x2 ] sobre el intervalo de
tiempo [t1 , t2 ] y en consecuencia se tiene
(x, t) +
(x, t)v(x, t) = 0
t
x
(3.32)
(x, t) +
f ((x, t)) = 0
t
x
(3.33)
(3.34)
Et + [v(E + P )x ] = 0
(3.35)
(x, t)
E(x, t)
entonces el sistema de ecuaciones (3.32), (3.34) y (3.35) puede ser escrito como
ut + [f (u)]x = 0
donde
(3.36)
f (u) = v 2 + P
v(E + P )
3.6.
Leyes de Conservaci
on Escalar
Ecuaci
on Lineal de Convecci
on o Advection
Asuma que la velocidad es constante, es decir, a = v(x, t) en la ecuacion (3.32),
entonces con una variable arbitraria u en lugar de se tiene
ut + aux = 0
(3.37)
x IR,
t>0
(3.38)
La solucion del Problema de Cauchy (3.37) y (3.38) es una traslacion de la condicion inicial u0 (x), esto es,
u(x, t) = u0 (x at),
t0
u(x(t), t) =
u(x(t), t) +
u(x(t), t).x0 (t)
dt
t
x
= ut + aux
= 0
Entonces, confirmamos que u es constante a lo largo de la curvas caractersticas
cuando u es suave.
Si suponemos que la velocidad v(x, t) = a(x), donde a(x) es una funcion suave.
Luego la ecuacion (3.32) tiene la forma
ut + (a(x)u)x = 0
que es una Ecuaci
on de Advection, la cual se puede escribir en la forma
ut + a(x)ux + a0 (x)u = 0
(3.39)
En este caso, las curvas caractersticas son solucion de la ecuacion diferencial or-
dinaria
x0 (t) = a(x(t))
x(0) = x0
Luego la solucion u(x, t) no es constante a lo largo de estas curvas, pues
d
u(x(t), t)) = a0 (x(t))u(x(t), t) 6= 0
dt
Dominio de Dependencia: Si consideramos el problema de valor inicial (3.37) (3.38). La solucion en el punto (x, t) depende del valor de u0 en el punto x0 , donde
x0 = x at. El punto x0 es llamado dominio de dependencia del punto (x, t). El
Dominio de Dependencia del punto (x, t) es el conjunto de todos los puntos en el
cual la solucion del Problema de Valor Inicial (3.37) - (3.38) en el punto (x, t) esta
dependiendo.
3.7.
Esquemas de Diferencias
Esquemas Explcitos
F orward time f orward space scheme
(U pwind hacia delante, a < 0)
n
vjn+1 vjn
vjn
vj+1
+a
= 0
k
h
(3.40)
(3.41)
(3.42)
(3.43)
(3.44)
(3.45)
Esquemas Implcitos
Backward time central space scheme
n+1
n+1
vjn+1 vjn
vj+1
vj1
+a
= 0
k
2h
(3.46)
(3.47)
(3.48)
(3.49)
Cada uno de los seis esquemas (3.40) y (3.45) puede ser escrito como una
combinacion lineal vjn+1 de los valores de v en los niveles n y n 1. Por ejemplo,
el esquema (3.40) puede ser escrito como
n
vjn+1 = (1 + Cr )vjn Cr vj+1
de v sobre estos niveles iniciales. Por ejemplo, para usar el esquema de Leapfrog
podemos especificar los valores de vj0 y vj1 para todo j, o podemos especificar alg
un
esquema para calcular los valores de vj1 de los valores de vj0 . En cualquier caso el
esquema Leapfrog debe ser usado para calcular vjn , n > 1.
Cuando nos referimos al esquema de Leapfrog no siempre hacemos distincion
entre los dos tiempos iniciales de calculo. Muchas de las propiedades del esquema
de Leapfrog son independientes del metodo usado para inicializar la solucion. Es
de practica usar un esquema de un paso para inicializar el esquema. Para tener
un panorama de los topicos que seguiremos discutiendo, hagamos un comentario
de los resultados obtenidos con los esquemas propuestos anteriormente, para ello
consideremos el (P V I), donde 0 < x < 2, t > 0 teniendo como condicion inicial
la funcion 2-periodica seccionalmente continua
0, 0 x < 2
,
u(x, 0) =
1, 2
x 4
,
3
3
0, 4 < x 2,
3
(3.50)
A continuacion se muestra una seccion de pseudocodigo, el cual nos permitira observar el comportamiento del periodo principal de esta funcion considerado como
una onda, para esto elejimos el esquema de Lax-Friedrichs, los demas tienen similar estructura
Do n = 0, Nmax
Do j = 19, 29
End Do
End Do
(3.51)
Como sabemos la solucion de este problema esta dado por u(x, t) = f (x at) y
es constante a lo largo de las caractersticas x at = cte. En las Figuras siguientes
presentamos las soluciones numericas que son aproximaciones de la solucion del
problema (??) y (3.50), con a = 1, h = 2/240, k = 2h/3 en t = 500k y en
t = 1000k.
upwind
LeapFrog
1.5
0.9
0.8
0.5
0.7
0
0.6
0.5
0.5
0.4
0.3
1.5
0.2
2
0.1
6
n=500,
8
n=1000
10
12
14
2.5
6
n=500,
8
n=1000
10
12
14
LaxWendroff
LaxFriedrichs
1.4
1.2
0.9
0.8
0.7
0.8
0.6
0.6
0.5
0.4
0.4
0.2
0.3
0.2
0.2
0.4
0.1
6
n=500,
8
n=1000
10
12
14
6
n=500,
8
n=1000
10
12
En las Figuras (3.7) - (3.10) la solucion exacta de la ecuacion diferencial esta dada por la linea punteada y la solucion numerica para el esquema respectivo, esta dado por la linea contnua en cada caso.
En las Figuras podemos observar que los esquemas de diferencias calculan la
solucion adecuadamente lejos de las discontinuidades, pero cerca de estas no se
muestra una buena aproximacion a la solucion.
Observar que el esquema Leapfrog genera muchas oscilaciones cerca de las discontinuidades. El metodo Lax-Wendroff tiene sobreoscilaciones y suboscilaciones
cerca de las discontinuidades, pero no se propagan immediantamente. Los metodos de Lax-Friedrichs y Upwind, no tienen suboscilaciones ni sobre oscilaciones,
sin embargo cerca de las discontinuidades se separan mucho de la soluacion exacta,
podramos decir que sus discontnuidades ocupan regiones mas extensas.
Para tiempos bajos (primeros niveles de t) la solucion tiene un comportamiento
aceptable, pero a medida que t crece la separacion de la solucion exacta es mas
acentuada, especialmente en los esquemas Upwind y Lax-Fredrichs, porque estos
14
3.7.1.
La Condici
on de Courant-Friedrichs-Lewy CFL
El n
umero de Courant se define por Cr = a, donde = hk . El comportamiento de este n
umero es muy importante, porque de el depende, la estabilidad
del esquema de diferencias finitas; ya que una condicion sobre este n
umero da
argumentos sobre la propagacion de la informacion.
La interpretacion de este n
umero se puede expresar en la forma
Cr =
a
1
a
x
t
(3.52)
3.8.
Esquemas Explcitos
3.8.1.
Esquema FTCS
(3.53)
(3.54)
t
.
x
(3.55)
(3.56)
= 1 iCr sen
asi que tenemos
|g|2 = 1 + Cr2 sen2
obviamente, |g|2 1 para todo , asi que el esquema (3.53) incondicionalmente
inestable.
Es posible hacer una analisis de estabilidad definiendo la funcion cuadratica en
funcion de sen.
Debido a la carencia de estabilidad del esquema FTCS, es que no es muy usado
para problema de conveccion pura.
Para tener una idea grafica del comportamienteo del factor de amplificacion
(3.56) efectuamos una parametrizacion en la forma
x=1
y = Cr sen
3.8.2.
Esquema Upwind
(3.57)
(3.58)
(3.59)
Ahora para tener una idea grafica del comportamiento del factor de amplificacion
complejo (3.59) lo parametrizamos en la forma
x = 1 (1 cos )Cr
y = Cr sen
donde [0, 2].
Graficando tenemos una de elipse con centro en (1 Cr, 0) y de radio Cr de
la forma
(x (1 Cr))2
y2
+
=1
Cr2
Cr2
(3.60)
t
1.
x
(3.61)
(3.62)
1.5
region inestabilidad
imaginaria
1.5
imaginario
g=1iCr sen(theta)
Cr=1.5
g=1
0.5
0.5
eje Y
eje y
real
0
real
0
Cr=0.5
region estabilidad
region estabilidad
0.5
Cr=1
0.5
region inestabilidad
1
1
1.5
1.5
1.5
0.5
0
eje x
0.5
1.5
2
3
2.5
1.5
0.5
0.5
1.5
eje X
(3.63)
3.8.3.
Esquemas Leapfrog
Siguiendo el mismo procedimiento como para la ecuacion de difusion buscamos una representacion mas precisa de (3.11) introduciendo diferencias centrales
en el espacio y el tiempo. Esto da el esquema Leapfrog
n
n
vjn+1 vjn1
vj1
vj+1
+a
= 0,
2t
2x
(3.64)
(3.65)
(3.66)
(3.67)
p
g = (1 Cr2 sen2 ) iCr sen,
(3.68)
cuya solucion es
(3.69)
si A = 2 i Cr sen entonces
vjn+1 = Avjn + vjn1
(3.70)
(3.71)
vjn+1
vjn
vjn
A 1
n1
vj
1 0
A 1
G=
1 0
(3.72)
1 Cr2 sen2
(3.73)
Cr2 sen2 1
= i (Cr sen
Cr2 sen2 1)
elevando al cuadrado
|1,2 |2 = 2 Cr2 sen2 2 Cr sen
Cr2 sen2 1 1
(3.74)
(3.75)
o lo que es lo mismo
Cr2 sen2 Cr sen
(3.76)
(3.77)
0 a2
u
= , A =
b1 0
v
Entonces tenemos
+ [A]
=0
t
x
(3.78)
(3.79)
t
(n nj1 )
x j
(3.80)
t
(1 ei )
x
(3.81)
t
|1
x
t
|1
x
p
g = iCr sen 1 Cr2 sen2
|
b1 a2
(3.82)
(3.83)
Cr2 sen2 1)
elevando al cuadrado
|g |2 = 2 Cr2 sen2 2 Cr sen
Cr2 sen2 1 1
p
Cr2 sen2 1 < 1
(3.84)
lo cual nunca es satisfecho puesto que Cr2 sen2 > 1 si arc sen(/2)
arc sen(1/Cr ) y arc sen(1/Cr ) arc sen(/2)). Por lo tanto concluimos que
el esquema Leapfrog es inestable para cualquier n
umero de Courant Cr como se
observa en la Figura (3.14).
x = 1 Cr2 sen2
y = C sen
r
Leapfrog
1.5
7
y
y
6
1
g+
0.5
Cr=1.5
x
0
2
0.5
1
|g| = 1
x
1
1.5
1.5
0.5
0.5
1.5
3.8.4.
1.5
0.5
0.5
Esquema Lax-Wendroff
1.5
sigue
un+1
= unj + t
j
u (t)2 2 u
+
+ O(t)3
t
2! t2
sustituyendo en un+1
tenemos
j
un+1
= unj + (a
j
u
(4t)2 2 2 u
)4t +
(a
) + O(t3 )
x
2
x2
unj
n
unj+1 unj1
unj+1 2unj + vj1
1 2
2
a4t(
) + a (4t) (
)
24x
2
(4x)2
n
n
vj+1
2vjn + vj1
x2
=0
(3.85)
(3.86)
(3.87)
(3.88)
x = 1 2Cr sen2
2
y = C sen
r
1
Cr=1.2
eje y
0.5
|g|= 1
Cr=0.8
Cr=0.5
real
0.5
1.5
2
2.5
1.5
0.5
0.5
1.5
eje x
3.8.5.
Esquema Lax-Friedrichs
(3.89)
y en forma de algoritmo
n
n
n
n
vjn+1 = 0.5(vj+1
+ vj1
) 0.5Cr (vj+1
vj1
).
(3.90)
Haciendo un analisis de consistencia se demuestra que el esquema (3.90) es consistente con (3.11) con un orden de truncamiento O(t) + O(x2 ).
Aplicando el analisis de estabilidad de von Neumann a (3.90) se llega a determinar el factor de amplificacion
g = cos i Cr sen
tal que |g|2 = cos2 + Cr2 sen2 < 1 si Cr 1.
(3.91)
x = cos
y = C sen
r
para todo [0, 2], cuya grafica se muestra en la Figura (3.16) para los n
umeros
de Courant Cr= 0.5, 0.8 y 1.5
3.9.
Esquems Implcitos
3.9.1.
Esquema BTCS
(3.92)
(3.93)
(3.94)
(3.95)
1
1 iCr sen
(3.96)
donde
|g|2 =
1
1 + Cr2 sen
(3.97)
3.9.2.
Esquema Crank-Nicolson
El esquema Crank-Nicolson es muy efectivo cuando lo aplicamos a la ecuacion de difusion unidimensional. En esta seccion daremos el esquema Crank-Nicolson
para la ecuacion de conveccion (3.11).
El esquema de diferencias finitas Crank-Nicolson puede ser escrito
vjn+1 vjn
+ a(0.5Lx vjn + 0.5Lx vjn+1 ) = 0
t
(3.98)
donde Lx vj = (vj+1 vj1 )/(2x). La ecuacion (3.98) produce el siguiente algoritmo tridiagonal
n+1
n+1
n
n
0.25Cr vj1
+ vjn+1 + 0.25Cr vj+1
= 0.25Cr vj1
+ vjn 0.25 Cr vj+1
(3.99)
el cual puede ser solucionado eficientemente usando eliminacion gaussiana y descomposicion LU.
El esquema de Crank-Nicolson (3.99) es consistente con (3.11) con un error de
truncamiento O(t)2 + O(x)2 .
Ahora aplicamos el analisis de estabilidad de von Neumann a (3.11)
0.25 Cr g(ei + ei ) + g = 0.25 Cr (ei ei ) + 1
0.25 g Cr(2isen) + g = 1 0.25 Cr (2isen)
g (i 0.5 Cr sen + 1) = 1 i 0.5 Cr sen
1 i 0.5Cr sen
1 + i 0.5Cr sen
(3.100)
g =
Entonces el modulo de g es
1 A2 2
2A 2
) +(
)
2
1+A
1 + A2
1 + 2A2 + A4
=
(1 + A2 )2
|g|2 = (
= 1
Por lo tanto el esquema de Crank-Nicolson es incondicionalmente estable y por el
Teorema de Lax es convergente de orden (O(t)2 , O(x)2 ).
Factor amplificacion, LaxFriedrichs
2
imaginario
1.5
1.5
Cr=1.5
imaginario
1
Cr=0.8
Cr=0.5
0.5
Cr=1.2
0.5
real
eje y
eje y
|g|=1.0
0
real
0
region estabilidad
0.5
0.5
Cr=0.5
|g|=1.0
1
region inestabilidad
1.5
1.5
2
1.5
0.5
0
eje x
0.5
1.5
2
2
1.5
0.5
0
eje x
0.5
1.5
La Figura (3.17) muestra los valores del factor de amplificacion para Cr = 1.0
en la que se observa que todos los valores son iguales a 1. Resultados similares se
obtienen para diferentes n
umeros de Courant mayores o menores que uno.
3.10.
Resultados Num
ericos
En esta seccion presentamos los esquemas de diferencias finitas para el problema de conveccion unidimensional analizadas en su estabilidad y consistencia
en la seccion anterior. De igual forma al problema de difusion hacemos un estudio comparativo de los esquemas resolviendo el mismo problema bajo condiciones
similares y dependencia de los parametros.
Ejemplo 28 Consideremos el problema siguiente
ut + aux = 0,
< x < ,
t>0
con condici
on inicial
u(x, 0) = f (x) =
1 , 0.05 x 0.15
0 , en otro lugar.
(3.101)
donde la constante a sera igual a 1. El rango para la variable x sera 1 < x < 2
y el espaciamiento en la malla para la variable x es x = 0.01.
Las Figuras (3.18), (3.19), (3.20) y (3.21) describen los resultados para los pasos n = 50, 100 y 150, y para Cr = 0.5 para los esquemas Lax-Wendroff, Upwind,
Leapfrog, Lax-Friedrichs donde las soluciones aproximadas son representantes por
lneas punteadas y la solucion exacta es representada por lneas continuas. En la
Figura (3.18) se muestra que el esquema de segundo orden dispersivo Lax-Wendroff
tiene oscilaciones cerca de las discontinuidades pero no se propagan inmediatamente. En la Figura (3.19) se muestra el esquema de primer orden disipativo Upwind
donde la curva de aproximacion se ajusta mas a la curva exacta, no presenta oscilaciones de ning
un tipo, cerca de las discontinuidades se separa mucho de la solucion
exacta. En la Figura (3.20) se muestra que el esquema de segundo orden dispersivo
tiene oscilaciones cerca de las discontinuidades que se propagan muy rapido en los
primeros niveles de tiempo.
tiempo = 0.250000
tiempo = 0.500000
0.8
0.8
0.6
0.6
u(x,t)
1.2
u(x,t)
1.2
0.4
0.4
0.2
0.2
0.2
0.2
0.1
0.2
0.3
0.4
0.5
eje x
0.6
0.7
0.8
0.9
0.1
0.2
0.3
0.4
0.5
eje x
0.6
tiempo = 0.750000
1.2
0.8
0.6
u(x,t)
0.4
0.2
0.2
0
0.1
0.2
0.3
0.4
0.5
eje x
0.6
0.7
0.8
0.9
0.7
0.8
0.9
tiempo = 0.250000
tiempo = 0.500000
0.8
0.8
0.6
0.6
u(x,t)
u(x,t)
0.4
0.4
0.2
0.2
0.1
0.2
0.3
0.4
0.5
eje x
0.6
0.7
0.8
0.9
0.1
0.2
0.3
0.4
0.5
eje x
tiempo = 0.750000
0.8
0.6
u(x,t)
0.4
0.2
0.1
0.2
0.3
0.4
0.5
eje x
0.6
0.7
0.8
0.9
0.6
0.7
0.8
0.9
tiempo = 0.250000
tiempo = 0.500000
1.4
1.2
1.2
1
1
0.8
eje u(x,t)
eje u(x,t)
0.8
0.6
0.6
0.4
0.4
0.2
0.2
0.2
0.2
0
0.1
0.2
0.3
0.4
0.5
eje x
0.6
0.7
0.8
0.9
0.1
0.2
0.3
0.4
0.5
eje x
tiempo = 0.750000
1.4
1.2
0.8
0.6
0.4
0.2
0.2
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
0.6
0.7
0.8
0.9
tiempo=0.250000
tiempo=0.500000
0.8
0.8
0.6
0.6
u(x,t)
u(x,t)
0.4
0.4
0.2
0.2
0.5
0.5
eje x
1.5
0.5
0.5
eje x
1.5
tiempo=0.750000
0.8
u(x,t)
0.6
0.4
0.2
0.5
0.5
eje x
1.5
sen(10x) ,
u(x, 0) =
0
,
0 x 0.1
0.1 < x 1.0
t > 0.
0
, 0 x at
u(x, t) =
sen(10(x at)) , at x at + 0.1
0
, at + 0.1 x 1.0
Para un valor a = 0.1 la solucion exacta en t = 0 y t = 8 son mostradas en
la Figura (3.22). La onda sinusoidal se propaga sin reduccir su amplitud a una
velocidad a = 0.1.
Soluciones computacionales fueron obtenidas con N = 40 para x en el intervalo
0 x 1.0 y con un n
umero de Courant Cr = 0.8. Utilizamos los esquemas de
primer orden disipativo Upwind: Figura (3.22), el esquema dispersivo de segundo
orden explcito Lax-Wendroff: (3.23), el esquema Lax-Friedrichs: (3.24) y el esquema FTCS: (3.25). La solucion computacional es suave, pero comparada con la
solucion exacta, tiene un ablandamiento (o difusividad); este es consistente con
la introduccion del termino de viscosidad artificial.
La solucion para el problema de Lax-Wendroff (3.85) es mostrada en la Figura
(3.23) para los valores de t = 8.0 y Cr = 0.8.
tiempo = 8.000000
tiempo = 8.200000
1.2
1.2
exacta t=0.0
0.8
exacta t=8.0
0.8
0.6
exacta t=0.0
u(x,t)
u(x,t)
0.6
exacta t=8.0
0.4
0.4
0.2
0.2
Upwind t=8.0
0.2
0.2
0
0.1
0.2
0.3
0.4
0.5
0.6
N=40, dx= 0.025, dt = 0.2
0.7
0.8
0.9
0.1
0.2
0.3
0.4
0.5
0.6
N=40, dx=0.025, dt=0.2
0.7
0.8
0.9
tiempo = 8.000000
tiempo = 3.000000
8
1.2
6
1
FTCS t=3.0
4
0.8
exacta t=0.0
2
exacta t=8.0
u(x,t)
u(x,t)
0.6
exacta t=8.0
exacta t=0.0
0
0.4
2
0.2
4
Lax Friedrichs t=8.0
0
6
0.2
0
0.1
0.2
0.3
0.4
0.5
0.6
N=40, dx=0.025, dt=0.2
0.7
0.8
0.9
0.1
0.2
0.3
0.4
0.5
0.6
N=40, dx=0.025 , dt=0.2
0.7
0.8
u(x, 0) =
1 ,
1 x 0
0 ,
0 < x 2.0
0.9
tiempo = 0.250000
tiempo = 0.500000
0.8
0.8
0.6
0.6
u(x,t)
1.2
u(x,t)
1.2
0.4
0.4
0.2
0.2
0.2
0.2
0
0.5
1.5
0.2
0.2
0.4
0.6
eje x
tiempo = 0.750000
1.2
0.8
0.6
0.4
0.2
0.2
0.2
0.8
1
eje x
u(x,t)
0.5
0.2
0.4
0.6
0.8
1.2
1.4
1.6
1.8
eje x
1.2
1.4
1.6
1.8
tiempo = 0.250000
tiempo = 0.500000
1.4
1.4
1.2
1.2
0.8
0.8
0.6
0.6
0.4
0.4
0.2
0.2
0.2
0.5
0.5
1.5
0.2
0.2
0.2
0.4
0.6
0.8
tiempo = 0.750000
1.4
1.2
0.8
0.6
0.4
0.2
0.2
0.2
0.2
0.4
0.6
0.8
1.2
1.4
1.6
1.8
1.2
1.4
1.6
1.8
tiempo = 0.250000
tiempo = 0.500000
0.8
0.8
0.6
0.6
u(x,t)
u(x,t)
0.4
0.4
0.2
0.2
0.2
0.2
0.4
0.6
0.8
1.2
1.4
1.6
1.8
0.2
0.2
0.4
0.6
0.8
eje x
1.2
1.4
1.6
eje x
tiempo = 0.750000
0.8
u(x,t)
0.6
0.4
0.2
0.2
0.2
0.4
0.6
0.8
1.2
1.4
1.6
1.8
eje x
u(x, 0) =
1 ,
0x1
0 ,
otro lugar
(3.102)
1.8
0.5
0.5
u(x,t)
u(x,t)
tiempo = 0.200000
1.5
0.5
0.5
0.5
0.5
1
1.5
2
2.5
FTCS, dx = 0.1, dt = 0.05, Tmax = 0.2 seg
3.5
0.5
0.5
1
1.5
2
2.5
FTCS, dx = 0.1, dt = 0.05, Tmax = 0.4 seg
3.5
Figura 3.29: Conveccion FTCS, Tmax=0.2 Figura 3.30: Conveccion FTCS, Tmax=0.4
tiempo = 2.000000
u(x,t)
0.8
0.6
0.4
0.2
0
1
0.5
0.5
1
1.5
2
2.5
Upwind, dx= 0.001, dt= 0.005, Tmax= 2.0 seg
3.5
Leapfrog
Leapfrog
1.4
1.4
t=0.25 seg
1.2
1.2
0.8
0.8
0.6
0.6
t = 0.375
u(x,t)
u(x,t)
0.4
0.4
0.2
0.2
0.2
0.2
0.4
1
0.5
0.5
1.5
2
Tmax = 0.25 seg
2.5
3.5
0.4
1
0.5
0.5
1.5
2
Tmax = 0.375 seg
2.5
3.5
Leapfrog
Leapfrog
1.4
1.4
t=1.0
1.2
t=2.0
1.2
t=0
0.8
0.8
0.6
0.6
u(x,t)
u(x,t)
0.4
0.4
0.2
0.2
0.2
0.2
0.4
1
0.5
0.5
1.5
2
Tmax = 1.0 seg
2.5
3.5
0.4
1
0.5
0.5
1.5
Tmax=2.0 seg
2.5
3.5
tiempo = 2.000000
u(x,t)
0.8
0.6
0.4
0.2
0
1
0.5
0.5
1
1.5
2
2.5
Laxfriedrics, dx = 0.01, dt = 0.05, Tmax = 2 seg
3.5
ut + aux = 0,
(3.103)
u(x, 0) =
y u(1, t) = 0.
, en otro lugar.
(3.104)
Upwind
1
t=0
0.8
Eje U
Cr=0.8
0.6
t = 2.4
0.4
0.2
0
1
0.5
0.5
1
Eje X,
1.5
2.5
h=0.1
Upwind
Upwind
Cr=0.8
Cr= 0.8
t=0
t=0
1
0.8
t = 2.4
0.8
0.6
Eje U
Eje U
t = 2.4
0.6
0.4
0.4
0.2
0
1
0.2
0.5
0.5
1
Eje X,
1.5
h=0.05
2.5
0
1
0.5
0.5
1
Eje X,
1.5
h=0.025
2.5
LaxFriedrichs
LaxFriedrichs
t=0
Cr = 0.8
t=0
Cr = 0.8
0.8
0.8
Eje U
0.6
t = 2.4
0.6
0.4
0.4
0.2
0.2
0
1
0.5
0.5
1
1.5
2
Eje X,
h = 0.1
2.5
3.5
0
1
0.5
0.5
1
Eje X,
t=0
Cr = 0.8
t = 2.4
0.8
0.6
0.4
0.2
0
1
0.5
0.5
1
1.5
2
Eje X,
h=0.025
1.5
2
h = 0.05
2.5
3.5
LaxFriedrichs
Eje U
Eje U
t = 2.4
2.5
3.5
LaxWendroff
1.2
Cr = 0.8
0.8
t = 2.4
Eje U
0.6
0.4
0.2
0.2
1
0.5
0.5
1
Eje X,
1.5
2.5
h = 0.1
LaxWendroff
LaxWendroff
1.2
1.2
1
Cr = 0.8
Cr = 0.8
t = 2.4
0.8
0.8
t = 2.4
Eje U
0.6
Eje U
0.6
0.4
0.4
0.2
0.2
0.2
1
0.5
0.5
Eje X,
1
h = 0.05
1.5
2.5
0.2
1
0.5
0.5
Eje X,
1
h=0.025
1.5
2.5
Leapfrog
1.2
1
Cr = 0.8
0.8
t = 2.4
eje U
0.6
0.4
0.2
0.2
1
0.5
0.5
1
Eje X, h =0.1
1.5
2.5
Leapfrog
Leapfrog
1.2
1.2
1
Cr = 0.8
Cr = 0.8
0.8
t = 2.4
0.8
t = 2.4
Eje U
0.6
Eje U
0.6
0.4
0.4
0.2
0.2
0.2
1
0.5
0.5
1
Eje X, h = 0.05
1.5
2.5
0.2
1
0.5
0.5
Eje X,
1
1.5
h = 0.025
2.5
Captulo 4
Ecuaci
ones Diferencales
Parab
olicas
En este captulo la ecuacion de difusion unidimensional sera usada como
un vehculo para desarrollar esquemas explcitos e implcitos.
Centraremos nuestra atencion en la estabilidad y consistencia de varios esquemas tanto explcitos como implcitos.
La forma canonica de la ecuacion parabolica de difusion unidimensional o ecuacion del calor esta dada por
U
2U
=K
,
T
X 2
0 X L,
T >0
(4.1)
U
T
X
, u=
, t=
L
U0
T0
181
(4.2)
(4.3)
puesto que
U
T
U dt
(U0 u) 1
=
t dT
t T
U0 u
=
T0 t
=
y
U0 2 u
2U
=
X 2
L2 x2
(4.4)
L2 u
2u
=
KT0 t
x2
Por consiguiente si elegimos
T0 =
L2
K
la ecuacion es
u
2u
=
donde 0 x 1
t
x2
(4.5)
Esta ecuacion adimensional puede ser usada para tipificar diversos procesos numericos y se puede deducir siguiendo las leyes fsicas siguientes
1. La cantidad de calor necesario para elevar la temperatura de un objeto de
masa m en una cantidad 4u es
m s 4u
2. La cantidad de calor que fluye a traves de una area por unidad de tiempo y
proporcional a la tasa de cambio de la temperatura con respecto a la distancia
perpendicular al area (normal)
Q = k A 4t
u
x
(4.6)
La cantidad de calor
u
= ku.
(1)
es la normal a la superficie S
* La Ley de Newton establece que la cantidad de calor que fluye durante
una unidad de tiempo a traves del area de la superficie del cuerpo al medio
(2)
u
x
< 0 esto
(4.7)
u
en x + 4x
x
(4.8)
QB = k A 4t
y
QC = k A 4t
QB QC = {k A 4t
u
u
|x+4x
|x } = m s 4u
x
x
= A 4x s 4u
u
|
x x+4x
4x
u
|
x x
=s
4u
4t
(4.9)
2u
u
=s
2
x
t
(4.10)
u
2u
= 2
t
x
(4.11)
donde =
k
s
t 6
]0, L[R+
valores de
valores de
frontera
frontera
valores iniciales
t 6
valores de
frontera tn
valores de
frontera
(j, n)
xj
valores iniciales
En los extremos podemos tener condiciones de frontera de Dirichlet o de Neumann ( Ver Figura 4.3).
u0(x)
u0 (x)=f(x)
A
uA=g (t)
B
aislamiento
uB=g (t)
2
uB =h (t)
uA x=h (t)
x=0
x=L
4.1.
Soluci
on de la Ecuaci
on del Calor
Construyamos la solucion exacta de la ecuacion del calor utilizando el Meto-
En efecto, sea
u(x, t) = T (t)X(x)
(4.12)
(4.13)
Tt (t)
Xxx (x)
=
=
T (t)
X(x)
(4.14)
Asi se obtiene
(4.15)
Xxx = X
(4.16)
X(0) = X(1) = 0
De donde concluimos que X = sen( x) y de las condiciones de contorno obtenemos que = k 2 2 , k Z. Asi
X(x) = sen(kx)
T (t) = T (0)ek
2 2 t
(4.17)
2 2 t
(4.18)
(4.19)
Como T (0), puede tomar cualquier valor, podemos asumir que depende tambien
de k, por tanto las soluciones particulares estan dadas por
uk (x, 0) = Tk (0)sen(kx)
(4.20)
Como la ecuacion (4.1) es lineal tendremos que las conbinaciones lineales de las
soluciones particulares (4.20) tambien satisfacen la ecuacion (4.1) y las condiciones
de contorno. Denotando por
m
u =
m
X
Tk (0)ek
2 2 t
sen(kx)
(4.21)
k=1
u (x, 0) =
m
X
Tk (0)sen(kx)
(4.22)
k=1
Se prueba en analisis que toda funcion C 1 (0, 1) que se anula en los extremos puede
ser expresada como una Serie de Fourier de la forma
u0 (x) =
Z
ck sen(kx),
ck =
u0 (x)sen(kx)dx
(4.23)
k=1
ck ek
2 2 t
sen(kx)
(3)
k=1
Formulaci
on matem
atica :
u
2u
= 0,16 2
t
x
(4.24)
u(0, t) = u(100, t) = 0
(4.25)
donde
f (x) =
60 ,
0 x 50
40 , 50 x 100
(4.26)
(4.27)
00
X.T = 0,16X T
0
00
T
X
=
= 2
0,16T
X
donde
0
T + 0,162 T = 0,
00
X + 2 X = 0
T = C1 e016 t ,
X = C2 cos x + C3 senx
Luego la solucion es
u = X.T
2
(4.28)
n
.
100
Entonces la
2 2
n
0,16 n100
x)
e
1002
(4.29)
n2 2
n=1
n
x)
100
(4.30)
bn sen(
n=1
n
x) = f (x)
100
Z 100
2
n
=
f (x)sen(
x) dx
100 0
100
Z 50
Z 100
n
n
2
2
f (x)sen(
f (x)sen(
x) dx +
x) dx
=
100 0
100
100 50
100
120
n
80
n
=
(1 cos
)+
(cos
cos n)
n
2
n
2
entonces tenemos
b1 =
200
,
b2 =
40
,
b3 =
200
, ...
3
(4.31)
200 16,102 t
40
2
2
e
sen
x + e64,10 t sen
x + ...
100
100
(4.32)
radiaci
on puede ocurrir hacia los alrededores. Asumiremos que la Ley de Newton
de enfriamiento es aplicable, la ecuaci
on que describe esto es
u
2u
= 2 C(u u0 )
t
x
(4.33)
(4.34)
Matem
aticamente las condiciones de aislamiento en los extremos de la barra colocados en x = 0 y x = L del eje X se expresan por
ux (0, t) = ux (L, t) = 0
(4.35)
Adem
as la temperatura inicial la indicamos por
u(x, 0) = u0 (x),
0<x<L
(4.36)
(4.37)
ux = e t (Asenx + B cos x)
evaluando en el extremo izquierdo x = 0
2
ux (0, t) = Be t = 0
B=0
Asi tenemos la solucion
2
u(x, t) = Ae t cos x
Derivamos (4.38) respecto a x
2
ux = Ae t senx
evaluando en el extremo derecho x = L
2
ux (L, t) = Ae t senL
senL = 0
L = n
L =
n
,
nZ
(4.38)
Entonces la solucion es
u(x, t) = Ae
n2 2 t
L2
(cos
n
x)
L
Para satisfacer la condicion inicial (4.36) debemos usar el Principio de Superposicion de Soluciones y asi llegamos a la solucion final
a0 X
n
2 2
2
u(x, t) =
+
an en t/L cos
x
2
L
n=1
entonces
a0 X
n
u(x, 0) =
+
an cos
x = u0 (x)
2
L
n=1
donde
2
an =
L
u0 (x) cos
0
n
x dx
L
Ejemplo 36 Conducci
on de calor en una barra con condiciones no cero
en los extremos
En el problema de Fourier: desde el punto de vista fsico no hay raz
on de no
haber escogido otras dos condiciones, tales como en el extremo x = 0 se mantienen
a 200 y el extremo en x = L a 600 .
Formulamos el problema de valor inicial
u
2u
= 2,
t
x
u(0, t) = 20,
(4.39)
u(L, t) = 60
(4.40)
u(x, 0) = 100
(4.41)
(4.42)
(4.43)
(4.44)
(4.45)
(4.46)
(4.47)
(x) = x +
(4.48)
(4.49)
Tambien sera bueno que los lados derechos de (4.45) y (4.46) sean ceros. Esto se
consigue si
(0) = 20,
Esto nos permite determinar que =
(L) = 60
40
L
(x) =
w(0, t) = 0,
w(L, t) = 0
(x) =
40
+ 20
L
(4.50)
n
, n = 1, 2, ...
L
2 2 t/L2
w(x, t) = Ben
sen
n
x
L
X
n=1
bn sen
n
x
L
(4.51)
n
x dx
L
bn en
2 2 t/L2
n=1
(sen
n
40
x) + ( x + 20)
L
L
(4.52)
(4.53)
u(0, y) = u(L, y) = 0
u(x, 0) = u0
Puesto que fiscamente es imposible que la temperatura llegue a ser infinita en todo
punto de la placa, asumiremos que u(x, y) es acotada, es decir,
|u(x, y)| < M,
(x, y)
frontera acotada
frontera aislada
frontera aislada
u(L,y)=0
u ( 0, y ) = 0
X
L
(4.54)
00
X Y + XY = 0
00
00
X
Y
=
= 2
X
Y
00
00
X
Y
2
= ,
= 2
X
Y
Resolviendo estas ecuaciones ordinarias obtenemos
X = C1 cos x + C2 senx,
Y = C3 ey + C4 ey
Entonces la solucion es
u = (C1 cos x + C2 senx)(C3 ey + C4 ey )
Asumiendo que > 0, esta solucion llega a ser no acotada cuando y . Para
evitar esto, escogemos C3 = 0,
u(x, y) = ey (A cos x + Bsenx)
Usamos la condicion de frontera izquierda u(0, y) = 0,
u(0, y) = Aey = 0
de donde se obtiene que A = 0. Entonces la solucion sera
u = Bey senx
Ahora usamos la condicion de frontera derecha u(L, y) = 0,
u(L, y) = Bey senL = 0
de donde se tiene que senL = 0 y =
n
,
L
nZ
u = Be L y sen
n
x
L
bn e L y sen
n=1
n
x
L
X
n=1
bn sen
n
x = u0
L
Z
n
2 L
u0 sen x dx
=
L 0
L
2u0 (1 cos n)
=
n
(4.55)
(4.56)
sen L x
2u0
tan1 [
]
senh L y
= a2 xu2 ,
u(0, t)
= u(L, t) = 0
u(x, 0)
= f (x),
(4.57)
0<x<L
ut (x, 0) = 0
0<x<L
00
= a2 X .T
00
T
= 2 = 2
aT
X
X
00
00
(4.58)
de donde
00
00
X + 2 X = 0,
T + a 2 2 T = 0
T = C3 cos at + C4 senat
Entonces la solucion es
u = (C1 cos x + C2 senx)(C3 cos at + C4 senat)
(4.59)
n
,
L
nZ
n
na
na
x (C3 cos
t + C4 sen
t)
L
L
L
(4.60)
n
na
na
na
na
x (C3
sen
t + C4
cos
t)
L
L
L
L
L
n
na
x. cos
t
L
L
(4.61)
bn sen
n=1
n
na
x. cos
t
L
L
Para el tiempo t = 0
u(x, 0) = u0 =
bn sen
n=1
n
x
L
u0 sen
0
n
x dx
L
Finalmente la solucion es
2X
u(x, t) =
(
L n=1
u0 sen
0
n
n
na
x dx) sen x. cos
t
L
L
L
(4.62)
=
a
t2
t
x2
(4.63)
0<x<L
u
(x, 0) = 0
t
(4.64)
(4.65)
(4.66)
(4.67)
00
X.T + X.T = a2 X .T
00
00
X
T
T
= 2 + 2 = 2
X
aT
aT
De donde obtenemos las ecuaciones diferenciales
00
00
X + 2 X = 0,
T + T + 2 a2 T = 0
X = C1 cos x + C2 senx;
donde w =
2 /4 2 a2
(4.68)
Incorporado la condicion de frontera izquierda u(0, t) = 0 a (4.68) se llega a determinar que C1 = 0. Del mismo modo con la condicion de frontera derecha u(L, t) = 0
llegamos a determinar que =
n
,
L
n Z.
n
x (C3 cos wt + C4 senwt)
L
(4.69)
La u
ltima condicion de frontera (4.66) nos da que C4 = 0. Mientras que para
satisfacer la condicion inicial (4.65) aplicamos el Prinicipio de Superposicion de
soluciones
u(x, t) = e
q
donde =
t/2
bn (sen
n=1
n2 2 a2
L2
2
.
4
n
x) cos
L
(4.70)
conduce a
u(x, 0) = f (x) =
bn =
4.2.
2
L
bn sen
n=1
f (x)sen
0
n
x dx
L
n
x dx
L
Principio del M
aximo, Decaimiento
Exponencial, Efecto Regularizante
Damos dos de las propiedades mas importantes de la ecuacion del calor, co-
(4.71)
de estas dos u
ltimas proposiciones tenemos que si
| x y |<
entonces
|
<
+ + =
3 3 3
fn (x)dx =
a
lm fn (x)dx =
a n
f (x)dx
(4.72)
, x [a, b]
ba
(4.73)
de donde obtenemos
Z
fn (x)dx
a
f (x)dx |
a
este resultado vale para toda sucecion de funciones integrables que converge uniforme.
4.2.1.
Principio del M
aximo
ut u = 0,
en R+
u |
= g
u(x, 0) = u (x),
en
0
(4.74)
Entonces, se cumple
n
M in inf g, inf u0
(4.75)
Prueba: Asumiremos que la funcion u0 y g son continuas y que existe una solucion
continua u dos veces derivable. Si uno de los extremos de la funcion u esta en la
frontera, entonces no tenemos nada a probar.
(4.76)
vt 0
y
vxx 0
=
vx
0
Por tanto, en el punto (x0 , t0 ) del interior del dominio, se cumple
vt (x0 , t0 ) v(x0 , t0 ) 0
Esta expresion contradice la desigualdad (4.76).
Por consiguiente, nuestra hipotesis de que el maximo de v esta en el interior del
dominio no es cierta. Por tanto el maximo de v debe de encontrarse en la frontera
del dominio, {0} R+ y se debe cumplir
2
2
= max sup (u(x, 0) + | x | ) , sup (u(x, t) + | x | )
2
2
= max sup u(x, 0) + L , sup u(x, t) + L
concluimos que
2
2
u(x, t) v(x, t) max sup u(x, 0) + L , sup u(x, t) + L
haciendo 0 obtenemos
(4.77)
4.2.2.
Decaimiento de Soluciones
E(t) :=
(4.78)
u =
(ux u) u2x
2 t
x
integrando sobre ]0, L[ obtenemos
d
E(t) = 2
dt
ux (x, t)2 dx
u dx
0
L
0
u2x dx
obtenemos
d
2
E(t) 2 E(t)
dt
L
resolviendo la desigualdad obtenemos
2
E(t) E(0)e L2 t
de donde sigue el teorema.
4.2.3.
Efecto Regularizante
Esta propiedad garantiza que, dependiendo de la condicion inicial, la solucion del problema de conduccion del calor puede ser suficientemente regular, se
resume en el siguiente teorema
Teorema 4.2.5 Sea u0 una funcion continua, para la cual exista una solucion u
para la ecuaci
on del calor
ut uxx = 0
en
]0, L[ ]0, [
(4.79)
(4.80)
4.3.
Esquemas Explcitos
Un esquema es explcito si la u
nica variable desconocida, por ejemplo vjn+1 ,
4.3.1.
Esquema FTCS
(4.81)
t 6
tn+1
(j, n + 1)
w
tn
(j 1, n)
(j, n)
(j + 1, n)
xj1
xj
xj+1
(4.82)
(4.83)
donde
g() = rei + (1 2r) + rei
(4.84)
(4.85)
tal que para cualquier valor de se tiene |g()| 1 si r 0.5. As (4.82) producira soluciones estables si r 0.5. Aplicando el Teorema de Lax tenemos que el
esquema es convergente de orden [O(t), O(x)2 ].
factor amplificcion ,FTCS
1.5
90
region inestabilidad
120
60
0.5
150
0.5
30
G(THETA)
r=0.2
0
r=0.3
r=0.2
180
0.5
region estabilidad
r=0.5
210
330
r=0.7
1.5
region inestabilidad
240
300
270
3
THETA
120
60
0.5
150
120
30
180
60
0.5
150
30
180
r=0.5
210
330
210
330
r=0.7
240
300
240
300
270
270
4.3.2.
Esquema Leapfrog
(4.86)
n+1
n
n1
j1
j+1
(4.87)
1 + 16r2 sen4
(4.88)
el cual es |g| > 1 para todo y todo r. Esto demuestra que el esquema es incondicionalmente inestable para todo r > 0. Por lo tanto, por el Teorema de Lax, no
puede ser convergente.
90
120
60
region estabilidad
0.5
150
30
G(THETA)
1
r=0.2
2
180
r=0.5
4
region inestabilidad
5
210
330
r=1.0
3
THETA
240
300
270
Coordenadas Polares, r = 1.0
Figura
4.11:
Leapfrog,
Polares
=
Figura 4.12: Leapfrog, Polares r = 1.0
120
60
0.5
150
30
180
330
240
60
0.5
150
210
120
300
30
180
210
330
240
300
270
270
La Figura (4.11) muestra que |g()| > 1 para cualquier r, en particular para
r = 0.2, 0.5, 1.0, vemos que los valores del factor de amplificacion caen fuera de la
region de estabilidad.
La Figura (4.12), (4.13) y (4.14) muestran el factor de amplificacion en coordenadas polares para radio = 0.2, 0.3, 1.0 donde se observa que los valores salen
fuera de la region de estabilidad.
4.3.3.
Esquema Du Fort-Frankel
El esquema Leapfrog (4.86) puede ser modificado para producir un algoritmo estable. Esto se puede lograr reemplazando vjn en (4.86) por 0,5(vjn+1 +vjn1 ).
El resultado es la ecuacion
n
n
vj+1
(vjn+1 + vjn1 ) + vj1
vjn+1 vjn1
=
.
2t
(x)2
(4.89)
t 6
tn+1
(j, n + 1)
w
tn
w
w
(j, n 1)
xj1
xj
xj+1
vjn+1
2r
1 + 2r
n
(vj1
n
vj+1
)
1 2r
1 + 2r
vjn1
(4.90)
tiempo. Este esquema esta dentro de los esquemas de paso multiple, pues usa una
diferencia central para aproximar la derivada temporal en el nodo (j, n), asi como
un promedio espacial en la parte difusa.
Igual que antes, consideremos una funcion diferenciable y utilizamos la serie
de Taylor en dimension uno.
El operador diferencial discreto Lnj aplicado a es dado por
Lnj
n+1
n1
nj+1 (n+1
+ jn1 ) + nj1
j
j
j
=
2t
x2
(4.91)
1 2
1
t tt (x, t) + t3 t (x, t) + O(t4 )
2!
3!
1 2
1
t tt + t3 ttt + O(t4 )
2!
3!
(4.92)
1 2
1
t tt (x, t) t3 ttt (x, t) + O(t4 )
2!
3!
1 2
1
t tt t3 ttt + O(t4 )
2!
3!
(4.93)
1
1
x2 xx (x, t) + x3 xxx (x, t) +
2!
3!
1
x4 xxxx (x, t) + O(x5 )
4!
y en la notacion de ndice
nj+1 = nj + xx +
1
1
1
x2 xx + x3 xxx + x4 xxxx (x, t) + O(x5 ) (4.94)
2!
3!
4!
1
1
x2 xx (x, t) x3 xxx (x, t) +
2!
3!
1
x4 xxxx (x, t) + O(x5 )
4!
nj1 = nj xx +
1
1
1
x2 xx x3 xxx + x4 xxxx + O(x5 )
2!
3!
4!
(4.95)
Sustituyendo estas espresiones (4.92), (4.93), (4.94) y (4.95) en las diferencias del
operador discreto Lnj dado en (4.91), teniendo en cuenta que derivadas t , x
, xx , xxx , xxxx , son evaluadas en el punto (xj , tn ) = (jx, nt), tenemos:
De (4.92) y (4.93) se tiene
n+1
n1
t2
j
j
= t +
ttt + O(t4 ) = t + O(t2 )
2t
3!
(4.96)
n+1
+ n1
= 2 + t2 tt + O(t4 )
j
j
(4.97)
(4.98)
Sustraendo (4.97) de (4.98), multiplicando esta difertencia por /x2 y adicionando con (4.96) obtenemos
Lnj = t xx + t2 x2 tt + O(t2 ) + O(x2 ) + O(t2 t2 )
y de ello se cumple
Lnj Lnj = t2 x2 tt + O(t2 ) + O(x2 ) + O(t2 x2 ) (4.99)
t
x
0,
cuando t 0,
x 0
Aplicamos el criterio de von Neumann para determinar la estabilidad del esquema (4.89). Para ello escribimos el esquema en la forma
n
n
+ vj1
+ (1 2r) vjn1
(1 + 2r) vjn+1 = 2r vj+1
(4.100)
(1 + 2r)g = 2r ei + ei + (1 2r)g 1
multiplicando por g a este ecuacion
(1 + 2r)g 2 2r ei + ei g (1 2r) = 0
(4.101)
4r(cos )
p
16r2 cos2 + 4(1 4r2 )
2(1 + 2r)
de esto
g =
2r cos
(4.102)
1. Si se cumple que
1 4r2 Sen2 0
(4.103)
es decir
0 1 4r2 sen2 1
por consiguiente
=1
1 + 2r
1 + 2r
Concluyendo, podemos decir que si se cumpliera (4.103), la condicion de von
Neumann es satisfecha
| g | 1
y el esquema de Du Fort-Frankel es estable.
2. Si cumpliera
1 4r2 Sen2 = b2 < 0
(4.104)
entonces
1 4r2 sen2 4r2 ,
es decir
2r cos ib
1 + 2r
(4.105)
de donde
| g |2 =
4r2 1
(2r + 1)(2r 1)
2r 1
=
=
<1
2
2
(1 + 2r)
(1 + 2r)
(1 + 2r)
en la u
ltima desigualdad utilizamos (4.105). Concluimos tambien que bajo
la restriccion (4.104), el esquema es estable.
Por tanto esquema de Du Fort-Frankel es estable sin restricciones, y por el
Teorema de Lax, el esquema es convergente, con la u
nica condicion exigida por la
consistencia, donde observamos que t vaya mas rapido a cero que el valor de x.
Para graficar el factor de amplificacion (4.102) para le caso r > 0.5 lo parametrizamos en la forma
x=
2r cos
1+2r
y=
4r2 sen2 1
1+2r
donde
4r2 cos2 4r2 sen2 1
+
(1 + 2r)2
(1 + 2r)2
4r2 (cos2 + sen2 ) 1
=
(1 + 2r)2
2r 1
=
1 + 2r
x2 + y 2 =
tenemos
x2 + y 2 =
2r 1
1 + 2r
(4.106)
q
2r1
1+2r
1
1
r=3
0.8
0.6
r=1
0.5
0.4
r=0.3
0.2
0
eje Y
G(THETA)
r=0.2
r=0.5
0
0.2
0.5
0.4
region estabilidad
0.6
1
0.8
1
1.5
3
2
0.8
0.6
0.4
0.2
0
eje X
0.2
0.4
0.6
0.8
1
0.8
r=1.0
r=0.2
0.6
r=0.5
r=0.3
0.4
eje Y
0.2
0
0.2
0.4
0.6
0.8
1
0.8
0.6
0.4
0.2
0
eje X
0.2
0.4
0.6
0.8
polares para radio = 0.2; 0.3; 0.5; 1.0 que se encuentra siempre dentro de la
circuferencia unitaria centrada en el origen(region de estabilidad).
Los esquemas de diferencias finitas dados hasta el momento tienen el inconveniente de ser condicionalmente estables. A continuacion presentamos una familia
de esquema que son incondicionalmente estables.
4.4.
Esquemas Implcitos
2u
Para esquemas implcitos el termino espacial
en (4.1) es evaluada, al
x2
4.4.1.
(4.107)
t 6
tn
(j 1, n)
(j, n)
(j + 1, n)
tn1
(j, n 1)
xj1
xj
xj+1
(4.108)
donde j = 0, 1, . . . , N , n = 0, 1, . . . , kmax y r = tx2 Verifiquemos la consistencia del esquema de Euler retrazado, para ello utilizamos la serie de Taylor. En
efecto: Observemos que el operador diferencial continuo esta dado por:
L=
2
2,
t
x
nj n1
nj+1 2nj + nj1
j
t
2
(4.109)
(4.110)
1
x2 xx (xj , tn )
2!
1
x3 xxx (xj , tn ) + O(x4 )
3!
En notacion indicial
nj1 = nj xx +
1
1
x2 xx x3 xxx + O(x4 )
2!
3!
(4.111)
1
1
x2 xx + x3 xxx + O(x4 )
2!
3!
(4.112)
1 n
j nj + tt + O(t2 )
t
1
1
n
2
3
4
n
j + xx + x xx + x xxx + O(x ) 2j
x2
2!
3!
1
1
n
2
3
4
j xx + x xx x xxx + O(x )
x2
2!
3!
Lnj =
1
2
tt + O(t4 )
x xx + O(x4 ) + O(x4 )
2
t
x
(4.113)
Observando que
Lnj Lnj 0,
cuando t, x 0
1
.
1 + 4rsen2 2
(4.114)
Analizando
0 sen2
1 entonces
2
1 1 + 4rsen2
1 + 4r
2
1
1
1
1 + 4r
1 + 4rsen2 2
es decir
0 < g 1 < 1 + t
Esta u
ltima desigualdad es la condicion de estabilidad de von Neumann, por
consiguiente el esquema de Euler retrasado es estable sin restricciones. Asi, por el
Teorema de Lax, el esquema es convergente.
120
60
r= 1.0
0.5
150
r= 0.3
r= 0.5
0.5
30
G(THETA)
r= 2.0
0
180
region estabilidad
0.5
210
330
240
1.5
3
THETA
300
270
Polares, r = 0.5
120
90
60
120
0.5
150
60
30
0.5
150
180
30
0
180
210
330
210
240
330
300
270
240
300
270
Polares, r = 1.0
Polares, r = 2.0
La Figura (4.20) muestra los valores del factor de amplificacion (4.114) para
r= 0.2, 0.5, 1.0, 2.0 que se encuentran dentro de la region de estabilidad. Las
Figuras (4.21), (4.22) y (4.23) muestran los valores del factor de amplificacion en
coordenadas polares para radio = 0.5, 1.0, 2.0 respectivamente que caen dentro
de la region de estabilidad |g| = 1.
4.4.2.
Esquema Crank-Nicolson
n+1
6
k
f(j, n + 1 )
2
j1
n+
1
2
w ?
j+1
n+1
n+1
n+1
n
n
vjn+1 vjn
2vjn + vj1
1 vj+1 2vj + vj1
1 vj+1
=
+
.
t
2
(x)2
2
(x)2
(4.115)
(4.116)
Este esquema utiliza 6 puntos de la malla, (xj1 , tn ), (xj , tn ), (xj+1 , tn ), (xj1 , tn+1 ),
En efecto: Por simplicidad denotamos por (x, t) el punto (xj , tn+ 1 ), tenemos
2
n+1
n+1
n+1
n+1
nj
1 j+1 2j + j1
1 nj+1 2nj + nj1
j
t
2
x2
2
x2
(4.117)
t
t
1
) = (x, t) xx (x, t)
t (x, t) + x2 xx (x, t) +
2
2
2!
2
t
1 t
x3
xxt (x, t) +
tt (x, t)
xxx (x, t)
2
2! 2
3!
t 2
t
t 3
2
3!
ttt (x, t)
2!
xxtt (x, t)
2!
x2 xxt (x, t) +
n+ 1
j 2
t
1
t
1
xx
t + x2 xx +
xxt +
2
2!
2
2!
t 3
t 2
t
x3
xxx
3!
3!
ttt
2!
xxtt
2!
t
2
x2 xxt
2
tt
(4.118)
t
1
t
1
+ xx
t + x2 xx
xxt +
2
2!
2
2!
t 3
t 2
t
n+ 1
j 2
x3
xxx
3!
3!
ttt +
2!
txtt
2!
t
2
2
tt
x2 xxt
(4.120)
t
1
t
1
xx +
t + x2 xx
xxt +
2
2!
2
2!
t 3
t 2
t
n+ 1
j 2
x3
xxx +
3!
3!
ttt
2!
xxtt +
2!
t
2
2
tt
x2 xxt
(4.121)
t
1
t
1
+ xx +
t + x2 xx +
txt +
2
2!
2
2!
t 3
t 2
t
n+ 1
j 2
x3
xxx +
3!
3!
ttt +
2!
xxtt +
2!
x
2
x2 xxt
2
tt
(4.123)
(4.124)
(4.125)
(4.126)
(4.127)
Esta u
ltima ecuacion puede ser escrita como
Lnj Lnj = O(t)2 + O(x)2
(4.128)
k
2
h2
2
h2
1 n+1 i(j+1)
r g e
2g n+1 eij + g n+1 ei(j1) +
2
1 n i(j+1)
r g e
2g n eij + g n ei(j1)
2
1
g 1 = r gei 2g + gei + r ei 2 + ei
2
2
(4.129)
g = 1 2grsen2 2rsen2
2
2
de donde obtenemos el factor de amplificacion
g=
1 2rsen2 2
1 + 2rsen2 2
(4.130)
4.4.3.
Esquemas Semi-Implcitos
(4.131)
0.5(x)2
(1 2)
(4.132)
90
120
1
60
r=0.5
0.5
150
r=1.0
30
0.5
G(THETA)
r=2.0
180
r=4.0
0.5
region estabilidad
210
330
1
240
1.5
300
270
3
THETA
6
Coordenadas polares, r=0.5
120
60
0.5
150
30
180
330
240
60
0.5
150
210
120
300
30
180
210
330
240
300
270
270
region de estabilidad.
En las Figuras (4.26), (4.27) y (4.28) se muestra el factor de amplificacion
del esquema (4.115) graficado en coordenadas polares para n
umeros de difusion
r = 0,5, 1,0, 2,0.
4.5.
(4.133)
ut = B uxx
(4.134)
(4.135)
Una gran simplificacion para ayudar a analizar (4.135) para sistemas hiperbolicos ocurren cuando el esquema tiene G como una funcion polinomial o racional de
la matriz A (e.g. los esquemas de Lax-Wendroff o Crank-Nicolson). Entonces la
misma matriz que diagonaliza la matriz A en (4.135) diagonaliza G, y la estabilidad
del esquema depende solo de la estabilidad de las ecuaciones escalares
ut + a t ux = 0
(4.136)
(4.137)
Esto implica que la estimacion (4.135) sera satisfecha para G si una estimativa
e
similar se tiene para G.
Para esquemas generales la situacion no es muy buena. Una condicion necesaria
para la estabilidad es
|g | = 1 + Kt
(4.138)
La matriz de amplificacion es
21
1 4rsen 2
G=
0
1
y asi
1
Gn =
0
4nrsen2 12
4.6.
Implementaci
on Num
erica
Esta seccion esta dedicada a la implementacion numerica de los esquemas
de diferencias finitas para el problema de conduccion del calor unidimensional analizados en la seccion anterior. Haciendo un estudio comparativo de los esquemas,
resolviendo el mismo problema bajo condiciones similares y para los mismos valores
de los parametros.
Para ello utilizamos el problema
ut
u(0, t)
(P )
u(L, t)
u(x, 0)
siguiente
= 2 uxx ,
= g1 (t),
t>0
= g2 (t),
t>0
= u0 (x),
0xL
(4.139)
= v1n
= v2n
..
.
n+1
n+1
rvi1
+ (1 + 2r)vin+1 rvi+1
= vin
..
.
n+1
n+1
n+1
n
rvN
= vN
1
2 + (1 + 2r)vN 1 rvN
al lado derecho de la ecuacion. Esto hace que el sistema resultante sea cerrado, es
decir solo tiene N 1 incognitas, el cual puede ser escrito, para cada n, en forma
matricial como
AU n+1 = bn
(4.140)
r
0 ...
0
0
1 + 2r
r
1 + 2r r . . .
0
0
..
A=
..
0
0
r 1 + 2r
r
0
0
r
1 + 2r
U n+1
v1n+1
v2n+1
..
.
n+1
vi
..
.
n+1
vN 2
n+1
vN 1
n+1
v1n + rv0
v2n
..
bn =
vin
..
n
vN
n+1
n
vN
1 + rvN
Para resolver el sistema (4.140), se puede utilizar eliminacion Gaussiana, asi como
tambien factorizacion lu. Para ello consideremos el sistema tridiagonal general de
la forma
0
0
d1 c1 0 . . .
a1 d2 c2 . . .
0
0
..
A=
..
0 0
aN 3 dN 2 cN 2
0 0
0
aN 2 dN 1
b1
b
2
.
..
= bi
..
.
bN 2
bN 1
v1n+1
v n+1
2
.
..
n+1
vi
..
.
n+1
vN 2
n+1
vN
1
d
1
c1
...
d2
c2
...
.
.
.
.
.
.
0
aN 2
dN 1
aN 1
l
0
0
cN 1
dN
0
0
...
...
.
.
.
.
.
.
0
lN 2
lN 1
u1
c1
...
u2
c2
...
.
.
.
.
.
.
0
uN 2
cN 1
uN
0
LET
u1 = d1
DO
i = 2 , N
l(i-1) = a(i-1)/u(i-1)
END DO
LET
DO
y1 = b1
i = 2 , N
yi = bi-l(i-1)*y(i-1)
END DO
LET
V(N) = y(N)/u(N)
DO
i = N-1 , 1
V(i) = (y(i)-c(i)*V(i+1))/u(i)
END DO
Este proceso es ejecutado para obtener U n+1 , y se repite para cada nivel el algoritmo para hallar la solucion del problema.
A continuacion implementamos el esquema implcito de Crank-Nicolson, para
ello, escribimos el esquema en la forma
n+1
n+1
n
n
= rvj+1
+ rvj1
+ (2 2r)vjn
+ (2 + 2r)vjn+1 rvj+1
rvj1
(4.141)
(4.142)
r
0 ...
0
0
2 + 2r
r
2 + 2r r . . .
0
0
..
A=
..
0
0
r 2 + 2r
r
0
0
r
2 + 2r
U n+1
v1n+1
v2n+1
..
.
n+1
vi
..
.
n+1
vN
2
n+1
vN
1
n
b =
n
rv0
+ (2 2r)v1n + rv2n + rv0n+1
..
.
n
n
n
+rvN
3 + (2 2r)vN 2 + rvN 1
n+1
n
n
n
+rvN
2 + (2 2r)vN 1 + rvN + rvN
Para resolver el sistema (4.142), se puede seguir los mismos pasos que se hizo al
resolver el sistema tridiagonal obtenido por el metodo de Euler retrasado
La implementacion del esquema de Du Fort-Frankel es simple porque es un
esquema explcito de paso m
ultiple, por consiguiente, este esquema se ejecutara a
partir del segundo nivel de tiempo, necesitando para ello los valores en los niveles
n = 0 (que es la condicion inicial) y n = 1, el cual puede ser obtenido, ejecutando
una u
nica vez, un esquema de un solo paso, como por ejemplo el esquema de Euler
progresivo o (FTCS)
n
n
), j = 1, . . . N 1,
+ vj1
vjn+1 = (1 2r)vjn + r(vj+1
n=0
(4.143)
vjn+1
2r
1 + 2r
n
vj+1
n
vj1
1 2r
1 + 2r
vjn1 ,
j = 1, . . . , N 1,
n1
(4.144)
Presentaremos resultados obtenidos por los tres esquemas en las mismas condiciones, luego se hace una comparacion con la solucion exacta, a fin de observar
sus propiedades presentadas en el analisis teorico, para ello, consideramos el valor
= 1, las condiciones de contorno
g1 = 0
g2 = 0
y la condicion inicial
u0 (x) = sen(x),
0<x<1
u(x, t) = e t sen(x)
=h
N
X
|vjn unj |2
j=0
donde unj = u(jh, nk). Hacemos enfasis, que solo en este caso, el primer paso del
esquema de Du Fort Frankel es calculado con la solucion exacta, de esta forma el
error que aparece en la siguiente tabla es un error del metodo mismo, y no de una
combinacion de dos metodos
Tambien presentaremos resultados para diferentes valores del parametro r y de
alli se obtendra el valor de k = h2 /r y si fijamos el valor que deseamos alcanzar
T , entonces el n
umero total de pasos en el tiempo se calcularan por la formula
Kmax = T /k
Por ejemplo los resultados que se presentan en las figuras (4.29), (4.30) y (4.31),
se obtuvieron con x = 0.1, T = 0.1, r = 1. Observemos que se presenta la
soluciones aproxiamadas comparadas ( lineas continuas ) con la solucion exacta
( lineas trazadas ), tambien se presenta el el error obtenido en el tiempo T con
la norma de l2 de las funciones discretas, es natural observar que el esquema de
Crank-Nicolson es mas preciso por tener una orden de precision (2,2), a diferencia
del esquema Euler retrasado que tiene orden de precision (1,2), y de Du FortFrankel que es de orden (2,2) pero con la condicion de k va para cero mucho mas
rapido que h.
Esquema
error=kek2
Euler Retrasado
1.0
0.01437
Crank-Nicolson
1.0
0.001933
Du Fort-Frankel
1.0
0.02303
Du FortFrankel
Euler retrasado
0.4
0.4
0.35
0.35
0.3
0.3
0.25
0.25
0.2
0.2
0.15
0.15
0.1
0.1
0.05
0.05
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
CrankNicolson
0.4
0.35
0.3
0.25
0.2
0.15
0.1
0.05
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
4.7.
Aplicaciones
Esta seccion esta dedicada a la presentacion de los esquemas de diferencias
estabilidad y consistencia en la seccion anterior. Hacemos un estudio comparativo de los esquemas resolviendo el mismo problema bajo condiciones similares y
dependencia de los parametros.
Presentamos resultados obtenidas con los esquemas de diferencias en las mismas condiciones, luego se hace una comparacion con la solucion exacta, a fin de
observar sus propiedades presentadas en el analisis teorico, para ello, consideramos
el siguiente problema
u
2u
= 2,
t
x
0 < x < 1,
u(0, t) = u(1, t) = 0,
u(x, 0) = sen(x),
t>0
t>0
0x1
u(x, t) = e t sen(x).
Tomamos los datos h = 0,1, = 1 y r = 1,0
FTCS Solucion exacta
0.4
0.35
0.35
0.3
0.3
0.25
eje y
eje y
0.25
0.2
0.2
0.15
0.15
0.1
0.1
0.05
0.05
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
Tmax=0.1, h=0.1, alfa=1, r=0.5
0.8
0.9
0.1
0.2
0.3
0.4
0.5
0.6
Tmax=0.1,h=0.1,alfa=1,r=1
0.7
0.8
0.9
CrnakNicolson exacta
0.4
0.35
exacta
0.3
cr
0.25
0.2
0.15
0.1
0.05
4
5
6
CrankNicolson, h=0.1, r=1
10
(4.145)
t=0 ,
u = U0 para y = 0
u = 0 para 0 < y < L
(b)Condiciones de frontera
t 0,
u = U0 para y = 0
u = 0 para y = L
40
40
35
35
30
30
25
25
U(M/SEG)
U(M/SEC)
20
t=1.08 seg
15
20
t=1.08
15
t=0.0 seg
10
10
10
0.01
0.01
0.02
Y(M)
0.03
0.04
0.05
10
0.01
t=0.0
0.01
0.02
Y(M)
0.03
0.04
0.05
En las Figuras (4.35) y (4.36)se muestra los perfiles de velocidad para la funcion
u.
Otro ejemplo es cuando se tienen fronteras aisladas, donde los extremos de
una barra se pueden aislar de modo que el calor no pueda escapar ni entrar por
los extremos. Sin embargo esto conduce a un problema trivial puesto que si la
barra inicialmente tiene una temperatura de 1000 C grados centgrados y si el
calor no puede escapar ni entrar por los extremos o por la superficie, entonces la
temperatura ciertamente permanecera a 1000 C grados centgrados todo el tiempo.
Para tener un problema no trivial tendramos que asumir que la temperatura
inicial no es constante. Asumamos por tanto que la barra esta aislada en ambos
extremos como tambien en la superficie pero que la distribucion de la temperatura
inicial esta especficada por alguna funcion f (x) donde 0 < x < L.
Formulaci
on matem
aticas: Usamos la ecuacion diferencial
u
2u
= 2
t
x
u
(0, t).
x
As,
u
(0, t)
x
u
(L, t)
x
= 0.
60 ,
0 x 50
f (x) =
40 , 50 < x 100
La solucion exacta obtenida por el metodo de separacion de variables para
= 0,16 es
u(x, t) = 50 +
n=1
40
n
sen
n
2
e16,10
6 n2 2 t
cos(
nx
)
100
(4.146)
Captulo 5
La Ecuaci
on del Transporte
5.1.
Ley de Conservaci
on de Masa
Sea > 0 la densidad de concentracion de un contaminante que se difunde
x2
(x, t)dx
(5.1)
x1
Suponiendo ademas que las paredes del tubo son impermeables y que la masa no se
crea ni se destruye, por lo tanto, la masa en el volumen de control puede cambiar
u
nicamente por el flujo de un extremo x1 y el otro extremo x2 . Por lo tanto el flujo
de la concentracion en el punto (x, t) esta dado por
f lujo de masa en (x, t) = a(x, t)(x, t)
(5.2)
x2
(5.3)
x1
x2
x1
(5.4)
x2
x1
de donde
d
(x, t)dx +
dt
x2
{
x1
x2
x1
(x, t) +
((x, t)a(x, t))}dx = 0
dt
x
(5.5)
(5.6)
la cual es valida para todo x en el intervalo [x1 , x2 ], para t > 0. Luego se concluye
que el integrando es identicamente cero, es decir,
d
(x, t) +
((x, t)a(x, t)) = 0
dt
x
(5.7)
(x, t) + a (x, t) = 0,
dt
x
(5.8)
Q
x
(5.9)
)
( x
2
Q
=
= 2
x
x
x
(5.10)
Si ademas consideramos una fuente generada por absorcion(proceso no estacionario, el cual hace crecer o decrecer la cantidad de sustancia disuelta y tiende a
1
El caudal drenado Q volumen por unidad de tiempo es proporcional a una seccion transversal
2
+a
= 2 +g
t
x
x
(5.11)
5.2.
Esquemas Explcitos
5.2.1.
Esquema FTCS
(5.12)
j = 0, ..., N ,
n = 0, ...T max,
y r = t/(x)2
(5.13)
Cr = at/x. donde la
1 2r
r 0.5Cr
0
r + 0.5Cr
0
r + 0.5Cr
1 2r
r 0.5Cr
0
0
r + 0.5Cr
1 2r
..
..
..
..
.
.
.
.
..
..
..
..
.
.
.
.
0
0
0
0
...
...
...
..
.
..
.
..
.
..
.
..
. r 0.5Cr
1 2r
r 0.5Cr
r 0.5Cr
0
1 2r
r + 0.5Cr
1 2r
r 0.5Cr
0
r + 0.5Cr
1 2r
..
..
..
.
.
.
..
..
..
.
.
.
0
0
0
...
...
...
..
.
..
..
.
.
1 2r
r + 0.5Cr
...
r 0.5Cr
1 2r
El desarrollo de la serie de Taylor alrededor del nodo (j, n)-esimo indica que
(5.13) es puntualmente consistente con (5.11) con un error de truncamiento de
orden O(t, x2 ). Para probar esto consideremos una funcion suficientemente
diferenciable y
desarrollemos la serie de Taylor de en el nodo (j, n)
(t)2
tt + O(t)3
2!
(x)2
(x)3
= nj + xx +
xx +
xxx + O(x)4
2!
3!
2
3
(x)
(x)
= nj xx +
xx
xxx + O(x)4 .
2!
3!
n+1
= nj + tt +
j
(5.14)
nj+1
(5.15)
nj1
(5.16)
(5.17)
x, t 0.
(5.18)
Z /h
1
=
ei j vbjn () d
2 /h
Z /h
1
=
ei (j1) vbjn () d
2 /h
Z /h
1
=
ei (j+1) vbjn () d
2 /h
vjn
n
vj1
n
vj+1
(5.19)
(5.20)
(5.21)
1
=
2
/h
/h
n+1
1
=
2
/h
ei j vbn+1 () d
(5.23)
/h
(5.24)
(5.25)
(5.26)
Observar que el metodo de Fourier es equivalente a reemplazar vjn por la expresion g n ei j lo cual es conocido con el nombre de An
alisis de Estabilidad de von
Neumann.
Para tener una idea grafica del comportamiento de los valores que toma el
factor de amplificacion (5.26) lo parametrizamos en la forma
x = 1 2r(1 cos )
y = C sen
r
(5.27)
(5.28)
(5.29)
por lo tanto
|g|2 = (4r2 Cr2 ) cos2 + 4r(1 2r) cos + 4r(r 1) + Cr2 + 1.
(5.30)
Esta u
ltima ecuacion es cuadratica en la variable cos y representa una curva
convexa con un mnimo o una curva concava con un maximo.
Se puede demostrar que para tener una solucion estable la funcion cuadratica
|g|2 no puede tener maximo global, es decir, la curva no puede ser concava. El
procedimiento matematico para demostrar esto es como sigue:
Calculamos la primera derivada:
d|g|2
= (4r2 Cr2 )(2 cos ) + 4r(1 2r)
d(cos )
(5.31)
y la segunda derivada
d2 |g|2
= 2(4r2 Cr2 )
d(cos )2
(5.32)
para que la funcion |g|2 tenga un maximo la segunda derivada tiene que ser negativa, es decir,
2(4r2 Cr2 ) < 0
Igualando a cero la ecuacion (5.31) obtenemos los puntos crticos
cos =
2r(1 2r)
4r2 Cr2
Cr2 (Cr2 4r + 1)
=
Cr2 4r2
Requerimos que el valor maximo de |g|2max debe ser menor que uno, es decir,
|g|2max 1
(5.33)
Puesto que (Cr2 2r)2 es siempre una cantidad positiva, la condicion de estabilidad
no puede cumplirse; y la funcion cuadratica representada por (5.30) la cual corresponde a una curva concava no es posible que sea alcanzada. Es decir, |g|2max > 1.
Concluimos que la funcion |g|2 no es concava pues si lo fuera tendra un valor
mayor que uno.
As se muestra que para una funcion concava donde
d2 |g|2
= 2(4r2 Cr2 ) < 0
2
d(cos )
Una solucion inestable se desarrollara debido al requerimiento que
(Cr2 2r)2 0
As la condicion para la inestabilidad como puede verse en la Figura (5.2), es
Cr2 > 2r
(5.34)
Inestabilidad FTCS
1.5
imaginario
C2 >2r
1
|g|=1
Eje Y
0.5
real
0
0.5
1.5
1.5
0.5
0
Eje X
0.5
1.5
(5.35)
factor de amplificacion g2
factor de amplificacion g2
g2
1.2
1.2
g2
0.8
0.8
r=0.5, C=0.5
0.6
0.6
0.4
r=0.4, C=0.3
0.4
2
estable
|g| <1
0.2
0
1
0.8
0.6
0.4
0.2
0
cos(theta)
0.2
0.4
| g |2 < 1
0.2
estable
0.6
0.8
0
1
0.8
0.6
1.2
r=0.5, C=1.1
0
cos(theta)
0.2
0.4
0.6
0.8
factor de amplificacion g2
1.2
g2
0.2
factor de amplificacion g2
0.4
g2
| g |2 > 1
r=0.4, C=1.0
Inestable
0.8
0.8
0.6
0.6
0.4
0.4
0.2
0.2
| g |2 > 1
0
1
0.8
0.6
0.4
0.2
0
cos(theta)
0.2
0.4
0.6
0.8
0
1
Inestable
0.8
0.6
0.4
0.2
0
cos(theta)
0.2
0.4
0.6
0.8
como
(a
t
t
)Cr 2
x
(x)2
1
aCr 2
x
es decir
a
x
2
Cr
2
Re
Cr
donde el n
umero de Reynolds es Re = a x
. Ademas, los extremos de cos es 1
(5.36)
8r(2r 1) 0
(5.37)
As, 2r 1 0 o r 0.5
Finalmente concluimos que las condiciones para la estabilidad son
r 0.5 o Re
2
Cr
(5.38)
t
x
2
0.5 y a
2
(x)
Cr
t
x
1 2
)(a
)
2
(x)
2 Cr
1
t
a
x
Cr
(5.39)
Por lo tanto Cr 1.
Finalmente podemos comparar el rango de estabilidad en el plano (r, Cr ) para
el esquema FTCS el cual es esquematizado en la Figura (5.7)
Comparacion r y Cr
1
2.0
0.9
0.8
0.7
Cr
0.6
1.0
0.5
5
0.4
0.3
0.5
10
0.2
0.2
0.1
0.05
0.1
0.15
0.2
0.25
r
0.3
0.35
0.4
0.45
0.5
El n
umero adimensional de Peclet, P e , que aparece en la Figura (5.7), se define
como P e = |a| x
o P e =
Cr
r
ser menor que 2. Observe que para r = 0, corresponde al caso de conveccion pura,
el esquema FTCS es inestable para todo Cr y .
En la Figura (5.8) se muestra el factor de amplificacion (5.26) para los parametros Cr y r. Para un Cr = 0,5 fijo, variamos r con valores de 0.2, 0.3, 0.4, 0.5
y observamos que las curvas generadas por el factor de amplificacion (5.26) caen
dentro de la region de estabilidad |g|2 1, esto es, la condicion (??) se cumple.
Para Cr = 0,5 y r = 0,6 la curva generado por (5.26) sale fuera de la region de
estabilidad, esto hace que el esquema FTCS se inestable para estos valores.
imaginario
1
|g|=1
Eje Y
0.5
real
0
r=0.5
r=0.6
0.5
r=0.3
r=0.4
r=0.2
region inestabilidad
region estabilidad
1.5
1.5
0.5
0.5
1.5
Eje X
5.2.2.
Esquema Upwind
Este esquema es la union del esquema Upwind para la ecuacion de conveccion con Cr 1 como condicion de estabilidad y el esquema FTCS para la
ecuacion de difusion con condicion de estabilidad r 0.5.
Si aplicamos el esquema Upwind a la ecuacion (5.11) obtenemos la ecuacion
algebraica
n
n
n
vjn+1 vjn
vjn vj1
vj1
2vjn + vj+1
+a
=
.
t
x
(x)2
(5.40)
(5.41)
(5.42)
(5.43)
t, x 0.
vjn
1
=
2
n
vj1
1
=
2
n
vj+1
1
=
2
/h
/h
/h
/h
/h
/h
ei j vbjn () d
(5.44)
ei (j1) vbjn () d
(5.45)
ei (j+1) vbjn () d
(5.46)
1
=
2
/h
/h
[(r + Cr ) ei + (1 2r Cr ) + r ei ] ei j vbjn () d
(5.47)
n+1
1
=
2
/h
/h
ei j vbn+1 () d
(5.48)
(5.49)
(5.50)
(5.51)
La formula (5.50) muestra que la ventaja de la solucion del esquema para un paso
de tiempo es equivalente a multiplicar la transformada de Fourier de la solucion por
un factor de amplificacion g(). El factor de amplificacion es llamado as porque su
magnitud es el valor de la amplitud de cada frecuencia en la solucion, dada por vbn
y es amplificado en un avance de la solucion sobre el paso del tiempo para obtener
vbn () = [g()]n vb0
(5.52)
(5.53)
En adelante podemos simplificar este procedimiento y solo reemplazar vjn por g n eij
lo cual se conoce con el nombre de analisis de estabilidad de von Neumann.
Para tener una idea geometrica del comportamiento del factor de amplificacion
(5.53) lo podemos parametrizar en la forma
y = C sen
r
(5.54)
(5.55)
de centro (1 2r Cr , 0) y semiejes 2r + Cr y Cr .
Se puede observar que para tener soluciones estables el modulo del factor de
amplificacion (5.53) debe ser menor que uno, lo cual implica que
Cr + 2r 1.
(5.56)
imaginario
1
g = 1.0
r=0.2
0.5
r=0.4
r=0.3
r=0.5
eje y
real
0
0.5
1.5
0.5
0
Cr = 0.5, r = 0.2; 0.3; 0.4; 0.5
0.5
(5.57)
5.2.3.
Esquema Du Fort-Frankel
(5.58)
o en forma de algoritmo
n
n
n
n
)
vjn+1 = vjn1 Cr (vj+1
vj1
) + 2r(vj+1
(vjn1 + vjn+1 ) + vj1
(5.59)
v (t)2 2 v (t)3 3 v
+
+
+ O(t)4
t
2! t2
3! t3
(5.60)
(5.61)
Desarrollamos la serie de Taylor alrededor del nodo (j, n) para la funcion v con un
incremento de x
n
vj+1
= vjn + (x)
v 1
2v
1
3v
+ (x)2 2 + (x)3 3 + O(x)4
x 2
x
3!
x
(5.62)
n
vj1
= vjn (x)
v 1
2v
1
3v
+ (x)2 2 (x)3 3 + O(x)4
x 2
x
3!
x
(5.63)
(5.64)
v (t)2 2 v (t)3 3 v
+
+ O(t)4
t
2! t2
3! t3
(5.65)
2v
) + O(t)4
t2
(5.66)
v
v
2v
t 2 v
t4
+a
2 + O(x)2 + O(t) + (
) 2 + O( 4 )
t
x
x
x t
x
(5.67)
v
v
2v
+a
2
t
x
x
(5.68)
tenemos
Lvjn Lnj v = O(x)2 + O(t) + (
t 2 v
t4
) 2 + O( 4 )
x t
x
(5.69)
(5.70)
(5.71)
[ 2 4(4r2 1)]1/2
2(1 + 2r)
(5.72)
Para analizar la estabilidad del esquema (5.59) se debe tener |g|2 < 1 la cual
por una experiencia computacional se demuestra que Cr 1, no existiendo restricciones para r.
En la Figura (5.10) tenemos el factor de amplificacion (5.72) que para los
valores de r = 0.5 fijo variamos Cr para 0.2 y 0.5. Sus valores caen dentro en la
region de estabilidad y para el n
umero de Courant Cr = 1.2 los valores caen fuera
de la region de estabilidad |g|2 < 1.
g, Cr=1.2
0.5
g, Cr=0.5
g+, Cr=1.2
g+, Cr=0.5
eje Y
g, Cr=0.2
g+, Cr=0.2
real
0.5
1.5
1.5
0.5
0
eje X
0.5
1.5
5.2.4.
Esquema Lax-Wendroff
n
n
2vjn + vj1
vj+1
x2
=0
(5.73)
(5.74)
Entonces el esquema Lax-Wendroff para la ecuacion del transporte (5.11) esta dado
por la union de los esquemas (5.73) y (??).
n
n
n
n
vjn+1 vjn
2vjn + vj1
vj1
vj+1
vj+1
+a
=
t
2x
x2
(5.75)
(5.76)
(5.77)
(5.78)
x = 1 2r (1 cos )
(5.79)
y = C sen
r
0.8
region estabilidad
0.6
0.4
eje y
0.2
r*=0.5
r*=0.6
r*=0.4
r*=0.3
r*=0.2
0.6
0.8
1
1.5
0.5
0
Cr=0.5, r*= 0.2; 0.3; 0.4; 0.5; 0.6
0.5
(5.80)
(5.81)
La Figura (5.11) se muestra el factor de amplificacion para el esquema LaxWendroff. Si tomamos los valores Cr = 0.5 y r = 0.2, 0.3, 0.4, 0.5 vemos que el
factor de amplificacion se encuentra dentro de la region de estabilidad |g|2 1 es
decir se cumple (5.81), mientras que para Cr = 0.5 y r = 0.6 la condicion (5.81)
no se cumple.
5.3.
Esquemas Implcitos
Para la ecuacion de difusion los esquemas implcitos son efectivos para
5.3.1.
Esquema de Crank-Nicolson
!
n+1
n+1
n
n
vjn+1 vjn
vj+1
vj1
vj+1
vj1
+ 0.5a
+
t
2x
2x
!
n+1
n+1
n
n
vj1
2vjn+1 + vj+1
vj1
2vjn + vj+1
= 0.5
+
x2
x2
(5.82)
el cual es la union de los esquemas de Crank-Nicolson para las ecuaciones de
difusion (4.115) y conveccion (3.98) vistas anteriormente.
El esquema (5.82) puede ser escrito como un algoritmo
n+1
n+1
(r + 0.5Cr )vj1
+ 2(1 + r)vjn+1 (r 0.5Cr )vj+1
n
n
= (r + 0.5Cr )vj1
+ 2(1 r)vjn + (r 0.5Cr )vj+1
.
(5.83)
Para probar la consistencia de (5.83) se considera una funcion , lo suficientemente diferenciable. Entonce hallamos los valores (xj1 , tn ), (xj , tn ), (xj+1 , tn ),
(xj1 , tn+1 ), (xj , tn+1 ) y (xj+1 , tn+1 ) en el nodo (xj , tn+1/2 ). Entonces utilizando
la serie de Taylor llegamos a determinar el operador diferencial discreto asociado
a (5.82)
Lnj = t + ax xx + O(t2 ) + O(x2 )
(5.84)
(5.85)
(5.86)
As que,
tal que
Lnj Lnj 0 ,
t, x 0.
(5.87)
x=
1A2 B 2
(1+A)2 +B 2
y=
2B
(1+A)2 +B 2
(5.88)
eje Y
0.5
r=0.8
0.5
r=0.5
r=0.2
1.5
1.5
0.5
0
eje X
0.5
1.5
5.4.
Resultados Num
ericos
En esta seccion presentamos los esquemas de diferencias descritos en las
secciones anteriores aplicados al problema de transmision de una fuente de temperatura con una velocidad a a traves de un fluido de difusividad termal .
Los esquemas de diferencias FTCS, Du Fort-Frankel, Upwind, Lax-Wendroff
y Crank-Nicolson son aplicados al problema de propagacion de una fuente de
temperatura, Figura (5.13)
2.0
0.0
2.0
u(2, t) = 0.
N
2X
(x at) exp[(2k 1)2 2 t/L2 ]
u(x, t) = 0.5
sen (2k 1)
k=1
L
2k 1
(5.89)
Hacemos una comparacion entre los esquemas discutidos en la seccion anterior resolviendo el mismo problema bajo condiciones similares y dependencia de
parametros.
0.8
0.8
0.6
0.6
u(x,t)
u(x,t)
exacta
0.4
0.4
0.2
LaxWendroff
exacta
0.2
ftcs
0
80
100
120
140
Cr=0.5, r=0.5 , t=110
160
180
200
0.2
30
40
50
90
100
110
120
0.9
0.9
0.8
0.8
0.7
dffr
exact
0.6
0.6
upwind
exacta
0.5
0.5
0.4
0.4
CN
0.3
0.3
0.2
0.2
0.1
0.1
70
80
Cr=0.4; r*=0.48
0.7
60
0
40
20
40
60
80
100
120
140
160
60
80
100
120
eje x
140
160
180
200
(5.90)
(5.91)
(5.92)
Se presentan las graficas de las soluciones aproximadas para los diferentes esquemas
estudiados para los valores L = 1 longitud del intervalo; a = 1 velocidad del
fluido efectiva, = 0.001 difusion molecular y k = 4.3 2 . Ademas el perfil de
concentracion del contaminante en el instante inicial t = 0 seg esta dado por
1 , 1 < x < 0
(x, 0) =
0 , otro lugar
k=
b
n
Captulo 6
Ecuaciones Elpticas
En este captulo hacemos un estudio en la forma mas simple de las ecuaciones elpticas, orientado principalmente al aspecto computacional antes que analisis.
6.1.
Ecuaci
on de Poisson
Consideremos como modelo de trabajo la ecuacion de Poisson
Au = f,
en
(6.1)
en
281
(6.2)
6.1.1.
Caso unidimensional
d2 u2
dx
= f,
en (0, 1)
u(0) = 0
(6.3)
u(1) = 0
j = 1, 2 . . . , N
(6.4)
uN +1 = 0
(6.5)
AU = F,
o en forma extensa como
AU =
x2
(6.6)
...
...
1 . . .
..
.
...
...
0 v1
0
0
v2
0
0
v3
..
.
2 1 vN 1
vN
1 2
0
f1
f2
f3
=
..
.
fN 1
fN
(6.7)
Observaci
on 1
Como era de esperar, el sistema continuo (6.1), es aproximado por el sistema de
ecuaciones algebraicas (6.6), es decir, el operador lineal continuo, infinito dimensional, A es sustituido por un operador lineal finito dimensional A
La matrz de coeficientes del sistema de ecuaciones (6.6) es una matrz tridiagonal, la cual esta dentro de la clase de matrices banda con ancho de banda w
como se presenta en la figura (6.1). La solucion del sistema anterior tiene la forma U = A1 F , pero la matriz A1 no es banda, por consiguiente es necesario
saber utilizar algoritmos adecuados para resolver este tipo de sistemas, en seguida
se describe un prodecimiento muy estandar para resolver este tipo de problemas.
1. Utilizar la descomposicion lu: que consiste en factorizar la matriz A en dos
matrices, una matriz triangular inferior y una matriz triangular superior y
escribir la matriz en la forma siguiente
A=
w- w
w-
w-
= lu
desacoplan-
do el sistema por cada fila. Para llevar a cabo esto, simplemente observamos
que, los autovectores son las columnas de la matriz S que tienen la forma
senjh
sen2jh
sen3jh
..
senN jh
(6.8)
...
1 . . .
2
...
..
.
...
...
0
0
0
0
2 1
1 2
senjh
sen2jh
sen3jh
..
..
senN jh
senjh
sen2jh
sen3jh
(6.9)
..
..
senN jh
Y =C
M U = (M A)U + F
(6.10)
Ahora resolvemos (6.10) iterativamente por sustitucion sucesiva. Comenzamos con un dado inicial X0 que puede ser cero. Luego si el valor actual es
Uk , en cualquier paso podemos obtener el nuevo vector Uk+1 por la formula
M Uk+1 = (M A)Uk + F
(6.11)
6.1.2.
Caso bidimensional
xu2 +
2u
y 2
= f,
en = (0, 1) (0, 1)
(6.12)
= g(x, y)
El dominio consta de un rectangulo unitario, la discretizacion del dominio lo realizamos en forma homogenea, x = y h = 1/(N + 1). Las derivadas 2 u/x2 y
2 u/y 2 tienen aproximacion de tres puntos, tanto en la direccion x como en la direccion y, combinando las dos direcciones obtenemos una molecula de cinco puntos
como se presenta en la figura (6.2). Las incognitas a determinarse corresponden
a los puntos de la malla, los cuales son un total de N 2 puntos enn el interior del
dominio.
j+1
z(i, j)x
j1
1
1
i1
i+1
(6.13)
i, j = 1, . . . N
vi0 = 0 = viN +1 ,
vj0 = 0 = vjN +1 ,
j = 0, . . . N + 1
Las incognitas vij pueden ordenarse en una manera natural, empezando con
v00 ,v10 , v20 ,. . ., vN 0 en la parte inferior del cuadrado. Este grupo sera denotado
por V0 . La proxima fila de puntos de la malla contribuyen a formar el grupo V1 .
i = 1, . . . N
j = 1, 2, . . . N
(6.14)
donde T es una matriz similar a la matriz A del sistema (6.6), con la diferencia
que en la diagonal en lugar de 2 tenemos el valor 4.
Utilizando los mismas tecnicas como en el caso unidimensional el sistema matricial (6.14), se puede escribir como:
KU = F
donde
K=
...
...
I . . .
..
.
...
...
0
0
0
0
T I
I T
6.2.
y
9
10
12
13
16
15
11
14 x
Por ejemplo, cuando se modela un reservorio con una frontera de roca impermeable, la condicion que no fluye agua a traves de la frontera es expresada por
u/n = 0, donde u es el potencial Hidraulico. En los generadores electricos la
corriente fluyendo en los hilos electricos, genera calor que es irradiada en el medio
ambiente. Este proceso es usualmente representado por una ecuacion, tal como:
u/n = g(x, y, u)
Para precisar consideremos la ecuacion de Laplace en el siguiente problema modelo:
2u
2u
= 0,
0 < x < 1,
0<y<1
2 + y 2
u
(6.15)
= 10,
y=1
y
u = 0,
en todos los otros lados
Es decir, en tres lados el rectangulo se aplica la condicion de la frontera de
Dirichlet, mientras que en un lado se aplica la condicion de Neumann.
4 1 1 0
1 4
0 1
1 0
4 1
0 1 1 4
analoga a
u1
u2
=
u3
u4
g6 + g16
g13 + g15
(6.16)
g7 + g9
g10 + g12
En la ecuacion (6.16) los valores nodales u9 y u10 son representados por g9 y g10 , los
cuales no son especificados, en lugar de ellos tenemos las derivadas de los valores
u
nodales u
,
tenemos varios metodos de tratar este problema.
y
y
9
6.2.1.
10
M
etodos de los puntos ficticios
u
y
rencias, como el esquema de diferencia central de 2do orden, por ejemplo en el nodo
9 escribimos
u
y
u9 u3
2y
u9 u3 + 2y
u
y
u1 + 4u3 u4 = g7 + 2y
u
y
u1 + 3u3 u4 = g7 + 2y
u
y
(6.17)
9
2
= 0 + (10)
3
De esta forma la ecuacion
6.2.2.
(6.18)
1 1 0 u1 0
u2 0
4
0 1
20
u3 3
0
3 1
20
u4
1 1 3
3
(6.19)
Aproximaci
on por series de Taylor (sin puntos ficticios)
Las aproximaciones de diferencias para puntos en la frontera, tambien puede hacerse sin recurrir a puntos ficticios usando series de Taylor como veremos
en el siguiente procedimiento, donde se presento una aproximacion de 2do orden.
Expresando los valores de u en los nodos 1 y 3 de la malla en terminos del valor
en el nodo 9 (ver figura(??))
u3 = u9 y
u
y
u1 = u9 (2y)
u
y
(y)2
+
2
9
(2y)2
+
2!
9
2u
y 2
(6.20)
9
2u
y 2
(6.21)
9
u
y
+ O(y 3 )
9
resolviendo para
u
y
u
y
=
9
u1 4u3 + 3u9
O(y 2 )
2y
y 9
2y
Es una aproximacion de segundo orden de la derivada en la frontera. Resolviendo
para u9 , tenemos
1
u
4u3 u1 + 2y
u9 =
3
y 9
1
10]
3
la cual, reescribiendo es
8
20
2
u1 + u3 u4 =
3
3
9
(6.22)
o
2u1 + 8u3 3u4 =
20
3
(6.23)
20
3
1 1 0 u1 0
u2 0
4
0 1
20
0
8
3
u3 3
20
u4
2 3 8
3
(6.24)
Este sistema, desde luego, es diferente del sistema (6.19) donde se usa puntos
ficticios.
6.2.3.
Si la frontera del dominio no coincide con las lineas de la malla. Las ecuaciones de diferencias finitas correspondientes a los nodos ubicados en la vecindad
inmediata de la frontera deben ser modificados ligeramente.
Consideremos el problema con una frontera curvada como se muestra en la
firgura (6.4), observe que los puntos A y B no pertencen a la malla, estan a una
distancia escalada por los n
umeros positivos a y b, menores que uno.
A
f
3
f
4
? f
y
?
by
ax
2
f
1
f
x -
u
(ax)2 2 u
= u4 + (ax)
+ O(x3 )
+
2
x
2!
x
4
2 4
2
u
(x ) u
= u4 (x)
+
+ O(x3 )
x 4
2!
x2 4
uB
u3
(6.25)
(6.26)
u
x
1
=
x
1
1a
a
uB
u4
u3
a(1 + a)
a
1+a
(6.27)
el cual tiene orden de precision 2. De igual manera, siguiendo los mismos argumentos, la eliminacion de (u/x)4 nos conduce a la siguiente formula
2u
x2
1
=
x2
2
2
2
uB +
u3 u4
a(1 + a)
1+a
a
(6.28)
2u
y 2
1
=
y 2
2
2
2
uA +
u2 u 4 ,
b(1 + b)
1+b
b
(6.29)
con precision del orden y. Por consiguiente, ( xu2 ) puede aproximarse en los puntos malla 1, 2 y 3 con una precision de orden (x)2 , mientras la precision global
es solo del orden (x) debido a las irregularidades en la frontera cerca del punto
malla 4.
Sustituyendo (6.28) y (6.29) apropiadamente en (6.16) obtenemos, para x =
y,
2
b(1+b)
2
a(1+a)
u1
u2
1
u3
1
2
2
u
+
4
a
b
0
g6 + g16
g13 + g15
g7 + g9
2
g
b(1+b) 10
2
g
a(1+a) 12
(6.30)
Captulo 7
Convergencia de problemas
lineales y no lineales
En este captulo haremos un analisis de esquemas de diferencias, en un nivel
un poco mas profundo, para problemas mas generales y problemas no lineales.
Basados en la experiencia de los captulos anteriores, por tanto se necesita un
poco mas de cuidado en el analsis de los mismos conceptos ya conocidos, haciendo
enfasis en los problemas no lineales, que por su naturaleza no tienen un tratamiento
igual al caso de los lineales.
296
7.1.
Estabilidad y Convergencia
Consideremos el (PVI) periodico para sistemas lineales de Ecuacines Dife-
x, t,
x
u,
(7.1)
u(x, 0) = u0 (x)
Suponga que u0 y los coeficientes son 2 periodicos en todas las direcciones del
espacio.
7.1.1.
Discretizaci
on del Dominio
j = 0. 1, 2, . . .
h = 2/(N + 1),
N un n
umero natural
7.1.2.
Discretizaci
on de la Variable
7.1.3.
Discretizaci
on de la Ecuaci
on
n+1
v0 =
q
X
Q v n ,
=0
u0 ,
n = q, q + 1, . . . ,
(7.2)
= 0, 1, . . . , q.
n+q
v
n+q1
v
vn =
..
vn
u0
q1
u0
,
u
=
0
..
u00
(7.3)
(7.4)
donde
Q(tn ) = Q(xj , tn ) =
Q1
1 Q0
Q1
1 Q1
...
Q1
1 Qq
I
I
I
(7.5)
es la matriz de Companion.
Por ejemplo, el esquema de Leapfrog puede ser escrito como
2kD0 I
n
n+1
v
=
v
I
0
(7.6)
(7.7)
Sh (tn , tn ) = I.
Si Q es independiente de t tenemos
Sh (tn , t ) = Qn ,
(7.8)
kvn kh =
q
X
!1/2
kv n+ k2h
=0
(7.9)
Definici
on 7.1.1 La Aproximaci
on de diferencias (7.4) es estable para 0 < h <
h0 , si existen constantes s , C, Ks tal que para todo h
kQ1
1 kh C,
(7.10)
K(tn ) = Ks es tn
(7.11)
Q1 v
n+1
v0 =
q
X
=0
u0 ,
Q v n + kF n ,
n = q, q + 1, . . . ,
(7.12)
= 0, 1, . . . , q.
n
T
F = (Q1
1 F , 0, . . . , 0) ,
(7.13)
v 0 = u0 ,
La version discreta del Principio de Duhamel es dado en el siguiente lema
Lema 2 La solucion de las ecuaciones (7.13) pueden ser escritas en la forma
v = Sh (tn , 0)u0 + k
n1
X
Sh (tn , t+1 )F .
(7.14)
=0
S t n
kv kh Ks e
ku0 kh + h (S , tn ) max kF kh
n
0n1
donde
h (S , tn )
n1
X
Z
S (tn t+1 )k
eS (t) d = (S , tn ).
=0
Prueba:
De la ecuacion (7.14) podemos tener
n
n1
X
=0
Q1 ve
q
X
Q ven + n ,
kn k
=0
(7.15)
n+1
q
X
Q wn + n .
(7.16)
=0
kwn kh C(tn ) .
k
(7.17)
= 1, 0, . . . , q
(7.18)
n+1
q
X
=
(Q + kR )v n .
(7.19)
=0
Ahora probaremos que la ecuacion (7.19) es estable si la aproximacion no perturbada (7.2) es estable. Para la prueba utilizamos el siguiente lema:
Lema 3 Sea Q un operador acotado. Para kkQkh < 1
k(I kQ)1 kh
1
1 kkQkh
kQkh
1 kkQkh
Por tanto usando el lema anterior, podemos escribir la ecucion (7.19) en la forma
en+1 = (Q + kR)e
v
vn
(7.20)
(7.21)
donde
n = Rwn .
F
El principio de Duhamel puede ser aplicado a la ecuacion (7.21). El operador
solucion Sh (tn , t ), que corresponde a ek Q se escribe como ek(n) QSh (tn , t ),
donde Sh es el operador solucion para la ecuacion (7.4).
As
kek(n) Sh (tn , t )kh KS e(S )(n) .
S + O(h).
7.2.
Definici
on 7.2.1 Sea u(x, t) una funcion suave de (7.1) entonces el error de truncaci
on local se define por
kjn
= Q1 u(xj , tn+1 )
q
X
=0
Q u(xj , tn ).
(7.22)
(7.23)
h2
k2
uxxx (xj , tn ) + uttt (xj , tn ) + O(h4 , k 4 )
3
3
(7.24)
h2
= ch2 utt uxxxx (xj , tn ) + O(h2+4 ).
2
12
(7.25)
Definici
on 7.2.2 Los datos iniciales tienen orden de precisi
on (r, p) si para toda
soluci
on suficientemente suave u(x, t) de (7.1)
kv u(., k)kh L1 (hr + k p ),
= 0, 1, . . . , q.
(7.26)
S t n
kv u(tn )kh KS e
kv u(0)kh +
kQ1
1 kh h (, tn )
max k kh
0jn1
= O(hr + k p ),
es decir, las soluciones del esquema de diferencias converge cuando h 0 a la
soluci
on de la ecuaci
on diferencial.
Prueba:
Si la solucion es sustituida en el esquema de diferencias (7.2) entonces obtenemos
Q1 u(xj , tn+1 ) =
q
X
Q u(xj , tn ) + kjn .
=0
q
X
Q wjn + kjn ,
=0
wj = u(xj , k) vj ,
= 0, 1, . . . , q
(7.27)
7.3.
(7.28)
(7.29)
Entonces para diferenciar llamaremos a la ecuacion (7.28) como ecuacion de Burgers disipativa.
7.3.1.
(7.30)
1
n
n
vj1
+ r vj+1
2vjn + vj1
vjn+1 = vjn vjn vj+1
2
donde = k/h2 y r = k/h2
7.3.2.
Consistencia
n+1
= nj + kt +
j
nj+1
nj1
(7.31)
7.3.3.
Estabilidad
= 1 i vsen 4r sen2
2
luego la condicion de estabilidad dice que
r
(1 4 r sen2 )2 + ( v sen)2 1
2
si y solamente s
|1 4rsen2 | 1 si y solamente s
2
1
2
h
(1 |1 2r|).
k
K > 0,
cR
X
h
|vj |
j
se tiene
1
n
|vjn+1 | |vjn ||1 2r| + |vjn ||vj+1
|
2
1
n
n
n
+ |vjn ||vj1
| + r|vj+1
| + r|vj1
|
2
sumando
1 X n n
n
|vj ||vj+1 | + |vjn ||vj1
| + 2rkvjn ks )
kvjn+1 ks kvjn ks |1 2r| +
2 j
1
(|1 2r| + 2r)kvjn ks + |v|max kvjn ks
2
|1 2r| + 2r + |v|max kv n ks
2
Como observamos en la u
ltima expresion la estabilidad existe, pero depende de
la solucion misma y se debe cumplir |1 2r| + 2r + 12 |v|max 1, esta expresion
indica que la condicion es mucho mas restrictiva que la condicion obtenida por el
criterio de von Neumann.
kvj k2 =
!1/2
|vj |2
kvjn+1 k22 = hvjn kvjn D0 vjn + k0 vjn , vjn kvjn D0 vjn + k0 vjn i
= kvjn k22 2khvjn Dv0 j n , vjn i + k 2 kvjn D0 vjn k22 2k 2 h0 vjn , vjn D0 vjn i
+2kh0 vjn , vjn i + k 2 kvjn k2
1
(|1 2r| + 2r)kvjn ks + |v|max kvjn ks
2
Prueba:
kvk2max
2
max |vj |
j
max |vj |2
j
|vj |2
|vj |2
X
h
h
|vj |2
h
j
1
kvk22
h
Lema 5
a) hu, vi kuk2 kvk2
b) hu, vi kuk2 kvk2
Desigualdad de Schwarz
kuk22 +kvk22
2
Prueba: Unicamente
probemos la parte b)
0 hv u, v ui = h
X
X
(vj uj )(vj ui ) = h
(vj2 2vj uj + u2j )
j
kvk22
2|hu, vi| +
kuk22
Prueba:
hu, vwi = h
uj vj wj hkvkmax h
uj wj
Lema 7
kD0 vk2 kD+ vk2 = kD vk2
donde
D+ v =
vj+1 vj
;
h
D v =
vj vj1
;
h
D0 =
D+ + D
2
Prueba:
kD0 vk2
1
1
kD+ vk2 + kD vk2
2
2
pero
X 1
|v vj |2
2 j+1
h
j
j
X
X 1
2
|v
v
|
=
h
|D vj |2 = kD vk22
= h
j1
2 j
h
j
j
kD+ vk22 = h
|D+ vj |2 = h
por tanto
kD+ vk22 = kD vk22
y de ello
1
1
kD0 vk2 = kD+ vk2 + kD+ vk2 = kD+ vk2
2
2
Lema 8
hv, 0 vi = hDv , Dv i
Prueba: Sabemos que 0 = D+ D , entonces
hv, 0 vi = h
vj D + D vj = h
vj D vj+1
X D vj+1 D vj
vj D vj =
(vj1 vj ) D vj = h
vj1 D vj
vj D vj
D vj D vj =
= hD v, D vi
Corolario 1
hv, 0 vi = kD vk22 = kD+ vk22
Ahora pasemos a analizar la estabilidad del esquema (FTCS) para la ecuacion de
Burgers. El esquema se puede escribir en la forma
vjn+1 = vjn kvjn D0 vjn + kD+ D vjn
tomando producto interno con vjn a la ecuacion anterior se tiene
hvjn+1 , vjn i kv n k22 hkvjn D0 vjn , vjn i + hkD+ D vjn , vjn i
kv n k22 + kkv n kmax kD0 v n k2 kv n k2 kkD v n k22
(7.32)
luego tenemos
kv n kmax
kv n k22 + 2k kD v n k2 kv n k2
2
1 kv n k2max
n 2
n 2
kv n k22
kv k2 + 2k
kD v k2 +
2
2 4
simplificando terminos se tiene
hvjn+1 , vjn i
k n 2
1+
kv kmax kv n k22
4
kv
n+1
k2
kC 2
1+
4
kv n k2
Observando la u
ltima desigualdad se puede iterar hasta tener
kv
n+1
k2
kC 2
1+
4
t0 C
kv 0 k2 e 4 kv 0 k
(7.33)
podemos observar la desigualdad (7.33) y decir que no hay estabilidad, para t0 sera
incontrolable. Sin embargo puede haber convergencia si hacemos w = u(tn , xj )vjn
y logramos que w satisfaga una ecuacion como la (7.33).
El siguiente lema presenta una estimativa para w
Lema 9 Si kC 2 /4 < 1, entonces se cumple
kwn+1 k2 (1 + C1 k)kwn k2 + kO(h2 )
con estos lemas podemos enunciar el teorema de converencia
Teorema 7.3.1 Bajo las hipotesis de los lemas precedentes el error w satisface
kwn k2 = O(h2 )
kwn kmax = O(h3 /2).
7.3.4.
Implementaci
on num
erica
1 2
u x=0
2
(7.34)
(7.35)
(7.36)
mas detalladamente
vjn+1 = vjn r
1 n 2 n 2
vj+1 vj1
,
2
(7.37)
u0 (x) =
1 , x<0
1x , 0x1
0 , x > 1.
(7.38)
No conservativo
Conservativo
1.3
1.3
1.2
1.2
1.1
1.1
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
-2
-1.5
-1
-0.5
0.5
1.5
0.3
-2
-1.5
-1
-0.5
0.5
1.5
Bibliografa
[1] Abbott, M. and Basco, D. Computational Fluid Dynamics. An Introduction for Engineers. Copublished in the United States with Jhon Wiley &
Sons, Inc. , New York, 1989.
[2] Angirasa, D., Eggles, J. and Nieuwstadt, F. Numerical Simulation of
Transient Natural Convection form an Isothermal Cavity Open on a Side,
Numerical Heat Transfer,Part A., 28,755-768, 1995.
[3] Bejan, A. Convection Heat Transfer. Jhon Wiley & Sons. Inc. USA. 1992.
[4] Claeyssen, J. Bravo, E. and Rubio, O. Thermally driven convetive flow
in the boundary condition for the pressure, Numerical Heat Transfer,34,
658-672, 2002
[5] Durran, R. Numerical Methods for Wave Equations in Geophysical Fluid
Dynamics.Edit. Springer-Verlag, New York 1998.
[6] Duchateau, P. and Zachmann, D. Ecuaciones Diferenciales Parciales,
Teora y Problemas, McGraw-Hill,Mexico 1998.
320
[7] Evans, G., Blackledge, J. and Yardley, P. Numerical Methods for Partial
Differential Equations. Edit. Springer-Verlag. London, 2000.
[8] Farlow, S. Partial Diferential Equations for Scientists and Engineers,
Jhon Wiley & Sons, New York, 1982.
[9] Fletcher, C. Computational Techniques for Fluid Dynamics 1. SpringerVerlag, Sydney, 1990.
[10] Flard, J., Manteuffel,T. and Mccormick,S. First-Order System Least
Squares(Fosls) For Convection-Difusion Problemas :Numerical Results,SIAM J. SCL.Comput. Vol.19, No.5,pp. 1958-1979, 1998.
[11] Gustafsson, B ; Kreiss, H. & Oliger, Joseph Time Dependent Problems and
Difference Methods Wiley-Interscience Publication. Jhon wiley & Sons,
INC; New York; 1995.
[12] Hoffmann, K. and Chiag S. Computational Fulid Dynamics for Engineers,
Vol I. Library of Congress Catalog Card Number 93-090310, 1993.
[13] Marshall, G. Solucion Numerica de Ecuaciones Diferenciales Vol. II,
Metodos Computacionales en Ciencia e Ingeniera. Edit. Reverte, Buenos
Aires, 1990.
[14] Murli M. A Survey of Some Second Order Difference Schemes for SteadyState Convection-Difussion Equation. SIAM, 1994, Vol.3 319-331.
[15] Mu
noz R. Introducci
on a los Metodos Numericos para Problemas de la
Fsica Matematica, Notas del Primer Seminario Internacional de Matematica Aplicada, Universidad de Lima, Lima, 1998.
[16] Lara, L. Metodo de Diferencias Finitas para la Solucion de las Ecuaciones
de Navier-Stokes. Tesis Maestria, UNT , 2000.
[17] Patankar, S. Numerical Heat Transfer and Fluid Flow. Edit. McGraw-Hill
Book Company, USA, 1990.
[18] Polanco, F. and May L. An Implicit Finite Difference Approximation
to the one-dimensional Transport Equation, ANZIAM J 42(E) pp.11791198, 2000.
[19] Quarteroni, V. A Numerical Approximation of Partial Differencial.
Springer-Verlag, U.S.A, 1994.
[20] Smith, G. Numerical Solution of Partial Diffrential Equations: Finite Difference Methods, Oxford University Press, 1990.
[21] Strange, G. Introductiona to Applied Mathematics , Wellesley-Cambridge
Press, Massachusetts, 1990.
[22] Striwerda, J. Finite Difference Schemes and Partial Differential Equations. University of Wisconsin-Madison. U.S.A, 1990.
[23] Sod, G. Numerical Methods in Fluid Dynamics: Initial and Initial
Boundary-Value Problems. Cambridge University Press. U.S.A, 1985.