MathGroup Archive 2003

[Date Index] [Thread Index] [Author Index]

Search the Archive

Re: guidance in mathematica programming in pde

  • To: mathgroup at
  • Subject: [mg43709] Re: guidance in mathematica programming in pde
  • From: Rob Knapp <rknapp at>
  • Date: Tue, 30 Sep 2003 16:43:14 -0400 (EDT)
  • References: <bl0unq$3vv$>
  • Sender: owner-wri-mathgroup at

farkhanda wrote:
> respected sir i have just started the mathematica and facing lot of
> problems. during my search on internet I found your group email
> address. you repliled to one student like me about the problem in
> mathematica.I am sending the problem with my approach of
> programming,if possible point out me my mistake and at the same time
> if write the programm of this problem again,it will better for me in
> the sense that it give me new idea to think on the same line.
> waiting for a very positive response.thank you very much
> PROBLEM "solve the heat equation using explicit difference
> approximation with the help of programming in mathematica"
>  let u= u(x,t) and  heat equation  is 
> derivative of u(x,t) wrt t =double derivative of u(x,t) wrt x 
>  with the given boundary conditions  
> u(0, t) =0, u(1, t) = 0  for  t greater than and equal to 0
> and the initial condition    u(x,0)=x    when  x belongs to [0,1/2]
>                              u(x,0)=1-x  when  x belongs to [1/2,1]
> The explicit difference scheme is given as 
> [u(x,t+delta t)-u(x,t)]/(delta t) =[u(x+delta x,t)-2u(x,t)+u(x-delta
> x,t)]/(delta x )^2
> with  a condition   
>            (delta t) / (delta x )^2 is les than equal to (1/2)
> my requirement is I need the graph in 3 dimension of u, x and t.
> ClearAll[u, x, t];
> deltat = 0.01;
> deltax = 0.2;
> (*we check the stability condition :*)
> Print["stability test = ", (deltat)/(deltax)^2 <= 1/2];
> kmax = Floor[1/deltax];
> nmax = Floor[1/deltat];
> (*boundary conditions :*)
> u[k_Integer /; k <= 0 || k >= kmax, n_Integer] = 0;
> (*initial conditions :*)
> u[k_Integer /; 0 <= k <= kmax/2, 0] := k deltax;
> u[k_Integer /; kmax/2 <= k <= kmax, 0] := 1 - k deltax;
> (*difference scheme :*)
>  (u[k, n + 1] - u[k, n])/(deltat) == (u[k + 1, n] - 2u[k, n] + 
>             u[k - 1, n])/(deltax)^2;
> eq = 0.25 u[-1 + k, n] + 0.50 u[k, n] + 0.25 u[1 + k, n]
> u[k_Integer /; 0 <= k <= kmax, n_Integer /; n > 0] :=
>     u[k, n] = 
>       0.25 u[-1 + k, n - 1] + 0.50 u[k, n - 1] + 0.25 u[1 + k, n - 1];
> abc = Table[u[k, n], {k, 0, kmax}, {n, 0, nmax - 1}] // Transpose // 
>     MatrixForm
> ListPlot3D[abc]

You do not have a "mistake" per se in that your program works and 
produces a result with the settings you have provided.  However, it is 
very limited in that the approach does not scale well to finer 
discretizations (which give more accurate solutions or might be required 
for a more complex equation or initial condition.)  If you change your 
parameters to deltax = 0.01 and deltat = 0.0005, you will begin to see 
that your approach takes noticeable time: on my machine the computation 
takes over 7 minutes (435 seconds).

One of the advantages of Mathematica is that there are usually many ways 
to approach a problem.  (Though as a beginner, you may not see this as 
an advantage yet!)  I will show several approaches to solving the 
problem, starting with the simplest.

First, define the initial condition as a function.

In[1]:= finit[x_] := If[x < 1/2, x, 1 - x]

The simplest thing to do is to simply use a built in command when 
possible.  To use NDSolve, you just need to give the equation, initial 
condition, and boundary conditions.  The one extra thing that is needed 
here is to specify a MaxSteps.  This is needed because the initial 
condition is not smooth and without it, NDSolve would end up using a 
very fine grid in the x direction, which, at least on my machine ends up 
exhausing the memory I have available.  Such a fine grid is not 
necessary because of the nature of the heat equation, but the numerical 
method does not know that the equation is purely dissipative.

In[2]:= Timing[sol = NDSolve[{
           D[u[t, x], t] == D[u[t, x], x, x],
           u[0, x] == finit[x],
           u[t, 0] == 0, u[t, 1] == 0},
        u, {t, 0, 1}, {x, 0, 1},
        MaxSteps -> {10000, 100}]]

    Using maximum number of grid points 100
      allowed by the MaxPoints or MinStepSize options for independent 

Out[2]= {0.18 Second, {{u ->

 >       InterpolatingFunction[{{0., 1.}, {0., 1.}}, <>]}}}

One advantage of using NDSolve is that you get the solution in an easily 
useable form.  You can evaluate at any point in the interval, or plot it 
using, for example Plot3D[First[u[t,x] /. sol],{t,0,1},{x,0,1}, 

The above discretization uses fourth order differences in space and and 
adaptive time step with method appropriate for stiff equations, which 
allows time steps much greater than the deltat/deltax^2 stability 
restriction you needed for your approach.

However, it is quite possible that this comes from an assignment where 
it was expected that you use the "standard" centered second order 
discretization in space and explicit Euler stepping method in time. 
Using more options, NDSolve will do this too:

Set up parameters which will be stable for a grid of 100 points in the 
spatial direction.

In[3]:= nx = 100; deltax = 1/nx; nt = 20000; deltat = 1/nt;

This tells NDSolve to use second order differences in space with the 
explicit Euler method in time (works with a fixed step size equal to
the value given by the StartingStepSize option).  The error message is 
because 100 spatial points is insufficient to satisfy spatial error 
tolerances (no number of points would because the initial condition is 
not smooth).

In[4]:= Timing[NDSolve[{D[u[t, x], t] == D[u[t, x], x, x], u[0, x] == 
       u[t, 0] == 0, u[t, 1] == 0}, u, {t, 0, 1}, {x, 0, 1},
     Method -> {"MethodOfLines",
         Method -> "ExplicitEuler",
         "SpatialDiscretization" -> {"TensorProductGrid",
             "DifferenceOrder" -> 2,
             "MaxStepSize" -> deltax, "MinStepSize" -> deltax}},
     StartingStepSize -> {deltat, Automatic},
     MaxSteps -> Infinity]]

    Warning: Estimated initial error on the specified spatial grid in the
      direction of independent variable x exceeds prescribed error 

Out[4]= {1.52 Second, {{u ->

 >       InterpolatingFunction[{{0., 1.}, {0., 1.}}, <>]}}}

This takes a lot longer because it has to take 20000 time steps.  The 
previous solution took far fewer time steps because it was using a more 
appropriate time stepping method.

This may still be quite a bit too simple for what the assignment 
intended, since after all, NDSolve is still doing all of the 
discretization automatically.

If you were expected to "do" the discretization yourself, there are 
still many approaches to to it in Mathematica.  In general, if you can 
work with lists in Mathematica rather than functions or individual 
elements, the resulting code will be much faster.  Thus, the approach I 
suggest is to define the solution at any given time step to be a list, 
U.  The solution at the next time will be computed by a function f[U] 
which returns a list.

This sets up a grid X, an initial value U0, which is the list of values 
representing the initial condition, and a parameter for convenience. 
Note that the grid is chosen to include the endpoints, so that U0 (and 
all instances of the solution U) will have the first and last parts 
equal to zero to match the boundary conditions.

In[6]:= X = N[Range[0, nx]]*deltax;U0 = Map[finit, X]; lambda = 

Now all that is necessary is to define the function f[U].  I will show 
two ways to do this efficiently.

The first approach is to consider the finite difference operator as a 
list correlation.  The trickiest part of using ListCorrelate for 
applications like this is getting the boundaries right.  Here I choose 
the boundary parameters so that it wraps around, which is possible since 
the endpoints are always 0.

In[7]:= f[U_] := ListCorrelate[{lambda, 1 - 2 lambda, lambda}, U, {2,2}];

In[8]:= Timing[data = Join[{U = U0}, Table[U = f[U], {nt}]];]

Out[8]= {2.9 Second, Null}

Note that when you get the matrix of values, data, it is easy to 
construct an InterpolatingFunction which can be used for plotting, and 
evaluation by using ListInterpolation[data, {{0,1},{0,1}}]

Another approach to handling finite differences is to use the fact that 
derivatives, and so finite differences are linear, and thus 
representable by a matrix (called a differentiation matrix).   Since the 
matrix is very sparse, I will use a SparseArray to represent it.  Also, 
because the equation is linear, I can include the whole action of the 
linear operator in the matrix, which saves some computation time.

Only the first few rows and columns are shown for brevity.

In[9]:= MatrixForm[Take[
         A = SparseArray[{
            {i_, i_} /; (1 < i < nx + 1) -> 1 - 2 lambda,
            {i_, j_} /; ((1 < i < nx + 1) && Abs[i - j] == 1) -> lambda},
           {nx + 1, nx + 1}, 0.],
         5, 5]]

Out[9]//MatrixForm= 0.    0.    0.    0.    0.

                     0.5   0.    0.5   0.    0.

                     0.    0.5   0.    0.5   0.

                     0.    0.    0.5   0.    0.5

                     0.    0.    0.    0.5   0.

Then we can define

In[10]:= f[U_] := A.U

In[11]:= Timing[data = Join[{U = U0}, Table[U = f[U], {nt}]];]

Out[11]= {0.96 Second, Null}

Computing the same thing somewhat faster.

Since the matrix above takes some thought to construct, I will show a 
final way of defining f, that, while not as fast, generalizes very 
easily to other equations.

The command NDSolve`FiniteDifferenceDerivative is described in the 
advanced documentation for NDSolve is Mathematica 5 (open the 
documentation for NDSolve and look for Advanced Documentation near the 
bottom of the page).  It produces a function which will approximate the 
given derivative for a function with values given at the given grid 
points using a specified difference order.

In[12]:= fddf = NDSolve`FiniteDifferenceDerivative[Derivative[2], X, 
DifferenceOrder -> 2]

Out[12]= NDSolve`FiniteDifferenceDerivativeFunction[Derivative[2], <>]

This just approximates the spatial derivative, so the update is Euler's 

In[13]:= f[U_] := U + deltat fddf[U]

In[14]:= Timing[data = Join[{U = U0}, Table[U = f[U], {nt}]];]

Out[14]= {1.91 Second, Null}

This takes about twice as long because of the extra scalar 
multiplication and addition required.  However, the advantage is that it 
is very easy to change the difference order, and it is quite easy to 
adapt the method to different (including nonlinear) equations.

All of this is probably quite a bit more than you were looking for, but 
I hope that you will find it useful nonetheless.

Rob Knapp
Wolfram Research, Inc.

  • Prev by Date: Bug 2?
  • Previous by thread: Re: guidance in mathematica programming in pde
  • Next by thread: ScientificForm in Ticks