NDSolve returns solution with single point domain

Question

Thanks to helpful comments from Michael E2 and George2079, I was able to focus in on exactly the source of the issue. With some simplification, I can reduce the problem to:

NDSolve[{Derivative[1][y][t] == (-2 E^(1/5 y[t]) + 1)/(-600 E^(1/5 y[t]) +   
        20338 + (2/10) t ), y[101690] ==  0}, y, {t, 0, 101690}, 
        MaxSteps -> \[Infinity], WorkingPrecision -> 150]

It quickly throws:

NDSolve::ndcf: Repeated convergence test failure at t == <<185>>; unable to continue. >>

But this problem has an analytical solution:

y[t_] := 1/600 (2 t + C[1] + 3000 
         ProductLog[C[2], 
                   -((E^(-(t/1500) - C[1]/3000) (101690 + t))/3000)])/. 
         {C[1] -> 0, C[2] -> -1}

Which you can verify via:

In[1]:= N[y[101690]]
Out[1]= 0.

In[2]:= Simplify[y'[t] - 
                (-2 E^(1/5 y[t]) + 1)/(-600 E^(1/5 y[t]) + 20338 + (2/10) t)]
Out[2]= 0

This function is pretty well behaved.

I am hoping to understand the related issues with NDSolve, because I'm relying on that engine for more complex problems and this issue is a barrier to its application.

The documentation for the ndcf error says:

If you see this message and the cause is not evident from the nature of the example, please send your example to Technical Support so that it can be investigated.

I get a step-size-too-small error. You have to increase the precision to around 100 to get effectively nonzero step sizes. But they're still so small (~10^-80 - 10^-30 for PrecisionGoal -> 8) that integration will never finish (effectively). — Michael E2, May 18 '16 at 02:38
wrapping everything in SetPrecision is unintentionally also converting integers to reals. Now your equations are expressed in terms of variables such as Cee[4.000000000000000000....]. Since now the symbols in the vars list (eg Cee[4]) do not appear in the equations you get a trivial solution. — george2079, May 18 '16 at 20:27
You can fix that issue by instead doing InitialConditions = {Cee[4][t0] == -270917. ... } /. x_Real :> SetPrecision[x, prec]. (With that fix I get an out of memory error though) — george2079, May 18 '16 at 20:32

score 2 · Answer 1 · answered Jan 22 '17 at 20:42

The simplified system makes me forget why I suggested the issue concerned working precision. (The edit history shows the OP was much more complicated.) Aside from the fact that @xzczd's answer shows a high WorkingPrecision is not needed, the convergence problem arises from PrecisionGoal and AccuracyGoal being automatically raised (to half WorkingPrecision). That imposes a rather strict requirement with WorkingPrecision -> 150, which apparently NDSolve cannot satisfy. There are two principal reasons for raising WorkingPrecision (WP):

One cannot achieve the precision and accuracy goals at the current WP. In such a case, raise WP but also manually set PrecisionGoal and AccuracyGoal to their current values (don't raise them).
One wants to raise PrecisionGoal and AccuracyGoal.

The present case is the first one. We should set PrecisionGoal -> 8 and AccuracyGoal -> 8 (or other desired values). With the high WP, NDSolve can get quite close to the singularity, making quite small steps. I suggest setting MaxSteps to a finite value; 100000 takes a few seconds. You could set MaxSteps -> 10000 (the default) to get the idea.

time = 0; steps = 0;
Dynamic@{time, steps}  (* monitor progress; it's educational *)

steps = 0;
{sol} = NDSolve[{Derivative[1][y][t] == (-2 E^(1/5 y[t]) + 1)/(-600 E^(1/5 y[t]) + 
       20338 + (2/10) t), y[101690] == 0}, y, {t, 0, 101690}, 
  MaxSteps -> 100000, WorkingPrecision -> 150, PrecisionGoal -> 8, 
  AccuracyGoal -> 8, StartingStepSize -> 0.001, 
  StepMonitor :> (time = t; steps++)]


ListLinePlot[N[y /. sol], PlotRange -> All]

Mathematica graphics

score 1 · Answer 2 · answered Oct 24 '16 at 07:23

Without the 2 options, NDSolve manages to give the correct solution:

nsol = NDSolveValue[{Derivative[1][y][
     t] == (-2 E^(1/5 y[t]) + 1)/(-600 E^(1/5 y[t]) + 20338 + (2/10) t), y[101690] == 0},
   y, {t, 0, 101690}]

NDSolveValue::ndsz

NDSolve spits out ndsz warning and stops at about 6883.23, but this behavior is reasonable, because when t < 6883.23 there's no real solution:

sol[t_] = 1/
     600 (2 t + C[1] + 
      3000 ProductLog[C[2], -((E^(-(t/1500) - C[1]/3000) (101690 + t))/3000)]) /. {C[
      1] -> 0, C[2] -> -1};

{{tL, tR}} = nsol["Domain"];
Plot[{Re@#, Im@#} &@sol@t // Evaluate, {t, 0, tR}, PlotRange -> All, 
  Epilog -> {Dashed, Line[{{tL, -10}, {tL, 20}}]}]~Show~
 Plot[nsol[t], {t, tL, tR}, PlotStyle -> {Thick, Dashed, Red}, PlotRange -> All]

Mathematica graphics

NDSolve returns solution with single point domain

2 Answers2

Linked