Capítulo 10 Ecuaciones diferenciales en derivadas parciales 10.1.

Clasificación de las ecuaciones en derivadas parciales

Los problemas en derivadas parciales que aparecen en Física se suelen clasificar en tres tipos principales: problemas parabólicos, elípticos e hiperbólicos. Los problemas parabólicos corresponden a ecuaciones del tipo de la de difusión del calor ! 2 " ∂ T (x, y, z,t) ∂ T (x, y, z,t) ∂ 2 T (x, y, z,t) ∂ 2 T (x, y, z,t) =k + + ∂t ∂ x2 ∂ y2 ∂ z2 Los problemas hiperbólicos a ecuaciones similares a la ecuación de ondas ! 2 " ∂ 2 ψ(x, y, z,t) ∂ 2 ψ(x, y, z,t) ∂ 2 ψ(x, y, z,t) 2 ∂ ψ(x, y, z,t) =c + + ∂t 2 ∂ x2 ∂ y2 ∂ z2

mientras que los problemas elípticos corresponden a las ecuaciones de Poisson

∂ 2 !(x, y, z,t) ∂ 2 !(x, y, z,t) ∂ 2 !(x, y, z,t) ρ(x, y, z,t) + + =− 2 2 2 ∂x ∂y ∂z ε0 y Laplace

∂ 2 !(x, y, z,t) ∂ 2 !(x, y, z,t) ∂ 2 !(x, y, z,t) + + =0 ∂ x2 ∂ y2 ∂ z2 para el potencial eléctrico. En general, si consideramos una ecuación diferencial en derivadas parciales en dos dimensiones de la forma a

∂ 2 ψ(x, y) ∂ 2 ψ(x, y) ∂ 2 ψ(x, y) ∂ ψ(x, y) ∂ ψ(x, y) + c + b +d +e + f ∂ ψ(x, y) + g = 0 2 2 ∂x ∂y ∂ x∂ y ∂x ∂x

se denomina hiperbólica si b2 − 4ac > 0, parabólica si b2 − 4ac = 0 y elíptica si b2 − 4ac < 0. En Física, la dimensión temporal es la única esencialmente distinta de las tres dimensiones espaciales en cuanto a su comportamiento, por lo que en caso de varias variables el tipo se refiere a la variable temporal y una de las espaciales. 203

204

CAPÍTULO 10. ECUACIONES DIFERENCIALES EN DERIVADAS PARCIALES

Ejercicio: Comprobar que la anterior definición general coincide con los tipos de las ecuaciones de ondas, difusión y Poisson avanzados anteriormente.

10.2.

La ecuación de propagación del calor (ecuación de difusión)

Vamos a estudiar en primer lugar la ecuación de conducción del calor, ∂ T (x,t) ∂ 2 T (x,t) =κ ∂t ∂ x2 que es más ilustrativa y sencilla al mismo tiempo. Esta ecuación representa la evolución del perfil de temperatura T (x,t) en función del tiempo en una barra de coeficiente térmico κ . Una ecuación similar es la ecuación de difusión, que da la concentración C(x,t) en función del tiempo de un soluto con coeficiente de difusión D(x) en un disolvente : ! " ∂C(x,t) ∂ ∂C(x,t) = D(x) ∂t ∂x ∂x Esta ecuación se reduce a la ecuación de propagación del calor si el coeficiente de difusión D es constante. Consideraremos en lo que sigue sólo ejemplos en dos dimensiones por sencillez, para evitar supermatrices en tres y cuatro dimensiones y para mantenerse en un volumen de cálculo moderado. Como tenemos que discretizar el espacio y el tiempo, introducimos la abreviación T (xi ,tn ) = Tin Vamos a considerar condiciones de contorno de Dirichlet T (x = −L/2,t) = Ta T (x = L/2,t) = Tb así como condiciones iniciales: T (x, 0) = f (x) Para discretizar el espacio y el tiempo, tomamos un paso de integración espacial de h y un paso temporal τ, con lo que tn = nτ y xi = ih − L/2. En general, estos espaciados son distintos, y como veremos, hay relaciones entre ellos que optimizan la precisión numérica del resultado. Supondremos que los índices i y n comienzan en 0, de acuerdo con la convención en los lenguajes C y C++. Tomaremos condiciones de contorno nulas por simplicidad.

10.2. LA ECUACIÓN DE PROPAGACIÓN DEL CALOR (ECUACIÓN DE DIFUSIÓN) 205 El primer esquema de integración que vamos a considerar es el llamado de diferencias finitas progresivas en el tiempo y centradas en el espacio, abreviado como FCTS ( de Forward Time Centered Space). Tomamos como aproximación a la derivada temporal # T (tn + τ, xi ) − T (tn , xi ) Tin+1 − Tin ∂ T (x,t) ## = = ∂t #tn ,xi τ τ

mientras que para la derivada espacial tomamos diferencias centradas

n − 2T n + T n ∂ 2 T (x,t) T (tn , xi + h) − 2T (tn , xi ) + T (tn , xi − h) Ti+1 i i−1 = = 2 2 2 ∂x h h

con lo que la versión discretizada de la ecuación de difusión queda de la forma n − 2T n + T n Ti+1 Tin+1 − Tin i i−1 =κ 2 τ h

De aquí podemos despejar los valores de la temperatura en un punto del espacio en un instante dado en función de la temperatura de este punto y los vecinos en instantes anteriores: % κτ $ n n Tin+1 = Tin + 2 Ti+1 − 2Tin + Ti−1 h que en forma matricial queda como con

Tn+1 = ATn 

   Tn =   

T0n = Ta T1n .. .

      

n TN−1 TNn = Tb κτ A=I+ 2 D h donde I es la matriz identidad y D la siguiente matriz tridiagonal   0 0 0 ··· 0 0  1 −2 1 · · · 0 0     0 1 −2 1 · · · 0     .. . . . . . . . .  D= . . . . . 0     0 · · · 1 −2 1 0     0 ··· 0 1 −2 1  0 ··· 0 0 0 0

De hecho, las filas y columnas primera y última corresponden a las condiciones de contorno, y no hace falta calcularlas en cada iteración, ya que las componentes 0 y N de T permanecen invariables, en caso de que las condiciones de contorno sean constantes en el tiempo, que es el caso

206

CAPÍTULO 10. ECUACIONES DIFERENCIALES EN DERIVADAS PARCIALES

que consideraremos en adelante. Cada paso de integración en el tiempo equivale a multiplicar el vector de valores iniciales una vez por la matriz A, con lo que tenemos Tn+1 = An+1 T0 Si descomponemos T0 en vectores propios de A, T0 = "k ak vk , tenemos Tn+1 = " λkn+1 ak vk

(10.1)

k

De esta expresión podemos deducir dos conclusiones: 1. Para valores grandes del tiempo (n), la solución viene dominada por el vector propio correspondiente al valor propio de mayor valor absoluto λm : Tn+1 " λmn+1 am vm 2. Si |λm | > 1, la solución divergerá para tiempos elevados, por lo que decimos que la solución es inestable: Tin+1 " λmn+1 am vmi −→ # Los valores propios dependen del espaciado en el espacio y el tiempo. Su número y valor depende del tamaño de la red de integración-

10.2.1.

Estabilidad del método de diferencias finitas progresivas

De lo dicho en el apartado anterior, vemos que el método FTCS es estable sólamente si el mayor valor propio de A cumple |λm | ≤ 1. Al valor absoluto del mayor valor propio de una matriz se le denomina radio espectral ρ(A) = |λm |

De 10.1vemos que ρ(An ) = ρ(A)n . Una forma numérica de calcular el radio espectral es calcular An v para valores sucesivos de n y para el vector v = (1, 1, . . . , 1). Cuando el cociente entre dos potencias sucesivas converja a una constante, calculamos el vector propio correspondiente al radio espectral como An v vm = n ||A v|| y λm = Avm De hecho, en el caso del método FTCS el radio espectral se puede obtener analíticamente. Escribiendo la matriz A como   1 0 ··· 0  α 1 − 2α α ···     .. .. .. ..  A= . . . .  = I + αD    ··· α 1 − 2α α  ··· 0 0 1

10.2. LA ECUACIÓN DE PROPAGACIÓN DEL CALOR (ECUACIÓN DE DIFUSIÓN) 207 con

yα=



   D=  

 0 0 ··· 0 α −2α α ···   .. .. .. ..  . . . .   ··· α −2α α  ··· 0 0 0

κτ , se puede demostrar que los valores propios de A son h2 ! " 2 kπ k = 1, . . . , N − 1 λk = 1 − 4α sin 2N

por lo que la condición m´ax(λk ) ≤ 1 equivale a ! " 1 2 kπ 0 ≤ α sin ≤ k = 1, . . . , N − 1 2N 2 Como esto se debe de cumplir para cualquier valor de k y N , entonces tenemos la condición 0≤α ≤

1 2

o lo que es lo mismo

κτ 1 ≤ 2 h 2 Vemos que los pasos de integración espaciales y temporales deben de cumplir la relación τ 1 ≤ 2 h 2κ para que el método sea estable. Si κ = 12 y h = 0,1, se debe cumplir τ ≤ 0,01 para que el método sea estable. Estta condición puede requerir un volumen de cálculo muy considerable para la integración temporal. Vamos a ver que hay otros métodos de integración con condiciones de estabilidad menos restrictivas. Se pude demostrar que la condición ded estabilidad implica que en un paso temporal de integración sólo pasa una cantidad apreciable de calor a las celdillas vecinas inmediatas. Si κ es muy grande, el calor se difunde deprisa y el paso temporal debe de ser pequeño.

10.2.2.

Método de diferencias finitas regresivas

Un método alternativo es utilizar diferencias finitas regresivas, de forma que la derivada temporal utiliza un instante anterior en el tiempo. La ecuación en diferencias toma la forma n T n − 2Tin + Ti−1 Tin − Tin−1 = κ i+1 τ h2

208

CAPÍTULO 10. ECUACIONES DIFERENCIALES EN DERIVADAS PARCIALES

que, denominando α =

κτ queda como h2 n n −αTi+1 + (1 + 2α)Tin − αTi−1 = Tin−1

Esta ecuación proporciona una relación de recurrencia que liga los valores de la temperatura en un instante con los valores en el instante anterior. Podemos escribir esta relación en forma matricial como BTn = Tn−1 con B una matriz tridiagonal



   B=  

1 0 ··· 0 −α 1 + 2α −α ··· .. .. .. .. . . . . ··· −α 1 + 2α −α ··· 0 0 1



    = I − αD  

Vemos que la matriz B no es la inversa de la matriz A por lo que el método de diferencias regresivas es esencialmente distinto del de diferencias progresivas: AB = I − α 2 D2

En el caso cde diferencias regresivas, tenemos que resolver un sistema de ecuaciones para avanzar un paso en el tiempo. Este sistema se puede resolver por factorización LU , y, si N es muy grande, por métodos iterativos como el método SOR o el método del gradiente conjugado. El método de diferencias regresivas tiene un buen comportamiento en lo que concierne a la estabilidad. Se puede demostrar analíticamente que los valores propios de B son ! " 2 kπ λk = 1 + 4α sin k = 1, . . . , N − 1 2N

Como α > 0, λk > 1, con lo que la matriz inversa existe, y además, como los valores propios de la inversa son los inversos de los valores propios ρ(B−1 ) < 1

por lo que el método siempre es estable. Sin embargo, estabilidad no es sinónimo de precisión numérica, y el hecho de que la derivada temporal no esté centrada en el tiempo, hace que el método tenga una precisión del orden de O(τ + h2 ), al igual que el método de diferencias progresivas. Sería deseable un método que convergiese más rápidamente en el tiempo, al menos como τ 2 .

10.2.3.

Algoritmo de Crank-Nicolson

La dependencia lineal del error con τ en los dos métodos anteriores proviene de la discretización de la derivada temporal mediante diferencias finitas no centradas en el tiempo. Tenemos ∂T ∂t ∂T ∂t

T (t + τ, x) − T (t, x) 1 ∂ 2 T + τ 2 + O(τ 2 ) τ 2 ∂t T (t, x) − T (t − τ, x) 1 ∂ 2 T = − τ 2 + O(τ 2 ) τ 2 ∂t =

10.2. LA ECUACIÓN DE PROPAGACIÓN DEL CALOR (ECUACIÓN DE DIFUSIÓN) 209 Si promediamos ambas expresiones, tendremos un error en O(τ 2 ) para la derivada primera temporal. Esta es la idea del algoritmo de Crank-Nicolson, que equivale a aproximar la derivada temporal mediante diferencias centradas en el tiempo. Si promediamos ambas expresiones n − 2T n + T n Ti+1 Tin+1 − Tin i i−1 = κ 2 τ h n+1 n+1 n+1 n T − 2Tin+1 + Ti−1 − Ti Ti i+1 = κ τ h2

obtenemos

Tin+1 − Tin κ = τ 2

,

n+1 n+1 n − 2T n + T n Ti+1 − 2Tin+1 + Ti−1 Ti+1 i i−1 + h2 h2

-

Esta ecuación se puede redisponer como , ! n n+1 n+1 n " − 2Tin+1 + Ti−1 − 2Tin + Ti−1 Tin+1 κ Ti+1 Tin κ Ti+1 − = + τ 2 h2 τ 2 h2 que se puede escribir en la forma matricial B1 Tn+1 = A1 Tn con

con



 1 0 0 ··· α  α  1−α ···    2  2 α  .. .. ..  A1 =  . . . . . . =I+ D  2 α  α  ···  1 − α   2 2 ··· 0 0 1   1 0 ··· 0 α  α  1+α − ···   −  2  2 α  ..  . . . . . . B1 =  . . . . =I− D   2  ··· −α 1 + α −α    2 2 ··· 0 0 1   0 0 ··· 0  1 −2 1 · · ·     .. .. .. ..  D= . . . .     · · · 1 −2 1  ··· 0 0 0

κτ si las condiciones de contorno son constantes. En A1 y B1 hemos definido α = 2 . La matriz h B1 es definida positiva con diagonal dominante, por lo que es invertible. Como es tridiagonal,

210

CAPÍTULO 10. ECUACIONES DIFERENCIALES EN DERIVADAS PARCIALES

el algoritmo de Crout es relativamente sencillo de aplicar, aunque para valores grandes de N es preferible utilizar los métodos de superrelajación. El método siempre es estable y converge como O(τ 2 + h2 ). Podemos escribir el método como T

n+1

.

α /−1 . α / n = I− D I+ D T 2 2

Desde el punto de vista numérico conviene evitar la multiplicación de grandes matrices. En la ecuación anterior, si definimos la matriz P por B−1 1 . α /−1 P= I− D 2

tenemos

A1 = I +

α α D = 2I − I + D = 2I − P−1 2 2

con lo que obtenemos T

n+1

0 . 1 α /−1 = P(2I − P )T = (2P − I)T = 2 I − D − I Tn 2 −1

n

n

que implica sólo la multiplicación de una matriz por un vector.

10.2.4.

Criterio de estabilidad de von Neumann

Existen diversos criterios para determinar la estabilidad de una ecuación diferencial en derivadas parciales. Entre ellos, el criterio de von Neumann es el más fácil de aplicar para determinar si un método de integración es estable. Consiste en suponer una solución compleja, de la forma T (x,t) = A(t)eikx con lo que discretizando queda

T jn = An eik jh

Se define el coeficiente de amplificación ξ como An+1 ξ= n A Si |ξ | > 1 la solución crecerá con el tiempo y el método es inestable. En caso contrario, la solución disminuye con el tiempo y el sistema es estable. Vamos a aplicar este método a los tres esquemas de integración anteriores. Para el esquema FCTS % κτ $ n n Ti+1 − 2Tin + Ti−1 2 h 3 κτ 2 An+1 eik jh = An eik jh + 2 An eik( j+1)h − 2An eik jh + An eik( j−1)h h Tin+1 = Tin +

tenemos

10.2. LA ECUACIÓN DE PROPAGACIÓN DEL CALOR (ECUACIÓN DE DIFUSIÓN) 211 y simplificando obtenemos ! " 3 κτ 2 ikh 2κτ 4κτ 2 kh −ikh ξ = 1+ 2 e +e − 2 = 1 + 2 [cos(kh) − 1] = 1 − 2 sin h h h 2 por lo que la condición |ξ | ≤ 1 sólo se cumple si ! " 4κτ 2 kh ≤2 sin h2 2 que a su vez se cumple siempre si

2κτ ≤1 h2 Esta condición, en el caso de la ecuación de difusión, implica que la longitud de difusión en un intervalo temporal τ, √ d = 2Dτ sea inferior a la longitud de una celdilla espacial h. Si aplicamos el criterio de von Neumann al método de diferencias finitas regresivas 5 4 n+1 n+1 Tin+1 = Tin + α Ti+1 − 2Tin+1 + Ti−1 obtenemos

3 κτ 2 n+1 ik( j+1)h n+1 ik jh n+1 ik( j−1)h A e − 2A e + A e h2 que, procediendo de forma análoga, obtenemos An+1 eik jh = An eik jh +

ξ=

1

! " 4κτ 2 kh 1 + 2 sin h 2

Por lo tanto, siempre se cumple |ξ | ≤ 1, y el método es incondicionalmente estable para cualquier paso de integración temporal τ. Análogamente, en el algoritmo de Crank-Nicolson , n+1 n+1 n − 2T n + T n − 2Tin+1 + Ti−1 Ti+1 Tin+1 Tin κ Ti+1 i i−1 = + + τ τ 2 h2 h2 tenemos An+1 eik jh = An eik2jh + 3 κτ n+1 ik( j+1)h n+1 ik jh n+1 ik( j−1)h n ik( j+1)h n ik jh n ik( j−1)h A e − 2A e + A e +A e − 2A e + A e 2h2 con lo que obtenemos ! " 2κτ 2 kh 1 − 2 sin h 2 ! " ξ= 2κτ kh 1 + 2 sin2 h 2 y el método es también incondicionalmente estable.

212

10.2.5.

CAPÍTULO 10. ECUACIONES DIFERENCIALES EN DERIVADAS PARCIALES

Condiciones de contorno

Las condiciones de contorno se clasifican en condiciones de Dirichlet, cuando se imponen sobre las funciones (la temperatura en los extremos de la barra en el caso de la ecuación del calor); T (0,t) = Ta T (L,t) = Tb y de Neumann cuando se imponen sobre las derivadas de las funciones # # ∂ T (x,t) ## ∂ T (x,t) ## = ja = jb ∂ x #x=0 ∂ x #x=L

que equivale a imponer condiciones sobre los flujos de calor. Otro tipo de condiciones, útiles para simular sistemas, son las condiciones de contorno periódicas, que consisten en igualar los valores de las funciones en los extremos T (0,t) = T (L,t) # # ∂ T (x,t) ## ∂ T (x,t) ## = ∂ x #x=0 ∂ x #x=L

Junto con las condiciones de contorno, válidas para cualquier valor del tiempo, hay que proporcionar condiciones iniciales, que proporcionan la distribución espacial del sistema en un instante inicial. Esto equivale a dar el perfil de temperaturas en el instante inicial T (x, 0) = f (x)

10.2.6.

Densidad de neutrones en un reactor nuclear

Un caso interesante es el de la determinación de la densidad de neutrones en un reactor nuclear. La densidad en un punto, n(r,t), es debida a la difusión de neutrones provinientes de otros puntos y a la creación de neutrones mediante fisión en el punto. La producción neta de neutrones resulta del balance de los neutrones producidos por fisión y los absorbidos por captura. Si la tasa neta de creación de neutrones por cada neutrón es C(r,t), tenemos ∂ n(r,t) = D$2 n(r,t) +C(r,t)n(r,t) ∂t Vamos a considerar el caso monodimensional ∂ n(x,t) ∂ 2 n(x,t) =D +C(x,t)n(x,t) ∂t ∂ x2 en una caja de longitud L, −L/2 ≤ L ≤ L/2, con condiciones de contorno de Dirichlet n(L/2,t) = n(−L/2,t) = 0 En el caso de 235 U, D " 105 m2 /s y C " 108 s−1 . Cuando la longitud L es mayor que una longitud crítica Lc , entonces se alcanza la masa crítica y se produce la reacción en cadena, con

10.2. LA ECUACIÓN DE PROPAGACIÓN DEL CALOR (ECUACIÓN DE DIFUSIÓN) 213 aumento exponencial de la densidad de neutrones. Se puede demostrar, resolviendo la ecuación por separación de variables, que 6 Lc = π D/C

En los ejemplos numéricos, tomaremos D = C = 1 con lo que Lc = π. El esquema de integración FCTS proporciona la ecuación en diferencias j

j+1

ni que resulta en j+1

ni

j

j

j

j n − 2ni + ni−1 − ni j = D i+1 +Cni 2 τ h j

= ni + j

/ Dτ . j j j j − 2n + n n i i−1 +Cτni i+1 h2

para i = 1, . . . , N − 1, con n0 = nN = 0.

10.2.7.

La ecuación de Schrödinger dependiente del tiempo

Un caso interesante de ecuación parabólica, es la ecuación de Schrödinger dependiente del tiempo: ∂ ψ(x,t) h¯ 2 ∂ 2 ψ(x,t) i¯h +V (x)ψ(x,t) =− ∂t 2m ∂ x2 Definiendo el operador Hamiltoniano como H =−

h¯ 2 ∂ 2 +V (x) 2m ∂ x2

podemos escribir una solución formal de la ecuación como 7 8 i ψ(x,t) = exp − Ht ψ(x, 0) h¯

(10.2)

Para resolver numéricamente esta ecuación discretizaremos el espacio y el tiempo, con pasos h y τ, respectivamente. En el método FTCS tenemos n − 2ψ n + ψ n ψin+1 − ψin h¯ 2 ψi+1 i i−1 i¯h =− +Vi ψin 2 τ 2m h

que podemos escribir como i¯h con la definición

h¯ 2 δi+1, j − 2δi, j + δi−1, j +Vi δi, j 2m h2 representa la delta de Kronecker. Podemos escribir la ecuación anterior como ! " iτ n+1 % = I − H %n (10.3) h¯ Hi j = −

donde δi, j

N ψin+1 − ψin = " Hi j ψ nj τ i=0

214

CAPÍTULO 10. ECUACIONES DIFERENCIALES EN DERIVADAS PARCIALES

donde



  % =  n

ψ1n ψ2n .. .

ψNn

    

La solución 10.3 se puede interpretar como el primer orden del desarrollo en serie de Taylor de la solución de la ecuación 10.2. Si utilizamos el método de diferencias regresivas tenemos i¯h

N ψin+1 − ψin = " Hi j ψ n+1 j τ i=0

con lo que recolectando términos obtenemos la ecuación ! " iτ I + H %n+1 = %n h¯

(10.4)

que nos permite avanzar un paso en el tiempo resolviendo un sistema de ecuaciones. Tomando el promedio de 10.3 y 10.4 obtenemos el algoritmo de Crank-Nicolson " ! " ! iτ iτ n+1 I+ H % = I − H %n 2¯h 2¯h que se puede reescribir como %

n+1

!

iτ = I+ H 2¯h

"−1 !

" iτ I − H %n 2¯h

(10.5)

que corresponde a aproximar la exponencial en la ecuación10.2 como ex =

1−x 1+x

Esta expresión no es otra cosa que el aproximante de Padé más bajo de la exponencial, que tiene un radio de convergencia más elevado que la serie de Taylor. Además el operador ! es unitario, al igual que el operador

iτ I+ H 2¯h

"−1 !

iτ I− H 2¯h

"

7 8 i exp − Ht h¯

Por lo tanto, el algoritmo de Crank-Nicolson no sólo está centrado en el espacio y el tiempo simultáneamente, sino que además conserva unitariedad. El algoritmo 10.5 nos obliga a calcular

10.2. LA ECUACIÓN DE PROPAGACIÓN DEL CALOR (ECUACIÓN DE DIFUSIÓN) 215 inversas de matrices, lo que puede resultar un serio inconveniente para grandes dimensiones. Una forma alternativa de plantear los cálculos es reordenar 10.5de la forma %

n+1

!

"−1 0 ! "1 iτ iτ 2I − I + H %n = I+ H 2¯h 2¯h 9 ! "−1 : ! "−1 iτ iτ n = 2 I+ H −I % = 2 I + H %n − %n 2¯h 2¯h

El primer término del último miempro se calcula como la solución !n del sistema de ecuaciones ! " iτ 1 I + H !n = %n 2 2¯h con lo que obtenemos para cada iteración %n+1 = !n − %n Por lo tanto, cada paso de integración implica sólamente la solución de un sistema de ecuaciones lineales complejo.

10.2.8.

Algunas ecuaciones no lineales en derivadas parciales con derivadas temporales primeras

En Física aparecen otras ecuaciones que implican derivadas primeras con respecto del tiempo, pero con términos no lineales. La forma de proceder para construir esquemas en diferencias finitas es similar, salvo que en los términos no lineales se utiliza las diferencias progresivas. Por ejemplo, en acústica no lineal aparece la ecuación de Burger ∂φ ∂φ ∂ 2φ = −φ −κ 2 ∂t ∂x ∂x Un método adecuado es el promedio de diferencias progresivas n n φ n − φi−1 φ n − 2φin + φi−1 φin+1 − φin = −φin i+1 − κ i−1 τ 2h h2

y regresivas n+1 n φ n+1 − 2φin+1 + φi−1 φ n − φi−1 φin+1 − φin = −φin i+1 − κ i−1 τ 2h h2 que da como resultado , n+1 n+1 n+1 n −φn n − 2φ n + φ n φ − 2φ + φ φ φ φin+1 − φin κ i i i−1 i−1 i−1 i−1 = −φin i+1 − + i−1 τ 2h 2 h2 h2

216

CAPÍTULO 10. ECUACIONES DIFERENCIALES EN DERIVADAS PARCIALES

Notemos que el término no lineal se ha evaluado de forma regresiva tanto en el método progresivo como en el regresivo, pero con derivadas espaciales centradas. Esta relación se puede poner como . α / n+1 . α / n I+ D ! = I − D ! − δ &(!n )C!n 2 2 con

 0 1 0 ··· 0 0  −1 0 1 ··· 0 0      0 −1 0 1 · · · 0     C =  ... . . . . . . . . . . . . 0      0 · · · −1 0 1 0    0 · · · 0 −1 0 1  0 ··· 0 0 −1 0   φ0n 0 0 · · · 0  0 φn 0 ··· 0    1 n   &(!n ) =  0 0 φ2 . . . 0   .. . . . . . .   . . . . 0  0 · · · 0 0 φNn   φ0n  φn   1  !n =  ..   .  φNn 

α=

κτ τ δ= 2 h 2h

Tenemos finalmente !

n+1

. . α /−1 . α / n α /−1 nT = I+ D I − D ! −δ I + D ! C!n 21 2 0 . 2 . /−1 α /−1 α = 2 I+ D − I !n − δ I + D &(!n )C!n 2 2

Resolviendo y

tenemos la relación.

.

. α / I + D % n = !n 2 I+

α / n D χ = &(!n )C!n 2

!n+1 = 2%n − !n − δ χ n Vemos que tenemos que resolver dos sistemas lineales por paso de integración.

10.2. LA ECUACIÓN DE PROPAGACIÓN DEL CALOR (ECUACIÓN DE DIFUSIÓN) 217 Ejercicio: Resolver la ecuación de Burger con las condiciones iniciales φ (x, 0) = −1 si x > 0 y φ (x, 0) = 1 si x < 0 y las condiciones de contorno φ (±L/2,t) = ∓1. Tomar L = 10 y κ = 1. Comparar con el resultado exacto para L −→ # φ (x,t) = donde

F(x,t) − F(−x,t) F(x,t) + F(−x,t)

! 7 "8 1 t−x x − 2t √ F(x,t) = e 1 − erf 2 2 t

con erf (x)la función de error. Otra ecuación no lineal de interés es la ecuación de Korteweg-DeVries(KdV), que describe el comportamiento de solitones: ∂ρ ∂ ρ ∂ 3ρ = −6ρ − ∂t ∂ x ∂ x3 Al igual que en caso anterior tomamos el promedio de diferencias finitas progresivas y regresivas, calculando el término no lineal de forma regresiva. Las diferencias espaciales se toman centradas. Para la derivada tercera tenemos d 3 f (x2 ) f4 − 2 f3 + 2 f1 − f0 " 3 dx 2h3 con lo que obtenemos para las diferencias progresivas n n + 2ρ n − ρ n ρ n − ρi−1 ρ n − 2ρi+1 ρin+1 − ρin i−1 i−2 = −6ρin i+1 − i+2 τ 2h 2h3

y regresivas n+1 n+1 n+1 n ρ n+1 − 2ρi+1 + 2ρi−1 − ρi−2 ρ n − ρi−1 ρin+1 − ρin = −6ρin i+1 − i+2 τ 2h 2h3

Tomando el promedio, tenemos 9 : n+1 n+1 n+1 n+1 n n n − 2ρ n + 2ρ n − ρ n ρi+2 − 2ρi+1 + 2ρi−1 − ρi−2 ρin+1 − ρin 1 ρi+2 i+1 i−1 i−2 n ρi+1 − ρi−1 = −6ρi − + τ 2h 2 2h3 2h3 Tenemos el sistema

(I + β F1 ) ρ n+1 = (I − β F1 ) ρ n − 6γ&(ρ n )Gρ n

que se puede reorganizar como ρ n+1 = (I + β F1 )−1 (I − β F1 ) ρ n − 6γ (I + β F1 )−1 &(ρ n )Gρ n

218

CAPÍTULO 10. ECUACIONES DIFERENCIALES EN DERIVADAS PARCIALES

Al igual que en caso anterior, poniendo I − β F1 = 2I − (I + β F1 ) tenemos 2 3 −1 n+1 ρ = 2 (I + β F1 ) − I ρ n − 6γ (I + β F1 )−1 &(ρ n )Gρ n con lo que tenemos la solución

ρ n+1 = 2!n − ρ n − 6γ χ n donde hemos definido

(I + β F1 ) !n = ρ n (I + β F1 ) χ n = &(ρ n )Gρ n

Vemos que cada paso de integración implica la resolución de dos sistemas de ecuaciones. Ejercicio:

" L Resolver la ecuación de KdV con las condiciones iniciales ρ x = ± = 0 y condiciones 2 de contorno periódicas. Comparar con la solución analítica ρ(x,t) = cosh−2 (x − 4t).

10.3.

!

La ecuación de ondas

En el caso de la ecuación de ondas 2 ∂ 2 u(x,t) 2 ∂ u(x,t) = c ∂t 2 ∂ x2

tenemos derivadas segundas tanto en el espacio como en el tiempo, por lo que la primera elección es tomar derivadas centradas tanto en el espacio como en el tiempo. Tomando un paso de integraj ción temporal de τ y un paso de integración espacial de h, y denominando u(x0 +ih,t0 + jτ) = ui , tenemos el esquema de integración j+1

ui

j

j−1

− 2ui + ui τ2

j j j 2 ui+1 − 2ui + ui−1 =c h2

2 2 Definiendo λ 2 = τ c2 , este esquema de integración se escribe como, h . / j+1 j j−1 j j j j j j j−1 2 ui = 2ui − ui − λ ui+1 − 2ui + ui−1 = 2(1 − λ 2 )ui + λ 2 (ui+1 + ui−1 ) − ui

(10.6)

Tomemos el dominio de integración 0 < x < l y t > 0 y consideremos condiciones de contorno constantes: u(0,t) = u0 , u(l,t) = ul para t > 0. Aparte de estas condiciones de contorno, hay que proporcionar dos condiciones iniciales, ya que la ecuación es de segundo orden en el tiempo. Tomemos estas condiciones como u(x, 0) = f (x)

10.3. LA ECUACIÓN DE ONDAS

219

y

∂ u(x, 0) = g(x) ∂t para 0 < x < l. En el caso de ondas en una cuerda, estas condiciones equivalen a proporcionar la posición y la velocidad de cada punto en el instante inicial. La relación dada por la Ec. 10.6 se puede escribir en forma matricial como u j+1 = Au j − u j−1 con

y



    A=   

2 0 0 0 ··· 0 2 2 2 λ 2(1 − λ ) λ 0 ··· 0 2 2 2 0 λ 2(1 − λ ) λ 0 0 .. ... ... ... . 0 0 .. .. ... . . λ 2 2(1 − λ 2 ) λ 2 0 0 ··· 0 0 2 

j

u0 j u1 .. .

   j u =   j  uN−1 j uN



        

      

Notemos que las filas primera y última de A lo único que hacen es mantener las condiciones j j de contorno de u0 y uN . En el caso de que estas condiciones de contorno sean nulas (u(0,t) = u(l,t) = 0 para t > 0), estas dos filas se pueden suprimir junto con las columnas primeras y última de A. El primer vector es   f (x0 )  f (x1 )      .. u0 =   .    f (xN−1)  f (xN ) Las componentes de u1 las podríamos tomar como u1i = f (xi ) + τg(xi ) pero esto sería una aproximación de sólo primer orden en el tiempo, mientras que las derivadas centradas son de segundo orden, lo que produciría un deterioro de la precisión de la solución en comparación con la naturaleza de segundo orden del método. Un método de proporcionar condiciones iniciales de segundo orden en τ es considerar el desarrollo en serie ∂ u(xi , 0) 1 ∂ 2 u(xi , 0) 2 τ+ τ + O(τ 3 ) ∂t 2 ∂t 2 De las condiciones iniciales tenemos ∂ u(xi , 0) = g(xi ) ∂t u1i = u0i +

220

CAPÍTULO 10. ECUACIONES DIFERENCIALES EN DERIVADAS PARCIALES

mientras que de la ecuación diferencial tenemos 2 ∂ 2 u(xi , 0) 2 ∂ u(xi , 0) = c = c2 f '' (xi ) 2 2 ∂t ∂x

Por lo tanto, obtenemos para u1i la siguiente expresión en segundo orden en τ 1 u1i = u0i + g(xi )τ + c2 τ 2 f '' (xi ) 2 En caso de que f (xi ) no esté disponible en forma continua sino sólo en la malla de puntos, o simplemente que la derivada analítica de f (x) sea compleja de calcular, podemos reemplazar f '' (xi ) por su derivada centrada, f '' (xi ) =

f (xi+1 ) − 2 f (xi ) + f (xi−1 ) h2

con lo que obtenemos λ2 [ f (xi+1 ) − 2 f (xi ) + f (xi−1 )] 2 λ2 = (1 − λ 2 ) f (xi ) + [ f (xi+1 ) + f (xi−1 )] + g(xi )τ 2

u1i = f (xi ) + g(xi )τ +

De esta forma, el error del algoritmo es O(τ 2 + h2 ).

10.3.1.

Criterio de estabilidad

Introduciendo en la ecuación 10.6 una solución unj = An eik jh tenemos ξ 2 − 2ξ (1 − λ 2 + λ 2 cos kh) + 1 = 0 que tiene como soluciones ; 2 ξ = 1 + λ (cos kh − 1) ± (1 + λ 2 (cos kh − 1)) − 1 2

El término fuera de la raíz está comprendido entre 1 y 1 −2λ 2 . Nunca puede hacerse mayor que 1 en valor absoluto si λ < 1, pero puede hacerse negativo y mayor que la unidad en valor absoluto si λ > 1. La raíz es siempre un número imaginario si λ < 1, ya que el primer término es siempre inferior a la unidad. Sin embargo, puede hacerse mayor que la unidad en valor absoluto si λ > 1. Por lo tanto, se cumple que |ξ | < 1 si λ < 1 o lo que es lo mismo cτ < h, que nos dice que para que la solución sea estable, la distancia que se propaga la onda durante un paso de integración temporal debe de ser menor que el paso de integración espacial.

10.4. LAS ECUACIONES DE POISSON Y DE LAPLACE

10.4.

221

Las ecuaciones de Poisson y de Laplace

En Electrostática, la ecuación de Poisson en tres dimensiones liga el potencial electrostático con la densidad de carga: $2 !(r) =

ρ(x, y, z) ∂ 2 !(x, y, z) ∂ 2 !(x, y, z) ∂ 2 !(x, y, z) + + = − ∂ x2 ∂ y2 ∂ z2 ε0

Cuando la densidad de carga es nula, esta ecuación se convierte en la ecuación de Laplace: $2 !(r) = 0 Vamos a considerar el problema en dos dimensiones ∂ 2 !(x, y) ∂ 2 !(x, y, z) + = f (x, y) ∂ x2 ∂ y2 definido en un dominio rectangular a ≤ x ≤ b, c ≤ y ≤ d. Tomamos condiciones de contorno de Dirichlet sobre el contorno !(x, y) = g(x, y) Para resolver el sistema discretizamos cada una de las coordenadas espaciales. Vamos a suponer en pricipio pasos de integración diferentes en cada una de las coordenadas. Denominemos h el paso de integración en la coordenada x y k al paso de integración en la coordenada y. Tendremos M subintervalos en la coordenada x (M + 1 puntos) h=

b−a M+1

y N subintervalos en la coordenada y (N + 1 puntos) k=

d −c N +1

Tomamos derivadas centradas en el espacio en ambas coordenadas, con lo que tendremos un método con error O(h2 + k2 ). La ecuación en diferencias queda como !(xi+1 , y j ) − 2!(xi , y j ) + !(xi−1 , y j ) !(xi , y j+1 ) − 2!(xi , y j ) + !(xi , y j−1 ) + = f (xi , y j ) h2 k2 (10.7) Las condiciones de contorno se escriben como !(x0 , y j ) !(xM , y j ) !(xi , y0 ) !(xi , yN )

= = = =

g(x0 , y j ) g(xM , y j ) g(xI , y0 ) g(xI , yN )

222

CAPÍTULO 10. ECUACIONES DIFERENCIALES EN DERIVADAS PARCIALES

para i = 0, 1, . . . , M, j = 0, 1, . . . , N.Podemos reescribir la ec. 10.7 como : 9! " ! "2 2 $ % % h h $ + 1 !(xi , y j )− !(xi+1 , y j ) + !(xi−1 , y j ) − !(xi , y j+1 ) + !(xi , y j−1 ) = h2 f (xi , y j ) 2 k k (10.8) Este es un sistema lineal con (N − 1) × (M − 1) incógnitas. Es conveniente reordenar los pares (i, j) de forma que la matriz de orden (N − 1)2 × (M − 1)2 sea una matriz de 5 bandas. Esto se consigue redefiniendo los índices (i, j) → s = i + (M − 1 − j)(N − 1) para i entre 1 y M − 1 y j entre 1 y N − 1. Alternativamente s = i + ( j − 1)(N − 1) (y ambas expresiones rremplazando i por j). La cuestión es que el índice varíe más rápidamente a lo largo del eje y que del eje x (o del eje x que del eje y) y que los índices sean consecutivos. Vamos a tomar en adelante la primera convención. De esta forma, ecuación ?? queda como 9! " : ! "2 % h 2 h $ 2 + 1 !s − (!s+1 + !s−1 ) − !s−(N−1) + !s+(N−1) = h2 fs k k con lo que tenemos una diagonal principal con todos los elementos iguales a 9! " : h 2 2 +1 k

con las dos subdiagonales superior e inferior con todos los elementos iguales a−1 y dos bandas separadas N − 1 elementos arriba y abajo también con todos sus elementos iguales a -1. Ejemplo: Consideremos un cuadrado dieléctrico de 4 cm de lado donde dos bordes están a tierra ( 0 V) y el vértice opuesto está a 100 V. Calcular La distribución de potencial, suponiendo que la densidad de carga es nula. Tenemos la ecuación de Laplace y como condiciones de contorno !(0, x) = !(y, 0) = 0. !(4, x) = 25x, !(y, 4) = 25y, suponiendo que satisface la ley de Ohm. Tomamos h = k = 1~cm, con lo que M = N = 4 y la dimensión del sistema es 9. Tenemos que calcular el potencial en 9 puntos (xi , y j ). Con la convención anterior tenemos (x1 , y1 ) = P7 , (x2 , y1 ) = P8 ,(x3 , y1 ) = P9 ,(x1 , y2 ) = P4 , (x2 , y2 ) = P5 ,(x3 , y3 ) = P6 ,(x1 , y3 ) = P1 , (x2 , y1 ) = P2 ,(x3 , y1 ) = P3 , y las ecuaciones en cada uno de estos puntos son P1 : 4!1 − !2 − !4 = P2 : 4!2 − !3 − !1 = P3 : 4!3 − !2 − !6 = P4 : 4!4 − !5 − !1 − !7 = P5 : 4!5 − !6 − !4 − !2 − !8 = P6 : 4!6 − !5 − !3 − !9 = P7 : 4!7 − !8 − !4 = P8 : 4!8 − !59 − !7 − !5 = P9 : 4!9 − !8 − !6 =

!0,3 + !1,4 !2,4 !4,3 + !3,4 !0,2 0 !4,2 !0,1 + !1,0 !2,0 !3,0 + !4,1

10.4. LAS ECUACIONES DE POISSON Y DE LAPLACE

223

Las condiciones de contorno dan !1,0 = !2,0 = !3,0 = !0,1 = !0,2 = !0,3 = 0, !2,4 = !4,2 = 50,!3,4 = !4,3 = 75. El sistema lineal asociado es      25 !1 4 −1 0 −1 0 0 0 0 0      −1 4 −1 0 −1 0 0 0 0    !2   50        0 −1 4 −1 0 −1 0 0 0   !3   150      !4   0   −1 0 −1 4 −1 0 −1 0 0       0 −1 0 −1 4 −1 0 −1 0   !5  =  0         !6   50   0 0 −1 0 −1 4 −1 0 −1           0 0 0 −1 0 −1 4 −1 0    !7   0        0 0  !8 0 0 0 −1 0 −1 4 −1 25 !9 0 0 0 0 0 −1 0 −1 4

Este sistema se puede resolver por el método LU ya que la matriz está bien condicionada. Cuando el sistema es muy grande es mejor utilizar el método iterativo de Gauss-Seidel. En este método ponemos 9 : ! "2 $ % % 1 h $ 3 !(xi+1 , y j ) + !(xi−1 , y j ) + !(xi , y j ) = 2$ %2 !(xi , y j+1 ) + !(xi , y j−1 ) + h2 f (xi , y j ) h k 2 +1 k

y ponemos todos los potenciales a cero salvo las condiciones de contorno. Calculamos cada !(xi , y j ) con las condiciones iniciales, que llamamos !(0) (xi , y j ) comenzando por las capas exteriores. En cada iteración recalculamos !(k) (xi , y j ) con los #valores obtenidos en la iteración # # (k+1) # (k−1) (k) anterior ! (xi , y j ) hasta que en dos iteraciones sucesivas#! (xi , y j ) − ! (xi , y j )# < ε, donde ε es una tolerancia prestablecida .

224

CAPÍTULO 10. ECUACIONES DIFERENCIALES EN DERIVADAS PARCIALES