• Author / Uploaded
• dora

Citation preview

Metodos Numericos para Ecuaciones Diferenciales Parciales y Dinamica de Fluidos Computacional Obidio Rubio Mercedes 28 de octubre de 1999

Pr´ ologo En la actualidad se est´a resolviendo problemas cuyos modelos matem´aticos son cada vez mas complejos, esta necesidad de resolverlas, lleva consigo la b´ usqueda de implementar alguna metodolog´ıa, en este caso los m´etodos num´ericos, resultaron ser los m´as u ´tiles, puesto que siempre generan un c´odigo o programa para ser implementado en un computador. La complejidad del problema exigir´a muchas veces un computador de alto desempe˜ no, pero cuyos resultados son indudablemente muy impactantes, en la actualidad por que permiten simular los problemas reales sin tener que hacer gasto de dinero o tiempo en simulaciones experimentales, generando una a´rea muy amplia, como es la matem´atica num´erica, entendi´endose que en ella se trata tanto los m´etodos num´ericos como el an´alisis num´erico de estos m´etodos , los cuales siempre estar´an enfocados en analizar el orden de consistencia y la velocidad de convergencia de las soluciones aproxmadas que generan los m´etodos. En el presente trabajo se hace una introducci´on al m´etodo de diferencas finitas, puesto que es uno de los m´etodos muy usados, por su simplicidad, puesto que puede ser accesado por un principiante con bastante ´exito, y por su capacidad de abordar una amplia variedad de problemas, esta presentaci´on se hace utilizando como modelos b´asiocs a ecuaciones diferenciales parciales lineales de primer orden, como la ecuaci´on de la onda y de segudno orden com la ecuaci´on del calor y la ecuci´on de Poisson. Teniendo en cuenta que un problema de ecuaciones diferenciales parciales, el cual es una abstracci´on de los modelos de par´ametros distribuidos, para ser formulado siempre necesita de un espacio infinito dimensional, entonces este problema debe ser aproximado por otro problema formulado en un espacio finito dimensional. El proceso principal de esta tarea es la discretizaci´on de un problema, que consiste en aproximar un problema de ecuaciones diferen1

2

ciales parciales, por un problema discreto, para ello se usa el proceso de discretizaci´on, el cual se ha sistematizado en los siguientes pasos: Discretizaci´on del dominio del problema, discretizaci´on de la variable o inc´ognita del problema y finalmente la discretizaci´on de la ecuaci´on diferencial del problema que consiste en usar un esquema de diferencias apropiado que aproxime la ecuaci´on. Cada vez que se hable de discretizaci´on de un problema se hace ´enfasis en estos tres pasos. El trabajo est´a organizado de la siguiente manera: se presenta una dscusi´on de los esquemas de diferencias finitas para problemas de ecuaciones diferenciales de primer orden, cuyo modelo patr´on es la ecuaci´on hiperb´olica de la onda escalar unidimensional, aqu´ı se introduce los conceptos fundamentales del an´alsis de los esquemas de diferencias finitas, como es el de convergencia, consistencia y estabilidad, luego se utlzan estos conceptos en un estudio detallado para las ecuaciones parab´olicas, cuyo modelo es la ecuaci´on del calor. Utilizamos los problemas el´ıpticas para hacer una presentaci´on simple de ellas con ´enfasis en el tratamiento de las condiciones de contorno, finalmente se hace una ligera introducci´on a temas m´as espec´ıficos, donde se necesita de mas agudeza en el estudio de la convergencia, tanto en problemas lineales como en problemas no lineales.

´Indice General 1 Introducci´ on 1.1 La naturaleza de los m´etodos num´ericos . . . . . . . 1.2 El concepto de discretizaci´on . . . . . . . . . . . . . . 1.3 Formulaci´on para derivar m´etodos de discretizaci´on . 1.3.1 Formulaci´on Serie de Taylor . . . . . . . . . . 1.3.2 Formulaci´on Variacional . . . . . . . . . . . . 1.3.3 Formulaci´on de Pesos Residuales . . . . . . . 1.3.4 Formulaci´on de Volumen de Control . . . . . . 1.3.5 Formulaci´on Integrales de Contorno y Teor´ıa Potencial . . . . . . . . . . . . . . . . . . . . 1.4 M´etodo de Diferencias Finitas . . . . . . . . . . . . . 1.4.1 Discretizaci´on del Dominio . . . . . . . . . . . 1.4.2 Discretizaci´on de la Variable . . . . . . . . . . 1.4.3 Discretizaci´on de la Ecuaci´on Diferencial . . . 1.4.4 Operadores de Diferencias . . . . . . . . . . . 1.5 Precisi´on de Aproximaciones . . . . . . . . . . . . . . 2 Ecuaciones Diferenciales Parciales 2.1 Ecuaci´on de la Onda Unidimensional . . . . 2.1.1 Sistemas de Ecuaciones Hiperb´olicas 2.1.2 Condiciones de contorno . . . . . . . 2.2 Esquemas de Diferencias Finitas . . . . . . . 2.2.1 Discretizaci´on del Dominio . . . . . . 2.2.2 Discretizaci´on de la Variable . . . . 2.2.3 Discretizaci´on de la Ecuaci´on . . . . 2.3 An´alisis de Esquemas de Diferencias Finitas 2.3.1 Convergencia y Consistencia . . . . . 2.3.2 Estabilidad . . . . . . . . . . . . . . 3

. . . . . . . . . .

. . . . . . . . . .

. . . . . . . . . .

. . . . . . . . . .

. . . . . . . . . .

7 7 7 8 8 8 9 9

. . . . . . . . . . . . . . de . . . . . . . . . . . . . .

10 11 11 12 13 16 16

. . . . . . . . . .

19 19 21 22 23 23 24 24 29 30 34

. . . . . . . . . .

´INDICE GENERAL

4

2.4

2.3.3 El Teorema de Equivalencia de Lax - Richtmyer . 2.3.4 La Condici´on de Courant-Friedrichs-Lewy . . . . An´alisis de Fourier . . . . . . . . . . . . . . . . . . . . .

3 Ecuaciones Diferenciales Parab´ olicas 3.1 Problema de Conducci´on del Calor . . . . . . . . . . 3.2 An´alisis de los Esquemas de Diferencias Finitas . . . 3.2.1 Teorema de Lax . . . . . . . . . . . . . . . . . 3.2.2 El esquema de Euler retrasado . . . . . . . . . 3.2.3 Consistencia y orden de precisi´on . . . . . . . 3.2.4 Estabilidad . . . . . . . . . . . . . . . . . . . 3.2.5 Esquema de Crank-Nicolson . . . . . . . . . . 3.2.6 Consistencia del esquema de Crank-Nicolson . 3.2.7 Estabilidad del esquema de Crank-Nicolson . . 3.2.8 Esquema de Du Fort-Frankel . . . . . . . . . . 3.2.9 Consistencia del esquema de Du Fort-Frankel 3.2.10 Estabilidad del Esquema de Du Fort-Frankel . 3.3 Implementaci´on Num´erica . . . . . . . . . . . . . . .

36 37 39

. . . . . . . . . . . . .

49 49 58 58 59 60 62 63 64 66 67 68 70 72

4 Disipaci´ on y Dispersi´ on 4.1 Disipaci´on . . . . . . . . . . . . . . . . . . . . . . . . . . 4.2 Dispersi´on . . . . . . . . . . . . . . . . . . . . . . . . . . 4.3 Velocidad de grupo y propagaci´on de paquetes de ondas .

81 81 83 85

5 Ecuaciones El´ıpticas: Condiciones de contorno 5.1 Ecuaci´on de Poisson . . . . . . . . . . . . . . . . . . . . 5.1.1 Caso unidimensional . . . . . . . . . . . . . . . . 5.1.2 Caso bidimensional . . . . . . . . . . . . . . . . . 5.2 Tratamiento de las condiciones de contorno . . . . . . . . 5.2.1 M´etodos de los puntos ficticios . . . . . . . . . . . 5.2.2 Aproximaci´on por series de Taylor (sin puntos ficticios) . . . . . . . . . . . . . . . . . . . . . . . . 5.2.3 Fronteras irregulares y curvadas . . . . . . . . . .

87 87 87 91 93 95

6 Convergencia: problemas lineales y no lineales 6.1 Estabilidad y Convergencia . . . . . . . . . . . . . . . . . 6.1.1 Discretizaci´on del Dominio . . . . . . . . . . . . .

99 99 99

. . . . . . . . . . . . .

96 97

´INDICE GENERAL

6.1.2 Discretizaci´on de la Varialble . . . . . 6.1.3 Discretizaci´on de la Ecuaci´on . . . . . 6.2 Consistencia y orden de precisi´on . . . . . . . 6.3 Problemas no lineales: Ecuaci´on de Burger’s . 6.3.1 Esquemas de diferencias finitas . . . . 6.3.2 Consistencia . . . . . . . . . . . . . . . 6.3.3 Estabilidad . . . . . . . . . . . . . . . 6.3.4 Implementaci´on num´erica . . . . . . .

5

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

100 100 106 108 108 109 109 115

6

´INDICE GENERAL

Cap´ıtulo 1 Introducci´ on 1.1

La naturaleza de los m´ etodos num´ ericos

La soluci´on num´erica de una ecuaci´on diferencial consiste de un conjunto de n´ umeros a partir de los cuales la distribuci´on de la variable dependiente puede ser construida. En este sentido, un m´etodo num´erico es semejante a un experimento de laboratorio, en el cual el conjunto de instrumentos leidos nos permiten establecer la distribuci´on de la cantidad medible en el dominio de investigaci´on. Tanto el analista num´erico como el experimentador de laboratorio deben quedarse con unicamente un n´ umero finito de valores num´ericos como el resultado. Finalmente podemos decir que un m´etodo num´erico trata como sus incognitas b´asicas los valores de las variables dependientes en un n´ umero finito de localizaciones (llamados los puntos de la malla) en el dominio de c´alculo. El m´etodo incluye la tarea de proveer un conjunto de ecuaciones algebraicas para estas incognitas y prescribir un algoritmo para resolver las ecuaciones.

1.2

El concepto de discretizaci´ on

Teniendo en cuenta que un problema de ecuaciones diferenciales parciales, el cual es una abstracci´on de los modelos de par´ametros distribuidos, para ser formulado siempre necesita de un espacio infinito demesional, entonces este problema debe ser aproximado por otro problema, discreto, y ser escrito en un espacio finito dimensional. Para ello se usa el proceso de discretizaci´on, que consiste en trasladar un problema infinito dimensional a otro finito dmensional, el cual se 7

´ CAP´ITULO 1. INTRODUCCION

8

ha sistematizado en los siguientes pasos: Discretizaci´on del dominio del problema, discretizaci´on de la variable o inc´ognita del problema y finalmente la discretizaci´on de la ecuaci´on diferencial del problema que consiste en usar un esquema de diferencias apropiado que aproxime la ecuaci´on. Cada vez que se hable de discretizaci´on de un problema se hace ´enfasis en estos tres pasos. Las tres etapas de este proceso lo comprenderemos como sigue: En primer lugar el dominio se subdivide en subdomios o elementos seg´ un una determinada prescripci´on, eligiendo un conjunto finito de puntos o nodos, que constituyen la malla, esta es la discretizaci´on del dominio. Focalizando la atenci´on en los valores de los puntos de la malla, el hecho de remplazar la informaci´on continua contenida en la soluci´on exacta de la ecuaci´on diferencial con valores discretos, asi hemos discretizado la variable dependiente Φ. Las ecuaciones algebraicas que involucran los valores de las variables de Φ en los puntos de la malla elegida se llaman ecuaciones de discretizaci´ on, las mismas que son derivadas de las ecuaciones diferenciales. La forma de derivar las ecuaciones algebraicas depende de como la varible dependiente varie en los puntos de la malla, generalmente se especifica una variaci´on de un determinado perfil.

1.3

Formulaci´ on para derivar m´ etodos de discretizaci´ on

Para una ecuaci´onn diferencial dada, los m´etodos de discretizaci´on requeridas pueden ser derivadas en muchas formas, a continuaci´on mencionamos algunos de los m´etodos mas comunes: 1.3.1

Formulaci´ on Serie de Taylor

El procedimiento usual para derivar las ecuaciones de diferencias finitas consiste en la aproximaci´on de las derivadas en la ecuaci´on diferencial via la serie de Taylor truncada. 1.3.2

Formulaci´ on Variacional

Otro m´etodo de obtener las ecuaciones de discretizaci´on estan basadas en el c´alculo de variciones, donde se verifica que, resolver una ecuaci´on

´ PARA DERIVAR METODOS ´ ´ 1.3. FORMULACION DE DISCRETIZACION

9

diferencial es equivalente a minimizar una funcional. La funcional es minimizada con respecto a los valores en los puntos de la malla de la variable dependiente, las condiciones resultantes generan las ecuaciones algebraicas requeridas. La formulaci´on variacional es muy comunmente empleada en m´etodos de elementos finitos. La principal dificultad de esta formulaci´on es su aplicabilidad ya que no todos las ecuaciones diferenciales de inter´es son interpretadas como minimizadoras de una funcional. 1.3.3

Formulaci´ on de Pesos Residuales

Un m´etodo muy potente de resolver las ecuaicones diferenciales es el m´etodo de los pesos residuales, El concepto b´asico es el siguiente: supo¯ ( la variable ner que la soluci´on aproximada de la ecuaci´on diferencial Φ dependiente) contiene n par´ametros por determinar, si sustituimos esta soluci´on en la Ecuaci´on diferencial obtendremos un residual. La idea es hacer este residual peque˜ no en algun sentido, generalmente esta u ´ltima tarea se hace multiplicando al residual por una funci´on peso W e la integraci´on sobre el dominio entero es puesta a cero. Por tanto eligiendon una sucesi´on de funciones peso, podremos generar todas las ecuaciones que se requieren para ecaluar los par´ametros. Muchas de las recientes desarrollos de las t´ecnicas de los elementos finitos estan basadas sobre perfiles por secciones usadas en conjunci´on con una clase de pesos residuales particulares conocidos como M´etodo de Galerkin 1.3.4

Formulaci´ on de Volumen de Control

El domino de c´alculo es dividido en un n´ umero de volumenes de control no traslapados de manera que existe un volumen de control encerrando un punto de la malla, la ecuaci´on diferencial es integrada sobre el volumen de control. Perfiles por secciones expresando la varici´on de la variable dependiente entre los puntos de la malla son usados para las integrales requeridas, el resultado es la ecuaci´on discretizada conteniendo los valores de la variable para un grupo de puntos de la malla. La ecuaci´on discretizada obtenida en esta manera expresa el principio de conservaci´on para Φ para el volumen de control finito tan igual como la ecuaci´on diferencial lo hace para un volumen infinitesimal.

10

1.3.5

´ CAP´ITULO 1. INTRODUCCION

Formulaci´ on Integrales de Contorno y Teor´ıa de Potencial

Este m´etodo, el mas moderno, como es natural, est´a basado en una teor´ıa mucho mas compleja, como es la Teor´ıa de Potencial, que consiste en formular las soluciones de ecuaciones dieferenciales como operadores de integrales de contorno, muchas veces singulares, de modo que la discretizaci´on ya no es en todo el dominio, sin´o unicamente en el contorno y la ecuaci´on se discretiza siguiendo la t´ecnica de Galerkin y de los elementos finitos, pero con una funci´on 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 m´etodos num´ericos para las ecuaciones diferenciales parciales esta dedicado al estudio de los esquemas de diferencias finitas. Nuestra intenci´on es hacer nuestro ingreso al campo de la Din´amica de Fluidos Computacional(DFC) de una manera accesible, desde luego que no dejamos de mencionar otros m´etodos que son importantes y que ser´an tratados en su debido tiempo, tales como los elementos finitos, m´etodos 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 fen´omenos naturales, mec´anicos y otros tipos de procesos que ocurren en las ciencias e Ingenier´ıa. Suponemos un conocimiento b´asico de las ecuaciones diferenciales parciales asi como algunos conocimientos del an´alisis, 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 [7], Farlow[1] y Tijonov [6] Presentamos primeramente las ecuaciones diferenciales parciales de tipo hiperb´olico, por considerar que dentro de ellas se encuentran ecuaciones simples, sin embargo modelan fon´omenos f´ısicos concretos. Despues de esta r´apida presentaci´on comezamos a presentar los conceptos b´asicos de los esquemas de diferencias finitas(EDF). Tambi´en presentamos los conceptos importantes de convergencia, consistencia y estabilidad los cuales son relacionados por el teorema de equivalencia de

´ 1.4. METODO DE DIFERENCIAS FINITAS

11

Lax-Richtmyer, utilizamos el an´alisis de Von Neumann para establecer la estabilidad de esquemas, y terminamos con un tratamiento de los esquemas de diferencias finitas para las ecuaciones de tipo parab´olico.

1.4

M´ etodo de Diferencias Finitas

En esta secci´on presentamos en forma general, y b´asica el concepto del m´etodo de diferencias finitas (MDF), por simplicidad lo presentamos en dimensi´on uno. Para esto , supongamos que tenemos un dominio, el cual ser´a un intervalo Ω = [a, b] ⊂ R donde se tiene que resolver una ecuaci´on diferencial. A continuaci´on hacemos el proceso de discretizaci´on en las tres etapas fundamentales.

1.4.1

Discretizaci´ on del Dominio

Por simplicidad consideremos el dominio como el intervalo Ω = [0, 1]. Este dominio es dividido en un n´ umero finito de subdominios o subintervalos igualmente espaciados de longitud Δx, el cual es un n´ umero positivo. Ahora damos la definici´on de malla Definici´ on 1.4.1. Sean Δx un n´ umero positivo, llamado la diferencia finita; una malla en el dom´ınio Ω ⊂ IR es un conjunto de puntos τ = {xn } ⊂ Ω, tal que cada xn esta en el interior o en la frontera de un subdominio y satisface xn+1 = xn + Δx xn−1 = xn − Δx es decir se cumple Δx = xn − xn−1 Esto puntos tambi´en se llaman nodos y el forma de elecci´on no es u ´nica, la obtenci´on de los subdominios junto con la elecci´on de los nodos se denomina discretizaci´on del dom´ınio Ω. En la figura 1.1 se presenta al punto xn de la malla, localizado entre los puntos xn−1 y xn+1 .

´ CAP´ITULO 1. INTRODUCCION

12

 f

xn−1

g

xn

Δx



 g

xn+1

Figura 1.1: malla unidimensional

Una de las elecciones mas comunes para los nodos son aquellos ubicados en la frontera del subdominio, es decir, sus coordenadas son de la forma n = 0, 1, 2, . . . , N xn = nΔx, , siendo N el n´ umero de subdominios. Pero tambi´en se puede elegir en el interior del subdominio, por ejemplo, si elegimos en el centro de cada subintervalo se tendr´ıa   1 Δx, n = 0, 1, 2, . . . , N − 1. xn = n + 2 1.4.2

Discretizaci´ on de la Variable

Definici´ on 1.4.2. Una funci´on discreta v, es aquella que est´a definida en una malla τ , tal que a cada punto xn le asocia un n´ umero vn . Observaci´ on 1.4.1. 1 .- Si tenemos una funci´ on continua u(t) definida en Ω, se le puede asociar una funci´ on discreta, es decir u puede ser discretizada sobre la malla, definiendo un = u(xn ). 2 .- Cuando la subdivisi´on del dominio no es uniforme, lo cual ocurre en aplicaciones pr´ acticas, se tiene una sucesi´ on de diferencias Δxn = xn+1 − xn para cada n. 3 .- Para efectos de un an´ alisis general, suponemos que vn = 0 para / Ω, denotando |n| > N , donde N es suficientmente grande tal que xn ∈ en este caso a la funci´on discreta por v = {vn } y n ∈ Z, es decir: v = {vn } = (· · · , v−2 , v−1 , v0 , v1 , v2 , · · · )

´ 1.4. METODO DE DIFERENCIAS FINITAS

13

Definici´ on 1.4.3. Sean f = {fn }, g = {gn } dos funciones discretas definidas sobre una malla τ , el producto interno de f y g se define por

 f, g h = hfn gn (1.1) n

donde h es un factor de escala, generalmente asociado al tama˜ no de la malla h = maxn {xn+1 − xn }. De la definici´on anterior se obtiene la L2 -norma o norma euclideana: 1/2

 1/2 f  = f, f h = hfi2 . (1.2) i

Tambi´en existen otras normas para funciones discretas, tales como: La norma de la suma

f  == h|fi | (1.3) i

La norma del m´aximo |f |∞ = max h|fi | i

1.4.3

(1.4)

Discretizaci´ on de la Ecuaci´ on Diferencial

El m´etodo de diferencias finitas, tambi´en llamado m´etodo taylor, porque resulta cunado se usa la formulaci´on serie de Taylor para realizar las aproximaciones de las derivadas en la ecuaci´on diferencial por diferencias finitas. Antes de determinar la aproximaci´on de las derivadas, es util revisar el significado de derivada. Las derivadas representan tasas de cambio, una derivada contiene informaci´on acerca de la variaci´on local de una funci´on. Por tanto, una aproximaci´on en diferencia a una derivada deberia atender a estos significados acoplando la informaci´on en los puntos vecinos de una malla, este acoplamiento se puede realizar con la serie de Taylor de la funci´on en una vecindad del punto. Asi la expansi´on de u(x + Δx) alrededor del punto x es d2 u (Δx)2 d3 u (Δx)3 du + 3 + ... u(x + Δx) = u(x) + Δx + 2 dx dx 2! dx 3!

(1.5)

´ CAP´ITULO 1. INTRODUCCION

14

utilizando la nomenclatura para la versi´on discreta, esta serie de Taylor puede ser escrita como 

un+1

du = un + dx





d2 u Δx + dx2 n



 3  du (Δx)2 (Δx)3 + + . . . (1.6) 3 2! dx 3! n n

De la misma forma, el valor u(x − Δx) puede ser escrito como un−1 y expandido como  un−1 = un −

du dx



 Δx +

n

d2 u dx2



(Δx)2 − 2! n



d3 u dx3



(Δx)3 + . . . (1.7) 3! n

Cuando la ecuaci´on (1.6) se resuelve para la primera derivada , obtenemos  2   3    du du (Δx) (Δx)2 un+1 − un du − + + ... = dx n Δx dx2 n 2! dx3 n 3! un+1 − un = + O(Δx) (1.8) Δx donde O(Δx) denota el t´ermino que contiene la primera y las potencias mayores de Δx. Si esta serie es truncada despues del primer t´ermino, lo que tenemos es una diferencia finita que aproxima la primera derivada en el punto n, esto es 

du dx

 ≈ n

un+1 − un Δx

(1.9)

El error (de truncamiento) que resulta de esta aproximaci´on se dice que es del orden de (Δx), porque (Δx) aparece en el primer t´ermino de la serie truncada. Por lo tanto el error de truncaci´on en la ecuaci´on (1.9), se dice que es primer orden, y por tanto, la (1.9) misma se dice que es exacta de primer orden. Observamos que la precedente aproximaci´on, se realiza en el punto n y utiliza el punto n + 1, por lo que se le llama una aproximaci´on de diferencia hacia adelante. De la misma manera, podemos derivar la aproximaci´on de la diferencia hacia atras de la ecuaci´on (1.7) como

´ 1.4. METODO DE DIFERENCIAS FINITAS



du dx

 n

15

 2   3  un − un−1 (Δx) (Δx)2 du du = + − + ... Δx dx2 n 2! dx3 n 3! un − un−1 + O(Δx) = Δx un − un−1 ≈ (1.10) Δx

la cual es tambi´en exacta de primer orden. Una aproximaci´on de diferencia mas precisa que la de primer orden, es la obtenida substrayendo la ecuaci´on (1.7) de la ecuaci´on (1.6 ) y du/dx, entonces   un+1 − un−1 du + +O(Δx)2 = dx n 2Δx un+1 − un−1 ≈ (1.11) 2Δx ´esta aproximaxci´on es llamada diferencia central porque esta centrada en el punto n. Las diferencias hacia adelante y hacia atras tambi´en pueden ser pensadas como diferencias centrales, pero no centradas en el punto n, sino en los puntos (n + 1/2) y en (n − 1/2) respectivamente. Para obtener una aproximaci´on de diferencia finita de la segunda derivada, las ecuaciones (1.7) y (1.6 ) son sumadas, y el resultado es resuelta para (d2 u/dx2 ), generando 

d2 u dx2



un+1 − 2un + un−1 + O(Δx)2 2 Δx

= n

(1.12)

El primer t´ermino del lado derecho, es usado para una aproximaci´on de diferencia finita para la segunda derivada, es decir; 

d2 u dx2

 ≈ n

un+1 − 2un + un−1 , Δx2

(1.13)

Esta aproximaci´on de la segunda derivada se llama diferencia central segunda, como observamso de (3.4), el error de truncamiento de esta aproximaci´on es de segundo orden.

´ CAP´ITULO 1. INTRODUCCION

16

1.4.4

Operadores de Diferencias

En base a las diferencias presewntadas anteriormente se definimos algunos operadores mas utilizados en lenguaje de los esquemas de diferencias finitas, que se pueden aplicar a funciones discretas 1. Operadores de desplazamiento (derecha, izquierda) S+ Φn = φn+1 ,

S− Φn = Φn−1 ,

2. Operador de diferencia hacia adelante D + Φn =

Φn+1 − Φn (S+ − I)Φn = Δx Δx

3. Operador de diferencia hacia atras D − Φn =

(I − S− )Φn Φn − Φn−1 = Δx Δx

4. Operador de diferencia central D 0 Φn =

(S+ − S− )Φn Φn+1 − Φn−1 = 2Δx 2Δx

Tambi´en existen otros tipos de operadores de mayor orden, que se pueden construir usando los operadores anteriores, tales como: Operador de diferencia central D− Φn+1 − D− Φn Δx Φn+1 − 2Φn + Φn−1 = Δx2

δ02 Φn = D− D+ Φn =

y se cumple D− D+ Φn = D+ D− Φn

1.5

Precisi´ on de Aproximaciones

Jhon Von Neumann y H.H. Goldstein han identificado cuatro principales fuentes de errores que son casi inevitables cuando estamos describiendo sistemas f´ısicos.

´ DE APROXIMACIONES 1.5. PRECISION

17

1. Errores de modelaci´ on: Los modelos matem´aticos de por s´ı son generalmente dise˜ nados con varias hip´otesis. Un sistema que es no lineal y dependiente del tiempo, es talv´ez representado por una ecuaci´on diferencial lineal con coeficientes constantes. Un sistema de par´ametros distribuidos (representado por una ecuaci´on diferencial parcial) puede ser aproximado por un modelo de par´ametros agrupados(sistema de ecuaciones diferenciales ordinarias). 2. Errores de medici´ on: La mayor´ıa de las descripciones de fen´omenos naturales requieren de datos medidos, estos datos de por s´ı tienen un error porque dependen de los equipos y otros factores que siempre aparecen en la obentci´on de ellos. 3. Errores de truncaci´ on: Casi todos las descripciones de problemas de campo involucran un continuo infinito. En el contexto del an´alisis num´erico, solo se puede tratar con un n´ umero finito de t´erminos de procesos limitantes los cuales son normalmente descritos por series infinitas. 4. Errores de redondeo: Errores introducidos por eliminaci´on de d´ıgitos; cuando los c´alculos son realizados, estos pueden ser hechos unicamente en un computador de precisi´on finita. Los errores aparecen en el tama˜ no limitado de los registros en la unidad de aritm´etica del computador. De los cuatro errores presentados, los errores de truncaci´on y redondeo caen en el campo del analista num´erico, estos errores son acumulados y tiene que controlarse para no dejarlos crecer. Aqu´ı se observa un hecho muy importante, Cuando los c´alculos usan una malla mas fina, el error de truncamiento decrece pero el error de redondeo crece debido a que los datos num´ericos son mas peque˜ nos y por el n´ umero creciente de c´alculos efectuados, por tanto hay que saber controlar ambos errores a fin de que el error total permanezca en rangos permisibles.

18

´ CAP´ITULO 1. INTRODUCCION

Cap´ıtulo 2 Ecuaciones Diferenciales Parciales 2.1

Ecuaci´ on de la Onda Unidimensional

Por su simplicidad, el prototipo de todas las ecuaciones diferenciales parciales es la ecuaci´on hiperb´olica de la onda unidimensional, ut + aux = 0,

(2.1)

donde a es una constante, t representa el tiempo, x la variable espacial, u(x, t) es la incognita de la ecuaci´on. El problema de valor inicial asociado a esta ecuaci´on es el siguiente: Dada u en el tiempo inicial, es decir u(x, 0) es igual a una funci´on dada u0 (x), para todo x ∈ R, deseamos determinar el valor u(x, t) para cada t > 0, tal que:  −∞ < x < ∞, t > 0 ut + aux = 0, (2.2) (P V I) −∞ < x < ∞. u(x, 0) = u0 (x), La soluci´on de este problema es facil determinar, por observaci´on podemos decir que la soluci´on de (PVI) es: u(x, t) = u0 (x − at),

(2.3)

La funci´on (2.3), soluci´on de (PVI) puede interpretarse de la siguiente manera: 1. En cualquier tiempo t0 > 0, la soluci´on es una r´eplica de la funci´on inicial, pero desplazada, a la derecha si a > 0 y a la izquierda si a < 0, hasta una posici´on |a|t0 . 2. La soluci´on en (x, t) depende unicamente del valor ζ = x − at. 19

CAP´ITULO 2. ECUACIONES DIFERENCIALES PARCIALES

20

3. Las rectas x − at = cte en el plano (x, t) se llaman caracter´ısticas de la ecuac´on. 4. El par´ametro a tiene dimensi´on de distancia dividido por el tiempo, y se llama velocidad de propagaci´on a lo largo de la caracter´ıstica. 5. Si la condici´on inicial es una onda, la soluci´on expresada en (2.3) dice que en cualquier tiempo la onda inicial se propaga con velocidad a sin cambiar de forma. Observaci´ on 2.1.1. La ecuaci´on (2.1) tiene sentido si u es diferenciable, la soluci´ on (2.3) no necesita de la diferenciabilidad de u0 . De acuerdo a esto, es permisible tener soluciones discontinuas para problemas hiperb´ olicos. Un ejemplo de una soluci´on discontinua es una onda de Shock(ver Farlow [1]), la cual es una clase de soluciones de ecuaciones hiperb´ olicas no lienales. Una generalizaci´on del problema hiperb´olico es el siguiente  ut + aux + bu = f (x, t), −∞ < x < ∞, t > 0 (P V I) −∞ < x < ∞. u(x, 0) = u0 (x),

(2.4)

donde a, b son constantes f es una funci´on. Basados en la observaciones anteriores, hacemos el cambio de variables   τ = t, t = τ, ´o (2.5) ζ = x − at, x = ζ + aτ. luego definiendo u¯(ζ, τ ) = u(x, t), donde (ζ, τ ) y (x, t) est´an relacionadas por el cambio de variable definida previamente, entonces la ecuaci´on en (2.4), se transforma en ∂t ∂x ∂ u¯ = ut + ux = ut + aux = −bu + f (ζ + aτ, τ ) ∂τ ∂τ ∂τ es decir

∂ u¯ = −b¯ u + f (ζ + aτ, τ ) ∂τ la cual es una ecuaci´on diferencial ordinaria en τ y cuya soluci´on depende de ζ y es de la forma  τ f (ζ + as, s)e−b(τ −s) ds, (2.6) u(x, t) = u0 (x − at)e−bτ + 0

´ DE LA ONDA UNIDIMENSIONAL 2.1. ECUACION

21

Usando la inversa del cambio de variable se tiene  t f (x − a(t − s), s) + as, s)e−b(τ −s) ds, (2.7) u(x, t) = u0 (x − at)e−bτ + 0

Notar que u(x, t) depende solamente de los valores de (x1 , t1 ) tal que x − at = x1 − at1 , es decir, la soluci´on depende unicamente de los valores de u y f sobre las caracter´ısticas que pasan por (x, t), para 0 ≤ t1 ≤ t. Este m´etodo se soluci´on puede extenderse facilmente a ecuaciones no lineales de la forma ut + aux = f (x, t, u). Tambi´en el m´etodo de las caracter´ısticas puede ser util para interpretar, aunque no expresar explicitamente la soluci´on de del problema no lineal, ut + g(u)ux = 0 Ya que se puede conocer la soluci´on sobre cada curva particular, pero no en todas al mismo tiempo(ver [1]). 2.1.1

Sistemas de Ecuaciones Hiperb´ olicas

Ahora damos un ligero vistazo a los sistemas de ecuaciones hiperb´olicas de tipo parab´olico en una dimensi´on, aunque la incognita u es un vector n-dimensional. Definici´ on 2.1.1. Un sistema de la forma ut + Aux + Bu = F (x, t), es hiperb´ olico si la matr´ız A es diagonalizable, es decir existe una matr´ız no singular P tal que P AP −1 = D donde D = diag(a1 , a2 , . . . , an ). Los autovalores de A son las velocidades caracter´ısticas del sistema. Bajo el cambio de variables w = P u, y poniendo B = 0 en el sistema, tenemos wt + Dw : x = P F (x, t) = F˜ (x, t), ´o wti + ai wxi = f˜i (x, t), las cuales son de la forma (2.4). Podemos decir que si la matriz B = 0, el sistema hiperb´olico unidimensional se reduce a un conjunto de n ecuaciones hiperb´olicas independientes, si B = 0 el sistema permanece acoplado.

22

CAP´ITULO 2. ECUACIONES DIFERENCIALES PARCIALES

Ejemplo 2.1.1. El sistema  ut + 2ux + vx = 0 vt + ux + 2vx = 0 se puede escribir como      2 1 u u + = 0. 1 2 v x v t 2.1.2

Condiciones de contorno

Si consideramos las EDPs sobre un intervalo finito en lugar de la recta real, entonces estamos frente a un problema que necesita imponer algunas condiciones correctas en los extremos del intervalo. La mayoria de las aplicaciones de las EDPs estan definidas en dominios con frontera. Las condiciones que relacionan la soluci´on de la ecuaci´on diferencial a datos de frontera lo llamamos, condiciones de contorno. El problema de determinar una soluci´on a la ecuaci´on diferencial con condici´on inicial y condici´on de contorno se llama problema de valor inicial y de contorno. La discusi´on de problemas de valores de contorno sirve para ilustrar de nuevo la importancia del concepto de caracter´ısticas, por ejemplo si consideramos el problema ⎧ 0 ≤< x < 1, t > 0 ⎨ ut + aux = 0, −∞ < x < ∞. (2.8) u(x, 0) = u0 (x), (P V IC) ⎩ u(0, t) = g(tx), 0 < t ≥ 0 < x < 1. Si a > 0, las caracter´ısticas en esta regi´on se propaga de izquierda a derecha segun la figura, examinando las mismas vemos que, es suficiente que se especifique, los datos iniciales y en el extremo x = 0, para tener la soluci´on en cada tiempo t > 0. Mas a´ un no es necesario tener datos en el otro extremo x = 1, lo cual convertir´a al problema en sobre determinado. Por ejemplo consideremos la condici´on inicial u(x, 0) = u0 (x) y la condici´on de contorno u(0, t) = g(t), la soluci´on del problema se puede escribir en la forma  u0 (x − at) si x − at > 0, u(x, t) = g(t − a−1 x) si x − at < 0

2.2. ESQUEMAS DE DIFERENCIAS FINITAS

23

t 6 >



> >

0

-

1

x

A lo largo de la caracter´ıstica x − at = 0 existir´a un salto de discontinuidad en u si u0 (0) = g(0). Si a < 0 el rol de las fronteras son intercambiados.

2.2

Esquemas de Diferencias Finitas

Para comenzar una discusi´on de esquemas de diferencia finitas consideramos el problema a tratar como un problema de valor inicial escalar  P(x, t, ∂t , ∂x )u = 0, x ∈ R, t > 0 (P V I) u(0) = u0 . (2.9) donde P(ξ1 , ξ2 , ξ3 , ξ4 ) es un polinomio en las variables ξ1 , ξ2 , ξ3 , ξ4 . Como sabemos la discretizaci´on lo hacemos en tres pasos 2.2.1

Discretizaci´ on del Dominio

Sabiendo que el dominio del problema es el semiplano superior del en el sistema (x, t), por tanto, la discretizaci´on del dominio consiste en dividirlo en subdominios y definir una malla de puntos en el semiplano (x, t), t > 0. En nuestro caso, por razones pr´acticas, los subdominios ser´an rect´angulos homog´eneos y en el extremo de cada rect´angulo se define un nodo como sigue Definici´ on 2.2.1. Sean h,k n´ umeros positivos; una malla es un conjunto

24

CAP´ITULO 2. ECUACIONES DIFERENCIALES PARCIALES

de puntos de la forma (xj , tn ) = (jh, nk), donde j, n son n´ umeros enteros arbitrarios con n > 0. 2.2.2

Discretizaci´ on de la Variable

Definici´ on 2.2.2. Una funci´on discreta v, es aquella que esta definida en una malla, tal que a cada punto (xj , tn ) le asocia un valor vjn . Observaci´ on 2.2.1. 1 .- Una funci´on continua u(x, t) definida en −∞ < x < ∞, puede ser discretizada sobre la malla, definiendo

t > 0,

unj = u(xj , tn ) 2 .- El conjunto de puntos de la malla (xj , tn ), con n fijo se llama nivel n de la malla. La discretizaci´on de la variable del problema consiste en utilizar un funci´on discreta vjn que aproxime la variable del problema, es decir a la soluci´on exacta u(x, t) en cada punto de la malla (jh, nk) vjn ≈ u(xj , tn )

(2.10)

En esta misma etapa tambi´en se discretiza las condciones de cotnorno as´ı com la condici´on inicial, por ejemplo en nuestro problema se tiene vj0 = u0 (xj ) 2.2.3

Discretizaci´ on de la Ecuaci´ on

En esta etapa, los operadores diferenciales son aproximados por operadores de diferencias finitas utilizando la serie de taylor, estos operadores son combinados adecuadamente formando un esquema de diferencias finitas, tal que el problema de valor inicial (P V I) es sustitu´ıdo por el problema de valor inicial discreto.  j ∈ Z, n ∈ Z+ Ph,k vjn = 0, (2.11) (P V ID) vj0 = u0 (jh), j ∈ Z. Con la utilizaci´on del operador dscreto de Ph,k el problema de valor inicial continuo definido para −∞ < x < +∞, t > 0, ha sido sustituido

2.2. ESQUEMAS DE DIFERENCIAS FINITAS

25

por un problema discreto definido en puntos discretos x = jh, t = nk, j ∈ Z, n ∈ Z+ . ´nico y dependiendo de la De esta manera el operador Ph,k , no es u forma en que se combinan los operadores de diferencias determina un esquema de diferencias finitas. La manera de conseguir estos esquemas es aproximando las derivadas por diferencias finitas, es fundamental en este proceso, la f´ormula de Taylor, esta t´ecnica tiene una justificaci´on que vemos a continuci´on. Suponiendo que u es suficientemente regular, la derivada de u con respecto a x en (x, t), se puede calcular con f´ormulas diferentes, como:  ∂u = lim−→0 u(x+,t)−u(x,t)  (x, t) (2.12) u(x+,t)−u(x−,t) ∂x = lim−→0 2 Si (x, t) es un punto de la malla esta derivada puede ser aproximada en los puntos de la malla por diferencias finitas de la forma  ∂u ≈ u((j+1)h,nk)−u(jh,nk) h (2.13) (jh, nk) u((j+1)h,nk)−u((j−1)h,nk) ∂x ≈ 2h respectivamente. Consideremos ahora, la ecuaci´on modelo ut + aux = 0,

(2.14)

Entre los Esquemas mas comunmente utilizados para discretizar esta ecuaci´on, tenemos: Upwind n vjn+1 − vjn − vjn vj+1 +a = 0, k h

a0

(2.16)

Leap-Frog, n n vjn+1 − vjn−1 − vj−1 vj+1 +a = 0, 2k 2h

(2.17)

CAP´ITULO 2. ECUACIONES DIFERENCIALES PARCIALES

26

Lax-Wendroff n n n vjn+1 − vjn − vj−1 v n − 2vjn + vj−1 vj+1 2 j+1 , +a = ka k 2h 2h2

(2.18)

Lax-Friedrichs n n n n + vj−1 ) vjn+1 − 12 (vj+1 − vj−1 vj+1 +a = 0. k 2h

(2.19)

En el esuqema Upwind , se presnetan dos esuqemas que depende del signo de a, el primero denotado por (FTFS) de Ingl´es ”forward time forward space”y el segundo se denota por (FTBS)”forward time Backward space”. En general el m´etodo para derivar esquemas de diferencias finitas es simple; sin embargo el an´alisis de que si un esquema de diferencias finitas es una buena aproximaci´on a la ecuaci´on diferencial necesita de una tarea matem´atica, a veces fuerte. Por otro lado, dado una lista de esquemas, es com´ un preguntarse, ¿cu´al de ellos son u ´tiles y cuales no? Para dar respuesta a esta pregunta b´asica, necesita algun tiempo y cuidado, y su respuesta se logra en varios aspectos. Comenzamos respondiendo en un primer aspecto, el cual ser´ıa , determinar cu´al esquema genera soluciones que aproximen a las soluciones de la ecuaci´on diferencial. En otro aspecto es determinar cu´ales esquemas son mas precisos que otros y finalmente investigamos la eficiencia de los esquemas. Cada uno de los esquemas antes mencionados pueden ser escritos expresando vjn+1 como combinaci´on lineal de los valores de v en los niveles n y n − 1. Por ejemplo el esquema (FTFS) se puede se escribe como: k . h La cantidad λ aparecera a menudo en el estudio de los esquemas para ecuaciones hiperb´olicas. Los esquemas permiten calcular la funci´on discreta v en el nivel n + 1 conociendo el valor de la misma en niveles inferiores, aquellos esquemas que solo necesitan un solo nivel inferior, por ejemplo del nivel n se llaman esquemas de un solo paso, entanto que los esquemas que involucran mas de un nivel inferior se llaman esuqemas de m´ ultiple paso, por ejemplo en la lista de esquemas anteriores, el esquema de Leap-Frog es de m´ ultiple n , vjn+1 = (1 + aλ)vjn − aλvj+1

con

λ=

2.2. ESQUEMAS DE DIFERENCIAS FINITAS

27

paso, usa el nivel n y el nivel n − 1 y todos los otros son esquemas de un solo paso. Dado los datos iniciales vj0 , un esquema de un solo paso puede utilizarse sin dfcultades para calcular vjn , para todos los valores positivos de n. Para los esquemas de m´ ultiple paso, como el esquema de Leapfrog no es suficiente tener los datos vj0 para obtener los demas valores en los otros niveles, sino que se tiene especificar v en los niveles suficientes a fin de que el esquema sea empleado, una forma de resolver esto es por ejemplo, usar un esquema de un solo paso hasta tener los niveles necesarios donde se comience a usar esquema de m´ ultiple paso. Para tener un panorama de los t´opicos a seguir discutiendo, hagamso un comentario de las resultados obtenidos con los esquemas propuestos anteriormente, para ello consideremos el (P V I), donde 0 < x < 2π, t > 0, teniendo como condici´on inicial la funci´on 2π-peri´odica seccionalmente continua ⎧ ⎨ 0, 0 ≤ x < 2π 3 , 2π (2.20) 1, 3 ≤ x ≤ 4π u(x, 0) = f (x) = 3 , ⎩ 0, 4π 3 < x ≤ 2π, A continuaci´on se muestra una secci´on de pseudoc´odigo, el cual nos permitir´a observar el comportamiento del periodo principal de esta funci´on considerado como una onda, para esto elijamos el esquema de LaxFriedrichs, los dem´as tienen similar estructura ⎧ Do n = 0, nmax ⎪ ⎪ ⎪ ⎪ ⎨ Do j = −19, 29 vj,n+1 ←− (vj+1,n + vj−1,n ) /2 − λ (vj+1,n − vj−1,n ) /2 ⎪ ⎪ End Do ⎪ ⎪ ⎩ End Do Como sabemos la soluci´on de este problema esta dado por u(x, t) = f (x − at) y es constante a lo largo de las caracter´ısticas x − at = cte en la figura 2.1, presentamos las soluci´ones num´ericas que son aproximaciones de la soluci´on del problema (2.14)-(2.20), con a = 1, h = 2π/240, k = 2h/3 en t = 500k y en t = 1000k.

CAP´ITULO 2. ECUACIONES DIFERENCIALES PARCIALES

28

upwind

Leap-Frog

1

1.4

0.9

1.2

0.8

1

0.7 0.8 0.6 0.6 0.5 0.4 0.4 0.2 0.3 0

0.2

-0.2

0.1 0 0

1

2

3

4

5

6

7

-0.4 0

1

2

Lax-Wendroff

3

4

5

6

7

5

6

7

Lax-Friedrich

1.2

1 0.9

1 0.8 0.8

0.7 0.6

0.6

0.5 0.4

0.4 0.3

0.2

0.2 0 0.1 -0.2 0

1

2

3

4

5

6

7

0 0

1

2

3

Figura 2.1: Soluciones num´ericas con los esquemas, en t = 500k,

4

´ 2.3. ANALISIS DE ESQUEMAS DE DIFERENCIAS FINITAS

29

La soluci´on exacta de la ecuaci´on diferencial est´a dada por la linea trazada y la soluci´on num´erica por el esquema respectivo, est´a dado por la linea cont´ınua en cada caso. En la figura podemos observar que el esquema de diferencias calcula la soluci´on adecuadamente bi´en, lejos de las discontinuidades, pero cerca de estas no se muestra una buena aproximaci´on a la soluci´on. Todas estas aproximaciones no son capaces de ser consideradas como buenas; porque, el m´etodo de Leap-frog genera muchas oscilaciones cerca de las discontinuidades. El m´etodo Lax-Wendroff tiene sobreoscilaciones y suboscilaciones cerca de las discontinuidades, pero no se propagan immediantamente. Los M´etodos de Lax-Friedrichs y Upwind, no tienen oscilaciones ni sobre oscilaciones, sin embargo cerca de las discontinuidades se separan mucho de la soluaci´on exacta, podr´ıamos decir que sus discontnuidades ocupan regiones mas extensas. Para tiempos bajos (primeros nivels de t) la soluci´on tiene un comportamiento aceptable, pero a medida que t crece la separaci´on de la soluci´on exacta es mas acentuada, especialmente en los esquemas Upwind y Lax-Fredrichs, porque estos esquemas de diferencias incorporan algunas propiedades que el problema original no lo tiene, como es el csao de la disipaci´on num´erica, como veremos en la pr´oxima secci´on.

2.3

An´ alisis de Esquemas de Diferencias Finitas

La serie de Taylor es una herramienta fundamental en los m´etodos num´ericos tanto para determinar los operadores de diferencias como para determinar los operadores de truncamiento de cualquier esquema de diferencias que aproxime una ecuaci´on diferencial. Convergencia es un concepto matem´atico familiar conocido para el caso de sucesiones de n´ umeros, sin embargo aqu´ı se refiere al hecho de que sucesiones de soluciones discretas obtenidas de las ecuaciones de diferencias, cada vez que la malla es refinada, se aproximan a la soluci´on exacta de la ecuaci´on diferencial continua. Por tanto hablaremos de convergencia, como el acercamiento, en algun sentido, de las soluciones vjn = vjn (h, k), cuando k, , h → 0. Desde luego que estrictamente hablando acerca de la convergencia de sucesiones de funciones, se nece-

30

CAP´ITULO 2. ECUACIONES DIFERENCIALES PARCIALES

sita definir el tipo de m´etrica con la que se puede medir la distancia entre funciones; sin embargo no necesitamos entrar en mucho detalle, es suficiente utilizar los puntos de la malla, por lo que hablaremos de convergencia puntual o convergencia uniforme y tambi´en por precisar algunos conceptos presentaremos esa misma definici´on en otro contexto. El t´ermino consistencia se refiere a una aproximaci´on entre el esquema mismo y la ecuaci´on diferencial. En un an´alisis de ecuaciones en diferenicas, la consistencia junto la establidad garantizan la convergencia, ´esto ser´a visto mas adelante en el teorema de Lax-Richtmyer. 2.3.1

Convergencia y Consistencia

Presentamos en primer lugar una ecuaci´on diferencial general, a fin de que estos conceptos sean presentados en un contexto te´orico y as´ı sean usados para en una amplia clase de aproximaciones de ecuaciones diferenciales. Consideremos la EDP lineal de la forma P (∂t , ∂x )u = f (x, t)

(2.21)

Donde P (ξ1 , ξ2 ) es un polinomio lineal en ξ1 , ξ2 , por razones de simplicidad, consideramos que la ecuaci´on anterior es de primer orden en la derivada con respecto al tiempo. Ejemplos de este tipo de ecuaciones son: ut − bux x + aux = 0, ut + butx + auxx + cuxxx = sin x. Supongamos, a fin de hacer la teor´ıa viable, que una ecuaci´on o sistemas de ecuaciones, de la forma (2.22), determina una u ´nica soluci´on si se especifica los datos iniciales u(x, 0). Para efectos de notaci´on, consideremos que un Esquema de Diferencias Finitas (EDF) tiene la forma Ph,k u = Bh,k f

(2.22)

Definici´ on 2.3.1. Dada la condici´on inicial u0 (x), para la cual existe una soluci´on u(x, t) de la EDP. Un esquema de diferencias finitas de un solo paso que aproxima la EDP se dice que es un esquema convergente

´ 2.3. ANALISIS DE ESQUEMAS DE DIFERENCIAS FINITAS

31

sii todas las soluciones vjn , dadas por el esquema de diferencias finitas, tal que vj0 converge a u0 (x) cuando jh −→ x; se cumple vjn −→ u(x, t) siempre que (jh, nk) −→ (x, t) cuando h, k −→ 0. Observar que: 1. Para el caso de esquemas de m´ ultiple paso, la definici´on es modificada precisando que. vji ,

converge a

u0 (xj ) para

0≤i≤J

2. Se debe precisar que la naturaleza de la convergencia de vjn a u(x, t) es en los puntos de la malla, ya que u est´a definida continuamente en el dominio que contiene a la malla, mientras que vjn est´a definida solo en los puntos de la misma. Probar que un esquema dado sea convergente, en general no es f´acil si intentamos hacerlo directamente, puesto que en general no se conoce la soluci´on de la EDP cont´ınua, sin embargo existen los conceptos t´ecnicos de consstencia y estabilidad, que son relativamente f´aciles de chequear que ayudar´an a determinar la convergencia sin tener la soluci´on antes mencionada; a contincuaci´on veremso la consistencia. Definici´ on 2.3.2. Dada una ecuaci´on diferencial parcial P u = f y un esquema de diferencias finitas Ph,k u = f ; decimos que el EDF es consistente con la EDP, sii para cualquier funci´on suave Φ(x, t) se cumple P Φ − Pn,k Φ −→ 0

cuando

h, k −→ 0

donde la convergencia es puntual, en cada punto de la malla. Observaci´on: 1. En algunos esquemas se necesita precisar la forma en que h y k tienden a cero a fin de obtener la consistencia. 2. Cuando nos referimos a una funci´on suave, se entiende que es una funci´on suficientemente diferenciable. 3. En el esquema de diferencias finitas, el operador de diferencias Ph,k , es el que caracteriza tal esquema, el cual puede estar definido en funci´on de los operadores discretos conocidos.

32

CAP´ITULO 2. ECUACIONES DIFERENCIALES PARCIALES

Presentamos, las definiciones anteriores en un contexto mas general, con el prop´osito de motivar a un estudio mas matem´atico de los esquemas de diferencias finitas. Definici´ on 2.3.3. Ph,k es consistente con P hasta un tiempo T en una cuando norma ||.||, sii ||Ph,k ϕ − P ϕ|| = kτ (h); y τ (h) → 0 k→0 τ (h) es el error de truncamiento local en tiempo nk . Definici´ on 2.3.4. Un esquema de diferencias finitas es convergente en la norma ||.||, sii ||unj − vjn || → 0 cuando h, k → 0 Es convergente de orden (p,q) en la norma ||.|| sii ||unj − vjn || = O(hp ) + O(k q ) Ejemplo: En el esquema de Lax-Friedrich, para la ecuaci´on de la onda, ∂ ∂ + a , o P Φ = Φt + aΦx P = ∂t ∂x − 12 (Φnj+1 + Φnj−1 ) Φn+1 Φnj+1 − Φnj−1 j +a , Ph,k Φ = k 2h

Φnj = Φ(jh, nk),

utilizando la serie de Taylor en el punto (xj , tn ) + + 1 1 Φn+ = Φnj − hΦx + h2 Φxx − h3 Φxxx + O(h4 ) j −1 2 6

de estas expresiones tenemos 1 1 n (Φj+1 + Φnj−1 ) = Φnj + h2 Φxx + O(h4 ) 2 2 Φnj+1 − Φnj−1 1 2 = Φx + h Φxxx + O(h4 ) 2h 6 k2 n+1 n Φj = Φj + kφt + Φtt + O(k 3 ). 2 Sustituyendo estas expresiones en el esquema tenemos 1 h2 1 2 h4 4 Ph,k Φ = Φt + aΦx + kΦtt − Φxx + ah Φxxx + O(h + + k2) 2 2k 6 k

´ 2.3. ANALISIS DE ESQUEMAS DE DIFERENCIAS FINITAS

33

podemos decir que, el esquema es consistente tan pronto como cuando h, k −→ 0, y se cumple

h2 k

−→ 0,

1 h2 1 h4 Ph,k Φ − P Φ = kΦtt − Φxx + ah2 Φxxx + O(h4 + + k 2 ) −→ 0, 2 2k 6 k lo cual significa que h2 debe ir mas rapidamente a cero que k. Ejemplo.- En la ecuaci´on ut + uk utilizamos el esquema FTFS n vjn+1 − vjn vj+1 − vjn + =0 k h que se puede escribir como

k n − vjn ) vjn+1 = vjn − (vj+1 h

···

(∗)

k h Como hicimos en el ejemplo anterior se puede verificar que n , vjn+1 = (1 + λ)vjn − λvj+1

λ=

1 1 Ph,k ϕ = ϕt + aϕx + kϕxt + ahϕxx + O(k 2 ) + O(h2 ) 2 2 Por lo tanto P ϕ − Ph,k ϕ = − 12 kϕtt + 12 ahϕxx + O(k 2 ) + O(h2 ) → 0, si (h, k) → 0 Por lo que el esquema es consistente. Es muy natural preguntarse, ¿cu´ando consistencia es suficiente para que que un esquema sea convergente?. En verdad, consistencia es necesaria para convergencia mas no suficiente, como veremos en el esquema del u ´ltimo ejemplo, el esquema puede ser consistente pero no convergente. Tomando la condici´on inicial  1, si −1 < x < 0 u0 (x) = 0, en otro lugar La soluci´on de la EDP u(x, t) = u0 (x − t) es un desplazamiento de u0 a la derecha por t.

CAP´ITULO 2. ECUACIONES DIFERENCIALES PARCIALES

34

En particular para t suficientemente grande existen x > 0 para los cuales u(x, t) = 0 Para el esquema de diferencias, los datos se obtienen por la discretizaci´on de u0 (x)  1, si −1 ≤ jh ≤ 0 vj0 = 0, en otro valor, en particular ∀j > 0 Seg´ un la ecuaci´on (∗), la soluci´on en (xj , tn ) depende solamente del valor en xj y xj+1 ; en el tiempo anterior, entonces, si j > 0,

n ≥ 0,

vjn = 0

∀j > 0,

n≥0

Por tanto vjn no puede converger a u(x, t), ya que para t > 0, x > 0, u no es id´enticamente cero, mientras que vjn lo ´es. 2.3.2

Estabilidad

En el ejemplo anterior vimos que para obtener convergencia, los esquemas deben satisfacer otras condiciones adem´as de consistencia, en esta secci´on nos referimos al importante concepto de estabilidad el cual es el paralelo de lo que, en l problema anal´ıtico es el concepto de problema bien puesto. Definici´ on 2.3.5. Un problema de valor inicial (P V 1) se dice que es bien puesto en una norma ||.|| si y solamente s´ı existe soluci´ on, es u ´nica y depende continuamente de los datos iniciales, es decir existe una constan-te C y α tal que ||u(x, t)|| ≤ Ceαt ||u0 (x)||

(a)

El concepto de estabilidad est´a relacionado con la u ´ltima desigualdad. Considerando el conjunto de todas las soluciones en un tiempo t, obtenidos mediante alg´ un esquema a partir de un mismo dato inicial (en t = 0). Este conjunto se obtuvo involucrando todas las relaciones que aparecen entre niveles de tiempo e intervalos de espacio, con ciertas restricciones. Entonces si los valores de las soluciones permanecen finitas sobre cada restricci´on que se consider´o, de manera que el conjunto de todas las posibles soluciones permanece acotado, se dice que el esquema es estable, dentro de las restricciones dadas. M´as precisamente, podemos enunciar la siguiente definici´on.

´ 2.3. ANALISIS DE ESQUEMAS DE DIFERENCIAS FINITAS

35

Definici´ on 2.3.6. Un esquema de diferencias finitas Ph,k vjn = 0 es estable si y solamente s´ı existen constantes K y β y alguna norma ||.|| tal que ||v n || ≤ Keβt ||v 0 ||,

t = nk,

K, β

independientes de

h, k

(a)

Observaci´ on 2.3.1. : Si el esquema fuera de paso m´ ultiple, la definici´on anterior podr´ıa modificarse en el sentido de que existe un entero J, y constantes K, β tal que

||v || ≤ Ke n

βt

J

||v i ||

i=0

Una de las normas comunmente usada es la norma euclideana o l2 norma ||.||2 . Con esta norma la desigualdad (a) puede escribirse en una manera m´as pr´actica como

h

∞

j=−∞

|vjn |2

≤ Ke h βt

∞

|vj0 |2

(b)

j=−∞

Para ciertos esquemas simples podemos establecer la estabilidad yendo directamente a la desigualdad de la forma (a) o´ (b) sin embargo en general es un poco complicado, por eso es necesario utilizar algunas t´ecnicas como la t´ecnica de Von Neumann, que veremos mas adelante. Ejemplo.- Probaremos una condicion suficiente de estabilidad, para los esquemas de paso hacia adelante en tiempo y paso hacia adelante en espacio(FTFS). El esquema general se debe escribir en la forma siguiente: n vjn+1 = α vjn + β vj+1

mostraremos que el esquema es estable si | α | + | β| ≤ 1, sabiendo que

CAP´ITULO 2. ECUACIONES DIFERENCIALES PARCIALES

36

2xy ≤ x2 + y 2 , entonces ∞

 n+1 2 v  = j j=−∞

≤ ≤

∞

 n  n 2 α vj + β vj+1 j=−∞ ∞ 

j=−∞ ∞





2 |α|2  vjn 

 n  n    2  n 2     + 2 |α| |β| vj vj+1 + |β| vj+1

 2  n 2    n   + |β|2 vj+1  |α|2  vjn  + 2 |α| |β| vjn  vj+1

j=−∞

≤



∞ 

 n 2  vj  |α| + 2 |α| |β| + |β| 2

2

j=−∞ ∞

 n 2  vj  ≤ (|α| + |β|) 2

j=−∞

De esta desigualdad tenemos que ∞ ∞



 n 2  0 2 2n  vj  ≤ (|α| + |β|)  vj  j=−∞

j=−∞

as´ı el esquema ser´a estable si | α | + | β| ≤ 1. En particular el esquema (FTFS) dado por n vjn+1 = (1 + a r) vjn − ar vj+1

ser´a estable si | 1 + a r |+| a r| ≤ 1. Lo cual significa que, de cada sumando ar < 0 y aλ > −1 ⇔ −1 ≤ λ a ≤ 0 umero de Courant, el n´ umero de El n´ umero Cr = r a se llama n´ t Courant indica que el esquema anterior es estable si a < 0 y r =  x = k 1 h < − a la cual genera una dependencia entre k y h a fin de mantener estabilidad. 2.3.3

El Teorema de Equivalencia de Lax - Richtmyer

Teorema 2.3.1. Un esquema de diferencias finitas consistente en una ecuaci´ on diferencial parcial, para la cual el problema de valor inicial es bien puesto, es convergente ⇔ es estable.

´ 2.3. ANALISIS DE ESQUEMAS DE DIFERENCIAS FINITAS

37

La prueba del teorema la omitimos por el momento(consultar a Strikwerda [5]), sin embargo dibi´endose observar la importancia de mismo teorema puesto que nos nos da una forma simple de verificar si un esquema es convergente, ya que si intentamos verificar la convergencia directamente por la definici´on no siempre es posible, sin embargo chequear consistencia y estabilidad pueden ser mucho mas f´acil ya que involucran generalmente c´alculos algebr´aicos solamente. Esto puede ser hecho en un lenguaje de computaci´on simb´olico como Maple, Mathematica. El Teorema de Lax - Richtmyer es un ejemplo del mejor tipo de teoremas matem a´ticos, pues relaciona un concepto importante que es d´ificil de establecer directamente con otros conceptos que son relativamente f´aciles de verificar y establece esta relaci´on muy adecuadamente. Notar que si unicamente tuvieramos la mitad del teorema, diciendo que consistencia y estabilidad implican convergencia, entonces podr´ia ocurrir que existan esquemas que siendo inestables sean convergentes. Pero si tuvieramos unicamente la otra mitad del teorema, diciendo que un esquema consistente y convergente es estable; entonces no podremos afirmar que un esquema consistente y estable es convergente. La importancia del uso del teorema de Lax-Richtmyer no radica tanto, desde el punto de vista de verificar consistencia y estabilidad, si no mas bi´en por la relaci´on precisa establecida entre estos conceptos y el concepto de convergencia. 2.3.4

La Condici´ on de Courant-Friedrichs-Lewy

Hemos definido el n´ umero de Courant Cr = ar, donde r = hk . El comportamiento de este n´ umero es muy importante, porque de ´el depende mucho, la estabilidad de esquemas de diferencias finitas; ya que una condici´on sobre este n´ umero da argumentos sobre la propagaci´on de informaci´on. La interpretaci´on de este n´ umero se puede expresar de la siguiente manera Cr =

a 1 r

=

velocidad de propagaci´on en la soluci´on anal´ıtica velocidad de propagaci´on en la soluci´on num´erica

En primer lugar veamos que una condici´on sobre Cr es una condici´on necesaria para la estabilidad de Esquemas Expl´ıcitos.

CAP´ITULO 2. ECUACIONES DIFERENCIALES PARCIALES

38

Esta condici´on, llamada condici´on de Courant-Friedrichs-Lewy (CFL) es |Cr | ≤ 1 Definici´ on 2.3.7. Un esquema de DF expl´ıcito es cualquier esquema que puede escribirse en la forma vjn+1 = Qvjn

donde Q es operador de diferencias lineal

en general, cuando el esquema es de m´ ultple paso, se puede escribir n

n+1 Qσ vjn−σ , σ≤n vj = σ=0

Los esquemas impl´ıcitos ser´an estudiados mas tarde, veamos un primer resultado, el cual contiene los esquemas de un solo paso. n n + βvjn + γvj+1 , un esquema expl´ıcito Teorema 2.3.2. Sea vjn+1 = αvj−1 para la ecuaci´ on hiperb´ olica con r constante; una condici´on necesaria para estabilidad es la condici´ on (CFL)

|Cr | ≤ 1 Para sistemas de ecuaciones, donde v es un vector, y α, β, γ son matri∀ autovalor ai de la matriz A. ces, debemos tener |ai r| ≤ 1 Prueba Para el caso escalar, supongamos |Cr | = |ar| > 1; entonces en el punto (x, t) = (0, 1) vemos que la soluci´on de la EDP u(x, t) = u0 (x − at) depende de los valores de u0 (x) en x = +a y x = −a (en realidad solamente de uno de ellos a o −a). Pero el esquema EDF tendr´a v0n dependiendo u ´nicamente de vj0 para j ≤ n, por la forma del esquemacon h = r−1 k, tenemos jh ≤ r−1 kn = r−1 , pues kn = 1 Esto dice que v0n depende de x solamente para|x| ≤ r−1 < |a|. Por tanto v0n no converge a u(0, 1) cuando h→0

con

h k

=1

por tanto |Cr | ≤ 1 De manera similar se puede probar que no existe esquemas consistentes expl´ıcitos para EDP hiperb´olicos que sean estables para todos los valores de r (con r constante cuando h, k → 0), presentamos el siguiente teorema; probado primeramente por Courant,Friedrichs y Lewy.

´ 2.4. ANALISIS DE FOURIER

39

Teorema 2.3.3. No existe EDF consistentes, incondicionalmente estables expl´ıcitos para sistemas hiperb´ olicos de EDP. A continuaci´on veamos un esquema impl´ıcito para la ecuaci´on de la onda (1), el mismo es consistente y estable para todo valor de r ilustrando que el teorema anterior no se puede extender a esquemas impl´ıcitos. a) Esquema de paso atr´as en tiempo y paso atr´as en el espacio (BTBS) n+1 vjn+1 − vj−1 vjn+1 − vjn +a =0 k h

(+)

En efecto; la consistencia es f´acil de ver, veamos la estabilidad: Escribiendo (+) en la forma (con λ = r) n+1 (1 + aλ)vjn+1 = vjn + aλvj−1

elevando al cuadrado ambos miembros,utilizando la desgualdad 2xy ≤ x2 + y 2 tenemos n+1 n+1 2 | + (aλ)2 |vj−1 | (1 + aλ)2 |vjn+1 |2 ≤ |vjn |2 + 2aλ|vjn ||vj−1 n+1 2 ≤ (1 + aλ)|vjn |2 + (aλ + (aλ)2 )|vj−1 |

Tomando la suma sobre j tenemos 2

(1 + aλ)

∞

|vjn+1 |2

≤ (1 + aλ)

j=−∞

∞

|vjn |2

2

+ (aλ + (aλ) )

j=−∞

∞

n+1 2 |vj−1 |

j=−∞

como (1 + aλ)2 − aλ − a(λ)2 = 1 + aλ, obtenemos: ∞

|vjn+1 |2

j=−∞

≤

∞

|vjn |2 .

j=−∞

Lo cual muestra que el esquema es estable para cualquier λ > 0.

2.4

An´ alisis de Fourier

En esta secci´on presentaremos esta herramienta importante, con la cual analizamos esquemas de diferencias finitas y sus soluciones del problema. En realidad usamos el An´alisis de Fourier para estudiar tanto los esquemas de diferencias finitas como las Ecuaciones Diferenciales Parciales.

40

CAP´ITULO 2. ECUACIONES DIFERENCIALES PARCIALES

Definici´ on 2.4.1. Para una funci´ on u(x) definida en la recta real R, su transformada de Fourier uˆ(w) se define por  ∞ 1 e−iwx u(x)dx, si la integral existe uˆ(w) = √ 2π −∞ Esta transformada de u es una funci´on de la variable real w(longitud de onda) y es u ´nicamente definida por u. La funci´on uˆ es una representaci´on alternativa de la funci´on u. Informaci´on de ciertas propiedades de u pueden ser deducidas de las propiedades de uˆ, por ejemplo, la raz´on en la cual uˆ decae para valores grandes de w est´a relacionado al n´ umero de derivadas que u tiene. La f´ormula de inversi´on de Fourier est´a dada por  ∞ 1 eiwx uˆ(w)dw u(x) = √ 2π −∞ Esta f´ormula expresa que u es una superposici´on de ondas eiwx con diferentes amplitudes uˆ(w). Definici´ on 2.4.2. Sea v una funci´ on discreta definida en Z; entonces su Transformada de Fourier se define por ∞ 1 −imζ e vm , vˆ(ζ) = √ 2π m=−∞

ζ ∈ [−π, π],

vˆ(π) = vˆ(−π)

Su f´ ormula de inversi´ on de Fourier est´ a dada por  π 1 eimζ vˆ(ζ)dζ vm = √ 2π −π Esta f´ormula de inversi´on tiene la misma interpretaci´on que la f´ormula de inversi´on en R. Observamos que el An´alisis de Fourier en Z es exactamente el estudio de las representaciones de las series de Fourier de funciones definidas en el intervalo [−π, π]. Si el espacio entre los puntos discretos es h, podemos establecer el an´alisis de Fourier en hZ por un cambio de variable   ∞ −π π 1 −imhζ , e vm h, ζ∈ − vˆ(ζ) = √ h h 2π m=−∞

´ 2.4. ANALISIS DE FOURIER

41

Con su f´ormula de inversi´on 1 vm = √ 2π



π h

− πh

eimhζ vˆ(ζ)dζ

(2.23)

Las definiciones anteriores nos garantizan que en la L2 -norma se satisface u||L2 ||u||L2 = ||ˆ

y ||v||2h =

∞

|vm |2 h =

n=−∞



π h

− πh

|ˆ v (ζ)|2 dζ = ||ˆ v ||2h

Llamadas f´ormulas de Parseval, usando estas relaciones podemos decir que la transformada de fourier se define para todas las funciones en L2 (R) y en L2 (hZ). En el siguiente ejemplo vemos la utildad pr´actica de la transformada de Fourier en algunas situaciones Ejemplo Sea la funci´on discreta ⎧ ⎪ ⎨ 1, si |xm | < 1 1 sobre hZ vm = 2 , si |xm | = 1 ⎪ ⎩ 0, si |xm | > 1 Supongamos que h = m mh = M ⎡

1 M,

donde M es un entero, entonces, xm =

M −1

⎤

h 1 1 vˆ(ζ) = √ ⎣ eiM hζ + eimhζ + e−iM hζ ⎦ 2 2π 2 m=−(M −1)   sen(M − 12 )hζ h = √ cosζ + sen 12 hζ 2π   senM hζcos 12 hζ − cosM hζsen 12 hζ h cosζ + = √ sen 12 hζ 2π   h h cosζ + senζcot ζ − cosζ = √ 2 2π h 1 = √ senζcot hζ 2 2π

CAP´ITULO 2. ECUACIONES DIFERENCIALES PARCIALES

42

La f´ormula de Parseval asegura que el c´alculo de la siguiente integral es:  π M −1

h2 4 1 2 1 2 21 2 2 v ||h = ||v||h = h[( ) + sen ζcot hζdζ = ||ˆ 1 + ( )2 ] 2π − π4 2 2 2 m=−(M +1)

1 = h[ + M − 1 − (−M + 1) + 1] 2 1 = h[ + 2M − 1] 2 1 = h[2M − ] 2 1 = 2 − h. 2 An´ alisis de Von Neumann

Con el uso del an´alisis de Fourier se puede dar condiciones necesarias y suficientes para la estabilidad de esquemas de diferencias finitas, esto es lo que se llama an´alisis de Von Neumann. Recordemos que un esquema de diferencias finitas de un solo paso puede escribirse en la forma vjn+1 = Qvjn , Siendo Q = Q(S+ , S− ), donde S+ , S− son los operadores de desplazan , S− vjn = miento hacia adelante y hacia atr´as respectivamente S+ vjn = vj+1 n vj−1 Observemos que, tomando la Transformada de Fourier a S+ vj tenemos S +v =

∞

j=−∞

−ijhζ

vj+1 e

=

∞

−i(m−1)hζ

vm e

ihζ

=e

m=−∞

∞

vm e−imhζ = eihζ vˆ

m=−∞

de la misma forma se obtiene −ihζ vˆ S −v = e

De esto se sigue que si v n+1 = Q(S+ , S− )v n entonces vn vˆn+1 = Q(eihζ , e−ihζ )ˆ Si denotamos por θ = hζ, entonces el t´ermino g(θ, k, h) = Q(eihζ , e−ihζ ) se llama factor de amplificaci´on.

´ 2.4. ANALISIS DE FOURIER

43

Condici´ on de Estabilidad de Von Neumann El factor de amplificaci´on g(θ, k, h) se dice que satisface la condici´on de estabilidad de Von Neumann(CVN), si ∃ una constante C > 0 (independiente de θ, k, h) tal que |g(θ, k, h)| ≤ 1 + Ck,

k es el peso del tiempo

esta condici´on puede ser vista en el siguiente teorema. Teorema 2.4.1. Un esquema de diferencias finitas de un paso es estable en la l2 norma ⇐⇒ satisface la condici´on de Von Neumann. Si g(θ, k, h) = g(θ) entonces la CVN se puede sustituir por |g(θ)| ≤ 1 Este teorema muestra que para determinar la estabilidad de un EDF, solamente necesitamos considerar el factor de amplificaci´on g(θ, k, h) Demostraci´ on. ⇐ Si la CVN es satisfecha entonces el el EDF es estable Sea T > 0 suficientemente grande, y supongamos que existe C > 0 tal que ∀ θ, 0 < k ≤ k0 , 0 < h ≤ h0 .

g(θ, k, h) ≤ 1 + Ck,

por la definici´on del factor de amplificaci´on se tiene Vˆ n = g Vˆ n−1 = g n Vˆ 0 Luego usando la identidad de Parseval  −π  h V n 2h = Vˆ n 2h = |Vˆ n (ζ)|2 dζ = − πh

≤ (1 + Ck)2n = (1 +



− πh

− πh

− πh

|g|2n |Vˆ 0 (ζ)|2 dζ

|Vˆ 0 (ζ)|2 dζ = (1 + Ck)2n Vˆ 0 2h

− πh Ck)2n V 0 2h

Como nk ≤ T , entonces k ≤ T /n, y del hecho que 1 + x ≤ ex , tenemos T (1 + Ck)2n ≤ (1 + Ck)2 k ≤ e2CT

x > 0,

Entonces usando los extremos de la desigualdad anterior obtenemos V n 2h ≤ (1 + Ck)2n V 0 2h ≤ e2CT V 0 2h

44

CAP´ITULO 2. ECUACIONES DIFERENCIALES PARCIALES

podemos decir que existe CT = e2CT > 0, asi como para las constantes h0 > 0, k0 > 0 elegidas en la CVN se cumple V n 2h ≤ CT V 0 2h ,

∀ h ≤ h0 , k ≤ k0 , nk ≤ T.

⇒ Para esto utilizamos la contraposici´on, es decir si la CVN no es satisfecha entonces el EDF no es estable En efecto: Supongamos que para cada C > 0, existe θc ∈ [−π, π] tal que g(θc , k, h) > 1 + Ck,

, 0 < k ≤ k0 , 0 < h ≤ h0 .

Entonces |g| tiene un comportamiento similar al mostrado en la figura (2.2) 6

|g| s

1 + Ck

s

θ1 θc

θ2

π

-

Figura 2.2:

Por la continuidad de g(θ) existe un intervalo [θ1 , θ2 ] que contiene a θc y adem´as |g(θ)| ≥ 1 + Ck, ∀ θ ∈ [θ1 , θ2 ] Elijamos una funci´on Vˆ 0 en el espacio de longitud de o nda(frecuencias) que sea nula en el exterior del intervalo [θ1 , θ2 ], por ejemplo   0, θ∈ / [θ1 , θ2 ] 0, ζ∈ / [ θh1 , θh2 ] 0 Vˆ (ζ) = = h h / [ θh1 , θh2 ] θ2 −θ1 , θ ∈ [θ1 , θ2 ] θ2 −θ1 , ζ ∈ Observamos que Vˆ 0  = 1.

´ 2.4. ANALISIS DE FOURIER

45

Entonces utilizando nuevamente la identidad de parseval tenemos  −π h |Vˆ n (ζ)|2 dζ V n 2h = Vˆ n 2h =  =

− πh

 =

− πh

−

θ2 h

θ − h1

− πh

|g|2n |Vˆ 0 (ζ)|2 dζ |g|2n |Vˆ 0 (ζ)|2 dζ 2n

≥ (1 + Ck)



−

−

θ2 h

θ1 h

ˆ|V 0 (ζ)|2 dζ = (1 + Ck)2n

Si elegimos C para el cual existe T > 0 y k con kn ∼ T tal que se cumpla T 2 + 2CT ≥ e2CT ∼ (1 + kC)2 k Podemos concluir que T

V n 2h ≥ (1 + Ck)2n = (1 + Ck)2 k T 1 1 ≥ (1 + Ck)2 k ∼ e2CT .1 2 2 1 2CT 0 2 = e V h 2 Por tanto para cada CT = 12 e2CT existe T > 0 tal que 1 + CT ≥ CT y se cumple V n 2h ≥ CT V 0 2h que es la no estabilidad del esquema. Ejemplo Consideremos el esquema de paso hacia adelante en tiempo y paso hacia adelante en espacio(forward time-forward space FTFS) n vjn+1 − vjn − vjn vj+1 +a = 0, k h el cual se escribe en la forma

vjn+1 = = = =

a > 0,

n vjn + aλvjn − aλvj+1 (1 + aλ)vjn − aλS+ vjn [(1 + λ)I − aλS+ ]vjn Q(S+ , S− )vjn

46

CAP´ITULO 2. ECUACIONES DIFERENCIALES PARCIALES

tomando transformada de Fourier vˆn+1 = g(θ, h, k)ˆ vn g(θ, h, k) = 1 + aλ − aλeihζ , con a > 0, λ =cte. |g|2 = (1 + aλ − aλcosθ)2 + (aλ)2 sen2 θ = 1 − 2(aλ)2 cosθ − 2aλcosθ + 2aλ + 2(aλ)2 = (1 + aλ)2 − 2aλ(1 + aλ)cosθ + (aλ)2 = (1 + aλ)[1 + aλ − 2aλcosθ] + (aλ)2 θ = (1 + aλ)[1 + aλ2sen2 − aλcosθ] + (aλ)2 2 θ θ θ = (1 + aλ)[1 + 2aλsen2 − aλ(cos2 − sen2 ) + (aλ)2 2 2 2 θ θ = (1 + aλ)[1 + 2aλsen2 − aλ(1 − 2sen2 )] + (aλ)2 2 2 θ = (1 + aλ)[1 + 4aλsen2 − aλ] + (aλ)2 2 θ = 1 − a2 λ2 + 4aλ(1 + aλ)sen2 + (aλ)2 2 θ = 1 + 4aλ(1 + aλ)sen2 2 Observemos que si θ = 0 ⇒ |g| > 1, por tanto el esquema es inestable, recordamos que ya hemos visto que este esquema no es convergente. Observaci´ on 2.4.1. No es necesario escribir las integrales (2.23) para tener la ecuaci´ on para el factor amplificador g. Un procedimiento equivalente y simple es remplazar vjn por g n eijθ en el esquema, para cada valor de j y n. Ejemplo Ahora utilizando el esquema paso adelante en el tiempo y diferencia central en el espacio (FTCS) n n vjn+1 − vjn − vj−1 vj+1 +a =0 k 2h ilustramos este procedimiento. Remplazando vjn por g n eijθ el esquema anterior se transforma en

g n ei(j+1)θ − g n ei(j−1)θ g n+1 eijθ − g n eijθ +a k 2h   iθ e − e−iθ n ijθ g − 1 =0 =g e +a k 2h

´ 2.4. ANALISIS DE FOURIER

47

luego de simplificar esta ecuaci´on obtenemos el factor de amplificaci´on g = 1 − iaλ sin θ con λ = k/h. Esta t´ecnica de obtener el factor de amplificaci´on es obviamente mas facil que el an´alisis presentado anteriormente. Si λ es constante, entonces g es independiente de h y k y |g(θ)|2 = 1 + a2 λ2 sin2 θ. Desde que |g(θ)| > 1 para θ diferente de 0 y de π, este esquema es inestable. Se debe tener en cuenta que, el hecho de sustituir vjn por g n eijθ no quiere decir que la soluci´on del esquema de diferencas tenga necesariamente la forma vjn = g n eijθ , En verdad esta sustituci´on es una forma corta de interpretar el m´etodo de la transformada de Fourier utilizado al inicio de la secci´on, donde se verifica que todas las soluciones de los esquemas de un solo paso son dados por la f´ormula vjn+1 = Qvjn la cual est´a lista para obtenr en forma pr´actica el factor de amplificaci´on. Un reordenamiento de los c´alculos para determinar el factor de amplificaci´on muestra que los dos procedimientos son equivalentes en la obtenci´on del factor de amplificaci´on.

48

CAP´ITULO 2. ECUACIONES DIFERENCIALES PARCIALES

Cap´ıtulo 3 Ecuaciones Diferenciales Parab´ olicas 3.1

Problema de Conducci´ on del Calor

En el presente cap´ıtulo se presenta esquemas de diferencias finitas para la ecuaci´on del calor, se pretende presentar de una manera detallada los aspectos num´ericos y computacionales en una clase de esquemas num´ericos para las ecucaciones difrenciales parciales de tipo parab´olico, como es la ecuaci´on del calor. Par ser mas especificos, tomamos el problema modelo de la siguiente forma ⎧ = α2 uxx , 0 < x < L ,t > 0 ut ⎪ ⎪ ⎨ t>0 u(0, t) = g1 (t), (P ) u(L, t) = g2 (t), t>0 ⎪ ⎪ ⎩ 0≤x≤L u(x, 0) = u0 (x), el cual es un problema t´ıpico de valor inicial y de contorno, α2 es el par´ametro de difusividad t´ermica y representa una proporci´on entre la conductividad t´ermica y la capacidad de almacenamiento t´ermico del material donde se realiza el estudio de difusi´on del calor. Como observamos, el problema es unidimensional, el cual puede ser una abstracci´on matem´atica de la f´ısica de la conducci´on del calor en una barra delgada homogenea de secci´on transversal constante y aislada por los costados, cuyos extremos se mantienen a temperaturas conocidas g1 (t) y g2 (t), e inicialmente se supone que tiene una distribuci´on de temperatura u0 (x) que depende unicamente de la posici´on de la secci´on 49

50

´ CAP´ITULO 3. ECUACIONES DIFERENCIALES PARABOLICAS

transversal x. El hecho de que este problema modelo, tiene soluci´on anal´ıtica expl´ıcita, es un buen modelo para probar la eficiencia de los esquemas num´ericos que se pretenden estudiar. Entre los esquemas de diferencias finitas que sirven para solucionar este problema tenemos:

• Esquemas explicitos; entre de los que se encuentran: El esquema de pso adelante en tiempo y diferencia central en el espacio tambi´en llamado esquema de Euler, el esquema de Du For-Frankel, esquema de Leaprog,etc.

• Esquemas impl´ıcitos, como el esquema de Euler Retrasado, el esquema de Crank-Nicolson. Para el desarrrollo del presente trabajo solamente trataremos con los esquemas Impl´ıcitos; Euler Retrazado y Crank-Nicolson, y el esquema explicito de Du Fort-Frankel, dejando los dem´as como un ejercicio al lector.

Discretizaci´ on del Dominio Observamos que nuestro dominio es el conjunto Ω × R+ =]0, L[×(0, ∞) como se muestra en la figura 3.1. Denotemos por h = Δx la longitud de la partici´on del dominio Ω ubicado en el eje de las abscisas, y por k = Δt la longitud de la partici´on en el eje de las ordenadas. Con estos n´ umeros dividimos el dominio en subdominios con lados de longitud h y k,como se muestra en la figura 3.2

´ DEL CALOR 3.1. PROBLEMA DE CONDUCCION

t

51

6

]0, L[×R+

valores de frontera

valores de frontera

-

0

valores iniciales

L

x

Figura 3.1: Dominio del problema

t

6

valores de frontera tn

e

(j, n)

valores de frontera

-

0

xj L valores iniciales Figura 3.2: Discretizaci´on del Dominio

x

De esta forma definamos los puntos de la malla como (xj , tn ), para j = 0, . . . , M + 1,

n = 0, . . . Nmax ;

los puntos en el eje de las abscisas toman la forma xj = jh, mientras que los puntos en el eje de las ordenadas toman la forma tn = nk

52

´ CAP´ITULO 3. ECUACIONES DIFERENCIALES PARABOLICAS

Discretizaci´ on de la variable Definici´ on 3.1.1. Cualquier funci´ on V definida sobre una red o malla se llama funci´on discreta, y sus valores se denotan por Vjn o (Vj,n ) Observaci´ on 3.1.1. . 1. Si u es una funci´on continua en (x, t) ( definida en todo el doan denotados por minio), Entonces los valores de u en (xj , tn ) ser´ n uj = u(xj , tn ) 2. El conjunto de puntos (xj , tn ) para j fijo son llamados nivel n de la malla 3. Muchas veces, para simplificar notaciones, se usa los ´ındices (j, n) para representar los puntos de la malla (xj , tn ). Discretizaci´ on de la Ecuaci´ on Definici´ on 3.1.2. Un esquema de diferencias finitas es una combinaci´ on de operadores discretos que permiten aproximar la ecuaci´ on diferencial parcial por una ecuaci´ on en diferencias finitas. Los valores de u definidos sobre las fronteras laterales y la frontera inferior son valores conocidos mientras que en los puntos interiores deben ser calculados. La aproximaci´on por diferencias finitas a las derivadas de una funci´on se obtiene utilizando el Teorema de Taylor, segun el an´alisis hecho en la secci´on anterior. La variable espacial x ser´a aproximada por la diferencia central, es decir utilizamos la siguiente aproximaci´on: uxx (x, t) 

u(x + h, t) − 2u(x, t) + u(x − h, t) h2

(3.1)

con error de truncamiento τ (x0 , t0 , h) =

h2 ∂ 4 u(η0 , t0 ). 12 ∂x4

Generalmente se utilizan estas diferencias porque los valores de u son conocidos en los extremos del intervalo [0, L] u(0, t) = g1 (t),

u(L, t) = g2 (t),

t>0

´ DEL CALOR 3.1. PROBLEMA DE CONDUCCION

53

En la variable temporal la situaci´on es diferente. Solamente son conocidos los valores de la funci´on en el tiempo t = 0 esto es: u(x, 0) = u0 (x),

x ∈ [0, L]

Por tanto no es posible utilizar diferencias centrales. Para este caso usualmente se utiliza la aproximaci´on de diferencia progresiva. ut (x, t) 

u(x, t + h) − u(x, t) k

(3.2)

con el error de truncamiento k ∂2 (t0 , ν0 ) τ (x0 , t0 , k) = − 2 ∂t2 Usando la notaci´on Vjn para la aproximaci´on en (xj , tn ) obtenemos Vjn+1 − Vjn ∂u (xj , tn ) ≈ ∂t k

(3.3)

y

n n − 2Vjn + Vj−1 Vj+1 ∂ 2u (xj , tn ) ≈ (3.4) ∂x2 h2 De las f´ormulas (3.3) y (3.4) obtenemos el esquema de diferencias progresiva para la ecuaci´on diferencial parcial del problema n Vjn+1 − Vjn V n − 2Vjn + Vj−1 2 j+1 =α k h2

(FTCS)

llamado comunmente esquema (FTCS), del Ingles forward time-central space, cuyo error de truncamiento es: τj,j

2 4 k ∂2 2h ∂ = u(xj , ξn ) − α U (ηj , tn ) ∈ O(k) + O(h2 ) 2 4 2 ∂t 12 ∂x

En el esquema (FTCS), la variable discreta Vjn , involucra los puntos de malla (xj−1 , tn ), (xj , tn+1 ), (xj+1 , tn ), (xj , tn ), que aproximan a la soluci´on u(x, t) del problema (P ), en el punto (xj , tn ). Este esquema se puede escribir de la siguiente manera  n n (3.5) + Vj−1 Vjn+1 = (1 − 2α2 μ)Vjn + α2 μ Vj+1 donde μ = kh−2

54

´ CAP´ITULO 3. ECUACIONES DIFERENCIALES PARABOLICAS

t

6

w

tn+1

(j, n + 1) w

tn

(j − 1, n)

w

(j, n)

w

(j + 1, n)

-

0

xj−1

xj

xj+1

L

x

Figura 3.3: Mol´ecula de c´alculo para FTCS

Observamos que los valores en el nivel n + 1 se obtiene unicamente sabiendo los valores en el nivel n, como se muestra en la figura 3.3, por eso, este esquema pertenece al conjunto de esquemas expl´ıcitos. Tambi´en podemos obtener otro esquema utilizando diferencia finita regresiva en el tiempo, es decir haciendo la siguiente aproximaci´on u(xj , tn+1 ) − u(xj , tn ) ∂u (xj , tn+1 ) ≈ ∂t k con error de truncamiento τj,j =

k ∂2 h2 ∂ 4 u(x , t ) − u(ηj , tn ) ∈ (O(k) + O(h2 )) j n 2 4 2 ∂t 12 ∂x

usando la notaci´on Vjn en la aproximaci´on num´erica anterior tenemos Vjn − Vjn−1 ∂u (xj , tn ) ≈ ∂t k

(3.6)

manteniendo la diferencia central en el espacio tenemos el esquema num´erico (BTCS) (del ingles ”backward time-central space”) para el problema de conducci´on del calor, tambi´en se le conoce como Euler retrasado . n Vjn − Vjn−1 V n − 2Vjn + Vj−1 2 j+1 =α k h2

(BTCS)

´ DEL CALOR 3.1. PROBLEMA DE CONDUCCION

55

En el esquema (BTCS), la variable discreta Vjn , involucra los puntos de malla (xj−1 , tn ), (xj , tn ), (xj+1 , tn ), (xj , tn−1 ), que aproximan a la soluci´on u(x, t) en el punto (xj , tn ) Este esquema se puede escribir de la siguiente manera  n n (1 − 2α2 μ)Vjn − α2 μ Vj+1 = Vjn−1 + Vj−1 (3.7) Observamos que los valores en el nivel n + 1 deben ser calculados simultaneamente sabiendo un solo valor en el nivel n, como se muestra en la figura 3.4, por eso, este esquema pertenece al conjunto de esquemas impl´ıcitos. t

6

w

tn

(j − 1, n)

w

(j, n)

w

(j + 1, n)

w

tn−1

(j, n − 1)

-

0

xj−1

xj

xj+1

L

x

Figura 3.4: Mol´ecula de c´alculo para Euler retrazado o´ BTCS

A continuaci´on presentamos los otros esquemas que ser´an estudiados en el presente trabajo. El esquema de Crank-Nicolson se obtiene promediando los esquemas Euler retrasado y Euler progresivo, de la siguiente manera: Sea el esquema de Euler progresivo en el punto (xj , tn ), utilizando el punto (xj , tn+ 21 ) obtenemos: n+1/2

Vj

k 2

− Vjn

=α

Vn 2 j+1

n − 2Vjn + Vj−1 h2

(FTCS)

Sea el esquema de Euler retrasado en el punto (xj , tn+1 ) utilizando el

´ CAP´ITULO 3. ECUACIONES DIFERENCIALES PARABOLICAS

56

punto (xj , tn+ 21 ) n+1/2

Vjn+1 − Vj k 2

=α

V n+1 2 j+1

n+1 − 2Vjn+1 + Vj−1 . h2

(BTCS)

Promediando los dos esquemas anteriores, (FTCS) y (BTCS), tenemos el esquema de Crank-Nicolson n+1 n+1 n+1 n n Vjn+1 − Vjn − 2Vjn + Vj−1 1 2 Vj+1 − 2Vj + Vj−1 1 2 Vj+1 = α + α (3.8) k 2 h2 2 h2 El nuevo esquema al ser promediado, se convierte en un esquema de diferencia central en la variable t, si lo observamos como discretizado en el punto (j, n + 1/2). Veremos que este esquema es tambi´en impl´ıcito, consistente, incondicionalmente estable y orden de presici´on (2,2). Observamos que este esquema utiliza seis puntos de malla (xj−1 , tn ), (xj , tn ), (xj+1 , tn ), (xj−1 , tn+1 ), (xj , tn+1 ), (xj+1 , tn+1 ) como se observa en la figura 3.5 

n+1

n

6 k 2 ?

w

h

-

w

w

f n + 1) (j, 2 w

j−1

w

j

6

k

n+

1 2

w ?

j+1

Figura 3.5: Mol´ecula de c´alculo para Crank-Nicolson

Usando diferencias centrales en el tiempo para la funci´on u(x, t) en el punto (xj , tn ) obtenemos la formula de diferencias finitas u(xj , tn+1 ) − u(xj , tn−1 ) ∂u (xj , tn ) = + τ (t1 , k) ∂t 2k donde τ (t1 , k) ∈ O(k 2 ), despreciando este error de truncamiento obtenemos la aproximacion de diferencias finitas centradas respecto a (xj , tn ) u(xj , tn+1 ) − u(xj , tn−1 ) ∂u (xj , tn ) ≈ ∂t 2k

´ DEL CALOR 3.1. PROBLEMA DE CONDUCCION

57

Usando la notaci´on Vjn tenemos el esquema de diferencias finitas de Leapfrog. n Vjn+1 − Vjn−1 V n − 2Vjn + Vj−1 2 j+1 , =α 2k h2

j = 1, ...N − 1 , n = 1, .... (3.9)

con error de truncamiento τj,j ∈ (O(k 2 ) + O(h2 ))

 En el esquema (3.9), podemos sustituir −2Vjn por − Vjn+1 − Vjn−1 , obtenemos;  n+1 n−1 n n Vjn+1 − Vjn−1 + Vj−1 V − V − V j j 2 j+1 =α (15) 2k h2 llamado esquema de Du Fort-Frankel, cuyo error de truncamiento es dado por τ (h, k) = O(h2 ) + O(k 2 ) + O(k −2 h2 ) como verificaremos en el cap´ıtulo siguiente. t

6

w

tn+1

(j, n + 1) w

tn

w w

(j, n − 1)

0

-

xj−1

xj

xj+1

L

x

Figura 3.6: Mol´ecula de c´alculo para Du Fort-Frankel

Este esquema es en realidad un esquema de cinco puntos aunque en la f´ormula y en la figura 3.6 aparecen explicitamente solo cuatro puntos, pero la discretizaci´on se hace en el punto denotado por (j, n).

58

3.2 3.2.1

´ CAP´ITULO 3. ECUACIONES DIFERENCIALES PARABOLICAS

An´ alisis de los Esquemas de Diferencias Finitas Teorema de Lax

El teorema de Lax que enunciamos aqu´ı es una variaci´on del teorema de Lax-Ritchmeyer presentado en el cap´ıtulo anterior, muy util, en el an´alisis de esquemas de diferencias finitas, puesto que, por medio de el y para una amplia clase de ecuaciones diferenciales, nos da la respuesta de la convergencia del esquema, utilizando unicamente los conceptos t´ecnicos de consistencia y estabilidad. Teorema 3.2.1. Sea un esquema de diferencias finitas lineal, estable y consistente con orden de precisi´ on (p, q). Entonces ´el es convergente de orden (p, q) Demostraci´on. Un esquema lineal, de un solo paso, con coeficientes constantes se puede escribir en la siguiente forma: Vjn+1 = QVjn ,

n≥0

(3.10)

donde Q = Qh,k . De acuerdo a esto se tiene V n = QV n−1 = Q(QV n−2 ) = . . . Qn−1 (QV 0 ) = Qn V 0 donde Qn es el operador de diferencias finitas Q aplicado n veces. Desde que el esquema es estable. Para todo T > 0, Existe CT > 0 tal que V n 2h ≤ CT V 0 2h Usando (3.10), tenemos que QV n−1 2h ≤ CT V 0 2h ←→ Qn V 0 2h ≤ CT V 0 2h Sea u(x, t) una soluci´on de la ecuaci´on P u = 0. Como el esquema es consistente con orden de precisi´on (p, q) se cumple, con Ph,k u = (Iun+1 − Qun ) Ph,k u = Ph,k u − P u = O(hp ) + O(k q ) es decir un = Qn−1 u + O(hp ) + O(k q ).

´ 3.2. ANALISIS DE LOS ESQUEMAS DE DIFERENCIAS FINITAS

59

Sea W n = un − V n el error en el n-´esimo paso, con W 0 = u0 − V 0 = 0, entonces, se cumple W n = QW n−1 + O(hp ) + O(k q ) = Q2 W n−2 + Q(O(hp ) + O(k q )) + O(hp ) + O(k q ) = ... n−1

n 0 Qj (O(hp ) + O(k q )) = Q W + j=0

=

n−1

Qj (O(hp ) + O(k q ))

j=0

De esto se obtiene, W  ≤ n

n−1

Qj (O(hp ) + O(k q ))

j=0

≤

n−1

CTj (O(hp ) + O(k q )),

Tj ∼ tj = jk

j=0

= O(hp ) + O(k q ) Por tanto el m´etodo es convergente de orden (p, q). A continuaci´on pasamos a analizar cada uno de los esquemas propuestos, cabe resaltar que para cada uno de los esquemas, estudiaremos la consistencia, orden de precisi´on, estabilidad y convergencia 3.2.2

El esquema de Euler retrasado

n Vjn − Vjn−1 V n − 2Vjn + Vj−1 2 j+1 =α (3.11) k h2 Como vemos, este esquema esta dentro del grupo de los esquemas impl´ıcitos. Probaremos que es consistente con orden de presici´on (1,2) e incondicionalmente estable. Denotemos por μ = kh−2 y λ = α2 μ. El esquema (3.11) se puede escribir en la forma n n +Vj−1 ), Vjn−1 = (1+2λ)Vjn −λ(Vj+1

j = 0, 1, 2, ..., N,

n = 0, 1, ..., kmax (3.12)

´ CAP´ITULO 3. ECUACIONES DIFERENCIALES PARABOLICAS

60

Si discretizamos las condiciones de contorno (Dirichlet) y condiciones iniciales del problema (P), tenemos: • Condiciones de Contorno 

V0n = u(x0 , tn ) = u(0, tn ) = T1 (tn ) VNn = U (xN , tn ) = u(L, tn ) = T2 (tn )

n = 0, 1, .., kmax

• Condiciones Iniciales Vi0 = u(xj , t0 ) = u(xj , 0) = u0 (xj ),

j = 0, 1, ..., N

Obteniendo el problema de diferencias finitas: ⎧ n n−1 n n n Vj −Vj 2 Vj+1 −2Vj +Vj−1 ⎪ ⎪ = α , j = 1, . . . , N − 1 n = 1, . . . , kmax ⎪ k h2 ⎨ n = T1 (tn ) V0 (P D1) n ⎪ = T2 (tn ) , n = 0, 1, .., k max V ⎪ ⎪ ⎩ VN0 = u0 (xj ), j = 1, 2, ..., N − 1 j 3.2.3

Consistencia y orden de precisi´ on

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: 2 ∂ 2 ∂ , P = −α ∂t ∂x2

el mismo que aplicado a una funci´on suave φ se tiene P φ = φt − α2 φxx . En tanto que el operador discreto Pk,h , aplicado a una variable discreta es el siguiente Ph,k Vjn

n Vjn − Vjn−1 V n − 2Vjn + Vj−1 2 j+1 = −α k h2

(3.13)

Utilicemos una funci´on suave φ y desarrollemos por Taylor en el punto (xj , tn ), las expresiones φn−1 = φ(xj , tn − k), φnj+1 = φ(xj + h, tn ) y j φnj−1 = φ(xj − h, tn )

´ 3.2. ANALISIS DE LOS ESQUEMAS DE DIFERENCIAS FINITAS

61

φ(xj , tn − k) = φ(xj , tn ) − kφt (xj , tn ) + O(k 2 ), en la notaci´on de ´ındices, tenemos φn−1 = φn−1 − kφt + O(k 2 ), j j

(3.14)

Luego la serie de Taylor de φ(xj − h, tn ) en (xj , tn ), es φ(xj − h, tn ) = φ(xj , tn ) −

1 1 hφx (xj , tn ) + h2 φxx (xj , tn ) − 1! 2!

1 3 h φxxx (xj , tn ) + O(h4 ) 3! En notaci´on indicial 1 2 1 h φxx − h3 φxxx + O(h4 ) 2! 3! De manera similar con φ(xj + h, tn ) en (xj , tn ), obtenemos: φnj−1 = φnj − hφx +

(3.15)

1 2 1 (3.16) h φxx + h3 φxxx + O(h4 ) 2! 3! las derivadas φt , φx , φxx , φxxx son evaluadas en el punto (xj , tn ) = (jh, nk). Sustituyendo las series (3.14),(3.15),(3.16) en el operador discretizado Ph,k φnj dado en (3.13), obtenemos: φnj+1 = φnj + hφx +

" 1! n φj − φnj + kφt + O(k 2 ) k   α2 n 1 2 1 3 − 2 φj + hφx + h φxx + h φxxx + O(h4 ) − 2φnj h 2! 3!   2 1 2 1 3 α n 4 − 2 φj − hφx + h φxx − h φxxx + O(h ) h 2! 3!

Ph,k φ =

la cual se puede escrbir como: α2  2 1 2 Ph,k φ = kφt + O(k ) − 2 h φxx + O(h4 ) + O(h4 ) k h Simplificando, y teniendo en cuenta la notaci´on del operador diferencial P y el hecho que O(h2 ) + O(h2 ) = O(h2 ), tenemos Ph,k φ = φt − α2 φxx + O(k) + O(h2 ) = P φ + O(k) + O(h2 )

62

´ CAP´ITULO 3. ECUACIONES DIFERENCIALES PARABOLICAS

De donde tenemos la forma P φ − Ph,k φ = O(k) + O(h2 )

(3.17)

Observando que P φ − Ph,k φ −→ 0,

cuando h, k −→ 0

Decimos que el esquema de diferencias finitas Euler retrasado es consistente con la ecuaci´on diferencial parcial ut = α2 uxx del problema (P ). De la ecuaci´on (3.17) vemos que el orden de precisi´on es (1, 2), 1 en el tiempo y 2 en el espacio. 3.2.4

Estabilidad

Para analizar la estabilidad del esquema de Euler retrasado, utilizamos el criterio de Von Neumann. Tomando el esquema en la forma expresada en (3.12) n+1 n+1 + (1 + 2λ)Vjn+1 − λVj−1 = Vjn −λVj+1

sustituyendo Vjn por g n eijθ , tenemos −λg n+1 ei(j−1)θ + (1 + 2λ)g n+1 eiθ − λg n+1 ei(j+1)θ = g n eijθ dividiendo por g n tenemos: −λgeijθ e−iθ + (1 + 2λ)geijθ − λgeijθ eiθ = eijθ multiplicando por e−ijθ −λge−iθ + (1 + 2λ)g − λgeiθ = 1 De esta forma tenemos el factor de amplificaci´on de la forma 1 g = iθ 1 − λ [e − 2 + e−iθ ] 1 = # iθ $2 iθ − 1−λ e2 −e 2 1 1 − λ [2i]2 sen2 2θ 1 = . 1 + 4λsen2 2θ =

´ 3.2. ANALISIS DE LOS ESQUEMAS DE DIFERENCIAS FINITAS

63

Observemos que |g| = g ≥ 0, por otro lado tenemos θ 0 ≤ sen2 ≤ 1 entonces 2

θ 1 ≤ 1 + 4λsen2 ≤ 1 + 4λ 2

De esto se tiene la desigualdad 0
0 entonces se tiene 1 − 2λ ≤ 0

(3.48)

luego el factor de amplificaci´on es complejo g± =

2λ cos θ ± ib 1 + 2λ

entonces el cuadrado del m´odulo es (2λ cos θ)2 + 4λ2 Sen2 θ − 1 |g| = (1 + 2λ)2 2

de donde 4λ2 − 1 (2λ + 1)(2λ − 1) 2λ − 1 = = |g| = 0 ut ⎪ ⎪ ⎨ t>0 u(0, t) = g1 (t), (P ) t>0 u(L, t) = g2 (t), ⎪ ⎪ ⎩ 0≤x≤L u(x, 0) = u0 (x), Donde la condici´on de compatibilidad entre la condici´on inicial y la condici´on de frontera esta dada por u0 (0) = u(0, 0) = g1 (0), u0 (L) = u(L, 0) = g2 (0),

Sea L = 1 y sea h = 1/N el espacio entre puntos de la malla. Como en algunos casos la estabilidad depende de la dependencia de h y k (el espacio entre los puntos temporales) la implementaci´on se realizar´a para algunos valores del par´ametro λ = α2 k/h2 Eligiendo uno de los esquemas de diferencias, el problema anterior puede escribirse en forma discretizada como ⎧ ⎪ ⎪ ⎪ Ph,k Vjn ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ V0n (P D) ⎪ ⎪ ⎪ ⎪ ⎪ VNn ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ V0 j

 = 0,

j = 1, 2, . . . , N − 1 n = 1, 2, . . . Kmax

= u(0, tn ) = g1 (tn ),

n = 0, 1, . . .

= u(L, tn ) = g2 (tn ),

n = 0, 2, . . .

= u(xj , 0) = u0 (xj ),

j = 0, 1 . . . , N

Comenzaremos implementando el esquema impl´ıcito de Euler retrazado, para ello, lo escribimos en la forma (3.12) n+1 n+1 −λVj−1 + (1 + 2λ)Vjn+1 − λVj+1 = Vjn

(3.49)

Escribiendo (3.49) para j = 1, . . . N − 1 aparece un sistema de N − 1

74

´ CAP´ITULO 3. ECUACIONES DIFERENCIALES PARABOLICAS

ecuaciones con N + 1 inc´ognitas −λV0n+1 + (1 + 2λ)V1n+1 − λV2n+1 −λV1n+1 + (1 + 2λ)V2n+1 − λV3n+1 .. . n+1 n+1 −λVi−1 + (1 + 2λ)Vin+1 − λVi+1 .. . n+1 n+1 −λVN −2 + (1 + 2λ)VNn+1 −1 − λVN

= V1n = V2n = Vin = VNn−1

Observamos que en este sistema de ecuaciones, los valores −λV0n+1 en la ecuaci´on j = 1 y −λVNn+1 en la ecuaci´on j = N − 1 son conocidas, por lo que deberian pasar al lado derecho de la ecuaci´on. Esto hace que el sistema resultante sea cerrado, es decir solo tiene N − 1 inc´ognitas, el cual puede ser escrito, para cada n, en forma matricial como AU n+1 = bn Donde A es una matriz tridiagonal ⎡ 1 + 2λ −λ 0 ⎢ −λ 1 + 2λ −λ ⎢ ⎢ ⎢ A=⎢ ⎢ ⎢ ⎣ 0 0 0 0

(3.50)

de orden (N − 1) × (N − 1) ⎤ ... 0 0 ... 0 0 ⎥ ⎥ .. ⎥ . ⎥ ⎥ .. ⎥ . ⎥ −λ 1 + 2λ −λ ⎦ −λ 1 + 2λ

mientras que U n+1 y bn son vectores de orden (N − 1) × 1 ⎡

U n+1

⎢ ⎢ ⎢ ⎢ ⎢ =⎢ ⎢ ⎢ ⎢ ⎣

V1n+1 V2n+1 .. . n+1 Vi .. . n+1 VN −2 VNn+1 −1

⎤

⎡

⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦

⎢ ⎢ ⎢ ⎢ ⎢ n b =⎢ ⎢ ⎢ ⎢ ⎣

y

V1n + λV0n+1 V2n .. . Vin .. . n VN −2 n VN −1 + λVNn+1

⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦

Para resolver el sistema (3.50), se puede utilizar eliminaci´on Gaussiana, asi como tambi´en factorizaci´on lu. Para ello consideremos el sistema

´ NUMERICA ´ 3.3. IMPLEMENTACION

75

tridiagonal general de la forma ⎡

d1 ⎢ a1 ⎢ ⎢ ⎢ A=⎢ ⎢ ⎢ ⎣ 0 0

c1 0 d2 c2

... 0 0 ... 0 0 .. . .. . aN −3 dN −2 cN −2 0 aN −2 dN −1

0 0

⎤

⎡

⎢ ⎥⎢ ⎥⎢ ⎥⎢ ⎥⎢ ⎥⎢ ⎥⎢ ⎥⎢ ⎦⎢ ⎣

V1n+1 V2n+1 .. . n+1 Vi .. . n+1 VN −2 VNn+1 −1

⎡

⎤

b1 b2 .. . bi .. .

⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥=⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎦ ⎣ bN −2 bN −1

⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦

de tal manera que sea diagonal dominante, es decir se cumple |di | > |ci | + |ai−1 |, El hecho de que una matriz sea diagonal dominante, es suficiente para usar eliminaci´on Gaussiana sin pivotamiento y factorizaci´on lu, el sistema (3.50), es un sistema diagonal dominante. La factorizaci´on lu de la matriz tridiagonal A se puede esquematizar en la siguiente forma ⎡ d 1 ⎢ a1 ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ 0 0

c1 d2

0 0

0 c2

... ... . . . . . .

0 0

aN −2 0

dN −1 aN −1

⎡ 1 ⎢ l1 ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎥ =⎢ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎣ ⎦ cN −1 0 0 dN 0 0

⎤

0 1

0 0

0 0

... ... . . . . . . lN −2 0

0 0

1 lN −1

0 ⎤ ⎡ u1 0 ⎥⎢ 0 ⎥⎢ ⎥⎢ ⎥⎢ ⎥⎢ ⎥⎢ ⎥⎢ ⎥⎢ ⎥⎢ ⎥⎢ ⎦⎣ 0 0 1 0

c1 u2

0 0

0 c2

... ... . . . . . . 0 0

0 0

0 0

uN −2 0

cN −1 uN

⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦

Supongamos que los elementos de la matriz tridiagonal se guardan en tres vectores columnas d, c y a, con los que obtendremos los vectores que conforman las matrices l y u, entonces el algoritmo para la factorizaci´on seria el siguiente: LET

u1 = d1

DO

i = 2 , N l(i-1) = a(i-1)/u(i-1) u(i) = d(i) - l(i-1)*c(i-1)

END DO

´ CAP´ITULO 3. ECUACIONES DIFERENCIALES PARABOLICAS

76

Luego se resuelve el sistema, calculando U a partir de l y de u, primeramente se resuelve el sistema triangular lY = B por sustituci´on progresiva

LET DO

y1 = b1 i = 2 , N yi = bi-l(i-1)*y(i-1)

END DO Finalmente se resuelve el sistema uU = Y por sustituci´on hacia atras

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 , repetido para cada nivel el algoritmo para la soluci´on del problema (PD). A continuaci´on implementamos el esquema impl´ıcito de Crank-Nicolson, para ello, escribimos el esquema en la forma (3.19) n+1 n+1 n n + (2 + 2λ)Vjn+1 − λVj+1 = λVj+1 + λVj−1 + (2 − 2λ)Vjn (3.51) −λVj−1

´ NUMERICA ´ 3.3. IMPLEMENTACION

77

Observamos que este esquema de diferencias, al igual que el anterior, genera un sistema de N − 1 ecuaciones con N − 1 inc´ognitas que puede ser escrito en forma matricial como: AU n+1 = bn Donde A es una matriz tridiagonal ⎡ 2 + 2λ −λ 0 ⎢ −λ 2 + 2λ −λ ⎢ ⎢ ⎢ A=⎢ ⎢ ⎢ ⎣ 0 0 0 0

(3.52)

de orden (N − 1) × (N − 1) ⎤ ... 0 0 ... 0 0 ⎥ ⎥ .. ⎥ . ⎥ ⎥ .. ⎥ . ⎥ −λ 2 + 2λ −λ ⎦ −λ 2 + 2λ

mientras que U n+1 y bn son vectores de orden (N − 1) × 1 ⎡ n+1 ⎤ V1 ⎢ V n+1 ⎥ ⎢ 2 ⎥ ⎢ .. ⎥ ⎢ . ⎥ ⎢ ⎥ n+1 y = ⎢ Vin+1 ⎥ U ⎢ . ⎥ ⎢ .. ⎥ ⎢ n+1 ⎥ ⎣ VN −2 ⎦ VNn+1 −1 ⎡ n λV−0 + (2 − 2λ)V1n + λV2n + λV0n+1 ⎢ λV1n + (2 − 2λ)V2n + λV3n ⎢ .. ⎢ . ⎢ ⎢ n n n b = ⎢ λVj−1 + (2 − 2λ)Vjn + λVj+1 ⎢ .. ⎢ . ⎢ ⎣ +λVNn−3 + (2 − 2λ)VNn−2 + λVNn−1 +λVNn−2 + (2 − 2λ)VNn−1 + λVNn + λVNn+1

⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦

Para resolver el sistema (3.52), se puede seguir los mismos pasos que se hizo al resolver el sistema tridiagonal obtenido por el m´etodo de Euler retrasado La implementaci´on del esquema de Du Fort-Frankel es simple porque es un esquema expl´ıcito de paso m´ ultiple, por consiguiente, este esquema se ejecutar´a a partir del segundo nivel de tiempo, necesitando para ello

78

´ CAP´ITULO 3. ECUACIONES DIFERENCIALES PARABOLICAS

los valores en los niveles 0 (que es la condici´on inicial) y 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 + Vj−1 ), j = 1, . . . N − 1, Vjn+1 = (1 − 2λ)Vjn + λ(Vj+1

n = 0 (3.53)

A continucaci´on , los valores en el nivel n + 1 se calcula por la f´ormula    λ 1 − 2λ n n Vjn+1 = 2 Vj+1 Vjn−1 , j = 1, . . . , N −1, + Vj−1 n≥1 + 1 + 2λ 1 + 2λ (3.54) Presentaremos resultados obtenidos por los tres esquemas en las mismas condiciones, luego se hace una comparaci´on con la soluci´on exacta, a fin de observar sus propiedades presentadas en el an´alisis te´orico, para ello, consideramos el valor α = 1, las condiciones de contorno g1 = 0

g2 = 0

y la codici´on inicial u0 (x) = sen(πx),

0