Práctica sobre el Método de Diferencias Finitas para E.D.P.'s

Realizada por: 		Pedro González Rodelas
Primera versión:	01-06-2005
Última modificación: 	30-05-2007

Ecuación de Poisson

Se trata de la resolución aproximada, mediante el método de diferencias finitas del siguiente problema de contorno de tipo Dirichlet para la ecuación de Poisson en un dominio rectángular:

		Ω≡{(x,y) /  a<x<b,  c<y<d}

	\(TraditionalForm
      \`Δu(x, y) ≡ 
          \(u\_xx\)(x, y) + \(u\_yy\)(x, y)\  = \ f(x, y)\),   ∀(x,y)∈Ω
	(Subscript[u, |∂Ω])≡φ, 
	(con \(TraditionalForm
      \`φ\_\(|{y = c}\) ≡ φ\_1, \ 
      φ\_\(|{x = b}\) ≡ φ\_2, \ 
      φ\_\(|{y = d}\) ≡ φ\_3, \ 
      φ\_\(|{x = a}\) ≡ φ\_4\))

empleando diferencias centradas para cada una de las derivadas parciales de segundo orden:

	\(TraditionalForm
      \`\(u\_xx\)(x\_i, y\_j) ≃ 
        \(u\_\(i - 1, j\) - 2  u\_\(i, j\) + 
            u\_\(i + 1, j\)\)\/Δx\^2\),
	
	\(TraditionalForm
      \`\(u\_yy\)(x\_i, y\_j) ≃ 
        \(u\_\(i, j - 1\) - 2  u\_\(i, j\) + 
            u\_\(i, j + 1\)\)\/Δy\^2\)

Resolveremos de hecho el siguiente ejemplo concreto:

	Δu(x, y)  =  x + y,    ∀(x,y)∈Ω≡{(x,y) /  0<x<2,  0<y<3}
	(Subscript[u, |∂Ω])≡φ,  
 (con \(TraditionalForm
      \`φ\_\(|{y = 0}\) ≡ sen(\(π\/2\) x), \ 
      φ\_\(|{x = 2}\) ≡ 0, \ 
      φ\_\(|{y = 3}\) ≡ sen(\(π\/2\) x), , \ 
      φ\_\(|{x = 0}\) ≡ 0\))

Definición de los datos del problema

Clear["Global`*"]

a := 0b := 2c := 0\(d := 3
\)\(f[x_, y_] := x + y
\)(Subscript[φ, 1])[x_] := Sin[(π/2) x](Subscript[φ, 2])[y_] := 0(Subscript[φ, 3])[x_] := Sin[(π/2) x](Subscript[φ, 4])[y_] := 0

Definición de la discretización numérica

Clear[u, x, y, Δx, Δy, i, imax, j, jmax]
      imax := 10\(jmax := 10
\)Δx := (b - a)/imax // NΔy := (d - c)/jmax // N\(


    Do[x\_i = a + i\ Δx, {i, 0, imax}]\)Do[(Subscript[y, j]) = c + j Δy, {j, 0, jmax}]

Table[(Subscript[x, i]), {i, 0, imax}]

Table[(Subscript[y, j]), {j, 0, jmax}]

Aplicación de la condición de contorno

Do[u[i, 0] = (Subscript[φ, 1])[(Subscript[x, i])], {i, 0, imax}]Do[u[imax, j] = (Subscript[φ, 2])[(Subscript[y, j])], {j, 0, jmax}]Do[u[i, jmax] = (Subscript[φ, 3])[(Subscript[x, i])], {i, 0, imax}]Do[u[0, j] = (Subscript[φ, 4])[(Subscript[y, j])], {j, 0, jmax}]

Table[u[i, j], {i, 0, imax}, {j, 0, jmax}] // TableForm

Planteamiento y resolución del sistema de ecuaciones

\(ecs = {}; 

    Do[
AppendTo[ecs, 
        \(u[i - 1, j] - 2  u[i, j] + u[i + 1, j]\)\/Δx\^2 + 
            \(u[i, j - 1] - 2  u[i, j] + u[i, j + 1]\)\/Δy\^2 
          == f[x\_i, y\_j]]
, {i, 1, imax - 1}, {j, 
        1, jmax - 1}]\)

ecs

Table[u[i, j], {i, imax - 1}, {j, jmax - 1}] // Flatten

?Solve

sol = 
      Solve[ecs, Table[u[i, j], {i, imax - 1}, {j, jmax - 1}] // Flatten]

Table[u[i, j], {i, 0, imax}, {j, 0, jmax}] /. sol[[1]]

?ListPlot3D

ListPlot3D[
      Table[u[i, j], {i, 0, imax}, {j, 0, jmax}] /. sol[[1]], 
      Ticks -> None];

ListDensityPlot[
      Table[u[i, j], {i, 0, imax}, {j, 0, jmax}] /. sol[[1]], 
      Ticks -> None];

Resolución de un problema mixto para la ecuación del calor (o de difusión lineal)

Se trata de la resolución aproximada, mediante el método de diferencias finitas, del siguiente problema mixto para la ecuación de ondas lineal:

	(Subscript[u, t]) = (α^2)(Subscript[u, xx]),   a<x<b,  t>0
	u(a,t) = (Subscript[h, 1])(t),   u(b,t) = (Subscript[h, 2])(t),     t>0
	u(x,0) = f(x),  a<x<b	
	
empleando como antes una diferencia dividida centrada para la derivada segunda  y otra "adelantada" para la aproximación de la derivada temporal 
 
			\(TraditionalForm
      \`\(u\_t\)(x\_i, t\_n) ≃ 
        \(u\_i\^\(n + 1\) - u\_i\^n\)\/Δt\)
			
(donde ahora estamos notando  ((Subscript[u, i])^n)≃ u((Subscript[x, i]), (Subscript[t, n])) la correspondiente aproximación de la solución exacta en los nodos de nuestra partición, tanto temporal como espacial).

	Este plantemiento da lugar al siguiente método explícito (notando  r ≡ (α^2)Δ(t/Δ)(x^2)), estable si y sólo si  r<(1/2):

	    ((Subscript[u, i])^(n + 1)) =  ((Subscript[u, i])^n) + r (\(TraditionalForm
      \`\(u\_\(i - 1\)\^n - 2  u\_i\^n + u\_\(i + 1\)\^n\) )\)
		= r ((Subscript[u, i - 1])^n) + (1 - 2  r) ((Subscript[u, i])^n) +  r ((Subscript[u, i + 1])^n)
	
Mientras que una aproximación análoga, pero "retrasada", evaluada en el nodo ((Subscript[x, i]), (Subscript[t, n + 1])), se podría escribir 

	 ((((Subscript[u, i])^(n + 1)) - ((Subscript[u, i])^n))/Δ)t≃ (Subscript[u, t])((Subscript[x, i]), (Subscript[t, n + 1])) 
	 	= (α^2) (Subscript[u, xx])((Subscript[x, i]), (Subscript[t, n + 1])) ≃ (α^2) \(TraditionalForm
      \`\(u\_\(i - 1\)\^\(n + 1\) - 2  u\_i\^\(n + 1\) + 
          u\_\(i + 1\)\^\(n + 1\)\)\/Δx\^2\)
	 
dando lugar ahora a un esquema de tipo implícito, como el que sigue:

 		-r ((Subscript[u, i - 1])^(n + 1))+(1+2r ) ((Subscript[u, i])^(n + 1)) - r ((Subscript[u, i + 1])^(n + 1)) =  ((Subscript[u, i])^n) 
 		
 	Otro posible método implícito que se puede plantear es el denominado método de Crank-Nicolson, promediando las diferencias centradas en los instantes (Subscript[t, n]) y (Subscript[t, n + 1]) para la aproximación de la derivada segunda en la variable espacial, 
 	
  ((((Subscript[u, i])^(n + 1)) - ((Subscript[u, i])^n))/Δ)t≃ (Subscript[u, t])((Subscript[x, i]), (Subscript[t, n])) =
  = (α^2) (Subscript[u, xx])((Subscript[x, i]), (Subscript[t, n]))  ≃ \(TraditionalForm
      \`\(\(α\^2\/2\) 
        \((\(u\_\(i - 1\)\^n - 2  u\_i\^n + 
                u\_\(i + 1\)\^n\)\/Δx\^2 + \ 
            \(u\_\(i - 1\)\^\(n + 1\) - 2  u\_i\^\(n + 1\) + 
                u\_\(i + 1\)\^\(n + 1\)\)\/Δx\^2)\)\ \)\)
  
  dando lugar al siguiente esquema numérico:
  
  -(r/2) ((Subscript[u, i - 1])^(n + 1))+(1+r )((Subscript[u, i])^(n + 1)) - (r/2) ((Subscript[u, i + 1])^(n + 1)) =  (r/2) ((Subscript[u, i - 1])^n)+(1-r )((Subscript[u, i])^n) + ((r )/2)((Subscript[u, i + 1])^n)
  
Estos esquemas de tipo implícito tienen la propiedad fundamental de resultar incondicionalmente estables, sea cual sea el valor de r; lo que hace que podamos emplear valores de Δt algo superiores que con los métodos de tipo explícito. Aunque también hay que tener en cuenta que un valor del paso temporal demasiado alto puede conllevar una pérdida de precisión temporal nada aconsejable.

Definición de los datos del problema

a := 0. \(b := 1
\)T := 1. \(α := 1
\)\(f[x_] := x\^2
\)(Subscript[h, 1])[x_] := 0(Subscript[h, 2])[x_] := 0

Definición de la discretización numérica

imax := 4nmax := 32

Clear[i, n, x, t, u, r]Δx := (b - a)/imax\(Δt := T/nmax
\)\(r := \(α\^2\) Δt\/Δx\^2
\)Do[(Subscript[x, i]) = a + i Δx, {i, 0, imax}]Do[(Subscript[t, n]) = n Δt, {n, 0, nmax}]

Table[(Subscript[x, i]), {i, 0, imax}]Table[(Subscript[t, n]), {n, 0, nmax}]

Aplicación de la condición de contorno

Do[u[0, n] = (Subscript[h, 1])[(Subscript[t, n])], {n, 0, nmax}]Do[u[imax, n] = (Subscript[h, 2])[(Subscript[t, n])], {n, 0, nmax}]

Aplicación de la condición inicial

Do[u[i, 0] = f[(Subscript[x, i])] // N, {i, 1, imax - 1}]

Table[u[i, n], {i, 0, imax}, {n, 0, nmax}] // TableForm

Evolución temporal mediante el esquema explícito

Do[Do[
        
u[i, n + 1] = 
          r u[i - 1, n] + (1 - 2  r) u[i, n] + 
            r u[i + 1, n]
, {i, 1, imax - 1}]
        
, {n, 0, nmax - 1}]

Table[u[i, n], {i, 0, imax}, {n, 0, nmax}] // TableForm

Evolución temporal mediante el esquema implícito retrógrado

\(matA = Table[0, {i, imax - 1}, {i, imax - 1}]; 

    Do[matA[\([i, i]\)] = 1 + 2  r, {i, imax - 1}]\)Do[matA[[i, i + 1]] = -r;  
      matA[[i + 1, i]] = -r, {i, imax - 2}]matA // MatrixForm

\(Do[
vecb = Table[u[i, n], {i, imax - 1}]; 
      
vecb[\([1]\)] = vecb[\([1]\)] + r\ u[0, n + 1]; 
      
vecb[\([imax - 1]\)] = 
        vecb[\([imax - 1]\)] + r\ u[imax, n + 1]; 
      
aux = LinearSolve[matA, vecb]; 
\ \ \ 
      Do[u[i, n + 1] = aux[\([i]\)], {i, imax - 1}]
, {n, 
        0, nmax - 1}]\)

Table[u[i, n], {i, 0, imax}, {n, 0, nmax}] // TableForm

Evolución temporal mediante el esquema implícito de Crank-Nicolson (aprovechando la estructura tridiagonal de la matriz de coeficientes)

Si queremos ahora aprovechar la estructura tridiagonal de la matriz de coeficientes del sistema lineal que tendremos que resolver a cada iteración de tiempo, vamos a cargar previamente un paquete de Mathematica que permite hacerlo sin mayor complicación mediante la orden  TridiagonalSolve

<< LinearAlgebra`Tridiagonal`

?TridiagonalSolve

\(aa = Table[\(-\(r\/2\)\), {i, imax - 2}]; \ cc = aa; 

    bb = Table[1 + r, {i, imax - 1}]; \)

\(Do[
vecb = 
        Table[r\/2\ u[i - 1, n] + \((1 - r)\) u[i, n] + r\/2\ u[i + 1, n], {
            i, imax - 1}]; 
      
vecb[\([1]\)] = vecb[\([1]\)] + r\/2\ u[0, n + 1]; 
      
vecb[\([imax - 1]\)] = 
        vecb[\([imax - 1]\)] + r\/2\ u[imax, n + 1]; 
      
aux = TridiagonalSolve[aa, bb, cc, vecb]; 
\ \ \ 
      Do[u[i, n + 1] = aux[\([i]\)], {i, imax - 1}]
, {n, 
        0, nmax - 1}]\)

Table[u[i, n], {i, 0, imax}, {n, 0, nmax}] // TableForm

Representación gráfica de la solución

ListPlot3D[Table[u[i, n], {i, 0, imax}, {n, 0, nmax}], 
      Ticks -> None, PlotRange -> All];

Animación gráfica de la solución

Do[
      ListPlot[Table[{(Subscript[x, i]), u[i, n]}, {i, 0, imax}], PlotJoined -> True, 
        PlotRange -> {-0.2, 1. }], 
{n, 0, nmax}];

Resolución de un problema mixto para la ecuación de ondas lineal

Se trata de la resolución aproximada, mediante el método de diferencias finitas del siguiente problema mixto para la ecuación de ondas lineal (prototipo de las ecuaciones de tipo hiperbólico):

	(Subscript[u, tt]) = (c^2)(Subscript[u, xx]),   a<x<b,  t>0
	u(a,t) = (Subscript[h, 1])(t),   u(b,t) = (Subscript[h, 2])(t),     t≥0
	u(x,0) = f(x),  (Subscript[u, t])(x, 0) =  g(x),  a<x<b	
	
Si seguimos notando  ((Subscript[u, i])^n)≃ u((Subscript[x, i]), (Subscript[t, n])) a la correspondiente aproximación de la solución exacta en los nodos de nuestra malla (producto cartesiano de una partición temporal por la correspondiente del intervalo espacial), podemos seguir planteando esquemas numéricos totalmente equivalentes a los ya considerados para el caso parabólico (de los cuales la ec. de difusión es un prototipo).

	Así pues un plantemiento centrado para cada  una de las derivadas de segundo orden 
	
	\(TraditionalForm
      \`\(u\_i\^\(n - 1\) - 2  u\_i\^n + u\_i\^\(n + 1\)\)\/Δt
            \^2 = \ c\^2\ 
          \(u\_\(i - 1\)\^n - 2  u\_i\^n + u\_\(i + 1\)\^n\)\/Δx
              \^2\)
	
da lugar al siguiente método explícito (notando ahora r ≡ (c^2)Δ((t^2)/Δ)(x^2)):

	((Subscript[u, i])^(n + 1)) = r ((Subscript[u, i - 1])^n) + 2 (1 - r) ((Subscript[u, i])^n) +  r ((Subscript[u, i + 1])^n) - ((Subscript[u, i])^(n - 1)).
 		
 	En lo que respecta a métodos de tipo implícito, ahora el conocido método de Crank-Nicolson resultará de aproximar
 	
  \(TraditionalForm
      \`\(u\_i\^\(n - 1\) - 2  u\_i\^n + u\_i\^\(n + 1\)\)\/Δt
          \^2\)≃ (Subscript[u, tt])((Subscript[x, i]), (Subscript[t, n])) =
  = (Subscript[u, xx])((Subscript[x, i]), (Subscript[t, n]))  ≃ \(TraditionalForm
      \`\(\(1\/2\) 
        \((\(u\_\(i - 1\)\^n - 2  u\_i\^n + 
                u\_\(i + 1\)\^n\)\/Δx\^2 + \ 
            \(u\_\(i - 1\)\^\(n + 1\) - 2  u\_i\^\(n + 1\) + 
                u\_\(i + 1\)\^\(n + 1\)\)\/Δx\^2)\)\ \)\)
  
  dando lugar al siguiente esquema numérico:
  
  -(r/2) ((Subscript[u, i - 1])^(n + 1))+(1+r )((Subscript[u, i])^(n + 1)) - (r/2) ((Subscript[u, i + 1])^(n + 1)) 
  		=  (r/2) ((Subscript[u, i - 1])^n)+(2-r )((Subscript[u, i])^n) + ((r )/2)((Subscript[u, i + 1])^n) - ((Subscript[u, i])^(n - 1))
  
Este y otros esquemas de tipo implícito para la ecuación de ondas también tienen la propiedad fundamental de tener un rango de estabilidad muy superior al de los esquemas de tipo explícito, lo que también va a permitir que podamos emplear valores de Δt algo superiores que con los métodos de tipo explícito. No obstante, se puede comprobar numéricamente que también hay que tener en cuenta que un valor de Δt  demasiado alto puede conllevar una pérdida de precisión temporal nada aconsejable.

	En cuanto a la aplicación de las correspondientes condiciones iniciales de este problema, podemos resaltar también que lo más inmediato sería empezar con  ((Subscript[u, i])^0) =  f((Subscript[x, i])) ≡ (Subscript[f, i]),  i=0,1,...,(Subscript[i, max])
y utilizar una simple aproximaci'on numérica de primer orden para la derivada temporal, (Subscript[g, i])≡ g((Subscript[x, i]))=(Subscript[u, t])((Subscript[x, i]), 0)≃((((Subscript[u, i])^1) - ((Subscript[u, i])^0))/(Δt ))+O(Δt), para despejar finalmente:  ((Subscript[u, i])^1) =  ((Subscript[u, i])^0) + Δt (Subscript[g, i]),  i=0,1,...,(Subscript[i, max]) (siendo dicha aproximación de segundo orden, es decir salvo un infinitésimo O(Δ(t^2)) ).

	No obstante, una mejor aproximación se puede obtener si utilizamos la ecuación satisfecha, dentro del siguiente desarrollo de Taylor (con  (Subscript[t, 0])=0  y Δt = (Subscript[t, 1]) - (Subscript[t, 0]))
	
	u((Subscript[x, i]), (Subscript[t, 1])) = u((Subscript[x, i]), (Subscript[t, 0])) + Δt (Subscript[u, t])((Subscript[x, i]), (Subscript[t, 0])) + \(TraditionalForm
      \`\(Δt\^2\/2\) \(\(u\_tt\)(x\_i, t\_0)\)\) + O(Δ(t^3))
	             = ((Subscript[u, i])^0) + Δt (Subscript[g, i]) + \(TraditionalForm
      \`\(Δt\^2\/2\) \(c\^2\) \(\(u\_xx\)(x\_i, t\_0)\)\) + O(Δ(t^3))
	             = ((Subscript[u, i])^0) + Δt (Subscript[g, i]) + \(TraditionalForm
      \`\(Δt\^2\/2\) \(c\^2\) \(f'\)' \((x\_i)\)\) + O(Δ(t^3))
	             = ((Subscript[u, i])^0) + Δt (Subscript[g, i]) + \(TraditionalForm
      \`\(Δt\^2\/2\) \(c\^2\) 
        \(f\_\(i - 1\) - 2  f\_i + f\_\(i + 1\)\)\/Δx\^2\) 
	             		+ O(Δ(t^3)) + O(Δ(t^2) Δ(x^2)) 
	             
donde también se ha empleado la aproximación habitual para la derivada segunda de la función de una variable f(x) en el nodo (Subscript[x, i]) (que sabemos que es de segundo orden):  f'' ((Subscript[x, i]))≃\(TraditionalForm
      \`\(f\_\(i - 1\) - 2  f\_i + f\_\(i + 1\)\)\/Δx\^2\)

Ejemplo de aplicación

Resolveremos en primer lugar el siguiente problema mixto para la ecuación de ondas: 

	           (Subscript[u, tt]) =  (Subscript[u, xx]),   0<x<π,  t>0
	u(0,t) = 0 = u(1,t),     t≥ 0
	u(x,0) = sen(x),  (Subscript[u, t])(x, 0) =  0,  0<x<π

cuya solución exacta se puede calcular mediante separación de variables y resulta ser:
			u(x,t) = cos(t) sen(x)

Definición de los datos del problema

a := 0\(b := π
\)\(T := 1
\)\(c := 1
\)f[x_] := Sin[x]\(g[x_] := 0
\)(Subscript[h, 1])[x_] := 0(Subscript[h, 2])[x_] := 0

Comprobación de la solución exacta

Analíticamente

Clear[x, t, u]uu[x_, t_] := Cos[t] Sin[x]vv[x_, t_] := -Sin[t] Sin[x]

\(∂\_{t, 2}\ uu[x, t] == ∂\_{x, 2}\ uu[x, t]\)uu[0, t]uu[π, t]uu[x, 0](∂_t) uu[x, t] /. t -> 0

Gráficamente, veamos que verifica las condiciones iniciales

\(t = 0; 

    Plot[{f[x], uu[x, t]}, {x, 0, 1}, 
      PlotStyle -> {RGBColor[1, 0, 0], Dashing[{0.1, 0.2}]}]; \)

\(t = 0; 

    Plot[{g[x], vv[x, t]}, {x, 0, 1}, 
      PlotStyle -> {RGBColor[1, 0, 0], Dashing[{0.1, 0.2}]}]; \)

Definición de la discretización numérica

imax := 10nmax := 20

Clear[i, n, x, t, u, r]Δx := (b - a)/imax\(Δt := T/nmax
\)\(r := \(c\^2\) Δt\^2\/Δx\^2
\)Do[(Subscript[x, i]) = a + i Δx, {i, 0, imax}]Do[(Subscript[t, n]) = n Δt, {n, 0, nmax}]

Table[(Subscript[x, i]) = a + i Δx, {i, 0, imax}]Table[(Subscript[t, n]) = n Δt, {n, 0, nmax}]

Aplicación de la condición de contorno

Do[u[0, n] = (Subscript[h, 1])[(Subscript[t, n])], {n, 0, nmax}]Do[u[imax, n] = (Subscript[h, 2])[(Subscript[t, n])], {n, 0, nmax}]

Aplicación de la condición inicial

Do[u[i, 0] = f[(Subscript[x, i])] // N; 
      
u[i, 1] = u[i, 0] + Δt g[(Subscript[x, i])], {i, 
        1, imax - 1}]

Otra posibilidad más precisa para la condición inicial

Do[u[i, 0] = f[(Subscript[x, i])] // N; 
      
u[i, 1] = 
        u[i, 0] + Δt g[(Subscript[x, i])] + 
          (r/2) (f[(Subscript[x, i - 1])] - 2  f[(Subscript[x, i])] + f[(Subscript[x, i + 1])]), {i, 
        1, imax - 1}]

Evolución temporal mediante el esquema explícito

Do[Do[
        
u[i, n + 1] = 
          r u[i - 1, n] + 2 (1 - r) u[i, n] + r u[i + 1, n] - 
            u[i, n - 1]
, {i, 1, imax - 1}]
        
, {n, 1, nmax - 1}]

Table[SetPrecision[u[i, n], 4], {i, 0, imax}, {n, 0, nmax}] // TableForm

Comparación entre la solución exacta y la aproximada en  t=T

Table[u[i, nmax], {i, 0, imax}]Table[uu[(Subscript[x, i]), T], {i, 0, imax}] // N

grafexact = 
      Plot[uu[x, T], {x, a, b}, PlotStyle -> RGBColor[1, 0, 0], 
        DisplayFunction -> Identity]; 
\(grafaprox = 
      ListPlot[Table[{x\_i, u[i, nmax]}, {i, 0, imax}], 
        PlotStyle -> PointSize[0.02], DisplayFunction -> Identity]; 

    


    \(Show[grafexact, grafaprox, DisplayFunction -> $DisplayFunction]; 
    \)\)

Conviene aquí, en función de los resultados obtenidos, probar con distintos valores de Δx  y Δt  para calibrar con cuáles se obtienen mejores resultados, teniendo en cuenta además que los esquemas de tipo explícito siempre van a tener que satisfacer ciertas restricciones en cuanto al valor del llamado coeficiente CFL (de Courant-Friedrich-Levy), intimamente relacionado con el valor que nosotros hemos denominado r≡ (c^2)Δ((t^2)/Δ)(x^2).

Evolución temporal mediante el esquema implícito de Crank-Nicolson (aprovechando la estructura tridiagonal de la matriz de coeficientes)

Si queremos ahora aprovechar la estructura tridiagonal de la matriz de coeficientes del sistema lineal que tendremos que resolver a cada iteración de tiempo, vamos a cargar previamente un paquete de Mathematica que permite hacerlo sin mayor complicación mediante la orden  TridiagonalSolve

<< LinearAlgebra`Tridiagonal`

?TridiagonalSolve

\(aa = Table[\(-\(r\/2\)\), {i, imax - 2}]; \ cc = aa; 

    bb = Table[1 + r, {i, imax - 1}]; \)

\(Do[
vecb = 
        Table[r\/2\ u[i - 1, n] + \((1 - r)\) u[i, n] + r\/2\ u[i + 1, n], {
            i, imax - 1}]; 
      
vecb[\([1]\)] = vecb[\([1]\)] + r\/2\ u[0, n + 1]; 
      
vecb[\([imax - 1]\)] = 
        vecb[\([imax - 1]\)] + r\/2\ u[imax, n + 1]; 
      
aux = TridiagonalSolve[aa, bb, cc, vecb]; 
\ \ \ 
      Do[u[i, n + 1] = aux[\([i]\)], {i, imax - 1}]
, {n, 
        0, nmax - 1}]\)

Table[u[i, n], {i, 0, imax}, {n, 0, nmax}] // TableForm

Representación gráfica de la solución

ListPlot3D[Table[u[i, n], {i, 0, imax}, {n, 0, nmax}], 
      Ticks -> None];

Animación gráfica de la solución

Para realizar una animación gráfica de la evolución temporal de la solución bastará con ejecutar el bucle que hay a continuación (ajustando convenientemente los valores de la opción PlotRange  si fuese necesario) y posteriormente seleccionar con el ratón el grupo de celdillas gráficas que se generarán y activar la animación desde el Menú Cell: Animate Selected Graphics (o bien con la combinación de teclas Ctrl-Y).

Do[
      ListPlot[Table[{(Subscript[x, i]), u[i, n]}, {i, 0, imax}], PlotJoined -> True, 
        PlotRange -> {-2, 2}], 
{n, 0, nmax}];

Ejercicios

1.-

Resolver de ahora el siguiente ejemplo concreto para la ecuación de Laplace en un cuadrado de lado  π:

Δu(x, y)  =  0,   ∀(x,y)∈Ω≡{(x,y) /  0<x<π,  0<y<π}
(Subscript[u, |∂Ω])≡φ,   
( (Subscript[φ, |{y = 0}]) ≡ 0,  
      (Subscript[φ, |{x = π}]) ≡ (cos^2)(y),  
      (Subscript[φ, |{y = π}]) ≡ 0,  
      (Subscript[φ, |{x = 0}]) ≡ 0)
	
cuya solución exacta se puede calcular mediante el método de separación de variables y resulta ser (comprobarlo):

	u(x,y) = sinh(x) sinh(y) +  
	
	+ \(TraditionalForm
      \`∑\+\(n = 3\)\%∞\( 2\/\(π\ \(sinh(n\ π)\)\)\) 
          \((\(1 - \((\(-1\))\)\^n\)\/\(2  n\) - 
              \(\(\((\(-1\))\)\^n\) n\)\/\(n\^2 - 1\))\)\)

2.-

Sea Ω el siguiente dominio rectangular 

	{(x,y)∈(ℝ^2) / 0 < x, y < 1}≡[0,1]×[0,1], 
	
con frontera  Γ≡∂Ω; plantear y resolver los siguiente problemas elípticos definidos en dicho dominio

i)	(Subscript[u, xx]) + (Subscript[u, yy])+ γ u = f(x,y),   (x,y)∈Ω
	∂(u/∂)n((x,y)) = g(x,y),            (x,y)∈∂Ω
	
siendo  γ≡γ(x,y) una función suficientemente regular (continua al menos) y ∂(u/∂)n la correspondiente derivada normal en cada uno de los puntos de la frontera donde tenga sentido.

ii)	((Subscript[α (Subscript[u, x])), x]) + (Subscript[((β (Subscript[u, y]))), y])+ γ u = f(x,y),   (x,y)∈Ω
	u((x,y)) = g(x,y),                (x,y)∈∂Ω
	
siendo α≡α(x,y),  β≡β(x,y), y  γ≡γ(x,y)  funciones suficientemente regulares en Ω .

3.-

Resolver el siguiente problema mixto para la ecuación del calor: 

	    (Subscript[u, t]) = (Subscript[u, xx]),  		0<x<1,  t>0
	u(0,t) = 0 = u(1,t),     	t≥ 0
	u(x,0) = x,  	       	  0<x<1

cuya solución exacta se puede calcular mediante separación de variables y resulta ser (comprobarlo):

	u(x,t) = cos(4π t) sen(π x) + 
	
	+\(TraditionalForm
      \`∑\+\(k = 1
            \)\%∞\(\( 2\ \((\(-1\))\)\^k\)\/\(π\ k\)\) 
          \(sen(k\ π\ x)\)\ e\^\(\(-k\^2\)\ π\^2\ t\)\)

experimentando con distintos valores del paso en x y en t, de manera que se obtenga una aproximación estable y lo suficientemente buena de la solución exacta.

4.-

Hacer lo propio con la solución del problema mixto parabólico: 

	 (Subscript[u, t]) = (Subscript[u, xx])+1,   	0<x<1,  t>0
	u(0,t) = 0 = u(1,t),     	t≥ 0
	u(x,0) = 0,  	         	0<x<1

y estudiar a qué tiende la solución cuando T->∞.

5.-

Resolver los siguientes problemas mixtos para la ecuación de ondas: 

	   (Subscript[u, tt]) =  (Subscript[u, xx]),   		0<x<1,  t>0
	u(0,t) = 0 = u(1,t),     	t≥ 0
	u(x,0) = f(x),  (Subscript[u, t])(x, 0) =  g(x),  0<x<1

con las siguientes condiciones iniciales:

	i) 	  f(x)=0, g(x)=1,  	x∈[0,1]
	
	ii)	f(x)=sen(π x), g(x)=0,  x∈[0,1]

6.-

Resolver finalmente el siguiente problema mixto para la ecuación de ondas: 

	           (Subscript[u, tt]) = 16  (Subscript[u, xx]),   0<x<1,  t>0
	u(0,t) = 0 = u(1,t),     t≥ 0
	u(x,0) = sen(π x),  (Subscript[u, t])(x, 0) =  x - (x^2),  0<x<1

cuya solución exacta se puede calcular mediante separación de variables y resulta ser (comprobarlo):

u(x,t) = cos(4π t) sen(π x) + \(TraditionalForm
      \`∑\+\(k = 1\)\%∞\( 2 \( sen(4 \((2  k - 1)\) t)\)\ 
            \(sen(\((2  k - 1)\) x)\)\)\/\(\(π\^4\)(2  k - 1)\)\^4\)