I'm trying to solve a set of differential equations numerically to get a 3D plot, but I am getting multiple different warnings and errors. First of all, here is the code:

tdot[\[Lambda]_, q_] := Simplify[1/(1 - 1/r)] /. r -> r[\[Lambda]]
rdot[\[Lambda]_, q_] := Simplify[(1/r)*Sqrt[r^2 - q^2*(1 - 1/r)]] /. r -> r[\[Lambda]]
\[Phi]dot[\[Lambda]_, q_] := Simplify[-(q/r^2)] /. r -> r[\[Lambda]]
tasym[\[Lambda]_, q_] := \[Lambda] + Log[\[Lambda]] - 1/\[Lambda] + (q^2 - 2)/(4*\[Lambda]^2) + (3*q^2 - 4)/(12*\[Lambda]^3) + (-3*q^4 + 8*q^2 - 8)/(32*\[Lambda]^4) + 
   (-9*q^4 + 20*q^2 - 16)/(80*\[Lambda]^5)
rasym[\[Lambda]_, q_] := \[Lambda] + q^2/(2*\[Lambda]) - q^2/(4*\[Lambda]^2) - q^4/(8*\[Lambda]^3) + (3*q^4)/(16*\[Lambda]^4) + (q^6/16 - q^4/20)/\[Lambda]^5
\[Phi]asym[\[Lambda]_, q_] := q/\[Lambda] - q^3/(3*\[Lambda]^3) + q^3/(8*\[Lambda]^4) + q^5/(5*\[Lambda]^5)
\[Lambda]min = 0; 
\[Lambda]max = 20; 
\[Lambda]inf = 999; 
qlist = Array[N[#1] & , 100, {-20, 20}]; 
maxder = 999999; 
eq[q_] := {t[\[Lambda]], r[\[Lambda]], \[Phi][\[Lambda]]} /. NDSolve[{Derivative[1][t][\[Lambda]] == tdot[\[Lambda], q], Derivative[1][r][\[Lambda]] == rdot[\[Lambda], q], 
      Derivative[1][\[Phi]][\[Lambda]] == \[Phi]dot[\[Lambda], q], t[\[Lambda]inf] == tasym[\[Lambda]inf, q], r[\[Lambda]inf] == rasym[\[Lambda]inf, q], 
      \[Phi][\[Lambda]inf] == \[Phi]asym[\[Lambda]inf, q], WhenEvent[{Abs[Derivative[1][t][\[Lambda]]] > maxder || 
         Abs[Derivative[1][r][\[Lambda]]] > maxder || Abs[Derivative[1][\[Phi]][\[Lambda]]] > maxder}, "StopIntegration"]}, 
     {t[\[Lambda]], r[\[Lambda]], \[Phi][\[Lambda]]}, {\[Lambda], \[Lambda]min, \[Lambda]inf}, 
     {"ExtrapolationHandler" -> {Indeterminate & , "WarningMessage" -> False}}][[1]]
eqlist = (eq[#1] & ) /@ qlist; 
tlist = eqlist[[All,1]]; 
rlist = eqlist[[All,2]]; 
\[Phi]list = eqlist[[All,3]]; 
surface = MapThread[{#2*Sin[#3], #2*Cos[#3], If[Abs[#3] < Pi, #1, Indeterminate]} & , {tlist, rlist, \[Phi]list}]; 
Show[ParametricPlot3D[surface, {\[Lambda], \[Lambda]min, \[Lambda]max}, PlotRange -> {{-20, 20}, {-30, 20}, {-30, 20}}], ImageSize -> Large]

Running this code, with the parameters I chose after playing around with them, gives no warning message, and results in what i'm trying to get, but only half of it. For certain values of the parameters $\lambda_{inf}$, $maxder$, the warning message in the title of the question appears, however after reading other questions regarding this issue I don't think this is too much of a problem.

The main problems start when i set $\lambda_{min}=-30$ which would give me the other half of the plot. the first warning message that appears is

NDSolve::mxst: Maximum number of 129336 steps reached at the point [Lambda] == -0.53469.

same thing for other values of $\lambda$. At first I tried overcoming this by increasing MaxSteps, however this didn't work for me as Mathematica would just use up all my RAM and my computer would stop working.

To investigate I tried to solve only for the negative values, I set $\lambda_{min}=-30$ and $\lambda_{max}=-5$, and set NDSolve to solve only between $\lambda_{min}$ and $\lambda_{max}$. What happens with these settings is the following warning message:

NDSolve`ProcessSolutions::nodata: No solution data was computed between [Lambda] == -30. and [Lambda] == -5..

Which is weird since I should have a solution, unless I made some dumb mistake. Reading other questions, I saw this could maybe be a case of backslide where since the solution doesn't adhere to the new standards of NDSolve it doesn't give any solution. However this is just speculation.

I also tried adding some methods to NDSolve with no results, since I don't know much about them. What I hope is that by tweaking some NDSolve parameters or using some methods I can manage to get a result, so any suggestion in this direction is very welcome.


1 Answer 1


Rationalize[qlist] and ParametricNDSolveValue evaluates without error

 tr\[Phi] =ParametricNDSolveValue[{Derivative[1][t][\[Lambda]] ==tdot[\[Lambda], q],Derivative[1][r][\[Lambda]] == rdot[\[Lambda], q],Derivative[1][\[Phi]][\[Lambda]] == \[Phi]dot[\[Lambda], q],t[\[Lambda]inf] == tasym[\[Lambda]inf, q],r[\[Lambda]inf] ==rasym[\[Lambda]inf, q], \[Phi][\[Lambda]inf] == \[Phi]asym[\[Lambda]inf, q] ,WhenEvent[{Abs[Derivative[1][t][\[Lambda]]] > maxder ||Abs[Derivative[1][r][\[Lambda]]] > maxder ||Abs[Derivative[1][\[Phi]][\[Lambda]]] > maxder}, "StopIntegration"] }
, {t ,r , \[Phi] }, {\[Lambda], \[Lambda]min, \[Lambda]inf}  , {q}, \{"ExtrapolationHandler" -> {Indeterminate &, "WarningMessage" -> False}}] 

Plot[Table[Through[ tr\[Phi][q][\[Lambda]]],{q,Rationalize[qlist]}]  , {\[Lambda], \[Lambda]min, \[Lambda]inf},Evaluated -> True]

WhenEvent seems to be unnecessary.

  • $\begingroup$ This does not seem to work for me for the same reason my code doesn't work, when i put $\lambda_{min}=-30$ the domain of the solution is $[-0.33,\lambda_{inf}]$ but I'm expecting a solution also for lower values of $\lambda$. $\endgroup$
    – Valentin
    Jun 17, 2021 at 12:17
  • $\begingroup$ tr\[Phi][-30 ] gives a solution in the range 0<\[Lambda]<999 $\endgroup$ Jun 17, 2021 at 12:26
  • $\begingroup$ Yes, but i'm also looking for $-30<\lambda<0$, not q=-30. That is my main problem, my original code already gives a good plot for $0<\lambda<999$. But Mathematica says there is no solution for that domain, or reaches MaxSteps and stops. Maybe i could force NDSolve to be less precise somehow? $\endgroup$
    – Valentin
    Jun 17, 2021 at 12:34
  • $\begingroup$ Try Method->"StiffnessSwitching" $\endgroup$ Jun 17, 2021 at 12:51
  • $\begingroup$ Thank you, combining this with your code, and some tweaking on precision and accuracy seems to give a result. I'm going to try and get a plot out of it now $\endgroup$
    – Valentin
    Jun 17, 2021 at 13:00

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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