How to eliminate the divide by zero error in NDSolve by implicit methods on fixed grids?

Questions
, ,
1

Hello, it is not my first post about how to use NDSolve for solving PDEs, but unfortunately, I have been unable to make any progress for more than a month. I think I want to accomplish a simple task: solve a certain PDE using exclusively fixed, uniform grids, and observe how the errors diminish when I refine the grids and/or increase discretization orders, and simultaneously increase the working precision. My goal is to obtain the solution having ca. 19 accurate digits. Unfortunately, whatever I do, I get different kinds of errors, and it seems nobody can tell me how to proceed. Is there noone there, in this list, who knows how NDSolve works? Maybe some authors of the marvel called NDSolve?
I now try to solve a PDE of the nonlinear reaction-diffusion-convection type, in which there is a coefficient standing at the temporal derivative (y is the temporal variable), which becomes zero when y=0. For this reason, I decided to use ImplicitRungeKutta methods. As far as I understand implicit methods, they should sample the PDEs at y>1, but never at y=0, so that there should be no problem with obtaining the solutions. But I get error messages stating that there are divide by zero errors. Here is my code:

ClearAll;
y0=0;
wmax=10;
ystep=1/10;
difford=2;
wprec=8;
pdesol=NDSolve[{D[psi[w,y],{w,2}]+2*w*D[psi[w,y],{w,1}]-4*y*(1-y)*D[psi[w,y],{y,1}]-4*(psi[w,y]^2-psi[w,y])==0,psi[w,y0]==Erf[w],psi[0,y]==0,psi[wmax,y]==Erf[wmax]},psi,{w,0,wmax},{y,y0,1},WorkingPrecision->wprec ,Method->{"FixedStep", Method->{"ImplicitRungeKutta","DifferenceOrder"->difford}},
StartingStepSize->ystep,
Method->{"MethodOfLines","DifferentiateBoundaryConditions"->true,"SpatialDiscretization"->{"TensorProductGrid", "MaxPoints"->20, "MinPoints"->20,"DifferenceOrder"->2}}];

What should I do to get at last the results I need? How to eliminate the divide by zero errors?
I know, of course, that I can select y0 not identical to zero, but equal to some small value, but this causes other error messages, depending also on the choice of other parameters. Some of these error messages suggest that my choice of MinPoints and MaxPoints is ignored, and NDSolve uses
10000 grid points for unclear reason. So, I am also not sure whether I really have fixed uniform grids.
If ImplicitRungeKutta methods are not suitable, than what other kind of fixed grid implicit methods can be used for my purpose? I mean here built-in methods, as I expect devising my own procedures would be much more difficult.
Leslaw

2

It seems that there is a compatibility issue on the w axis:

eqs = {D[psi[w, y], {w, 2}] + 2*w*D[psi[w, y], {w, 1}] -
     4*y*(1 - y)*D[psi[w, y], {y, 1}] -
     4*(psi[w, y]^2 - psi[w, y]) == 0,
   psi[w, 0] == Erf[w], psi[0, y] == 0,
   psi[10, y] == Erf[10]};
eqs[[1]] /. y -> 0 /. psi -> Function[{w, y}, Erf[w]]
3

Thank you, indeed, I had an error in my PDE (one missing “y”), so that the corrected code is:

ClearAll;
y0=0;
wmax=10;
ystep=1/10;
difford=2;
wprec=8;
pdesol=NDSolve[{D[psi[w,y],{w,2}]+2*w*D[psi[w,y],{w,1}]-4*y*(1-y)*D[psi[w,y],{y,1}]-4*y*(psi[w,y]^2-psi[w,y])==0,psi[w,y0]==Erf[w],psi[0,y]==0,psi[wmax,y]==Erf[wmax]},psi,{w,0,wmax},{y,y0,1},WorkingPrecision->wprec ,Method->{"FixedStep", Method->{"ImplicitRungeKutta","DifferenceOrder"->difford}},
StartingStepSize->ystep,
Method->{"MethodOfLines","DifferentiateBoundaryConditions"->true,"SpatialDiscretization"->{"TensorProductGrid", "MaxPoints"->20, "MinPoints"->20,"DifferenceOrder"->2}}];

However, after correcting this mistake, I still get the same error message: Power: Infinite expression 1/0, etc. plus many following messages. As far as I can guess, they arise as a result of dividing the pde by the factor 4y(1-y) standing at the first derivative of psi with respect to y, at y=0. But I don’t understand why such a calculation is performed; it should not occur in implicit methods.
I don’t know whether my specification of “Method” is correct: is it OK to have two different such specifications? Maybe one of them is ignored? How else I could select methods for spatial and temporal discretizations?
Leslaw

4

I am no expert in numerical methods, but when y == 0 the differential equation give no information on what direction to go as y increases. Something like DSolve[{y f'[y] == 3, f[0] == 1}, f, y].

By the way, I suspect that ClearAll; achieves no purpose at all, in spite of its name.

5

Yes, the disappearance of the temporal derivative at y=0 worries me, too. But I have not invented this equation. It is taken from some literature and certainly possesses a solution which was verified numerically, albeit not with a high accuracy that I want to achieve.
Leslaw

6

My previous posts didn’t seem to solve your problem. This will be no different. My apologies. I have two remarks, which are the sources of two of the problems you have in the code for NDSolve.

  1. Only the first option setting is used (in any function). So the second Method setting is ignored. There are reasons that this behavior is convenient for programs, but none of them are relevant here.

  2. The basic data computed for each step by the ODE solvers in NDSolve (including those used in the method of lines) consists of, at a minimum, the current time, value(s) of the dependent variable(s), and the value(s) the temporal derivative(s). This last one means that the value of the derivative is computed for the initial condition, which is why you get the 1/0. error. It is also stored in the output (the interpolating functions returned). In fact, this will happen no matter what method you use, because it happens in the initialization phase before any steps are taken. (Try NDSolve`ProcessEquations[<<NDSolve arguments go here>>], and you will see the error.) Furthermore, you should expect a possible problem when (or if) the integration reaches the other singular point y == 1. I don’t have a solution, sorry.

  3. Maybe a third point. Not all implicit methods avoid evaluating the ODE at the end points of a time step. Many do avoid them, and I don’t think it’s an issue here.

I hope someone who can help you finds this discussion.

Appendix

The thingsbelow may be read about in the document (Tech Notes) for NDSolve.

The complete solution data structure components and names:

Table[
 component -> 
  NDSolve`SolutionDataComponent[{t, {x}, {u}, {Derivative[0, 1][
      u]}, {a}, {b}, {c}, {d}},
   component],
 {component, {"Time", "Space", "DependentVariables", 
   "TemporalDerivatives", "DiscreteVariables", "IndexedVariables", 
   "Parameters", "SensitivityParameters"}}]
(*
{"Time" -> t,
 "Space" -> {x},
 "DependentVariables" -> {u},
 "TemporalDerivatives" -> {Derivative[0, 1][u]},
 "DiscreteVariables" -> {a},
 "IndexedVariables" -> {b},
 "Parameters" -> {c},
 "SensitivityParameters" -> {d}}
*)

Not all solution variables are used in each problem:

state = First@
   NDSolve`ProcessEquations[{Derivative[2, 0][u][x, t] == 
      Derivative[0, 1][u][x, t], u[x, 0] == Exp[-x^2], u[0, t] == 1, 
     u[1, t] == 1/E}, u, {x, 0, 1}, {t, 0, 1}];
state@"VariableData"
(*  {t, {x}, {u}, {Derivative[0, 1][u]}, {}, {}, {}, {}}  *)

Executing the following for state above will show the actual data, at the beginning and then after a time-stepping through to t == 0.1:

state@"SolutionData"["Forward"]

NDSolve`Iterate[state, 0.1]

state@"SolutionData"["Forward"]
7

  1. If the second setting for the “Method” is ignored, then how does NDSolve decide about spatial discretization? Is there no way to choose and control simultaneously both spatial and temporal discretizations?

  2. I would argue that the obligatory calculation of the y-derivative from the PDE at y=0 is in disagreement with the usual mathematics of IBVPs. Normally, a PDE is to be obeyed at 0 < y < infinity, whereas at y=0 initial conditions are to be satisfied, not the PDE. And the y-derivative at y=0 is computed a posteriori, based on the numerical solution obtained, not from the PDE itself.
    (Of course I mean here PDEs involving only the first order y-derivative, not those of higher order).

  3. Which implicit methods sample the PDE at the beginnng of an integration step? In my opinion, the notion of “implicit” means that discretization is made at y values larger than the initial y value in any step. If the discretization is at the initial y value, we have an explicit method.

Lesław

8

  1. Here’s an example of the usage of Method adapted from the advice in the documentation for NDSolve. Note that it — um — “runs” but does not really solve your problem. And it only runs until it hits the other singular point y == 1, at which it blows up.

    y0 = 0 + 10^-100; (* you wish to avoid this, I know )
    wmax = 10;
    ystep = 1/10;
    difford = 2;
    wprec = 68;
    pdesol =
    NDSolve[{D[psi[w, y], {w, 2}] + 2    wD[psi[w, y], {w, 1}] -
    4    y*(1 - y)D[psi[w, y], {y, 1}] - 4(psi[w, y]^2 - psi[w, y]) == 0
    , psi[w, y0] == Erf[w]
    , psi[0, y] == 0
    , psi[wmax, y] == Erf[wmax]}
    , psi, {w, 0, wmax}, {y, y0, 1}
    , WorkingPrecision → wprec, PrecisionGoal → 10
    , Method → {“MethodOfLines”,
    “DifferentiateBoundaryConditions” → True,
    “SpatialDiscretization” → {“TensorProductGrid”
    , “Coordinates” → {N[Subdivide[0, wmax, 20]]}
    , “DifferenceOrder” → 2}
    , Method → {“FixedStep”, Method → {“ImplicitRungeKutta”
    , “DifferenceOrder” → difford
    , “ImplicitSolver” → {“Newton”, AccuracyGoal → 20,
    PrecisionGoal → 20, “IterationSafetyFactor” → 1}}}}
    , StartingStepSize → ystep]

  2. I know a lot about how NDSolve works, but not everything. Either I don’t know how to make it do what you describe, or NDSolve is not set up to do it.

  3. Any IRK method that has a leading 0 in the c-vector (of standard Butcher table) evaluates the ode at the beginning of the time step. And any RK method whose a-matrix has nonzero entries on the diagonal or above is called “implicit.” It means the algebraic equations of the RK method cannot be solved explicitly for the updates at each stage in terms of earlier stages. E.g. y_{k+1} = [f(x_k, y_k) + f(x_k+h,y_{k+1})]\,h/2 is the implicit trapezoidal rule and evaluates the ode at both endpoints of the time step. Or Lobatto methods such as

    {a, b, c} =
    NDSolve`ImplicitRungeKuttaLobattoIIIACoefficients[
    4, MachinePrecision]

The c-vector is {0., 0.5, 1.}.

9

Michael: Thank you. But, in view of former letters, what are actually the discretizations along y and w variables, keeping in mind that one of the settings for “Method” is ignored? Is there any way to find out which method NDSolve actually uses?

Leslaw

10

The w and y steps used are available with

psi["Coordinates"] /. First@pdesol

As for the method:

  1. For the Method option, it’s always the first occurrence that is used.

  2. For the PDE method, it’s the method of lines.

  3. I assume you’re asking about syntax. For the time integration method (the Method suboption of the main Method option), it’s the fixed-step strategy with the submethod specified by the Method sub-suboption of the Method suboption, which is "ImplicitRungeKutta". Note there is only one Method option specified for NDSolve, and it has Method options nested within it. The general syntax is

    Method → { name , methodoptions }

Some NDSolve methods are classified as “strategies” in the documentation. These may apply one or more step-integrators in some way. The step-integrators are specified with Method -> meth. Mathematica supplies default ones if no method is specified.

Here is the single method option, with one item per line, indented to reflect nesting:

Method -> {
  "MethodOfLines"
  , "DifferentiateBoundaryConditions" -> True
  , "SpatialDiscretization" -> {
    "TensorProductGrid"
    , MaxPoints -> 20 (* might be better than "Coordinates" *)
    , "DifferenceOrder" -> 2}
  , Method -> { (* time integration method for MethodOfLines *)
      "FixedStep"
      , Method -> { (* submethod for time steps in FixedStep method *)
          "ImplicitRungeKutta"
          , "DifferenceOrder" -> difford
          , "ImplicitSolver" -> {
              "Newton"
              , AccuracyGoal -> 20
              , PrecisionGoal -> 20
              , "IterationSafetyFactor" -> 1}}}}

If you were asking, how to be sure NDSolve did what it was told to do, that’s less clear. I assume it does. I have checked sometimes, usually for class demonstrations, that it does what I’m going to tell the class it does. I might have gotten messages saying in effect, “This method failed. Using method instead.” NIntegrate does that, but I can’t remember if NDSolve does or not.

The following are probably less helpful. You’ll have to read through “Components and Data Structures” and maybe the “Method Plugin Framework” tutorials to understand it.

This sometimes shows the call to initialize a method. For whatever reason, it doesn’t show all the methods, either because they are initialized in another way or because they can’t be traced.

Trace[
 NDSolve[ <NDSolve args> ],
 _NDSolve`InitializeMethod,
 TraceInternal -> True]

This returns the result of initialization, an NDSolve`StateData object, about which you will have to read.

state = First@NDSolve`ProcessEquations[ <NDSolve args> ]

This will return some information about the method:

NDSolve`Iterate[state, 1] (* you need to advance the integration at least one step *)
state@"MethodData"[]
11

In this case, I’m “Updating Name”. That happens when it takes too long to fill in the answer and I get logged out. It happens to others, too. You can’t always be sure whom “Updating Name” refers to.

12

One technical question: Is there any way to save in a file and/or print the messages from this forum? I cannot see any option for this.

Leslaw

13

I don’t know of one, other than to save or print the webpage, which can be done from the browser. Saving as HTML captures the text but not all the styling. Saving as or “printing” to PDF looks better. Saving as a “webarchive”, which I never tried before, seems to save to disk a faithful reproduction of the page (worked on this page; took up 3.5MB of space).

14

What do you exactly mean by saying that it “runs”? And how do you find that the solution blows up at y=1? When I try to run this code, I get lots of error messages, and although some Interpolating Function object is created, nothing can be plotted (the plotting box is empty).

One more question: I have an independent numerical procedure that calculates the derivative
d(psi[w,y])/ dy at y=0 and returns it as a function psi1[w] (this can be done with a pretty good accuracy, but not analytically). The derivative is finite. Is there any way to combine such calculated derivative with the pde to be solved, so that NDSolve does not encounter the 1/0 problem while testing the pde at y=0? I am trying to replace the pde by the statement

D[psi[w,y],{y,1}]== If[y==0, psi1[w], b]

where “b” is the symbolic formula for the derivative obtained from the pde. However, this approach also produces errors, namely there is a message that 0 is not a valid variable.

Leslaw

15

Hello again.

I now manged to force NDSolve to produce something which looks like a reasonable solution, albeit not without weird error messages, and not very accurately. I replaced the initial condition at y0=0 by an approximate initial condition at y0 > 0, equal to the two terms of the solution expansion into power
series of y. For this, the first y-derivative of the solution at y=0 is calculated separately by solving an appropriate two-point BVP by the shooting method (this takes a bit of time).
Unfortunately, when I start changing options in NDSolve, solving the pde fails with numerous error messages, and I cannot get solutions. I need to use exclusively non-adaptive methods, and fixed
y- and w-grids, so that I can perform an elementary convergence experiment by refining the grids. Sorry to say this guys, but it is incomprehensible for me why NDSolve does not allow one to perform such an experiment (as it seems). With this deficiency, it is just a fancy toy but not a scientific tool
for solving real-life problems. If anybody knows how to change my code (by setting appropriate options in NDSolve) so that I can test the convergence using fixed, uniform grids, and achieve 19-20 digits of accuracy, please help. My current code is included below. Please note that I am particularly interested in determining accurately the first w-derivative of the solution at w=0, which is plotted by the final commands in the code.

Leslaw

ClearAll;
(*-----------------------------------------------------*)
wmax=10;
(*-----------------------------------------------------*)
(* First order expansion term (in powers of y) for y>0 *)
(*  psi[w,t]=Erf[w]+Psi1[w]*y + ...  *)
wstep1=1/100;
difford1=25;
wprec1=80;
SS1[wm_,g0_,wp_]:=NDSolve[{psi1''[v]+2*v*psi1'[v]-4*psi1[v]==-4*Erf[v]*Erfc[v],psi1[0]==0,psi1'[0]==g0},psi1,{v,0,wm},WorkingPrecision->wp ,Method->{"FixedStep", Method->{"ImplicitRungeKutta","DifferenceOrder"->difford1}}, StartingStepSize->wstep1];
(* Assume the interval of g0 for the Brent method:  0.3 < g0 < 0.35  *)
(* Numerical solution for psi1 at wmax by the shooting method *)
psi1wmax[wm_,g0_?NumericQ,wp_]:=Module[{},sol1=SS1[wm,g0,wp];psi1[wm]/.sol1];
(* Determination of the gradient at w=0 that makes psi1wmax = 0 *)
g0acc1[wm_,g0min_,g0max_,wp_]:=Last[First[FindRoot[psi1wmax[wm,g0,wp]==0,{g0,g0min,g0max},Method->"Brent",WorkingPrecision->wp]]];
(* Calculation of the computed psi1 value for given w *)
Psi1[w_]:=First[(psi1[w]/.sol1)];
gradient0=g0acc1[wmax,30/100,35/100,wprec1]
Psi1[wmax]
Plot[Evaluate[Psi1[w]],{w,0,wmax},PlotRange->All]
(*----------------------------------------------------*)
(* Now pass to solving the pde *)
y0=10^-20;
ystep=1/10;
difford=2;
wprec=16;
(*-----------------------------------------------------*)
pde=D[psi[w,y],{w,2}]+2*w*D[psi[w,y],{w,1}]-4*y*(1-y)*D[psi[w,y],{y,1}]-4*y*(psi[w,y]^2-psi[w,y])==0;
pdesol=NDSolve[{pde,psi[w,y0]==Erf[w]+Psi1[w]*y0,psi[0,y]==0,psi[wmax,y]==Erf[wmax]},psi,{w,0,wmax},{y,y0,1}];
Plot3D[Evaluate[psi[w,y]/.pdesol],{w,0,wmax},{y,y0,1},PlotRange->All]
Psip[w_,y_]:=Evaluate[D[psi[w,y],{w,1}]/.pdesol];  (* First w-derivative *)
Psip[0,0]
Psip[0,1]
Plot[Psip[0,y],{y,y0,1},PlotRange->All]
Plot3D[Psip[w,y],{w,0,wmax},{y,y0,1},PlotRange->All]
16

By “runs”, I mean NDSolve[...] evaluates to its normal output, namely, rules with interpolating functions. This means NDSolve[...] accepted the input and ran the solver (a few steps only, in this case).

By “doesn’t run”, I would mean that NDSolve[...] returns the input as output. That happens, for instance, if there is an error in the input (equations, options, etc.).

I would consider the numerous messages as “warnings” not “errors”, indicating that there seems to be difficulty in computing the solution. Many of the warnings indicate possible severe numerical errors, so calling them “errors” is reasonable. I call them “warnings” to distinguish them from syntax/input errors that prevent NDSolve from running.

The solution returned is computed at y == 0, 0.1, 0.2,..., 0.9, after which point there is a numerical error in the solve that NDSolve cannot recover from. That must occur while computing the step to y == 1.

You might be able to use Piecewise[] to define separate derivative formulas for y == 0 and y != 0. I have a busy day at work, and I’m not sure how to quickly check that there aren’t problems with that suggestion. I’ve done it with success in ODEs, though.