Numerical Methods: Mathematica 3.0 Notebooks
(c) John H. Mathews, 1998
To accompany the text:
Numerical Methods:  for Mathematics, 
Science, and Engineering, 2nd Ed, 1992
Prentice Hall, Inc.
Simon & Schuster,  One Lake Street
Upper Saddle River, New Jersey 07458  U.S.A.
Prentice Hall, Inc.; USA, Canada, Mexico ISBN 0-13-624990-6
Prentice Hall, International Editions:      ISBN 0-13-625047-5
www.Prenhall.com
This free software is compliments of the author.
E-mail address:      in%"mathews@fullerton.edu"

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.

Clear["Global`*"];
Off[General::"spell1"];
Clear[a, a0, b, b0, c, c0, CoeffMat, d, d0, Dirichlet, f, g, g1, 
          g2, GridU, i, j, k, m, n, Neumann, q, qList, r, SolveU, t, TriMat, 
          u, w, x];
On[General::"spell1"];

Algorithm 10.1 Finite-Difference Solution for the Wave Equation.    To approximate the solution for the wave equation:
(Subscript[u, tt]) (x, t) =  (c^2) (Subscript[u, xx]) (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),  (Subscript[u, t]) (x, 0) = g (x)  for   0 ≤ x ≤ a.
Section  10.1   Hyperbolic Equations   Page  507

Step sizes:  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.

Execute this cell to activate Finite-Difference Solution for the Wave Equation.
GridU[n_,m_]
SolveU[n_,m_]
n_ is the number of subintervals along the x-axis,
m_ is the number of subintervals along the t-axis.

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] + 
              1\/2\ r\^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 (π x) +  sin (2  π x).
Compare the solution with the exact solution:
	u (x, t)  =  
        sin (π x) cos (π t) +  
          sin (2  π x) cos (4  π 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[c, f, g, h, i, k, m, n, r]; 

    f[i_] := N[
        Sin[\(π\ \((i - 1)\)\)\/\(n - 1\)] + 
          Sin[\(2\ π\ \((i - 1)\)\)\/\(n - 1\)]]; 
\(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[π x] Cos[2 π t] + Sin[2 π x] Cos[4 π t], {
        x, 0, 1}, {t, 0, 1}];

To see the numerical values enter the command:

TableForm[Transpose[Chop[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 
u (x, 0) = f (x) ={({{x                for   0 ≤ x ≤(( 3)/5)}, {(5/2)-((5 x)/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.
For the boundary conditons we need  11  function values for  f (x).

\(Clear[c, f, g, h, i, k, m, n, r]; 
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:

TableForm[Transpose[Chop[N[u, 3]]]]

Algorithm 10.2  Forward-Difference Method for the Heat Equation.  
To approximate the solution for the heat equation  
   (Subscript[u, t]) (x, t) =  (c^2) (Subscript[u, xx]) (x, t)  over   R = {(x, t) :  0 ≤ x ≤ a,   
          0 ≤ t ≤ b }  
with   
   u(x,0)  =  f(x)   for 0 ≤ x ≤ a    and   0 ≤ x ≤ a  and,
   u(0,t)  =  (Subscript[c, 1]),   u(a,t)  =  (Subscript[c, 2])      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  (Subscript[g, 1])[j]  furnishes the values  u[[1, j]]  and
the function   (Subscript[g, 2])[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].

Execute this cell to activate Forward-Difference Method for the Heat Equation.
GridU[n_,m_]
SolveU[n_,m_]
n_ is the number of subintervals along the x-axis,
m_ is the number of subintervals along the t-axis.

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 - 2\ r)\)\ 
                  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)  =  4  x  -  4  (x^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[c, f, g1, g2, h, i, k, m, n, r, v]; 
g1[j_] := 0.0; 

    \(g2[j_] := 0.0; \)\)f[i_] := N[4 (0.2 (i - 1)) (1 - 0.2 (i - 1))];

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}];

To see the numerical values enter the command:

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}];

To see the numerical values enter the command:

TableForm[N[Transpose[u], 6]]

Algorithm 10.3  Crank-Nicholson Method for the Heat Equation.  
To approximate the solution for the heat equation  
(Subscript[u, t]) (x, t) =  (c^2) (Subscript[u, xx]) (x, t)    over    R = {(x, t) :  0 ≤ x ≤ a,   
          0 ≤ t ≤ b }  
with   
u (x, 0) =  f (x)  for 0 ≤ x ≤ a    and    0 ≤ x ≤ a  and,
u (0, t) =  (Subscript[c, 1]),  u (a, t) =  (Subscript[c, 2])  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  (Subscript[g, 1])[j]  furnishes the values  u[[1, j]]  and
the function   (Subscript[g, 2])[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].

Execute this cell to activate Crank-Nicholson Method for the Heat Equation.
GridU[n_,m_]
CoeffMat[n_]
TriMat[a_,d_,c_,b_,n_]
SolveU[n_,m_]
a_,d_,c_,b_ are vectors needed in the solution,
n_ is the number of subintervals along the x-axis,
m_ is the number of subintervals along the t-axis.

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]]]; ]; \)]; 
		
        Print["The solution has been computed."]]; 

    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 (π 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 (π x) (e^(-t(π^2))) +   
          sin (3  π x) (e^(-9 t(π^2))).
Solution.  We will use  c = 1,  h = 0.2,  k = 0.04 .  This forces  This forces  r = (1^2) = (0.04/((((0.2))^2))) = 1.
First set up the boundary conditons:

\(Clear[c, f, g1, g2, h, i, k, m, n, r]; 
Clear[g1, g2, f]; 

    \(g1[j_] := 0.0; \)\)g2[j_] := 0.0; f[i_] := 
      N[Sin[π 0.1 (i - 1)] + Sin[3 π 0.1 (i - 1)]];

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[π x] (E^(-0.1  t (π^2))) +   
        Sin[3  π x] (E^(-0.1 9  t (π^2))), {t, 
        0, 1}, {x, 0, 1}, ViewPoint -> {4, 2, 3}];

Algorithm 10.4  Dirichlet Method for Laplace's Equation.  
To approximate the solution of  
(Subscript[u, xx]) (x, y) + (Subscript[u, yy]) (x, y) = 0    over    R = {(x, y) :  0 ≤ x ≤ a,   
          0 ≤ y ≤ b }  
with 
u (x, 0) = (Subscript[f, 1]) (x),   
      u (x, b) = (Subscript[f, 2]) (x)     for   0 ≤ x ≤ a   and,
u (0, y) = (Subscript[f, 3]) (y),   
      u (a, y) = (Subscript[f, 4]) (y)     for   0 ≤ y ≤ b.  
It is assumed that  Δx = h,  Δy = h  and that integers  n  and  m  exist so that  a = n h  and  b = m h.
Section  10.3   Elliptic Equations   Page  531

First, execute the following group of cells:

Execute this cell to activate Dirichlet Method for Laplace's Equation.  
Dirichlet[q_List,n_,m_,w_]
Neumann[q_List,n_,m_,w_]
q_List are the initial starting values,
w_ is a coefficient to be supplied,
n_ is the number of subintervals along the x-axis,
m_ is the number of subintervals along the y-axis.

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 = 
                      \(1\/4.0\) 
                        \((w\ \((
                              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 = 
                      \(1\/4.0\) 
                        \((w\ \((
                              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 = 
                  \(1\/4.0\) 
                    \((w\ \((
                          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[i, j, u, w]; 
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]] = 
      (1/2) (u[[1, 2]] + 
            u[[2, 1]]); u[[1, 9]] = 
      (1/2) (u[[1, 8]] + 
            u[[2, 9]]); u[[9, 1]] = 
      (1/2) (u[[8, 1]] + 
            u[[9, 2]]); u[[9, 9]] = 
      (1/2) (u[[8, 9]] + 
            u[[9, 8]]);

Next, solve it.

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

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

To see the numerical values enter the command:

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:    (Subscript[u, x]) = 0
Right:             0
First set up the boundary values.

\(Clear[i, j, u, w]; 
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]] = 
      (1/2) (u[[1, 2]] + 
            u[[2, 1]]); u[[1, 9]] = 
      (1/2) (u[[1, 8]] + 
            u[[2, 9]]); u[[9, 1]] = 80; u[[9, 9]] = 0; For[j = 1, j ≤ 9, j++, 
      u[[9, j]] = 
        u[[9, 1]] + 
          (1/8) (j - 1) 
            (u[[9, 9]] - 
                u[[9, 1]]); ];

Next, solve it.

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

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

To see the numerical values enter the command:

TableForm[N[u, 4]]