Re: guidance in mathematica programming in pde
- To: mathgroup at smc.vnet.net
- Subject: [mg43709] Re: guidance in mathematica programming in pde
- From: Rob Knapp <rknapp at wolfram.com>
- Date: Tue, 30 Sep 2003 16:43:14 -0400 (EDT)
- References: <bl0unq$3vv$1@smc.vnet.net>
- Sender: owner-wri-mathgroup at wolfram.com
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.
>
> MY APPROACH OF PROGRAMING
> 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}]]
NDSolve::mxsst:
Using maximum number of grid points 100
allowed by the MaxPoints or MinStepSize options for independent
variable
x.
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},
PlotRange->All]
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] ==
finit[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]]
NDSolve::eerri:
Warning: Estimated initial error on the specified spatial grid in the
direction of independent variable x exceeds prescribed error
tolerance.
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 =
deltat/deltax^2;
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
method:
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.