Ordinary and Partial Differential Equations

In the previous part I gave a Mathematica version of the material in Dick Clements article [Clements 1995]. The present article deals similarly with his subsequent papers [Clements 1995a] and [Clements 1995b]. I have found this a most informative and helpful exercise both because of the interesting material that he presents and because it has forced me to clarify my ideas about how to use Mathematica.

The main difficulty I have found has been, not with how to use system for doing mathematics, but with how and what to communicate with it — the educational problem.  This is as it should be: we need to see through the technology to the mathematics and the message

A proper mathematical education must involve students and teachers in routinely using the proper tools for the job. I hope that these exercises will help in deciding what these are.

As to style, I suggest that we might hide quit a lot of the routine input and programming so as to more clearly expose the essential material; but it should still be accessible (this is easily achieved on-screen using close cells and groups of cells). Of course, if we are teaching Mathematica then we will show more.

Ordinary Differential Equations

Here is an ordinary differential equation

A solution of it over a range [a,b] of t is a function, x, of t, such that

for all t in [a,b].
Replacing dx(t)/dt by its approximation (x(t+h) - x(t))/h  (for small h) we get

Which can be rearranged to give

x(a+h) = x(t) + hf(t,x(t))						(1)

The hope is that if, starting with a particular x[a], we use (1) to construct

x(a), x(a+h), x(a+2h), . . . , x(a+nh)

then this will be close to the sequence of values at a, a+h, a+2h ... of a solution of the equation.

Let's see how our idea works (from now on I shall usually use the Mathematica forms  f[x] ... in place of f(x)....)

The function step[h] below gives  {t+h, x+h f[t,x]} from {t,x}.

step[h_][{t_, x_}] := {t + h, x + h f[t, x]}

Let

f[t_, x_] := ((2 (x^2))/(1 + 2 t))

and try to find a solution over [0,3]  with x[0] = 1/2.

Using h = 3/10 we get

sol10 = NestList[step[(3.0/10)], {0, 0.5}, 10]

We expect to do better with smaller h's. Here are the pictures for h = 3/10, 3/20,... ,3/90.

These seem to be approaching a limit, hopefully a solution. In this case we can confirm that this is so because our differential equation is exactly solvable.

dsol=DSolve[{x'[t]==f[t, x[t]],x[0] == .5},x[t],t][[1, 1, 2]]

We can fill in the gaps for the last approximation by using an interpolating function.

sol90 = Interpolation[NestList[step[(3/90)], {0, 0.5}, 90]];

Plotting this together with dsol we get

This is not very good after  t=2.5: let's try h = 3/450.

sol450 = Interpolation[NestList[step[(3/450)], {0, 0.5}, 450]];

In fact, theory predicts that the errorof an approximate solution, soln, defined by
			error[soln][t]  = soln[t] - dsol[t], 
should be proportional to h. So for sol90 it should be about five times that for sol450. Let's see.

Unfortunately, as shown below, error[sol450] still increases rapidly as the solution bends upwards:

Of course, we expect this, since the steps follow the tangent rather than the curve:

An obvious check is to substitute the approximate solution in the equation and find the difference between the two sides. The interpolating function, which is differentiable, is crucial for this. For sol450 at t = 2 we get

equationerror[sol450][2]=f[t, sol450[t]]-sol450'[t]/.t -> 2

The plot of equationerror[sol450]over [0,3] is

Once more we see the bad effect of the upward bending.

The method that we have been using is called Euler's method; you may have some ideas about how it can be  improved, for example how to deal with the problem of the bend.
But let's see how NDSolve, the built-in numerical differential equation solver, measures up.

ndsol=NDSolve[{x'[t]==f[t, x[t]],x[0] == .5},x,{t, 0, 3}][[1, 1, 2
          ]];

In fact, NDSolve is a very clever program. To get an understanding of the way that it works you should start by finding out about the Runge-Kutta methods (see Appendix) An account of the specific techniques used can be found in [Hindmarsh 1983] and [Petzold 1983].

I need to clear  f of its old definitions before using it further.

Clear[f]

Partial Differential Equations

Steady temperature in a plate

The first two problems are to find steady state temperature distributions in a uniform plate; time is not involved.

Given temperatures on the edges

A uniform square plate has temperature  20 on three sides and 60 along the fourth side. Find the temperatures at points inside the plate.

Take an (x,y) coordinate system in the plate, as shown below, and let the temperature at (x,y) be  T(x,y).

T(x,y) must satisfy the Laplace equation.

However, as before, instead of trying to solve this by calculus we will turn it into a related problem in simple arithmetic — and then get the computer to do the arithmetic.

For convenience, define a function for getting the standard approximation to the derivative of an expression expr with respect to a variable x for a small change h in x.

s[expr_, x_, h_] := (((expr /. x -> x + h) - expr)/h)

For example

s[f[x, y], x, h]

Using this twice gives an approximation to the second derivative with respect to x

Simplify[s[s[f[x, y], x, h], x, -h]]

From this we get the following approximation to Laplace's equation for  T[x,y].

Simplify[
      s[s[T[x, y], x, h], x, -h] + s[s[T[x, y], y, h], y, -h] == 0]

Which is equivalent to

"T[x,y]" == ("T[x,y-h] + T[x,y+h] + T[x-h,y] + T[x+h,y]"/4)

This is amazingly simple: it just says that  T at (x,y) is the average of T at the four points at distance h,either side, above and below (x,y). 
Here's how I use it.

Take an 8 by 8 uniform array of points on the plate .

I'm going to try to find temperatures at these points which are as given in the problem along the edges and which, at each inside point are the average of the temperatures at the point's four neighbours (above, below,to the right and to the left). Thats all! I hope that these temperatures will be close to those of a solution to the problem.

I start with the following array of temperatures.

temp = Table[Which[n == 1, 60.0, True, 20.0], {m, 8}, {n, 8}]

This doesn't satisfy the condition; so I'll try to make it satisfy the condition by the brute force method of replacing values by what I want them to be. The catch is, that each replacement can spoil what I have already achieved. Still, let's see what happens.

Define a function step to make the replacement for each pair {m,n} with  m = 2,3,...7 and n = 2, 3, ...7. (this leaves the values on the edges unchanged).

step[temp_] := 
      Block[{T = temp}, 
        Do[T[[m, n]] = 
            (1/4) (
                T[[m - 1, n]] + 
                  T[[m, n - 1]] + 
                  T[[m, n + 1]] + 
                  T[[m + 1, n]]), {m, 
            2, 7}, {n, 2, 7}]; T]

The list tempslist8 ,below, is temp followed by the results of 8 successive applications of step. I show the  last member.

Last[tempslist8 = NestList[step, temp, 8]]

Let's look at what these tables give for the lines of temperatures 22, 25, 30, 35, 40, 45, 50, 55. We get Mathematica to estimate where these lines will run. Here is the sequence of plots (it is readily animated on the computer).

These seem to be approaching a stable state , hopefully a solution. But maybe we should keep going rather longer. We use the function FixedPointList and tell it to stop when the the last application of step has changed the temperatures by less than .0001.

ans = 
      FixedPoint[step, temp, 
        SameTest -> (Max[Abs[#1 - #2]] < .0001&)];

Here is the plot of ans for temperatures 22,25,30,35,40,45,50,55.

To find the temperature at intermediate points I interpolate;

ansf = 
      Interpolation[
        Flatten[Table[{m, n, 
              ans[[m, n]]}, {m, 8}, {n, 
              8}], 1]];

But remember that m refers to the row and n to the column. In addition, if we want to use the (x,y) coordinate system shown at the beginning then a bit of scaling is needed. For example, if the plate has side length 5 and there is no scaling then the appropriate function is

ansfxy[x_, y_] = ansf[8 - ((7 y)/5), 1 + ((7 x)/5)];

The plot obtained from this, shown below, is not quite the same as the one that we found above. This opens up another fascinating question for investigation — what  is a good interpolation?

The interpolation that we have used tries to fit a smooth function. This was useful when we needed to substitute into a differential equation, but it may not be best here; for example it can lead to overshooting, as shown below.

int = 
      Interpolation[{{0, 0}, {1, 0}, {2, 0}, {3, 0}, {4, 1}, {5, 1}, {6, 1}, {
            7, 1}}];

Plot[int[x], {x, 0, 7}];

It may be better to use linear interpolation

int1 = 
      Interpolation[{{0, 0}, {1, 0}, {2, 0}, {3, 0}, {4, 1}, {5, 1}, {6, 1}, {
            7, 1}}, InterpolationOrder -> 1];

Plot[int1[x], {x, 0, 7}];

Let's try this out on our plate temperature problem:

ansf1 = 
      Interpolation[
        Flatten[Table[{m, n, 
              ans[[m, n]]}, {m, 8}, {n, 
              8}], 1], InterpolationOrder -> 1]; ansf1xy[x_, y_] = ansf1[8 - ((7 y)/5), 1 + ((7 x)/5)];

Which gives the plot

One edge insulated

Modify the previous problem: instead of the right hand edge having temperature 20 , have it insulated, so that  no heat flows across  it.

The new condition can be expressed by the following equation along the insulated edge.

The approximation for this is

Simplify[s[T[x, y], x, -h] == 0]

Which can be rearranged to 
				 T[x,y] = T[x-h,y]                   (2)

Now to find a solution.
I start with the same points and temperature array ,temp, as before.

Condition (2) translates into temp[[8,n]]= temp[[7,n]]  for  n = 1,2,...,8; this is already satisfied by temp.

The appropriate step to create a new array of temperatures is now

step2[temp_] := 
      Block[{T = temp}, 
        Do[T[[m, n]] = 
            (1/4) (
                T[[m - 1, n]] + 
                  T[[m + 1, n]] + 
                  T[[m, n - 1]] + 
                  T[[m, n + 1]]), {m, 
            2, 7}, {n, 2, 7}]; 
        Do[T[[m, 8]] = 
            T[[m, 7]], {m, 1, 8}]; T]

Using this and the function FixedPoint as before we get the approximate solution

ans2 = 
      FixedPoint[step2, temp, 
        SameTest -> (Max[Abs[#1 - #2]] < .0001&)]

Here is the plot of ans2 for temperatures 22,25,30,35,40,45,50,55.

Diffusion of temperature along a rod

The next problem is to find how the temperature distribution in a rod changes with time.

A uniform rod begins at temperature 20 everywhere. It is kept at temperature 0 at one end and insulated at the other end. Find how the temperature at each point varies with time.

Let x be the distance along the rod from the end kept at temperature zero, and let the temperature at x at time t be T(x,t).

The differential equation systemfor this problem is

The approximate version is

Simplify[s[s[T[x, t], x, dx], x, -dx] == a s[T[x, t], t, dt]]

s[T[x, t], x, -dx] == 0

which, with  b = dt/(a dx2), is equivalent to

T[x, t + dt] == 
      T[x, t] + b (T[x - dx, t] - 2 T[x, t] + T[x + dx, y]) 3

T[x,t] == T[x-dx,t]					                            (4)

The method of solution is similar to that already used.

I take 8 equally spaced points along the rod and a list, temp3 of temperatures at these points.

temp3 = Table[Which[x == 1, 0., True, 20.], {x, 1, 8}]

Following conditions (3) and (4), I use the function step3, defined below, to construct a sequence of lists that I hope will be close to the lists of temperatures at the chosen points as time increases by increments  dt.

step3[temp_] := 
      Block[{newtemp = temp}, 
        Do[newtemp[[m]] = 
            temp[[m]] + 
              b (temp[[m - 1]] - 
                    2 temp[[m]] + 
                    temp[[m + 1]]), {
            m, 2, 7}]; 
        newtemp[[8]] = 
          newtemp[[7]]; newtemp]

(Notice that temp not newtemp is used on the right of the assignment 
	newtemp[[m]] = temp[[m]]+b(temp[[m-1]]-2 temp[[m]]+temp[[m+1]])
)

Before calculating I need to give b a numerical value.

b = 0.4;

I apply step3 50 times to generate ans3 below – the first eight members are shown.

(Take[#1, 8]&)[ans3 = NestList[step3, temp3, 50]]

The progressive cooling is vividly show by the plots lines of equal temperatures

We can also plot the temperature along the bar for various times. As before I use an interpolation function and scale the entries. Assume that the bar has length 6 and that each step gives the change over 0.1 minute. The appropriate function is found by

pf = 
      Flatten[Table[{(1/7) (m - 1) 6, (n - 1) .1, 
            ans3[[n, m]]}, {m, 8}, {n, 
            51}], 1]; ans3fsc = Interpolation[pf];

Here are the temperatures given by ans3fsc along the rod at times 1,2,3,4 and 5 minutes.

Appendix

Runge-Kutta

The standard second order Runge-Kutta scheme  with constant step size h for a a differential equation
	dy/dx = f(x,y)
 is to take the next point after {x0,y0} to be {x1,y1}
 	{x1,y1} = step[h,f] [{x0,y0}] 
 where

step[f_, h_][{x_, y_}] := 
      (k1 = f[x, y]; 
        k2 = f[x + h, y + h k1]; {x + h, y + (1/2) h (k1 + k2)});

But this is only the beginning.

References

Clements R R,  Spreadsheets as a Tool for Teaching Elementary Numerical 
	Analysis,  Maths&Stats Vol 6, No 1, 2-7, (1995).
Clements R R, Teaching the numerical solution of ODEs using spreadsheets, 	     
          Maths&Stats Vol 6, No 2, 10-14, (1995a).
Clements R R, Illustrating finite difference solutions of PDEs using spreadsheets,   
 	 Maths&Stats Vol 6, No3, 8-12, (1995b).
Hindmarsh A C, "ODEPACK: A Systematized Collection of ODE Solvers", in 	
	Scientific Computing, R. S. Stepleman et. al. (eds). North-Holland, Amsterdam 
	(1983), pp. 55-64.
Petzold  L R,  "Automatic Selection of Methods for Solving
    	 Stiff and Nonstiff Systems of Ordinary Differential Equations", SIAM Journal of 	 
    	 Scientific and Statistical Computing, 4 (1983),pp. 136-148.
MathSource, http://www.wolfram.com/mathsource/