3
$\begingroup$

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.

$\endgroup$
  • 1
    $\begingroup$ 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). $\endgroup$ – Michael E2 May 18 '16 at 2:38
  • $\begingroup$ 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. $\endgroup$ – george2079 May 18 '16 at 20:27
  • $\begingroup$ 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) $\endgroup$ – george2079 May 18 '16 at 20:32
2
$\begingroup$

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

| improve this answer | |
$\endgroup$
1
$\begingroup$

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

| improve this answer | |
$\endgroup$

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.