Message Boards Message Boards

Which parameters of NDSolve to change to get a more accurate PDE solution?

Posted 5 months ago

Hello, I have been trying to solve PDEs using NDSolve, in particular the following one:

xmax=10;
wmax=10;
xmin=10^-20;
psiinit[w_,x_]:=(Exp[-w^2]/Sqrt[π]-w*Erfc[w])*x;
sol=NDSolve[{2*x*D[psi[w,x],x]==D[psi[w,x],{w,2}]+2*w*D[psi[w,x],w]-(x*psi[w,x])^2,psi[w,xmin]==psiinit[w,xmin],Derivative[1,0][psi][0,x]==-x,psi[wmax,x]==0},psi,{w,0,wmax},{x,xmin,xmax},PrecisionGoal->20,AccuracyGoal->10,Method->{"MethodOfLines","SpatialDiscretization"->{"TensorProductGrid","MaxStepSize"->1/5000,"DifferenceOrder"->4}}]
N[psi[0,7]/.sol,20]
{1.14428}

However, as can be seen at the end of the above code, the solution obtained is provided with only 6 significant digits, whereas I would be happy to obtain more accurate solutions. Changing parameters such as PrecisionGoal and AccuracyGoal does not seem to help, but cause some error messages. I managed to obtain some improvement only by increasing MaxStepSize. It is less clear whether changing DifferenceOrder brings any improvement. Which other parameters one can change to get a better accuracy of the solution, without causing error messages? Manuals do not seem to provide any discussion of this essential issue. Manuals are also incomplete, as they do not explain which are default settings for various parameters, so that one cannnot even know which particular method or available option is used. Ideally, I would like to make use of the multiprecision available in MATHEMATICA, in order to obtain solutions with as many as 19 accurate significant digits. Is this possible at all? If yes, I would also be happy to learn how to split the PDE solving calculations into several runs, as they are likely to be time-consuming. Or, perhaps there are more accurate and efficient alternatives to NDSolve?
Leslaw

POSTED BY: Leslaw Bieniasz
9 Replies
Posted 5 months ago

There are lists and descriptions of methods, including implicit one-step methods:

LSODA is the default solver for ODE's. IDA for DAEs. Not sure about other forms of DEs.

According to the Overview, "FixedStep" is restricted to one-step methods, so you can't use BDF or LSODA with it. A list of such submethods is in the docs. for NDSolve. Strangely, the FixedStep section gives an example that claims to show how to use Method -> "BDF" using fixed steps, but the steps are in fact of different sizes. You can use BDF for time integration in a PDE, but I do not know how to fix the step size or to restrict the minimum step size.

I do not know of any implementation of Rosebrock methods in Mathematica/WL. How to implement your own methods is described in the Plugin Framework.

POSTED BY: Updating Name
Posted 5 months ago

Responding to the person "Updating Name": I tried your version of the code, but I get a message: "The initial conditions did not evaluate to an array of numbers of depth 1 on the spatial grid. Initial conditions for partial differential equations should be specified as scalar functions of the spatial variables". To make your code work I had to remove WorkingPrecision -> 40, but then I get an accuracy of no more than 4-11 digits, depending on the choice of MaxPoints and MinPoints, while increasing AccuracyGoal or PrecisionGoal resulted in an enormous slowing down. My feeling is that NDSolve does not tolerate multiprecision, at least not in adaptive calculations.

Concerning the pseudospectral discretization: isn't it practical for periodic boundary conditions only?

Responding to Michael Rogers: I'll try your suggestions, but I am not happy about "ExplicitRungeKutta". Is there any possibility to employ implicit methods, say BDF? Or, other robust implicit schemes, like for example Rosenbrock methods? Some of these are very handy for nonlinear problems, as they do not require Newton iterations. I noticed in some examples on the internet that LSODA was used in NDSolve, and LSODA is based on BDF, if I remember. Explicit methods are likely to cause stability problems. Unfortunately, I cannot find anywhere a complete list of available integrators and other options for NDSolve. Leslaw

POSTED BY: Leslaw Bieniasz
Posted 5 months ago

Suppose that I don't want to use automatic features of NDSolve, but just prefer fixed, uniform spatial and temporal grids, combined with my choice of the spatial and temporal discretization orders, plus my choice of WorkingPrecision. How to do this? I suspect I might choose SpatialDiscretization->MaxSteps == SpatialDiscretization->MinSteps == (my value) to eliminate automatic spatial grid selection. Is this correct? But how to set the other options so that PrecisionGoal and AccuracyGoal are effectively not used by NDSolve? Given the fragility of NDSolve adaptive calculations, the above approach might be the best one, if only no criptic error messages occur; I could test the accuracy of the procedure by comparing the actual errors with some known solutions, instead of relying on the adaptive algorithm. Leslaw

POSTED BY: Leslaw Bieniasz

Used "FixedStep" + StartingStepSize for a fixed time step and whichever method for the spatial grid. Ignore the error/warning message about StartingStepSize. You've already complained about error messages, and I have no defense for this one. (It seems StartingStepSize may be used for both temporal and spatial discretization, and it's okay for one but not the other in the examples below.)

Example 1:

\[Nu] = 0.01;
mysol = First[
  NDSolve[{D[u[x, t], t] == \[Nu] D[u[x, t], x, x] - 
      u[x, t] D[u[x, t], x], u[x, 0] == E^-x^2, u[-5, t] == u[5, t]}
   , u, {x, -5, 5}, {t, 0, 4}
   , Method -> {"MethodOfLines", 
     "SpatialDiscretization" -> {"TensorProductGrid", 
       "MinPoints" -> 257, "MaxPoints" -> 257,
        "DifferenceOrder" -> 4},
     Method -> {"FixedStep", Method -> "ExplicitRungeKutta"}}, 
   StartingStepSize -> 1/2^13]]

Example 2:

mygrid =  Join[-5. + 10  Range[0, 48]/80, 1. + Range[1, 4  70]/70];
\[Nu] = 0.01;
mysol = First[
  NDSolve[{D[u[x, t], t] == \[Nu] D[u[x, t], x, x] - 
      u[x, t] D[u[x, t], x], u[x, 0] == E^-x^2, u[-5, t] == u[5, t]}
   , u, {x, -5, 5}, {t, 0, 4}
   , Method -> {"MethodOfLines", 
     "SpatialDiscretization" -> {"TensorProductGrid", 
       "Coordinates" -> {mygrid}},
     Method -> {"FixedStep", Method -> "ExplicitRungeKutta"}}, 
   StartingStepSize -> 1/2^13]]

The following gives the distinct step sizes:

u["Coordinates"] /. mysol // Map@Differences // Map@DeleteDuplicates

Examples adapted from the [Method of Lines tutorial.]

POSTED BY: Michael Rogers
Posted 5 months ago

I don't understand the connection of your remarks about errors and error message to the problem at hand. I get no errors with this:

Clear[foo, sol, x, w, psi];
ClearSystemCache["Numeric"]; (* Why? IDK. Seems to fix a bug when re-run. *)
nstep = 0; neval = 0; (* steps, PDE evaluations *)
PrintTemporary@Dynamic@{Clock@Infinity, foo, nstep, neval}; (* monitors progress *)
xmax = 10;
wmax = 10;
xmin = 10^-20;
psiinit[w_, x_] := (Exp[-w^2]/Sqrt[\[Pi]] - w*Erfc[w])*x;
sol = NDSolve[{2*x*D[psi[w, x], x] == 
    D[psi[w, x], {w, 2}] + 2*w*D[psi[w, x], w] - (x*psi[w, x])^2, 
   psi[w, xmin] == psiinit[w, xmin], 
   Derivative[1, 0][psi][0, x] == -x, psi[wmax, x] == 0}, 
  psi, {w, 0, wmax}, {x, xmin, xmax},
  PrecisionGoal -> 20, AccuracyGoal -> 10,
  WorkingPrecision -> 40,
  StepMonitor :> (++nstep; foo = x),
  EvaluationMonitor :> (++neval),
  Method -> {"MethodOfLines", 
    "SpatialDiscretization" -> {"TensorProductGrid", 
      "MaxStepSize" -> 1/5000, "DifferenceOrder" -> 4}}]

Takes a few seconds to set up the first step. Takes several minutes to get to x > xmin. I have only 32GB RAM and NDSolve[] uses around 38.5GB, so much of the slowness is probably due to swapping.


Have you tried "Pseudospectral"?:

Method -> {"MethodOfLines",
  "SpatialDiscretization" -> {
    "TensorProductGrid",
    (*"MinPoints" -> 65,(* or higher *)*)
    "DifferenceOrder" -> "Pseudospectral"}}

A trial run shows that psi approaches zero very rapidly, going below $10^{-20}$. To get 20 digits of precision, the AccuracyGoal needs to be more than the maximum of $20-\log_{10} | \psi |$ throughout the solution domain. I haven't been able to determine the maximum. To achieve, say, AccuracyGoal -> 60, I need WorkingPrecision to be around 400 or more, or some internal numerical glitch happens and I get an error like the one you described (but did not provide code for!). It also takes a long time. Note: Compilers do not give run-time errors, which is what this is. I used to get an uncaught exception and a "core dump" that could be used with a debugger, if the program had been compiled with debugger metadata. Mathematica gives you a stack trace, but only of user-level expressions. I don't get to trace the numerical glitch, which probably occurs in an untraceable library anyway. It would be nice if it indicated the "non-numerical value" it encountered, instead of just giving the time. It would be a good hint to what went wrong.

POSTED BY: Updating Name
Posted 5 months ago

I tried some settings for WorkingPrecision, but then, if I remember, I get a series of error messages stating that something is not a number. The messages are totally incomprehensible, like most of error messages in MATHEMATICA, giving no clues to where and why the errors occurred. Forgive my irritation, but I really don't understand why the error handling is so bad in MATHEMATICA. Any compiler for programming languages is usually able to provide a precise location of errors, together with suggestions for eliminating the errors. Why this can't be done in MATHEMATICA? This should be particularly easy given the fact that compilation/execution is limited to a clear sequence of, or even to single commands. Leslaw

POSTED BY: Leslaw Bieniasz

Have you tried WorkingPrecision -> 40 or some such setting for it?

POSTED BY: Michael Rogers
Posted 5 months ago

Well, thanks, this should be useful for output purpose, but my main concern and question is how to force NDSolve to produce results that have a precision of number representation higher than the standard 16 digits, and how to activate an effective use of higher order accurate discretizations (hopefully) built-into the procedure. My goal is to obtain the solution having many accurate digits, ideally 19-20 digits. This cannot be normally achieved by using standard floating point numbers and discretizations that are only 2nd or even 4th order accurate. Leslaw

POSTED BY: Leslaw Bieniasz

Sorry, the computer ate my post. I don't have time to rewrite it. :/

Here is the full precision (53-bits, 15.95 digits):

Style[psi[0, 7] /. sol, PrintPrecision -> 20]
(*  {1.1442768817404132}  *)

Read http://reference.wolfram.com/language/tutorial/Numbers.html for more about kinds of numbers.

POSTED BY: Michael Rogers
Reply to this discussion
Community posts can be styled and formatted using the Markdown syntax.
Reply Preview
Attachments
Remove
or Discard

Group Abstract Group Abstract