Re: Mathematica NDSolve max number of grid points error
- To: mathgroup at smc.vnet.net
- Subject: [mg22952] Re: Mathematica NDSolve max number of grid points error
- From: Robert Knapp <rknapp at wolfram.com>
- Date: Fri, 7 Apr 2000 02:54:34 -0400 (EDT)
- Organization: Wolfram Research, Inc.
- References: <8chc3f$9fv@smc.vnet.net>
- Sender: owner-wri-mathgroup at wolfram.com
Scott Beckman wrote:
>
> Hello,
>
> I'm trying to solve a simple PDE with Mathematica's NDSolve function,
> but I'm getting the error "Using maximum number of grid points 1500
> allowed by the MaxSteps option.". I've used the MaxSteps function to
> increase the number of steps up to the point that I've consumed all the
> available memory. This isn't the right way to solve the problem.
>
> Has anyone else had this problem and is there a way around it? Are
> there books on using Mathematica for numerical problems such as this?
>
> Below are the lines I'm having problems with. Its a simple little
> program to calculate the heat conduction in x and t for one fixed T
> boundary condition and one insulated boundary condition and an initial
> temperature distribution that is constant.
>
> >alpha = 1.32*10^-5
>
> >solution =
> NDSolve[{D[u[x, t], t] == alpha D[u[x, t], {x, 2}],
> u[x, 0] == (UnitStep[(x - 0.000000001)]*(423 - 455) +
> 455), ((D[u[x, t], x]) /. x -> 0.00005) == 0, u[0, t] ==
> 455},
> u, {x, 0, 0.00005}, {t, 0, 0.00003}, MaxSteps -> 1500];
>
> >interf[tt_] := Evaluate[u[0.00005, tt] /. Last[solution]]
>
> >FindRoot[interf[tim] == 428, {tim, 0.00001}]
>
The problem NDSolve is running to is that the initial temperature
distribution is not quite constant: it has a discontinuity at x ==
0.000000001. NDSolve uses error estimates to try to decide how many
spatial points to use. The error estimates need to assume that the
underlying functions are smooth. Since NDSolve is using a purely
numerical solution method it does not know that the dissipative nature
of the heat equation will smooth out intial errors quite rapidly.
Because of this, the solution computed with the default MaxSteps will
likely be plenty accurate despite the warning message. In fact, that
many spatial points is excessive for the problem. Here are a couple of
ways you can do it faster:
alpha = 1.32*10^-5;
xmax = 0.00005;
tmax = 0.00003;
solution =
NDSolve[{D[u[x, t], t] == alpha D[u[x, t], {x, 2}],
u[x, 0] == (UnitStep[(x - 0.000000001)]*(423 - 455) +
455), ((D[u[x, t], x]) /. x -> xmax) == 0, u[0, t] ==
455},
u, {x, 0, xmax}, {t, 0, tmax},
StartingStepSize -> {xmax/50, Automatic}];
Note that I have given the StartingStepSize options as a list, one
component for the x and the Automatic for t.
or it is possible to use a nonuniform grid to better capture the
dissipation of the initial discontinuity:
alpha = 1.32*10^-5;
xmax = 0.00005;
tmax = 0.00003;
grid = Join[{x}, Range[0., xmax/5 - xmax/200, xmax/200],
Range[xmax/5, xmax, xmax/40]];
solution =
NDSolve[{D[u[x, t], t] == alpha D[u[x, t], {x, 2}],
u[x, 0] == (UnitStep[(x - 0.000000001)]*(423 - 455) +
455), ((D[u[x, t], x]) /. x -> xmax) == 0, u[0, t] ==
455}, u,
grid, {t, 0, tmax}];
In both solutions, I find the solution doesn't quite get to 428 at the
edge, so I would suggest you solve to a bigger time t before using
FindRoot.
Rob Knapp
Wolfram Research, Inc.