I am experiencing trouble when trying to solve a PDE for various parameter values within mathematica: The PDE in principle can easily be solved numerically, however the stepsize of the spatial discretization (method of lines) has to be altered depending on the choice of parameters for NDSolve to converge. The code is as follows,

Sol[a_, b_] := Sol[a, b] = NDSolve[{PDE[a, b], BC}, P, {x, xmin, xmax}, {\[Phi], 0, 2*Pi}, {t, 0, tmax}, MaxStepFraction -> {1/200, 1/20, 1/10}, Method -> {"BDF", "MaxDifferenceOrder" -> 2}, EvaluationMonitor :> (CurrTime = t;), InterpolationOrder -> All]

where $P(x,\phi,t)$ is the unknown function entering the PDE. Unfortunately, the xStepFraction of $\frac{1}{200}$ is not small enough for some values of $a$ and $b$, for some other values it must however not be chosen too small (otherwise NDSolve does not converge). I know that up to now mathematica does not provide adaptive grid refinement for the method of lines. In addition, the right xStepFraction does depend rather arbitrarily on $a$ and $b$, so I cannot define a function $\text{xStepFraction}(a,b)$.
So, my question is, whether it is possible to construct some kind of grid refinement, say to run NDSolve with a rather coarse grid spacing, and if NDSolve returns an error message like

NDSolve::eerr: Warning: scaled local spatial error estimate of 220.4370802734196` at t = 0.1646666008001777` in the direction of independent variable x is much greater than the prescribed error tolerance. Grid spacing with 201 points may be too large to achieve the desired accuracy or precision. A singularity may have formed or a smaller grid spacing can be specified using the MaxStepSize or MinPoints method options. >>

to automatically reduce the xStepFraction. Any help would be highly appreciated.

  • $\begingroup$ Have you seen reference.wolfram.com/language/tutorial/… ? $\endgroup$ – Dr. belisarius Feb 11 '15 at 17:16
  • $\begingroup$ Yes, but that gave me no answers to my problem. $\endgroup$ – Alex Feb 11 '15 at 17:26
  • $\begingroup$ I think it'll be better if you can add the specific equations to your question, maybe someone can find a combination of options that'll circumvent your problem. $\endgroup$ – xzczd Feb 12 '15 at 5:46

Here's one way, using Burgers' equation, to illustrate the general principle. One can use Check to react to any given message. Here we reduce the step size by 1/2. Another factor may be more appropriate in another case.

PDE[a_] := {D[u[x, t], t] == a D[u[x, t], x, x] - u[x, t] D[u[x, t], x]};
BC = {u[0, t] == u[1, t], u[x, 0] == Sin[2 Pi x]};

sol::restart = "Restart integration with step size ``";
sol[a_, stepFraction_] := 
 Check[NDSolve[{PDE[a], BC}, u, {x, 0, 1}, {t, 0, 0.5}, 
   MaxStepFraction -> {stepFraction, 1/10}, 
   Method -> {"BDF", "MaxDifferenceOrder" -> 2}, 
   EvaluationMonitor :> (CurrTime = t;), InterpolationOrder -> All],
  Message[sol::restart, stepFraction/2];
  sol[a, sol["stepFraction"] = stepFraction/2],

{usol} = sol[0.01, 1/50]
(* messages omitted *)

One can check the final max step fraction with

(*  1/200  *)

Plot3D[u[x, t] /. usol, {x, 0, 1}, {t, 0, 0.5}, AxesLabel -> Automatic]

Mathematica graphics

One probably ought to use a variable to limit the number of reductions of the step size, or set an absolute minimum step size at which the recursive calls to sol would stop.

  • $\begingroup$ Thank you very much, that was exactly what I was looking for. $\endgroup$ – Alex Feb 12 '15 at 10:57
  • $\begingroup$ @Alex You're welcome. If it ever seems unsatisfactory, I would echo xzczd's comment that there may be better approaches depending on the actual PDE. $\endgroup$ – Michael E2 Feb 12 '15 at 14:03

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.