Numerical Methods: Mathematica Notebooks
(c) John H. Mathews, 1996

To accompany the text:

Numerical Methods:  for Mathematics, 
Science, and Engineering, 2nd Ed, 1992
Prentice Hall, Inc.
Englewood Cliffs, New Jersey 07632  U.S.A.
Prentice Hall, Inc.; USA, Canada, Mexico ISBN 0-13-624990-6
Prentice Hall, International Editions:      ISBN 0-13-625047-5

This free software is compliments of the author.
E-mail address:  in%"mathews@fullerton.edu"
http://titan.fullerton.edu/~mathews/

CONTENTS

Chapter 10.	Solution of Partial Differential Equations
	Algorithm 10.1	Finite-Difference Solution for the Wave Equation	
	Algorithm 10.2	Forward-Difference Method for the Heat Equation	
	Algorithm 10.3	Crank-Nicholson Method for the Heat Equation	
	Algorithm 10.4	Dirichlet Method for Laplace's Equation

Algorithm 10.1 (Finite-Difference Solution for the Wave Equation).  To approximate the solution for the wave equation  

   utt(x,t) = c2uxx(x,t)    over  R = {(x,t): 0 <= x <= a, 0 <= t <= b}  with   

   u(0,t) = 0,     u(a,t) = 0              for 0 <= t <= b     and  
   u(x,0) = f(x), ut(x,0) = g(x)        for 0 <= x <= a.  
Section	10.1	Hyperbolic Equations   Page  507

h  is the step in the x-direction 
k  is the step in the y-direction.

We will solve the problem on the grid
u[[i,j]]  where  i = 1,2,...,n  is the extent along the x-axis.
               and j = 1,2,...,m  is the extent along the t-axis.
The method proceds by simple iteration.
It is known that  r  must satisfy   r <= 1  for the method to be stable.

First, execute the following group of cells:

Off[General::spell1];
Clear[GridU,n,m,u];
GridU[n_,m_] :=
  Module[{},
  Clear[u];
  u = Table[1,{n},{m}];
  For[i=1,i<=n,i++,
    u[[i,1]] = f[i];
    u[[i,2]] = (1-r^2) f[i] + k g[i] +
                r^2/2 (f[i+1]+f[i-1]); ];
  For[j=1,j<=m,j++,
    Module[{},
      u[[1,j]] = 0;
      u[[n,j]] = 0; ]];
  Print["The set up grid for u[[i,j]] has been formed."];];

Clear[SolveU,n,m,u];
SolveU[n_,m_] :=
Module[{i,j},
For[j=3,j<=m,j++,
  Module[{},
  For[i=2,i<=n-1,i++,
    u[[i,j]] = (2-2 r^2) u[[i,j-1]] +
                r^2 (u[[i+1,j-1]]+u[[i-1,j-1]]) -
                u[[i,j-2]]; ]]];                
  Print["The solution has been computed."];];
On[General::spell1];

Chapter 10, Example 10.1, Page 504.

Consider the wave equation where  c^2 = 1.  The length of the rod is  L = 1.
Assume that the initial position is   u(x,0) = f(x) = sin(Pi x) + sin(2Pi x).
Compare the solution with the exact solution:
     u(x,t) = sin(Pi x)cos(Pi t) + sin(2Pi x)cos(4Pi t)
Solution.  We will use  c = 1
                                    h = 0.2.
                                    k = 0.04
This forces  r = 1^2 0.04/(0.2)^2  =  1.

First enter the boundary conditons:

Clear[f,g];
f[i_] := Sin[Pi 1/(n-1) (i-1)] +
         Sin[2 Pi 1/(n-1) (i-1)]//N;
g[i_] := 0.0;

Now set up the table of solutions.

c = 2.0;
h = 0.1;
k = 0.05;
r = c k/h;
n=11;
m=21;
GridU[n,m]

SolveU[n,m]

ListPlot3D[u,
AxesLabel->{"t(i)","x(i)","u"},ViewPoint->{4,2,3}];

Compare with the analytic solution.

SetOptions[Plot3D, PlotPoints->{n,m}];
Plot3D[Sin[Pi x]Cos[2 Pi t] +
         Sin[2 Pi x]Cos[4 Pi t],
         {x,0,1},{t,0,1}];

To see the numerical values enter the command:

Transpose[N[u,3]]

Chapter 10, Example 10.2, Page 506

Consider the wave equation where  c^2 = 1.
The length of the rod is  L = 1.  Assume that the initial position is
u(x,0) = f(x) = 5x/3   for  0 < x < 3/5
                      5/2 - 5x/2  for  3/5 < x < 1
Solution.  We will use  c = 1
                                    h = 0.2.
                                    k = 0.04
This forces  r = 1^2 0.04/(0.2)^2  =  1.

Boundary conditons.
We need 11 function values for f.

f[1]  = 0.0;  f[2]  = 0.1;
f[3]  = 0.2;  f[4]  = 0.3;
f[5]  = 0.4;  f[6]  = 0.5;
f[7]  = 0.6;  f[8]  = 0.45;
f[9]  = 0.3;  f[10] = 0.15;
f[11] = 0.0;

g[i_] := 0.0;

Now set up the table of solutions.

c = 2.0;
h = 0.1;
k = 0.05;
r = c k/h;
n=11;
m=21;
GridU[n,m]

SolveU[n,m]

ListPlot3D[u,
AxesLabel->{"t(i)","x(i)","u"},ViewPoint->{4,2,3}];

To see the numerical values enter the command:

Transpose[N[u,3]]

Algorithm 10.2  (Forward-Difference Method for the Heat Equation).  To approximate the solution for the heat equation  

   ut(x,t) = c2uxx(x,t)   over  R = {(x,t): 0 <= x <= a , 0 <= t <= b}  with  

   u(x,0) = f(x)  for 0 <= x <= a    and   
   u(0,t)  =  c1,  u(a,t)  =  c2    for 0 <= t <= b.  
Section	10.2	Parabolic Equations   Page  516

We will solve the problem on the grid

u[[i,j]]  where  i = 1,2,...,n  is the extent along the x-axis.
               and j = 1,2,...,m  is the extent along the t-axis.
               
The function  g1[j]  furnishes the values  u[[1,j]]  and
the function  g2[j]  furnishes the values  u[[n,j]]  and
the function  f[i]  furnishes the initial condition u[[i,1]].
For some problems the corners must be modified.
The grid  u[[i,j]  is formed with the subroutine  GridU[n,m].

First, execute the following group of cells:

Off[General::spell1];
GridU[n_,m_] :=
  Module[{},
  Clear[u];
  u = Table[1,{n},{m}];
  For[i=1,i<=n,i++,
      u[[i,1]] = f[i] ];
  For[j=1,j<=m,j++,
    Module[{},
      u[[1,j]] = g1[j];
      u[[n,j]] = g2[j]; ]];
  Print["The grid u[[i,j]] has been formed."];];

SolveU[n_,m_] :=
Module[{i,j},
  For[j=2,j<=m,j++,
    For[i=2,i<=n-1,i++,
      u[[i,j]]=(1-2r)*u[[i,j-1]] +
                 r*(u[[i-1,j-1]] + u[[i+1,j-1]]) ];];
  Print["The solution has been computed."]; ];
On[General::spell1];

Chapter 10, Example 10.3, Page 511.

Consider the heat equation where  c^2 = 1.   The length of the rod is  L = 1.
Assume that the ends of the rod are held at the temperature u=0.
Assume that the initial temperature distribution is
      u(x,0) = f(x) = 4x - 4x^2.
Apply the forward difference method with r=1 and
obtain temperature distributions for t=0.04, 0.08, 0.12,...,0.40.

Solution.  We will use  c = 1
                                    h = 0.2.
                                    k = 0.04
This forces  r = 1^2 0.04/(0.2)^2  =  1.

Enter the boundary conditons:

Clear[g1,g2,f];
g1[j_] := 0.0;
g2[j_] := 0.0;
f[i_] := 4*(0.2*(i-1))*(1 - 0.2*(i-1))//N;

Now set up the table of solutions.

n=6;
m=11;
GridU[n,m]

Next, solve it.

v = u;
c = 1;
h = 0.2;
k = 0.02;
r = c^2 k/h^2;
SolveU[n,m]

ListPlot3D[u,
AxesLabel->{"t(i)","x(i)","u   "},ViewPoint->{4,2,3}];
TableForm[N[Transpose[u],6]]

Now investigate what happens when the step size is too large.

n=6;
m=11;
GridU[n,m]

Next, solve it.

v = u;
c = 1;
h = 0.2;
k = 1/30;
r = c^2 k/h^2;
SolveU[n,m]

ListPlot3D[u,
AxesLabel->{"t(i)","x(i)","u   "},ViewPoint->{4,2,3}];
TableForm[N[Transpose[u],6]]

Algorithm 10.3  (Crank-Nicholson Method for the Heat Equation).  To approximate the solution for the heat equation  

   ut(x,t) = c2uxx(x,t)    over  R = {(x,t): 0 <= x <= a , 0 <= t <= b}   with    

   u(x,0) = f(x)   for  0 <= x <= a   and   
   u(0,t)  =  c1,  u(a,t)  =  c2     for 0 <= t <= b.  
Section	10.2	Parabolic Equations   Page  517

We will solve the problem on the grid

u[[i,j]]  where  i = 1,2,...,n  is the extent along the x-axis.
               and j = 1,2,...,m  is the extent along the t-axis.
               
The function  g1[j]  furnishes the values  u[[1,j]]  and
the function  g2[j]  furnishes the values  u[[n,j]]  and
the function  f[i]  furnishes the initial condition u[[i,1]].
For some problems the corners must be modified.
The grid  u[[i,j]  is formed with the subroutine  GridU[n,m].

First, execute the following group of cells:

Off[General::spell1];
Clear[GridU,n,m,u];
GridU[n_,m_] :=
  Module[{i,j},
  Clear[u];
  u = Table[1,{n},{m}];
  For[i=1,i<=n,i++,
    u[[i,1]] = f[i] ];
  For[j=1,j<=m,j++,
    Module[{},
      u[[1,j]] = g1[j];
      u[[n,j]] = g2[j]; ]];
  Print["The grid u[[i,j]] has been formed."];];

Clear[CoeffMat,n];
CoeffMat[n_] :=
Module[{},
  Clear[r];
  Va = Table[-1,{n-1}];
  Va[[n-1]]=0;
  Vc = Table[-1,{n-1}];
  Vc[[1]]=0;
  Vd = Table[2 + 2/r,{n}];
  Vd[[1]]=1;
  Vd[[n]]=1;
  Print["The coefficient matrix has been formed."];]

Clear[TriMat,a,d,c,b,n];
TriMat[a_,d_,c_,b_,n_] :=
  Module[{a0,b0,c0,d0,r,t},
    a0=a;
    b0=b;
    c0=c;
    d0=d;
    Clear[x0];
    x0 = Table[0,{n}];
    For[r=2,r<=n,r++,
      Module[{},
        t = a0[[r-1]]/d0[[r-1]];
        d0[[r]] = d0[[r]] - t c0[[r-1]];
        b0[[r]] = b0[[r]] - t b0[[r-1]] ]; ];
    x0[[n]] =  b0[[n]]/d0[[n]];
    For[r=n-1,1<=r,r--,
       x0[[r]] = (b0[[r]] - c0[[r]] x0[[r+1]])/
                  d0[[r]] ];
    Return[x0] ];

Clear[SolveU,n,m];
SolveU[n_,m_] :=
Module[{i,j},
  Clear[b,x];
  b = Table[0,{n}];
  x = Table[0,{n}];
  For[j=2,j<=m,j++,
    Module[{},
      b[[1]] =  g1[j];
      b[[n]] =  g2[j];
      For[i=2,i<=n-1,i++,
        b[[i]] = u[[i-1,j-1]] +
          (2/r -2)*u[[i,j-1]] + u[[i+1,j-1]] ];
      x = TriMat[Va,Vd,Vc,b,n];
      For[i=1,i<=n,i++,
        u[[i,j]] = x[[i]]]; ]; ]; ];
On[General::spell1];

Chapter 10, Example 10.4, Page 515.

Consider the heat equation where  c^2 = 1.  The length of the rod is  L = 1.
Assume that the ends of the rod are held at the temperature u=0.
Assume that the initial temperature distribution is
     u(x,0) = f(x) = sin(Pi x).
Apply the Crank-Nicholson method with r=1 and
obtain temperature distributions for t=0.04, 0.08, 0.12,...,0.40.

Compare the solution with the exact solution:
     u(x,t) = sin(Pi x) exp(-t/Pi^2).

Solution.  We will use  c = 1
                                   h = 0.2.
                                   k = 0.04
This forces  r = 1^2 0.04/(0.2)^2  =  1.

First set up the boundary conditons:

Clear[g1,g2,f];
g1[j_] := 0.0;
g2[j_] := 0.0;
f[i_] := Sin[Pi 0.1 (i-1)]+Sin[3 Pi 0.1 (i-1)] //N;

Now set up the table of solutions.

n=11;
m=11;
GridU[n,m]

Setting up the tri-diagonal matrx with n rows.  Indeed, we could get away
with n-2 rows, but the implementation is nice this way.
The following matrix will usually use r = 1.

CoeffMat[n]

Next, solve it.

c = 1;
h = 0.1;
k = 0.01;
r = c^2 k/h^2;
SolveU[n,m]

ListPlot3D[u,
AxesLabel->{"t(i)","x(i)","u   "},ViewPoint->{4,2,3}];

Compare with the analytic solution.

SetOptions[Plot3D, PlotPoints->{n,m}];
Plot3D[Sin[Pi x]Exp[-0.1 Pi^2 t]+
          Sin[3 Pi x]Exp[-0.1 9 Pi^2 t],
         {t,0,1},{x,0,1},ViewPoint->{4,2,3}];

Algorithm 10.4  (Dirichlet Method for Laplace's Equation).  
To approximate the solution of  
   uxx(x,y) +  uyy(x,y) = 0             over  R = {(x,y): 0 <= x <= a , 0 <= y <= b}  with 
   u(x,0) = f1(x), u(x,b) = f2(x)     for 0 <= x <= a   and  
   u(0,y) = f3(y), u(a,y) = f4(y)     for 0 <= y <= b.  
It is assumed that Dx = Dy = h and that integers  n  and  m  exist so that
a = nh  and  b = mh. 
Section	10.3	Elliptic Equations   Page  531

First, execute the following group of cells:

Off[General::spell1];
Clear[Dirichlet,q,n,m,w];
Dirichlet[q_List,n_,m_,w_]:=
  Module[{i,j,k},
  Err = 1; k = 0;
  While[(Err > 0.001) && (k <= 70),
    Module[{},
      Err = 0.0;
      For[j=2,j<=m-1,j++,
        Module[{},
          For[i=2,i<=n-1,i++,
            Module[{},
              Relax = w/4.0 ( u[[i,j+1]] + u[[i,j-1]] + 
                u[[i+1,j]] + u[[i-1,j]] - 4.0 u[[i,j]]);
              u[[i,j]] = u[[i,j]] + Relax;
              Err = Max[Err,Abs[Relax]];
            ]]]];
      Print["Max grid change = ",Err];
      k = k+1]]];

Clear[Neumann,qList,n,m,w];
Neumann[q_List,n_,m_,w_]:=
  Module[{i,j,k},
  Err = 1; k = 0;
  While[(Err > 0.001) && (k <= 70),
    Module[{},
      Err = 0.0;
      For[j=2,j<=m-1,j++,
        Module[{},
          For[i=2,i<=n-1,i++,
            Module[{},
              Relax = w/4.0 ( u[[i,j+1]] + u[[i,j-1]] + 
                u[[i+1,j]] + u[[i-1,j]] - 4.0 u[[i,j]]);
              u[[i,j]] = u[[i,j]] + Relax;
              Err = Max[Err,Abs[Relax]];
            ]]]];
      For[j=2,j<=m-1,j++,
        Module[{},
          Relax = w/4.0 ( 2 u[[n-1,j]] + 
            u[[n,j-1]] + u[[n,j+1]] - 4.0 u[[n,j]]);
          u[[n,j]] = u[[n,j]] + Relax;
          Err = Max[Err,Abs[Relax]]]];
      Print["Max grid change = ",Err];
      k = k+1]]];
On[General::spell1];

Chapter 10, Example 10.7, Page 529.

Solve Laplace's equation over a 9 by 9 grid with boundary conditions
Top:       180
Left:        80
Bottom:   20
Right:        0

First set up the boundary values.

Clear[u]; u = Table[ 70, {9},{9}];
For[i=2,i<9,i++, u[[i,1]]=80]
For[i=2,i<9,i++, u[[i,9]]=0]
For[j=2,j<9,j++, u[[1,j]]=180]
For[j=2,j<9,j++, u[[9,j]]=20]
u[[1,1]] = (u[[1,2]]+u[[2,1]])/2;
u[[1,9]] = (u[[1,8]]+u[[2,9]])/2;
u[[9,1]] = (u[[8,1]]+u[[9,2]])/2;
u[[9,9]] = (u[[8,9]]+u[[9,8]])/2;

Next, solve it.

w = 4/(2+Sqrt[4-(Cos[Pi/8]+Cos[Pi/8])^2])//N;
Dirichlet[u,9,9,w]

ListPlot3D[Transpose[u],
AxesLabel->{"y(i)","x(i)","u   "},ViewPoint->{4,2,3}];
TableForm[N[u,4]]

Chapter 10, Example 10.7, Page 529.

Solve Laplace's equation with the Neuman boundary condition
over a 9 by 9 grid with boundary conditions
Top:            180
Left:             80
Bottom:   ux = 0
Right:             0

First set up the boundary values.

Clear[u]; u = Table[ 70, {9},{9}];
For[i=2,i<9,i++, u[[i,1]]=80]
For[i=2,i<9,i++, u[[i,9]]=0]
For[j=2,j<9,j++, u[[1,j]]=180]
For[j=2,j<9,j++, u[[9,j]]=20]
u[[1,1]] = (u[[1,2]]+u[[2,1]])/2;
u[[1,9]] = (u[[1,8]]+u[[2,9]])/2;
u[[9,1]]=80;
u[[9,9]]=0;
For[j=1,j<=9,j++,
    u[[9,j]] = u[[9,1]]+ (j-1)/8 (u[[9,9]]-u[[9,1]]); ];

Next, solve it.

w = 4/(2+Sqrt[4-(Cos[Pi/8]+Cos[Pi/8])^2])//N;
Neumann[u,9,9,w]

ListPlot3D[Transpose[u],
AxesLabel->{"y(i)","x(i)","u   "},ViewPoint->{4,2,3}];
TableForm[N[u,4]]