excerpt from Numerical and Analytical Methods for Scientists and Engineers

The material in this notebook is copyright John Wiley and Sons,2003. Bulk reproduction or resale is prohibited.

Multiple Spatial Dimensions

Explicit Methods

The explicit methods discussed previously generalize easily to PDEs in multiple spatial dimensions. Continuing with our heat equation model, we now consider the problem in two dimensions:

For simplicity, we take the boundary to be a rectangular box of sides (Subscript[L, x]) and (Subscript[L, y]), and create a rectangular spatial grid in x  and y according to

where Δx = (Subscript[L, x])/(Subscript[M, x]) and Δy = (Subscript[L, y])/(Subscript[M, y]). (Nonrectangular boundaries are treated in Sec. 6.2.2.)

The FTCS method on this grid is a trivial extension of the method in one-dimension:

Equation (6.2.28) is valid for the points interior to the rectangle, and must be supplemented by boundary conditions that determine the points on the boundary, as well as by initial conditions specifying (Subscript[T, j k]^( 0)).von Neumann stability analysis for this difference equation considers the stability of  a Fourier mode of the form \(TraditionalForm\`ξ\^n\ \(e\^\(i\ \(k\_x\) x\)\) e\^\(i\ k\_y\ y\)\), and leads to a Courant condition that

(see the exercises).  The other explicit methods discussed previously, CTCS and Lax, also can be easily extended to multiple dimensions.

Example: Wave Equation with Nonuniform Wave Speed

As an example of an explicit method in two dimensions,  we construct the numerical solution to a wave equation in two dimensions with nonuniform wave speed,

\(TraditionalForm\`\(∂\^2\/\((∂t)\)\^2\) \(z(x, y, 
            t)\) = \(\(c\^2\)(x, 
              y)\)\ [∂\^2\/∂x\^2 + ∂\^2\/∂y\^2] \(z(x, y, t)\)\),

with initial conditions \(TraditionalForm\`z(x, y, 0) = \(z\_o\)(x, y), \ \(z\& . \)(x, y, 
          0) = \(v\_o\)(x, y)\). Here we take the simple case (Subscript[z, o]) = (Subscript[v, o]) = 0. For boundary conditions, we choose a square domain with unit sides and take an oscillating boundary condition,

z(0, y, t) = z(1, y, t) = z(x, 1, t) = 0;  
      z(x, 0, t) = x(1 - x) sin 2  t.

For the wave speed, we choose

c(x, y)  = 
        (1/10) (1  - 
              0.7 exp[-2 (((x - 0.5))^2) - 5 (((y - 0.5))^2)]).

The wave speed, plotted in Cell 6.72,  has a minimum at the center of the box.

\(\(Clear["Global`*"];\)
\)
\(\(cw[x_, y_]\  = \ 
        1/10\ \((1\  -  .7\ Exp[\(-2\) \((x - 0.5)\)\^2 - 
                    5 \((y - 0.5)\)\^2])\);\)
\)
Plot3D[cw[x, y], {x, 0, 1}, {y, 0, 1}, 
        AxesLabel -> {"x", "y", ""}, 
        PlotLabel -> "wave speed"];

For the numerical method, we choose CTCS. This explicit method works well for the wave equation, being second-order in both the  space and time step sizes, and allowing relatively large steps according to the von Neumann stability criterion: for Δx=Δy=Δ,  (c^2)((Δt)^2)/(Δ^(2 ))≤ 1/2 for stability (see the exercises). Therefore, if  c ≤  1/10 and Δ =  1/30, Δt ≤ 1/3 (Sqrt[2])=0.24. Below, we choose Δt = 0.2:

(Subscript[M, x]) = 30;  (Subscript[M, y]) =  30;
(Subscript[L, x ]) =  (Subscript[L, y]) = 1;
Δx  =  (Subscript[L, x])/(Subscript[M, x]);  Δy  =  
          (Subscript[L, y])/(Subscript[M, y]);
Δt  =  0.2;
x[j_]  :=  j Δx;
y[k_]  :=  k Δy;

  
z0[j_, k_]  :=  0;
v0[j_, k_]  :=  0;

  
z[j_, k_, 0]  :=  z0[j, k];

  
c[j_, k_]  :=  cw[x[j], y[k]];

  
z[0, k_, n_]  :=  0;
z[(Subscript[M, x]), k_, n_]  :=  0;
z[j_, 0, n_]  :=  
            x[j] (1 - x[j]) Sin[2  n Δt];
z[j_, (Subscript[M, y]), n_]  :=  0;

The CTCS method for the wave equation is

(Subscript[z, j k]^(n + 1))=2  (Subscript[z, j k]^n)-(Subscript[z, j k]^(n - 1))+(((Subscript[c, j k]^( 2)) ((Δt)^2))/((Δx)^2))((Subscript[z, j + 1 k]^n) - 2  (Subscript[z, j k]^n) + 
                  (Subscript[z, j - 1 k]^( n)))+(((Subscript[c, j k]^( 2)) ((Δt)^2))/((Δy)^2))((Subscript[z, j k + 1]^n) - 2  (Subscript[z, j k]^n) + 
                  (Subscript[z, j k - 1]^( n)) ).

We construct the recursion relation for this method below:

z[j_, k_, 
        n_]  :=  (z[j, k, n]  =  
            2 z[j, k, n - 1]  - 
              z[j, k, n - 
                  2]  +   ((Δ(t^2) 
                      c[j, k]^2)/Δ)(x^2) ((z[j + 1, k, 
                          n - 1]  -  2 z[j, k, n - 1]  +  
                        z[j - 1, k, n - 1])  +  (z[j, k + 1, 
                          n - 1]  -  2 z[j, k, n - 1]  +  
                        z[j, k - 1, n - 1])))  /;  
        n ≥  2

The method can only be used for n ≥  2 because evaluation of the nth time step  refers back to the n - 2nd step. Therefore, for the first step, we use the trick employed in Chapter 1 for second-order ODEs [see Eq. (1.4.28) and the surrounding discussion]:

(Subscript[z, j k]^1)=(Subscript[z, j k]^0)+(Subscript[((Subscript[v, o])), j k])Δt+(((Δt)^2)/2)(Subscript[a, j  k]),
where (Subscript[((Subscript[v, o])), j k]) is the initial velocity (equal to zero for our choice of initial conditions), and (Subscript[a, j k]) is the initial acceleration, given by

(Subscript[a, j k])=(((Subscript[c, j k]^( 2)) )/((Δx)^2))((Subscript[z, j + 1 k]^0) - 2  (Subscript[z, j k]^0) + 
                  (Subscript[z, j - 1 k]^( 0)))+(((Subscript[c, j k]^( 2)) ((Δt)^2))/((Δy)^2))((Subscript[z, j k + 1]^0) - 2  (Subscript[z, j k]^0) + 
                  (Subscript[z, j k - 1]^( 0)) ).

(this is also equal to zero for our initial conditions). Therefore, we simply have (Subscript[z, j k]^1) = (Subscript[z, j k]^0)=0:

z[j_, k_, 1]  :=  0;

(Of course, this cell must be modified if other initial conditions are used.) At each time step, we construct an interpolation function of the solution in x and y:

zsol[n_]  :=  zsol[n]  =  
        Interpolation[
          Flatten[Table[{j Δx, k Δy, 
                z[j, k, n]}, {j, 0, (Subscript[M, x])}, {k, 0, (Subscript[M, y])}], 1]]

In Cell 6.77 we plot the solution.

Table[
        Plot3D[Evaluate[zsol[n][x, y]], {x, 0, 1}, {y, 0, 1}, 
          PlotPoints -> 60, PlotRange -> {- .5,  .5}, 
          Mesh -> False, 
          PlotLabel -> "t="  <> ToString[n Δt]], {n, 0, 
          140, 4}];

The wavefronts curve toward the region of lower c according to Snell's law, and as they do so the amplitude and wavelength decrease, as expected from our WKB analysis of this equation in Chapter 5.  However, note the formation of a focus at a point behind the region of low wave speed, where the amplitude increases as wave energy is concentrated. This region is acting as a lens for the waves [see Exercise (13) in Sec. 5.2, and the examples in Sec. 5.3.1 for  similar problems in which the ray trajectories for the waves are evaluated].