Introducción al Método de Diferencias Finitas y su Implementación Computacional Antonio Carrillo Ledesma y Omar Mendoza Bernal Facultad de Ciencias, UNAM http://www.mmc.geofisica.unam.mx/acl/ Una copia de este trabajo se puede descargar de http://www.mmc.geofisica.unam.mx/acl/Textos/ Invierno 2015, Versión 1.0α

Introducción al Método de Diferencias Finitas y su Implementación Computacional

Índice 1 Expansión en Series de Taylor 1.1 Aproximación de la Primera Derivada 1.1.1 Diferencias Progresivas . . . . . 1.1.2 Diferencias Regresivas . . . . . 1.1.3 Diferencias Centradas . . . . . 1.2 Derivadas de Ordenes Mayores . . . . 1.2.1 Derivada de Orden Dos . . . . 1.2.2 Derivadas de Ordenes Mayores 1.3 Derivadas en Dos Dimensiones . . . . 1.4 Derivadas en Tres Dimensiones . . . .

. . . . . . . . .

. . . . . . . . .

. . . . . . . . .

. . . . . . . . .

. . . . . . . . .

. . . . . . . . .

. . . . . . . . .

. . . . . . . . .

. . . . . . . . .

. . . . . . . . .

. . . . . . . . .

. . . . . . . . .

. . . . . . . . .

. . . . . . . . .

. . . . . . . . .

2 2 2 3 3 4 5 5 6 7

2 Método de Diferencias Finitas en Una Dimensión 10 2.1 Problema con Condiciones de Frontera Dirichlet . . . . . . . . . . 10 2.2 Problema con Condiciones de Frontera Neumann . . . . . . . . . 16 2.3 Problema con Condiciones de Frontera Robin . . . . . . . . . . . 17 2.4 Discretización del Tiempo . . . . . . . . . . . . . . . . . . . . . . 23 2.4.1 Ecuación con Primera Derivada Temporal . . . . . . . . . 23 2.4.2 Ecuación con Segunda Derivada Temporal . . . . . . . . . 27 2.5 Consistencia, Estabilidad, Convergencia y Error del Método de Diferencias Finitas . . . . . . . . . . . . . . . . . . . . . . . . . . 30 3 Consideraciones Sobre la Implementación de Métodos de Solución de Grandes Sistemas de Ecuaciones Lineales 3.1 Métodos Directos . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.1.1 Factorización LU . . . . . . . . . . . . . . . . . . . . . . . 3.1.2 Factorización Cholesky . . . . . . . . . . . . . . . . . . . . 3.1.3 Factorización LU para Matrices Tridiagonales . . . . . . . 3.2 Métodos Iterativos . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2.1 Método de Gradiente Conjugado . . . . . . . . . . . . . . 3.2.2 Método Residual Mínimo Generalizado . . . . . . . . . . . 3.3 Estructura Óptima de las Matrices en su Implementación Computacional . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.3.1 Matrices Bandadas . . . . . . . . . . . . . . . . . . . . . . 3.3.2 Matrices Dispersas . . . . . . . . . . . . . . . . . . . . . . 3.3.3 Multiplicación Matriz-Vector . . . . . . . . . . . . . . . .

33 33 33 35 35 36 38 40 41 43 44 45

4 Implementación Computacional del Método de Diferencias Finitas para la Resolución de Ecuaciones Diferenciales Parciales 47 4.1 Implementación en SciLab . . . . . . . . . . . . . . . . . . . . . . 47 4.2 Implementación en Octave (MatLab) . . . . . . . . . . . . . . . . 49 4.3 Implementación en C++ . . . . . . . . . . . . . . . . . . . . . . . 51 4.4 Implementación en Python . . . . . . . . . . . . . . . . . . . . . 53 5 Bibliografía

[email protected]

56

1

Antonio Carrillo Ledesma

Introducción al Método de Diferencias Finitas y su Implementación Computacional

1

Expansión en Series de Taylor

Sea f(x) una función definida en (a, b) que tiene hasta la k−ésima derivada, entonces la expansión de f(x) usando series de Taylor alrededor del punto xi contenido en el intervalo (a, b) será


(x − xi ) df

(x − xi )2 d2 f

(x − xi )k dk f

f (x) = f(xi ) + + + ... + (1.1) 1! dx
xi 2! dx2
xi k! dxk
ε donde ε = xi + θ(x − xi ) y 0 < θ < 1.

1.1

Aproximación de la Primera Derivada

Existen distintas formas de generar la aproximación a la primera derivada, nos interesa, una que nos de la mejor precisión posible con el menor esfuerzo computacional. 1.1.1

Diferencias Progresivas

Considerando la Ec.(1.1) con k = 2 y x = xi + ∆x, tenemos

df

∆x2 d2 f

f(xi + ∆x) = f(xi ) + ∆x + dx
xi 2! dx2
εp

(1.2)

de esta ecuación obtenemos la siguiente expresión para aproximación de la primera derivada

df

f(xi + ∆x) − f(xi ) ∆x2 d2 f

= − (1.3) dx
xi ∆x 2! dx2
εp

en este caso la aproximación de f ′ (x) mediante diferencias progresivas es de primer orden, o sea O(∆x). Siendo Op (∆x) el error local de truncamiento, definido como
∆x2 d2 f

Op (∆x) = − . (1.4) 2! dx2
εp Es común escribir la anterior expresión como
df

f (xi + ∆x) − f (xi ) = − Op (∆x) dx
xi ∆x

como

f ′ (xi ) =

fi+1 − fi ∆x

(1.5)

(1.6)

para simplificar notación.

[email protected]

2

Antonio Carrillo Ledesma

Introducción al Método de Diferencias Finitas y su Implementación Computacional

1.1.2

Diferencias Regresivas

Considerando la Ec.(1.1) con k = 2 y x = xi − ∆x, tenemos

df

∆x2 d2 f

f (xi − ∆x) = f (xi ) − ∆x + dx
xi 2! dx2
εr

(1.7)

de esta ecuación obtenemos la siguiente expresión para aproximación de la primera derivada

df

f (xi ) − f (xi − ∆x) ∆x2 d2 f

− = (1.8) dx
xi ∆x 2! dx2
εr

en este caso la aproximación de f ′ (x) mediante diferencias regresivas es de primer orden, o sea O(∆x). Siendo Or (∆x) el error local de truncamiento, definido como
∆x2 d2 f

. (1.9) Or (∆x) = 2! dx2
εr

Es común escribir la anterior expresión como
df

f (xi ) − f(xi − ∆x) = + Or (∆x) dx
xi ∆x

como

f ′ (xi ) =

fi − fi−1 ∆x

(1.10)

(1.11)

para simplificar notación. 1.1.3

Diferencias Centradas

Considerando la Ec.(1.1) con k = 3 y escribiendo f(x) en x = xi + ∆x y x = xi − ∆x, tenemos


df

∆x2 d2 f

∆x3 d3 f

f (xi + ∆x) = f (xi ) + ∆x + + (1.12) dx
xi 2! dx2
xi 3! dx3
εp y

f(xi − ∆x) = f(xi ) − ∆x




df

∆x2 d2 f

∆x3 d3 f

+ − dx
xi 2! dx2
xi 3! dx3
εr

restando la Ec.(1.12) de la Ec.(1.13), se tiene 


 df

∆x3 d3 f

d3 f

f (xi + ∆x) − f (xi − ∆x) = 2∆x + + dx
xi 3! dx3
εp dx3
εr

(1.13)

(1.14)

esta última expresión lleva a la siguiente aproximación de la primera derivada mediante diferencias centradas
df

f(xi + ∆x) − f(xi − ∆x) = + Oc (∆x2 ) (1.15)
dx 2∆x xi

[email protected]

3

Antonio Carrillo Ledesma

Introducción al Método de Diferencias Finitas y su Implementación Computacional

con un error local de truncamiento de segundo orden Oc (∆x2 ), es decir 

 2 3
3
∆x d f d f
+
Oc (∆x2 ) = (1.16) 3! dx3
εp dx3
εr

comparado el error local de truncamiento de la aproximación anterior Oc (∆x2 ), con los obtenidos previamente para diferencias progresivas y regresivas Op (∆x) y Or (∆x), se tiene que lim Oc (∆x2 ) < lim Op (∆x).

∆x→0

∆x→0

Es común encontrar expresada la derivada1
f (xi + ∆x) − f (xi − ∆x) df

= dx
xi 2∆x

como

f ′ (xi ) =

fi+1 − fi−1 2∆x

(1.17)

(1.18)

(1.19)

para simplificar notación. Derivadas Usando Más Puntos Utilizando el valor de la función en más puntos se construyen fórmulas más precisas para las derivadas2 , algunos ejemplos son f ′ (xi ) = f ′ (xi ) = f ′ (xi ) = f ′ (xi ) = f ′ (xi ) =

1.2

−3fi + 4fi+1 − fi+2 + O(∆x2 ) (1.20) 2∆x 3fi − 4fi−1 + fi−2 + O(∆x2 ) 2∆x 2fi+1 + 3fi − 6fi−1 + fi−2 + O(∆x3 ) 6∆x fi−2 − 8fi−1 + 8fi+1 − fi+2 + O(∆x4 ) 12∆x −25fi + 48fi+1 − 36fi+2 + 16fi+3 − 3fi+4 + O(∆x4 ) 12∆x

Derivadas de Ordenes Mayores

De forma análoga se construyen aproximaciones en diferencias finitas de orden mayor, aquí desarrollaremos la forma de calcular la derivada de orden dos en diferencias centradas. 1 En el caso de que la derivada sea usada en una malla no homogénea, es necesario incluir en la derivada a que ∆x se refiere, por ejemplo
en cada punto i, tenemos la ∆xi− (por la f (xi +∆xi− )−f (xi −∆xi+ ) df
izquierda) y la ∆xi+ (por la derecha), i.e. dx .
= (∆xi− )+(∆xi+ ) xi 2 Al usar estas derivadas en el método de diferencias finitas mostrado en la sección (2) las matrices generadas no serán tridiagonales.

[email protected]

4

Antonio Carrillo Ledesma

Introducción al Método de Diferencias Finitas y su Implementación Computacional

1.2.1

Derivada de Orden Dos

Partiendo del desarrollo de Taylor f (xi +∆x) = f(xi )+∆xf ′ (xi )+

∆x2 ′′ ∆x3 ′′′ ∆x4 (4) f (xi )+ f (xi )+ f (ξ p ) (1.21) 2! 3! 4!

y f (xi −∆x) = f(xi )−∆xf ′ (xi )+

∆x2 ′′ ∆x3 ′′′ ∆x4 (4) f (xi )− f (xi )+ f (ξ r ) (1.22) 2! 3! 4!

y eliminado las derivadas primeras, sumando las ecuaciones anteriores y despejando se encuentra que f ′′ (xi ) =

f(xi − ∆x) − 2f(xi ) + f(xi + ∆x) ∆x2 (4) − f (ξ c ) ∆x2 12

(1.23)

así, la aproximación a la segunda derivada usando diferencias centradas con un error de truncamiento Oc (∆x2 ) es f ′′ (xi ) =

f(xi − ∆x) − 2f (xi ) + f(xi + ∆x) ∆x2

(1.24)

Es común escribir la anterior expresión como f ′′ (xi ) =

fi−1 − 2fi + fi+1 ∆x2

para simplificar notación. Derivadas Usando Más Puntos Utilizando el valor de la función en más puntos se construyen fórmulas más precisas para las derivadas, algunos ejemplos son f ′′ (xi ) = f ′′ (xi ) = f ′′ (xi ) = 1.2.2

fi+2 − 2fi+1 + fi + O(∆x) ∆x2 −fi+3 + 4fi+2 − 5fi+1 + 2fi + O(∆x2 ) ∆x2 −fi+2 + 16fi+1 − 30fi + 16fi−1 − fi−2 + O(∆x4 ) 12∆x2

(1.25)

Derivadas de Ordenes Mayores

De forma análoga se construyen derivadas de ordenes mayores utilizando el valor de la función en más puntos, algunos ejemplos para derivadas terceras son f ′′′ (xi ) = f ′′′ (xi ) = f ′′′ (xi ) =

fi+3 − 3fi+2 + 3fi+1 − fi + O(∆x) (1.26) ∆x3 fi+2 − 2fi+1 + 2fi−1 − fi−2 + O(∆x2 ) 2∆x3 fi−3 − 8fi−2 + 13fi−1 − 13fi+1 + 8fi+2 − fi+3 + O(∆x4 ) 8∆x3

[email protected]

5

Antonio Carrillo Ledesma

Introducción al Método de Diferencias Finitas y su Implementación Computacional

Algunos ejemplos para derivadas cuartas son f ′′′′ (xi ) = f ′′′′ (xi ) = f ′′′′ (xi ) =

1.3

fi+4 − 4fi+3 + 6fi+2 − 4fi+1 + fi + O(∆x) (1.27) ∆x4 fi+2 − 4fi+1 + 6fi − 4fi−1 + fi−2 + O(∆x2 ) ∆x4 −fi+3 + 12fi+2 − 39fi+1 + 56fi − 39fi−1 + 12fi−2 − fi−3 + O(∆x4 ) 6∆x4

Derivadas en Dos Dimensiones

De forma análoga se construyen aproximaciones en diferencias finitas de primer y segundo orden en dos dimensiones. Usando el teorema de Teylor para funciones en dos variables x y y, es posible escribir de forma exacta para el punto xi y yj ∂f (xi , yj ) ∆x ∂ 2 f (xi + θ 1 ∆x, yj ) + ∂x 2 ∂x2 (1.28) ∂f (xi , yj ) ∆y ∂ 2 f (xi , yj + θ2 ∆y) f (xi , yj + ∆y) = f (xi , yj ) + ∆y . + ∂y 2 ∂y 2

f(xi + ∆x, yj ) = f (xi , yj ) + ∆x

Así, la aproximación en diferencias hacia delante de ∂f /∂x y ∂f /∂y es ∂f (xi , yj ) ∂x ∂f (xi , yj ) ∂y

≃ ≃

f (xi + ∆x, yj ) − f (xi , yj ) ∆x f (xi , yj + ∆y) − f (xi , yj ) ∆y

(1.29)

o en su forma simplificada (para simplificar la notación, asociamos ∆x = h y ∆y = k), tenemos ∂f (xi , yj ) ∂x ∂f (xi , yj ) ∂y

≃ ≃

fi+1,j − fi,j h fi,j+1 − fi,j . k

(1.30)

La aproximación en diferencias hacia atrás de ∂f /∂x y ∂f /∂y es ∂f (xi , yj ) ∂x ∂f (xi , yj ) ∂y

≃ ≃

f (xi , yj ) − f (xi − ∆x, yj ) ∆x f (xi , yj ) − f (xi , yj − ∆y) ∆y

(1.31)

o en su forma simplificada, tenemos ∂f (xi , yj ) ∂x ∂f (xi , yj ) ∂y [email protected]

≃ ≃ 6

fi,j − fi−1,j h fi,j − fi,j−1 . k

(1.32)

Antonio Carrillo Ledesma

Introducción al Método de Diferencias Finitas y su Implementación Computacional

La aproximación en diferencias centradas de ∂f /∂x y ∂f /∂y es ∂f (xi , yj ) ∂x ∂f (xi , yj ) ∂y

≃ ≃

f (xi + ∆x, yj ) − f (xi − ∆x, yj ) 2∆x f (xi , yj + ∆y) − f (xi , yj − ∆y) 2∆y

(1.33)

o en su forma simplificada, tenemos ∂f (xi , yj ) ∂x ∂f (xi , yj ) ∂y

fi+1,j − fi−1,j 2h fi,j+1 − fi,j−1 . 2k

≃ ≃

(1.34)

Por otro lado, la aproximación en diferencias centradas de ∂ 2 f/∂x2 y ∂ 2 f/∂y 2 es ∂ 2 f (xi , yj ) ∂x2 2 ∂ f (xi , yj ) ∂y2

≃ ≃

f (xi + ∆x, yj ) − 2f (xi , yj ) + f (xi − ∆x, yj ) (1.35) ∆x2 f (xi , yj + ∆y) − f (xi , yj ) + f (xi , yj − ∆y) ∆y2

o en su forma simplificada, tenemos ∂ 2 f (xi , yj ) ∂x2 2 ∂ f (xi , yj ) ∂y 2

1.4

≃ ≃

fi+1,j − 2fi,j + fi−1,j h2 fi,j+1 − 2fi,j + fi,j−1 . k2

(1.36)

Derivadas en Tres Dimensiones

De forma análoga se construyen aproximaciones en diferencias finitas de primer y segundo orden en tres dimensiones. Usando el teorema de Teylor para funciones de tres variables x, y y z, es posible escribir de forma exacta para el punto xi , yj y zk ∂f (xi , yj , zk ) ∆x ∂ 2 f (xi + θ1 ∆x, yj , zk ) + ∂x 2 ∂x2 (1.37) ∂f (xi , yj , zk ) ∆y ∂ 2 f (xi , yj + θ2 ∆y, zk ) f(xi , yj + ∆y, zk ) = f(xi , yj , zk ) + ∆y + ∂y 2 ∂y 2 2 ∂f (xi , yj , zk ) ∆z ∂ f (xi , yj , zk + θ 3 ∆z) f (xi , yj , zk + ∆z) = f(xi , yj , zk ) + ∆z + ∂z 2 ∂z 2

f(xi + ∆x, yj , zk ) = f(xi , yj , zk ) + ∆x

Así, la aproximación en diferencias hacia delante de ∂f /∂x, ∂f /∂y y ∂f /∂z

[email protected]

7

Antonio Carrillo Ledesma

Introducción al Método de Diferencias Finitas y su Implementación Computacional

es ∂f (xi , yj , zk ) ∂x ∂f (xi , yj , zk ) ∂y ∂f (xi , yj , zk ) ∂z

≃ ≃ ≃

f (xi + ∆x, yj , zk ) − f (xi , yj , zk ) ∆x f (xi , yj + ∆y, zk ) − f (xi , yj , zk ) ∆y f (xi , yj , zk + ∆z) − f (xi , yj , zk ) ∆z

(1.38)

o en su forma simplificada (para simplificar la notación, asociamos ∆x = h, ∆y = l y ∆z = m), tenemos ∂f (xi , yj , zk ) ∂x ∂f (xi , yj , zk ) ∂y ∂f (xi , yj , zk ) ∂z

≃ ≃ ≃

fi+1,j,k − fi,j,k h fi,j+1,k − fi,j,k l fi,j,k+1 − fi,j,k . m

(1.39)

La aproximación en diferencias hacia atrás de ∂f /∂x, ∂f /∂y y ∂f /∂z es ∂f (xi , yj , zk ) ∂x ∂f (xi , yj , zk ) ∂y ∂f (xi , yj , zk ) ∂z

≃ ≃ ≃

f (xi , yj , zk ) − f (xi − ∆x, yj , zk ) ∆x f (xi , yj , zk ) − f (xi , yj − ∆y, zk ) ∆y f (xi , yj , zk ) − f (xi , yj , zk − ∆z) ∆z

(1.40)

o en su forma simplificada, tenemos ∂f (xi , yj , zk ) ∂x ∂f (xi , yj , zk ) ∂y ∂f (xi , yj , zk ) ∂z

≃ ≃ ≃

fi,j,k − fi−1,j,k h fi,j,k − fi,j−1,k l fi,j,k − fi,j,k−1 . m

(1.41)

La aproximación en diferencias centradas de ∂f /∂x, ∂f/∂y y ∂f /∂z es ∂f (xi , yj , zk ) ∂x ∂f (xi , yj , zk ) ∂y ∂f (xi , yj , zk ) ∂z

≃ ≃ ≃

f (xi + ∆x, yj , zk ) − f (xi − ∆x, yj , zk ) 2∆x f (xi , yj + ∆y, zk ) − f (xi , yj − ∆y, zk ) 2∆y f (xi , yj , zk + ∆z) − f (xi , yj , zk − ∆z) 2∆z

[email protected]

8

(1.42)

Antonio Carrillo Ledesma

Introducción al Método de Diferencias Finitas y su Implementación Computacional

o en su forma simplificada, tenemos ∂f (xi , yj , zk ) ∂x ∂f (xi , yj , zk ) ∂y ∂f (xi , yj , zk ) ∂z

fi+1,j,k − fi−1,j,k 2h fi,j+1,k − fi,j−1,k 2l fi,j,k+1 − fi,j,k−1 . 2m

≃ ≃ ≃

(1.43)

Por otro lado, la aproximación en diferencias centradas de ∂ 2 f /∂x2 , ∂ 2 f /∂y 2 y ∂ f /∂z 2 es 2

∂ 2 f (xi , yj , zk ) ∂x2 ∂ 2 f (xi , yj , zk ) ∂y 2 2 ∂ f (xi , yj , zk ) ∂z 2

≃ ≃ ≃

f (xi + ∆x, yj , zk ) − 2f (xi , yj , zk ) + f (xi − ∆x, yj , zk ) ∆x2 (1.44) f (xi , yj + ∆y, zk ) − f (xi , yj , zk ) + f (xi , yj − ∆y, zk ) ∆y2 f (xi , yj , zk + ∆z) − f (xi , yj , zk ) + f (xi , yj , zk − ∆z) ∆z 2

o en su forma simplificada, tenemos ∂ 2 f (xi , yj , zk ) ∂x2 2 ∂ f (xi , yj , zk ) ∂y2 ∂ 2 f (xi , yj , zk ) ∂z 2

≃ ≃ ≃

[email protected]

fi+1,j,k − 2fi,j,k + fi−1,j,k h2 fi,j+1,k − 2fi,j,k + fi,j−1,k l2 fi,j,k+1 − 2fi,j,k + fi,j,k−1 . m2

9

(1.45)

Antonio Carrillo Ledesma

Introducción al Método de Diferencias Finitas y su Implementación Computacional

2

Método de Diferencias Finitas en Una Dimensión

Consideremos la ecuación diferencial parcial ′

(p (x) u′ (x)) + q (x) u′ (x) − r (x) u (x) = f (x)

(2.1)

donde: u (a) = uα y u (b) = uβ

en a ≤ x ≤ b

con condiciones de frontera Dirichlet o cualquier otro tipo de condiciones de frontera. Para usar el procedimiento general de solución numérica mediante el método de diferencias finitas, debemos de hacer: 1. Generar una malla del dominio, i.e. una malla es un conjunto finito de puntos en los cuales buscaremos la solución aproximada a la ecuación diferencial parcial. 2. Sustituir las derivadas correspondiente con alguna de las formulas de diferencias finitas centradas (véase secciones 1.1 y 1.2), en cada punto donde la solución es desconocida para obtener un sistema algebraico de ecuaciones Au = f. 3. Resolver el sistema de ecuaciones (véase capítulo 3), y así obtener la solución aproximada en cada punto de la malla.

2.1

Problema con Condiciones de Frontera Dirichlet

Consideremos un caso particular de la Ec.(2.1) definido por la ecuación u′′ (x) = f (x),

0 ≤ x ≤ 1,

u(0) = uα ,

u(1) = uβ

(2.2)

con condiciones de frontera Dirichlet. Para usar el procedimiento general de solución numérica mediante el método de diferencias finitas, debemos de hacer: 1. Generamos una malla homogénea del dominio3 xi = ih,

i = 0, 1, ..., n,

h=

1 = ∆x n

(2.3)

2. Sustituimos la derivada con Ec.(1.24) en cada punto donde la solución es desconocida para obtener un sistema algebraico de ecuaciones. Así, en cada punto xi de la malla aproximamos la ecuación diferencial por4 u′′ (xi ) ≈

u(xi − h) − 2u(x) + u(xi + h) h2

(2.4)

3 En el caso de que la malla no sea homogénea, es necesario incluir en la derivada a que h se refiere, por ejemplo en cada punto i, tenemos la hi− (por la izquierda) y la hi+ (por la

derecha), i.e. u′′ (xi ) ≈

u(xi −hi− )−2u(x)+u(xi +hi+ ) . (hi− )(hi+ )

4 Notemos que en cada punto de la malla, la aproximación por diferencias finitas supone la solución de tres puntos de la malla xi−1 , xi y xi+1 . Al conjunto de estos tres puntos de la malla son comúnmente llamados el esténcil de diferencias finitas.

[email protected]

10

Antonio Carrillo Ledesma

Introducción al Método de Diferencias Finitas y su Implementación Computacional

o en su forma simplificada u′′ (xi ) ≈

ui−1 − 2ui + ui+1 h2

(2.5)

definiendo la solución aproximada de u(x) en xi como ui como la solución del siguiente sistema lineal de ecuaciones uα − 2u1 + u2 h2 u1 − 2u2 + u3 h2 ui−1 − 2ui + ui+1 h2 un−3 − 2un−2 + un−1 h2 un−2 − 2un−1 + uβ h2

= f (x1 ) = f (x2 ) .. . (2.6)

= f (xi ) .. . = f (xn−2 ) = f (xn−1 ).

Este sistema de ecuaciones se puede escribir como un matriz A y los vectores u y fde la forma          

− h22 1

h2



1 h2

− h22 1

h2

1 h2

− h22 ..

.

1 h2

..

.

1 h2

..

. − h22 1

h2

1 h2

− h22

u1 u2 u3 .. .

         un−2 un−1





        =      



f(x1 ) − uhα2 f (x2 ) f (x3 ) .. . f (xn−2 ) u f(xn−1 ) − hβ2

    .   

factorizando 1/h2 del sistema lineal Au = f, tenemos 

−2 1  1 −2 1  1 −2 1    .. h2  .  

1 .. . 1



u1 u2 u3 .. .

     ..  .  −2 1   un−2 1 −2 un−1





        =      

f(x1 ) − uhα2 f (x2 ) f (x3 ) .. . f(xn−2 ) u f(xn−1 ) − hβ2



    .   

esta última forma de expresar el sistema lineal algebraico asociado es preferible para evitar problemas numéricos al momento de resolver el sistema lineal por métodos iterativos (véase capítulo 3.2) principalmente cuando h → 0. [email protected]

11

Antonio Carrillo Ledesma

Introducción al Método de Diferencias Finitas y su Implementación Computacional

3. Resolviendo el sistema de ecuaciones (véase capítulo 3), obtenemos la solución aproximada en cada punto interior de la malla. La solución completa al problema la obtenemos al formar el vector

 uα u1 u2 u3 · · · un−2 un−1 uβ . Para el problema general dado por la Ec.(2.1) ′

(p (x) u′ (x)) + q (x) u′ (x) − r (x) u (x) = f (x) en a ≤ x ≤ b

(2.7)

donde: u (a) = uα y u (b) = uβ

el término q (x) u (x) algunas veces es llamado el término advectivo5 si u es la velocidad. Existen dos diferentes técnicas de discretización que se detallan a continuación: ′

• Si las funciones p(x), q(x) y r(x) son constantes, el esquema de diferencias finitas centradas para todas las derivadas, este está dado por el esténcil pi

ui−1 − 2ui + ui+1 ui+1 − ui−1 − ri ui = fi + qi 2 h 2h

i = 1, 2, ..., n. (2.8)

donde la ventaja de esta discretización, es que el método es de segundo orden de exactitud. La desventaja es que los coeficientes de la matriz generada pueden no ser diagonal dominante si r (x) > 0 y p (x) > 0. Cuando la advección |p (x)| es grande, la ecuación se comporta como la ecuación de onda. • Tomando las funciones p(x, t), q(x, t) y r(x, t) más generales posibles, es necesario hacer una discretización que garantice que el método es de segundo orden de exactitud, esto se logra mediante la siguiente discretización para (p (x, t) u′ (x, t))′ mediante



∂ ∂u ∆x u (x + ∆x, t) − u (x, t) p (x, t) ≃ p x + ,t (2.9) ∂x ∂x 2 ∆x

 ∆x u (x, t) − u (x − ∆x, t) − p x− ,t /∆x 2 ∆x entonces se tiene que i−1 p(ui+1 , t) ui+1h+ui − p (ui−1,t ) ui −u ui+1 − ui−1 h + qi − ri ui = fi (2.10) h 2h

para i = 1, 2, ..., n, (véase [51] pág. 78 y 79). 5 Cuando

la advección es fuerte, esto es cuando |q (x)| es grande, la ecuación se comporta como si fuera una ecuación de onda.

[email protected]

12

Antonio Carrillo Ledesma

Introducción al Método de Diferencias Finitas y su Implementación Computacional

• El esquema mixto, en donde se usa el esquema de diferencias finitas centradas para el término de difusión y el esquema upwind para el término de advección pi

ui−1 − 2ui + ui+1 ui+1 − ui−1 + qi − ri ui 2 h h

= fi , si qi ≥ 0

ui−1 − 2ui + ui+1 ui − ui−1 + qi − ri ui h2 h

= fi , si qi < 0

pi

(2.11)

el propósito es incrementar el dominio de la diagonal. Este esquema es de orden uno de exactitud y es altamente recomendable su uso si |q (x)| ∼ 1/h, en caso de no usarse, se observará que la solución numérica oscila alrededor del cero. Veamos unos ejemplos desarrollados en SCILAB6 Ejemplo 1 Sea −u′′ (x) + u (x) = 0,

0 ≤ x ≤ 1,

u(0) = 0,

u(1) = 1

entonces el programa queda implementado como: a=0; // Inicio dominio c=1; // Fin dominio M=50; // Partición N=M-2; // Nodos interiores h=(c-a)/(M-1); // Incremento en la malla Y0=0; // Condición inicial en el inicio del dominio Y1=1; // Condición inicial en el fin del dominio A=zeros(N,N); // Matriz A b=zeros(N); // Vector b

P=2/(h^2); Q=-1/(h^2)+1/(2*h); R=-1/(h^2)-1/(2*h);

// Primer renglon de la matriz A y vector b A(1,1)=P; A(1,2)=Q; b(1)=-Y0*R; // Renglones intermedios de la matriz A y vector b for i=2:N-1 6 Scilab es un programa open source para el cálculo numérico el cual provee un poderoso ambiente de cálculo para aplicaciones Científicas y de Ingeniería [http://www.scilab.org].

[email protected]

13

Antonio Carrillo Ledesma

Introducción al Método de Diferencias Finitas y su Implementación Computacional

A(i,i-1)=R; A(i,i)=P; A(i,i+1)=Q; end // Relglon final de la matriz A y vector b A(N,N-1)=R; A(N,N)=P; b(N)=-Y1*Q; // Resuleve el sistema lineal Ax=b x=inv(A)*b; // Prepara la graficación xx=zeros(M,1); for i=1:M xx(i)=a+h*(i-1); end yy=zeros(M,1); yy(1)=Y0; // Condición inicial for i=1:N yy(i+1)=x(i); end yy(M)=Y1; // Condición inicial // Grafica la solución de la Ecuación Diferencial Parcial en 1D plot2d(xx,yy)

Ejemplo 2 Sea u′′ (x) = π2 cos (πx) ,

0 ≤ x ≤ 1,

u(0) = 1,

u(1) = −1

entonces el programa queda implementado como: function y=LadoDerecho(x) y=-%pi*%pi*cos(%pi*x); endfunction function y=SolucionAnalitica(x) y=cos(%pi*x); endfunction a=-1; // Inicio dominio c=2; // Fin dominio M=100; // Partición N=M-2; // Nodos interiores h=(c-a)/(M-1); // Incremento en la malla Y0=-1; // Condición inicial en el inicio del dominio [email protected]

14

Antonio Carrillo Ledesma

Introducción al Método de Diferencias Finitas y su Implementación Computacional

Y1=1; // Condición inicial en el fin del dominio A=zeros(N,N); // Matriz A b=zeros(N); // Vector b R=1/(h^2); P=-2/(h^2); Q=1/(h^2); // Primer renglon de la matriz A y vector b A(1,1)=P; A(1,2)=Q; b(1)=LadoDerecho(a)-Y0*R; // Renglones intermedios de la matriz A y vector b for i=2:N-1 A(i,i-1)=R; A(i,i)=P; A(i,i+1)=Q; b(i)=LadoDerecho(a+h*(i-1)); end // Relglon final de la matriz A y vector b A(N,N-1)=R; A(N,N)=P; b(N)=LadoDerecho(a+h*N)-Y1*Q; // Resuleve el sistema lineal Ax=b x=inv(A)*b; // Prepara la graficación xx=zeros(M,1); zz=zeros(M,1); for i=1:M xx(i)=a+h*(i-1); zz(i)=SolucionAnalitica(xx(i)); end yy=zeros(M,1); yy(1)=Y0; // Condición inicial for i=1:N yy(i+1)=x(i); end yy(M)=Y1; // Condición inicial // Grafica la solución de la Ecuación Diferencial Parcial en 1D plot2d(xx,[yy,zz])

[email protected]

15

Antonio Carrillo Ledesma

Introducción al Método de Diferencias Finitas y su Implementación Computacional

2.2

Problema con Condiciones de Frontera Neumann

Consideremos el problema u′′ (x) = f (x),

0≤x≤1

(2.12)

con condiciones de frontera Neumann du du = cte1 en u(0) y = cte2 en u(1) dx dx para usar el procedimiento general de solución numérica mediante el método de diferencias finitas, primeramente debemos de discretizar las condiciones de frontera, una manera seria usar para la primera condición de frontera una aproximación usando diferencias progresivas Ec.(1.5)
du

u(xi + h) − u(xi ) = dx
xi h

quedando

u1 − u0 = cte1 (2.13) h para la segunda condición de frontera una aproximación usando diferencias regresivas Ec.(1.10)
u(xi ) − u(xi − h) du

=
dx xi h

quedando

un − un−1 = cte2 (2.14) h pero el orden de aproximación no seria el adecuado pues estamos aproximando el dominio con diferencias centradas con un error local de truncamiento de segundo orden Oc (∆x2 ), en lugar de ello usaremos diferencias centradas Ec.(1.15) para tener todo el dominio con el mismo error local de truncamiento. Para usar diferencias centradas Ec.(1.15)
du

u(xi + h) − u(xi − h) =
dx xi 2h en el primer nodo necesitamos introducir un punto de la malla ficticio x−1 = (x0 − ∆x) con un valor asociado a u−1 , entonces u1 − u−1 = cte1 2h

(2.15)

así también, en el último nodo necesitamos introducir un punto de la malla ficticio xn+1 = (xn + ∆x) con un valor asociado a un+1 , obteniendo un+1 − un−1 = cte2 . 2h

(2.16)

Estos valores no tienen significado físico alguno, dado que esos puntos se encuentran fuera del dominio del problema. Entonces debemos de hacer: [email protected]

16

Antonio Carrillo Ledesma

Introducción al Método de Diferencias Finitas y su Implementación Computacional

1. Generamos una malla homogénea del dominio xi = ih,

i = 0, 1, ..., n,

h=

1 = ∆x. n

(2.17)

2. Sustituimos la derivada con Ec.(1.24) en cada punto donde la solución es desconocida para obtener un sistema algebraico de ecuaciones. Así, en cada punto xi de la malla aproximamos la ecuación diferencial por u′′ (xi ) ≈

u(xi − h) − 2u(x) + u(xi + h) h2

(2.18)

definiendo la solución aproximada de u(x) en xi como ui como la solución del siguiente sistema lineal de ecuaciones u1 − u−1 2h u0 − 2u1 + u2 h2 ui−1 − 2ui + ui+1 h2 un−2 − 2un−1 + un h2 un+1 − un−1 2h

= cte1 = f (x1 ) .. . = f (xi )

(2.19)

.. . = f (xn−1 ) = cte2 .

3. Resolviendo el sistema de ecuaciones (véase capítulo 3), obtenemos la solución aproximada en cada punto de la malla.

2.3

Problema con Condiciones de Frontera Robin

El método de un punto de la malla ficticio es usado para el manejo de las condiciones de frontera mixtas, también conocidas como condiciones de frontera Robin. Sin perdida de generalidad, supongamos que en x = a, tenemos αu′ (a) + βu (b) = γ donde α = 0. Entonces usando el punto de la malla ficticio, tenemos que α

u1 − u−1 + βun = γ 2h

o u−1 = u1 +

[email protected]

2β 2hγ un − α α

17

Antonio Carrillo Ledesma

Introducción al Método de Diferencias Finitas y su Implementación Computacional

introduciendo esto en términos de diferencias finitas centradas, en x = x0 , entonces se tiene que

2 2β 2 2γ − 2+ un + 2 u1 = f0 + h αh h αh o



1 β 1 f0 γ − 2+ un + 2 u1 = + h αh h 2 αh

lo que genera coeficientes simétricos en la matriz. Consideremos el problema u′′ (x) = f (x),

(2.20)

0≤x≤1

con condiciones de frontera Dirichlet y Neumman y

u(0) = uα

du = cte1 en u(1). dx

respectivamente. Para usar el procedimiento general de solución numérica mediante el método de diferencias finitas, primeramente debemos de expresar la condición de frontera Neumann mediante diferencias centradas Ec.(1.15)
du

u(xi + h) − u(xi − h) = dx
xi 2h

en el último nodo necesitamos introducir un punto de la malla ficticio xn+1 = (xn + ∆x) con un valor asociado a un+1 quedando un+1 − un−1 = cte2 . 2h

(2.21)

Este valor no tiene significado físico alguno, dado que este punto se encuentra fuera del dominio del problema. Entonces debemos de hacer: 1. Generamos una malla homogénea del dominio xi = ih,

i = 0, 1, ..., n,

h=

1 = ∆x. n

(2.22)

2. Sustituimos la derivada con Ec.(1.24) en cada punto donde la solución es desconocida para obtener un sistema algebraico de ecuaciones. Así, en cada punto xi de la malla aproximamos la ecuación diferencial por u′′ (xi ) ≈

u(xi − h) − 2u(x) + u(xi + h) h2

[email protected]

18

(2.23)

Antonio Carrillo Ledesma

Introducción al Método de Diferencias Finitas y su Implementación Computacional

definiendo la solución aproximada de u(x) en xi como ui como la solución del siguiente sistema lineal de ecuaciones uα − 2u1 + u2 h2 u1 − 2u2 + u3 h2 ui−1 − 2ui + ui+1 h2 un−2 − 2un−1 + un h2 un+1 − un−1 2h

= f (x1 ) = f (x2 ) .. . (2.24)

= f (xi ) .. . = f (xn−1 ) = cte1 .

3. Resolviendo el sistema de ecuaciones (véase capítulo 3), obtenemos la solución aproximada en cada punto de la malla. La solución completa al problema la obtenemos al formar el vector 

uα u1 u2 u3 · · · un−2 un−1 un . Veamos unos ejemplos desarrollados en SCILAB7 Ejemplo 3 Sea u′′ (x) = −π2 cos (πx) ,

0 ≤ x ≤ 0.5,

u(0) = 1,

u′ (0.5) = −π

entonces el programa queda implementado como: function y=LadoDerecho(x) y=-%pi*%pi*cos(%pi*x); endfunction function y=SolucionAnalitica(x) y=cos(%pi*x); endfunction a=0; // Inicio dominio c=0.5; // Fin dominio M=40; // Partición N=M-1; // Nodos interiores h=(c-a)/(M-1); // Incremento en la malla Y0=1; // Condición Dirchlet inicial en el inicio del dominio 7 Scilab es un programa open source para el cálculo numérico el cual provee un poderoso ambiente de cálculo para aplicaciones Científicas y de Ingeniería [http://www.scilab.org].

[email protected]

19

Antonio Carrillo Ledesma

Introducción al Método de Diferencias Finitas y su Implementación Computacional

Y1=-%pi; // Condición Neumann inicial en el fin del dominio A=zeros(N,N); // Matriz A b=zeros(N); // Vector b R=1/(h^2); P=-2/(h^2); Q=1/(h^2); // Primer renglon de la matriz A y vector b A(1,1)=P; A(1,2)=Q; b(1)=LadoDerecho(a)-Y0*R; // Frontera Dirichlet // Renglones intermedios de la matriz A y vector b for i=2:N-1 A(i,i-1)=R; A(i,i)=P; A(i,i+1)=Q; b(i)=LadoDerecho(a+h*(i-1)); end // Relglon final de la matriz A y vector b A(N,N-1)=-1/(h^2); A(N,N)=-1/(h^2); b(N)=Y1/h; // Frontera Neumann // Resuleve el sistema lineal Ax=b x=inv(A)*b; // Prepara la graficación xx=zeros(M,1); zz=zeros(M,1); for i=1:M xx(i)=a+h*(i-1); zz(i)=SolucionAnalitica(xx(i)); end yy=zeros(M,1); yy(1)=Y0; // Condición inicial for i=1:N yy(i+1)=x(i); end // Grafica la solución de la Ecuación Diferencial Parcial en 1D plot2d(xx,[yy,zz])

[email protected]

20

Antonio Carrillo Ledesma

Introducción al Método de Diferencias Finitas y su Implementación Computacional

Ejemplo 4 Sea −u′′ (x) − k2 u(x) = 0,

0 ≤ x ≤ 1,

u(0) = 1,

u′ (1) = iku(1)

con k = 150, entonces el programa queda implementado como8 : TEST = 1; // (0) Diferencias finitas, (1) Diferencias finitas exactas segun Yau Shu Wong y Guangrui Li function y=LadoDerecho(x) y=0.0; endfunction function y=SolucionAnalitica(x, k) //y=cos(k*x)+%i*sin(k*x); y=exp(%i*k*x); endfunction K = 150; KK = K*K; a=0; // Inicio dominio c=1; // Fin dominio M=300; // Partición N=M-1; // Nodos interiores h=(c-a)/(M-1); // Incremento en la malla Y0=1; // Condición Dirchlet inicial en el inicio del dominio Y1=%i*K; // Condición Neumann inicial en el fin del dominio A=zeros(N,N); // Matriz A b=zeros(N); // Vector b if TEST = 0 then R=-1/(h^2); P=2/(h^2)-KK; Q=-1/(h^2); else R=-1/(h^2); P=(2*cos(K*h)+(K*h)^2)/(h^2) - KK; Q=-1/(h^2); end 8

Esta ecuación se conoce como la ecuación de Helmholtz la cual es dificil de resolver si r (x) ≤ 0 y |r (x)| es grande, i.e. r (x) ∼ 1/h2 . Para resolver dicha ecuación, se usa diferencias finitas y diferencias finitas exactas, procedimiento desarrollado en: Exact Finite Difference Schemes for Solving Helmholtz equation at any wavenumber, Yau Shu Wong and Guangrui Lim International Journal of Numerical Analysis and Modeling, Volumen 2, Number 1, Pages 91-108, 2011.

[email protected]

21

Antonio Carrillo Ledesma

Introducción al Método de Diferencias Finitas y su Implementación Computacional

// Primer renglon de la matriz A y vector b A(1,1)=P; A(1,2)=Q; b(1)=LadoDerecho(a)-Y0*R; // Frontera dirichlet // Renglones intermedios de la matriz A y vector b for i=2:N-1 A(i,i-1)=R; A(i,i)=P; A(i,i+1)=Q; b(i)=LadoDerecho(a+h*(i-1)); end // Relglon final de la matriz A y vector b if TEST = 0 then A(N,N-1)=1/(h^2); A(N,N)=-1/(h^2)+ Y1/h; b(N)=LadoDerecho(c)/2; else A(N,N-1)=1/(h^2); A(N,N)=-1/(h^2)+ %i*sin(K*h)/(h^2); b(N)=LadoDerecho(c)/2; end // Resuleve el sistema lineal Ax=b x=inv(A)*b; ESC = 5; xxx=zeros(M*ESC,1); zzz=zeros(M*ESC,1); for i=1:M*ESC xxx(i)=a+h/ESC*(i-1); zzz(i)=SolucionAnalitica(xxx(i),K); end // Prepara la graficación xx=zeros(M,1); zz=zeros(M,1); for i=1:M xx(i)=a+h*(i-1); zz(i)=SolucionAnalitica(xx(i),K); end yy=zeros(M,1); yy(1)=Y0; // Condición inicial for i=1:N yy(i+1)=x(i); end // Gráfica la solución de la Ecuación Diferencial Parcial en 1D plot2d(xx,yy,15) [email protected]

22

Antonio Carrillo Ledesma

Introducción al Método de Diferencias Finitas y su Implementación Computacional

plot2d(xxx,zzz)

2.4

Discretización del Tiempo

Hasta ahora se ha visto como discretizar la parte espacial de las ecuaciones diferenciales parciales, lo cual nos permite encontrar la solución estática de los problemas del tipo elíptico. Sin embargo, para ecuaciones del tipo parabólico e hiperbólico dependen del tiempo, se necesita introducir una discretización a las derivadas con respecto del tiempo. Al igual que con las discretizaciones espaciales, podemos utilizar algún esquema de diferencias finitas en la discretización del tiempo. 2.4.1

Ecuación con Primera Derivada Temporal

Para la solución de la ecuaciones diferenciales con derivada temporal (ut ), se emplean diferentes esquemas en diferencias finitas para la discretización del tiempo. Estos esquemas se conocen de manera general como esquemas theta(θ). Definiendo la ecuación diferencial parcial general de segundo orden como

donde

ut = Lu

(2.25)

Lu = (p (x) u′ (x)) + q (x) u′ (x) − r (x) u (x) − f (x)

(2.26)

ut = (1 − θ) (Lu)j + θ (Lu)j+1 .

(2.27)

′

aquí, los coeficientes p, q y r pueden depender del espacio y del tiempo. Entonces el esquema theta está dado por

Existen diferentes casos del esquema theta a saber: • Para θ = 0, se obtiene un esquema de diferencias finitas hacia delante en el tiempo, conocido como esquema completamente explicito, ya que el paso n + 1 se obtiene de los términos del paso anterior n. Es un esquema 2 sencillo, el cual es condicionalmente estable cuando ∆t ≤ ∆x 2 . • Para θ = 1, se obtiene el esquema de diferencias finitas hacia atrás en el tiempo, conocido como esquema completamente implícito, el cual es incondicionalmente estable. • Para θ = 12 , se obtiene un esquema de diferencias finitas centradas en el tiempo, conocido como esquema Crank-Nicolson, este esquema es también incondicionalmente estable y es más usado por tal razón. Para la implementación del esquema Crank-Nicolson se toma una diferencia progresiva para el tiempo y se promedian las diferencias progresivas y regresivas

[email protected]

23

Antonio Carrillo Ledesma

Introducción al Método de Diferencias Finitas y su Implementación Computacional

en el tiempo para las derivadas espaciales. Entonces si tenemos la Ec. (2.26), las discretizaciones correspondientes son: ut ≃

uj+1 − uji i ∆t

(2.28)

 j  j j j+1 j+1 j+1 u − 2u + u − 2u + u u p i−1 i i+1 i−1 i i+1 (p (x) u′ (x)) ≃ + 2 ∆x2 ∆x2  j  j j+1 j+1 q ui−1 + ui+1 ui−1 + ui+1 ′ q (x) u (x) ≃ + 2 2∆x 2∆x ′

(2.29) (2.30)

además de r(x)uji y fij . Entonces, una vez que se sustituyen las derivadas por su forma en diferencias finitas, lo siguiente es formar el sistema Auj+1 = Buj + f j

(2.31)

esto se logra, colocando del lado izquierdo la igualdad de los términos que contengan el paso del tiempo correspondiente a j + 1 y del lado derecho a los correspondientes términos de j. A continuación, veamos un ejemplo del esquema Crank-Nicolson desarrollados en SCILAB9 Ejemplo 5 Sea ut − a(x)u′′ (x) − b(x)u′ (x) + c(x)u = f,

l0 ≤ x ≤ l,

0