Rev. Int. M´ et. Num. C´ alc. Dis. Ing. Vol. 19, 3, 331–340 (2003)

Revista Internacional de M´ etodos Num´ ericos para C´ alculo y Dise˜ no en Ingenier´ıa

Resoluci´ on de ecuaciones diferenciales en derivadas parciales dependientes del tiempo de segundo orden utilizando diferencias finitas generalizadas Francisco Ure˜ na y Juan Jos´e Benito ETS Ingenieros Industriales de la UNED Apdo. de Correos 60149, 28080 Madrid, Espa˜ na Tel.: 34-91-398 64 57; Fax: 34-91-398 60 46 e-mail: [email protected]

Luis Gavete ETS Ingenieros de Minas de Madrid UPM c/ R´ıos Rosas 21, 28003 Madrid, Espa˜ na Tel. 34-91-336 64 64; Fax: 34-91-336 32 30 e-mail: [email protected]

´ Ram´on Alvarez ETS Ingenieros Indsutriales de Madrid UPM c/ Jos´ e Guti´ errez Abascal 3, 28006 Madrid, Espa˜ na Tel.: 34-91-310 51 58; Fax: 34-91-398 60 46

Resumen En este art´ıculo se muestra la eficiencia del m´etodo de diferencias finitas generalizadas (MDFG) en la resoluci´ on, por el m´etodo expl´ıcito, de ecuaciones diferenciales en derivadas parciales dependientes del tiempo de segundo orden para una o dos dimensiones espaciales. La obtenci´ on de f´ ormulas expl´ıcitas en diferencias finitas generalizadas permite establecer un sencillo criterio de estabilidad, que viene expresado en funci´ on de los coeficientes de la ecuaci´ on de la estrella.

Palabras clave: Diferencias finitas generalizadas, m´etodos sin malla, m´ınimos cuadrados m´ oviles. RESOLUTION OF SECOND ORDER PARTIAL DIFFERENTIAL EQUATION TIME DEPENDENTS BY GENERALIZED FINITE DIFFERENCE METHOD

Summary This paper shows the efficience of generalized finite difference method (MDFG) in the solution, by the explicit method, of second order partial differential equation time dependents for one or two space dimensions. The explicit finite difference formulae obtained allows us to establish an easy criterion of stability, which is expressed in function of the coeficients of the star equation.

Key words: Generalized finite differences, meshless methods, moving least squares.

c
Universitat Polit` ecnica de Catalunya (Espa˜ na). ISSN: 0213–1315 Recibido: Abril 2002 Aceptado: Octubre 2002

´ F. Ure˜ na, J.J. Benito, L. Gavete y R. Alvarez

332

´ INTRODUCCION La aparici´ on de potentes computadoras ha permitido la utilizaci´on del m´etodo de diferencias finitas para resolver ecuaciones diferenciales en derivadas parciales en mallas irregulares. Al desarrollo de este m´etodo, incluido en los llamados M´etodos sin Malla1,4,10 , han contribuido los trabajos de Forsythe y Wasow3 , Jensen6 , Perrone y Kao12 , Kaczkowski7 , Liska8,9 y Orkisz11 . Las expresiones expl´ıcitas de las f´ormulas en diferencias finitas utilizadas en el M´etodo de Diferencias Finitas Generalizadas (MDFG), as´ı como la influencia de los par´ ametros que intervienen han sido estudiadas en la referencia 2. En este art´ıculo se resuelven ecuaciones diferenciales en derivadas parciales dependientes del tiempo, en una y dos dimensiones espaciales, utilizando el m´etodo de diferencias finitas generalizadas y el m´etodo expl´ıcito, desarroll´ andose, tambi´en un criterio de estabilidad. ´ FORMULAS EXPL´ ICITAS EN DIFERENCIAS FINITAS GENERALIZADAS En este apartado se presentan las expresiones lineales expl´ıcitas que aproximan las derivadas parciales en un punto del dominio (f´ ormulas en diferencias finitas) con el objeto de sustituir dichas expresiones en la siguiente ecuaci´ on diferencial en derivadas parciales de segundo orden Caso unidimensional ∂U ∂2U ∂U = A 2 +B ∂t ∂x ∂x con la condici´on de contorno




t > 0,

a≤x≤b

(1.a)

U (a, t) = g(t) t>0

(1.b)

U (b, t) = h(t) con la condici´on inicial U (x, 0) = f (x)

(1.c)

donde g(t), h(t) y f (x) son funciones conocidas. Caso bidimensional ∂ 2U ∂U ∂U ∂2U ∂2U ∂U = A 2 +B 2 +C +D +E ∂t ∂x ∂y ∂x∂y ∂x ∂x

t > 0,

(x, y) ∈ Ω ⊂ R2

(2.a)

con la condici´on de contorno α

∂U + βU = G(t) ∂n

en

Γ

(2.b)

con la condici´on inicial U (x, y, 0) = F (x, y) donde G(t) y F (x, y) son funciones conocidas y Γ es la frontera de Ω.

(2.c)

333

Resoluci´ on de ecuaciones diferenciales en derivadas parciales

A cada nodo de la malla se le asocia un n´ umero de nodos N , tambi´en de la malla, escogidos siguiendo un determinado criterio (distancia y cuadrante)2 . Al conjunto de estos N + 1 nodos se le denomina estrella y al nodo a partir del cual se define la estrella, nodo central. El criterio de la distancia consiste en seleccionar los N nodos m´ as pr´ oximos al central, mientras que en el criterio del cuadrante se eligen N/4 nodos m´ as pr´ oximos de cada uno de los cuadrantes definidos por unos ejes cartesianos situados con el origen en el nodo central (caso bidimensional) o N/2 en cada semieje a partir del nodo central (caso unidimensional). Si U0 es el valor de la funci´ on en el nodo central y Ui el valor de la funci´ on en el resto de los nodos, con i = 1, . . . , N , entonces, de acuerdo con el desarrollo del polinomio de Taylor se tiene para el caso unidimensional (3.a) y para el caso bidimensional (3.b) ∂U0 h2i ∂ 2 U0 + ... + ∂x 2 ∂x2   2 2 ∂U0 ∂U0 1 ∂ 2 U0 2 ∂ U0 2 ∂ U0 Ui = U0 + hi + ki + 2hi ki + ki + hi + ... ∂x ∂y 2 ∂x2 ∂y 2 ∂x∂y Ui = U0 + hi

(3.a)

(3.b)

donde para el caso unidimensional se tiene que x0 es la coordenada del nodo central y xi las coordenadas del nodo i de la estrella y hi = xi − x0 . Para el caso bidimensional (x0 , y0 ) son las coordenadas del nodo central, (xi , yi ) las coordenadas del nodo i de la estrella y hi = xi − x0 , ki = yi − y0 . Si en (3.a) y (3.b) se truncan los segundos miembros a partir de los t´erminos superiores al de segundo orden, se obtienen las aproximaciones de segundo orden del valor Ui que se indicar´a por ui (tanto en el caso unidimensional como en el bidimensional). De esta manera se pueden definir las funciones B2 (u) y B5 (u) de la forma siguiente 2  N   ∂u0 h2i ∂ 2 u0 + B2 (u) = u0 − ui + hi w(hi ) ∂x 2 ∂x2 i=1

B5 (u) =

N   i=1

∂u0 ∂u0 h2i ∂ 2 u0 ki2 ∂ 2 u0 ∂ 2 u0 + ki + + + h k u0 − ui + hi i i ∂x ∂y 2 ∂x2 2 ∂y 2 ∂x∂y

(4.a)

2

 w(hi , ki )

(4.b) donde w(hi ) y w(hi , ki ) son funciones de ponderaci´ on. Si se minimizan las normas (4.a) y (4.b) con respecto a las derivadas parciales, se obtienen los siguientes sistemas de ecuaciones lineales, que en forma matricial son A2 Du2 = b2

(5.a)

A5 Du5 = b5

(5.b)

Las matrices A2 y A5 tienen dimensiones de 2 × 2 y 5 × 5 respectivamente, de ah´ı los sub´ındices empleados y los vectores Du2 y Du5 son respectivamente  T ∂u0 ∂ 2 u0 , Du2 = (6.a) ∂x ∂x2  Du5 =

∂u0 ∂u0 ∂ 2 u0 ∂ 2 u0 ∂ 2 u0 , , , , ∂x ∂y ∂x2 ∂y 2 ∂x∂y

T (6.b)

´ F. Ure˜ na, J.J. Benito, L. Gavete y R. Alvarez

334

Puesto que las matrices A2 y A5 son sim´etricas, es posible utilizar el m´etodo de Cholesky para resolver los sistemas de ecuaciones (5.a) y (5.b), descomponi´endolas en producto de matrices triangulares A2 = L2LT 2

(7.a)

A5 = L5LT 5

(7.b)

donde los elementos de las matrices L2 y L5 son denotados por l(i, j). Una vez resueltos los sistemas (5.a) y (5.b) se obtienen las expresiones expl´ıcitas de las f´ ormulas en diferencias finitas generalizadas de segundo orden2 con P = 2 para el caso unidimensional y P = 5 para el caso bidimensional DuP

1 = l(k, k)

−u0

N 

M (k, i)ci +

i=1

N  j=1

uj

5 

M (k, i)dji

(k = 1, . . . , P ) (8)

i=1

1  l(i, k)M (k, j) l(i, i) k=j i−1

M (i, j) = (−1)1−δij

M (i, j) =

1 l(i, i)

M (i, j) = 0

con

j < i (i, j = 1, . . . , P )

con

j = i (i, j = 1, . . . , P )

con

j > i (i, j = 1, . . . , P )

con δij funci´ on delta de Kronecker y

ci =

N 

2

dji , dj1 = hj W , dj2

j=1

h2j 2 kj2 2 2 = W , dj3 = kj W , dj4 = W , dj5 = hj kj W 2 2 2

en los que W = w(hi ) para el caso unidimensional y W = w(hi , ki ) para el bidimensional. Si la expresi´ on (8) se sustituye en la ecuaci´on (1.a) o (2.a) seg´ un sea el caso unidimensional o bidimensional, se obtiene la expresi´ on siguiente N

 n+1 n n u0 = u0 (1 − ∆tm0 ) + ∆t mi ui (9) i=1

el valor de la soluci´ on num´erica en el nodo central en el donde se ha designado por un+1 0 instante n + 1; un0 el valor de la soluci´ on num´erica en el nodo central en el instante n; uni el valor de la soluci´ on num´erica de los nodos de la estrella en el instante n; m0 coeficiente del valor de la soluci´ on num´erica en el nodo central un0 , que para el caso de la ecuaci´ on (1.a) viene dado por la expresi´ on  N   M (1, i) M (2, i) +B (10) A m0 = ci l(2, 2) l(1, 1) i=1 y el de la ecuaci´ on (2.a) por  N   M (4, i) M (5, i) M (1, i) M (3, i) M (2, i) +B +C +D +E m0 = A ci l(2, 2) l(4, 4) l(5, 5) l(1, 1) l(3, 3) i=1

(11)

335

Resoluci´ on de ecuaciones diferenciales en derivadas parciales

donde mi es el coeficiente del valor de la soluci´on num´erica del nodo i en la estrella cuyo nodo central tiene por soluci´ on num´erica un0 . Para el caso de la ecuaci´on (1.a) viene dado por la expresi´ on  N   M (1, i) M (2, i) +B mi = (12) A dij l(2, 2) l(1, 1) i=1 y para el de la ecuaci´ on (2.a) por   N  M (4, i) M (5, i) M (1, i) M (3, i) M (2, i) mi = +B +C +D +E A dij l(2, 2) l(4, 4) l(5, 5) l(1, 1) l(3, 3) i=1

(13)

CRITERIO DE ESTABILIDAD EN DIFERENCIAS FINITAS GENERALIZADAS Un m´etodo de diferencias para un problema con valor inicial en la frontera, como se sabe, es estable en el sentido von Neumann, si cualquier soluci´ on para la ecuaci´ on en diferencias (9) de la forma unj = ξ n eiαlj h

(14)

|ξ| ≤ 1

(15)

tiene la propiedad Si se aplica el criterio de von Neumann a la ecuaci´ on (9) de la siguiente manera un0 = ξ n eiαlh , ξ n+1 eiαlh

unj = ξ n eiαlj h , lj = 1 + λj

N  = ξ n eiαlh (1 − ∆tm0 ) + ∆t ξ n eiαlh eiαλj h mj

(16) (17)

j=1

Al simplificar, efectuar operaciones y tomar valor absoluto, se tiene  N   iαλ h  |ξ| = 1 + ∆t e j mj − m0 

(18)

j=1

teniendo en cuenta (15) y que m0 =

N 

mj

(19)

j=1

se tiene 0 ≤ ∆t ≤

1 2  ⇒ 0 ≤ ∆t ≤  N   |m0 | mj (1 − eiαλj h )  

(20)

j=1

Puesto que las mallas son irregulares, los m0 son distintos para cada estrella por lo tanto se escoge como criterio de estabilidad el siguiente 1 0 ≤ ∆t ≤ (21) |m0 |max

´ F. Ure˜ na, J.J. Benito, L. Gavete y R. Alvarez

336

´ RESULTADOS NUMERICOS Se presenta en este apartado la resoluci´ on de ecuaciones diferenciales en derivadas parciales dependientes del tiempo tanto en el caso unidimensional como en el bidimensional

CASO UNIDIMENSIONAL Resoluci´on de la ecuaci´ on ∂ 2 u ∂u ∂u = 0, 0001 2 − ∂t ∂x ∂x

t>0 0