A Heat Flow Problem

Discussion of the Problem

Our second example has its origins in a physical problem of heat flow. 
Suppose we have a rod of material of unit length, say, situated on the
x-axis, for 0 <=  x <= 1$.  Suppose further that the rod is laterally
insulated, but has a known internal heat source.  Also suppose that
both the left and right ends of the rod are held at 0 degrees Celcius. 
When sufficient time passes, the temperature of the rod at each point
will settle down to ``steady state'' values, dependent only on position
x.  It can be shown from physical principles that if  u  is the
temperature function of the rod, i.e.,  u(x)   is the temperature of the
rod at   x  , then   u(x)   and   f(x)   are related by the following equations:

	-u''(x) = f(x),  
	  u(0)   =   0,  
          u(1)   =   0.
	
Here ' denotes derivatives with respect to   x   and   f(x)  is a known
function which depends on the material and the internal heat
distribution of the rod.  The type of problem just described is called
a   boundary value problem  and abbreviated to BVP.  Note that if we
knew how to integrate   f(x)   twice, it would be easy to determine
 u(x)  : simply use the boundary conditions at 0 and 1 to determine the
two constants of integration.  The interesting case is when   f(x)   is
not easily integrated.  In fact, a closed formula for  f(x)   may not
even be given.  How do we solve the problem then?

We can formulate a discrete approximation to the BVP above as follows. 
Choose an integer  n  and divide the interval [0, 1] into   n+1   equal
subintervals with endpoints  0 = x_0, x_1, ...  ,x_{n+1} = 1   (these
points are called  nodes  ).  Then the width of each subinterval is
  h = 1/(n+1)   .  Next let  u_i    be our approximation to  u(x_i) . 
Observe that we can approximate the first derivative by the forward difference formula 
	
(1)   	u'(x) ~(u(x+h) - u(x)) / h .

Likewise, we can approximate the second derivative by the
backward difference formula 

(2)   	u''(x)  ~ (u'(x-h) -u'(x)) / (-h) 

Now use  the formula (1)  with   x   and  x-h   replacing  x  and substitute into Formula  (2)  to obtain the second derivative formula

(3)   	u''(x) ~ (u(x+h) - 2u(x) + u(x-h)) / (h^2) 

At each unknown node   x_i   replace the second derivative
in the equation   -u''(x_i) = f(x_i) ,   i = 1,...,n  by the approximation
given by Formula (3) ; then multiply the equation by  h^2. 

Observe that   u(x_i  +  h) = u_{i + 1}   and one obtains a system of
  n   equations in the   n   unknowns   u_1, u_2,...,u_n,   where the ith
equation takes the form 
	
	 -u_{i-1} + 2u_i - u_{i+1} = h^2f(i h)   
	 
 Remember that   u_0   and   u_{n+1}   are known to be 0, so the terms   u_0  
and   u_{n+1}   disappear.

It is reasonable to expect that the smaller   h   is, the more accurately
 u_i     will approximate   u(x_i)  .  This is indeed the case.  But consider 
what we are confronted with when we take  n = 5  , i.e.,  h = 1/6 , 
which is hardly a small value of $h$.  The system of five equations in
five unknowns becomes


2x_1 -  x_2                                          =   f(1/6)   
 -x_1 + 2x_2  -  x_3                            =   f(2/6)  
             -x_2  + 2x_3 -  x_4                =   f(3/6) 
                          -x_3  + 2x_4  -  x_5  =  f(4/6)  
                                     -x_4  + 2x_5  =   f(5/6) 
 


This problem is already about as large as we would want to work by hand.  

Now let's try a specific example.  Use N[Sin[Pi*x]] for f[x].  Define the required equations yourself and use Mathematica to compute the approximate solution by filling in the  following cells with your equations.  Use x1, x2, x3, x4, and x5 for your variables.  BE SURE TO USE DECIMALS SOMEWHERE IN YOUR EQUATIONS TO FORCE MATHEMATICA TO VIEW EVERYTHING AS NUMBERS.  The last cell will show you a graph of your approximate heat function u(x).

vars = {x1, x2, x3, x4, x5}

equations = {    , 
   ,    
   , 
   ,
        }

rules = Solve[equations,vars]

ListPlot[{{0.,0},{N[1/6],x1},{N[2/6],x2},{N[3/6],x3},
{N[4/6],x4},{N[5/6],x5},{1.,0.}}/.rules [[1]], PlotJoined->True]

Exercises