# NDSolve fails for certain choices of parameters and solve range

I'm trying to solve a pair of coupled ODEs with NDSolve.

I know roughly what the solution should look like (both should give periodic functions, pi/2 out of phase, the amplitude of which damp towards a common limiting value), but I am trying to work out what the limiting amplitude is for a range of parameters.

For certain choices of parameters, the solution looks how I would expect, but for others the NDSolve sporadically takes very large steps and the output either diverges, becomes constant, or shows other noisy/glitchy behaviour. Changing the NDSolve range also makes these effects come and go, and reducing the MaxStepSize or changing the Working Precision to anything but 16 makes things worse - so I'm pretty sure this is a computational problem.

I'm pretty new to Mathematica so this is probably something simple, but every suggestion I've seen in other posts (increasing workingprecision, limiting step size etc) has only made things worse so I'm a bit stumped!


(Define some constants)
wc0 = 0.9645*^8*B/mass /. {mass -> 40, B -> 1.845};
fc0 = wc0/2/[Pi];
(Define wz, w1, wc and wm as functions \
of trap voltage V)
wz[V1_] = Sqrt[4*0.9645*^8*1*V1/40/0.0189^2];
w1[V1_] = Sqrt[wc0^2/4 - wz[V1]^2/2];
wc[V1_] = wc0/2 + w1[V1];
wm[V1_] = wc0/2 - w1[V1];
(Want to find limiting value as these parameters are varied)
V = 244;
d = 100000;
rate = 2*^6;
rm0 = 5*^-6;
rc0 = 5*^-6;
(Solve for dynamics)
mysolve =
NDSolve[{
rc'[t]==rateB/w1[V]^2/2/rc[t]-ratewc[V]Arc[t]/w1[V]/2+drm[t],
rm'[t]==rateB/w1[V]^2/2/rm[t]+rate*wm[V]Arm[t]/w1[V]/2-drc[t],
rm[0] == rm0, rc[0] == rc0} /. {B -> 3.4^-4, A -> 3*^-3},
{rc[t], rm[t]}, {t, 0, 0.009}, MaxSteps -> Infinity
];

range = 1*^-6;
(Plot dynamics)
Plot[Evaluate[rc[t] /. mysolve], {t, 0.0, 0.008},
PlotRange -> {-range, range}, AspectRatio -> 0.25, ImageSize -> 800]
Plot[Evaluate[rm[t] /. mysolve], {t, 0.0, 0.008},
PlotRange -> {-range, range}, AspectRatio -> 0.25, ImageSize -> 800]

(Find limiting value)
intmin = 0.004;
intmax = 0.008;
lrc = Sqrt[2*NIntegrate[Evaluate[rc[t]^2 /. mysolve], {t, intmin, intmax}]/(intmax-intmin)]
lrm = Sqrt[2*NIntegrate[Evaluate[rm[t]^2 /. mysolve], {t, intmin, intmax}]/(intmax-intmin)]

(**)


The parameters listed above give a sensible result (first image) - but if, for instance, I change the NDSolve window to {t,0,0.004} or {t,0.04}, or change 'd=100000' to 'd=80000' (second image) it starts behaving strangely.

Any help would be very much appreciated!

• Be sure to use enough plot points in your plots. Try eg. PlotPoints->1000. The artifacts at ~0.0035 and ~0.0057 in the second plot remind me of ones I encountered when using default PlotPoints for a larger range of oscillating functions like these. – DaveStrider Sep 8 '14 at 10:59

Second answer -- OK, the first answer was hogwash (the curious can inspect the edit history). It pays sometimes to write out the equation and think about it first. The code for this one looks so complicated, but the equations basically have the form (here d = 80000)

rc'[t] == 2.05594*10^-10/rc[t] - 8189.14 rc[t] + 80000. rm[t],
rm'[t] == -80000. rc[t] + 2.05594*10^-10/rm[t] + 2189.14 rm[t]


We see that there are singularities along both axes rc = 0 and rm = 0. So the solution should not be crossing the axes at all. Here is a simplified vector field that will illustrate the problem NDSolve faces.

StreamPlot[{1/rc - rc + 2 rm, -2 rc + 1/rm + rm},
{rc, -3, 3}, {rm, -4.5, 4.5}, Axes -> True]


As a trajectory moves clockwise away from an axis, it approaches the next axis before turning toward an equilibrium point or out toward infinity, depending on the quadrant. The turn can be rather sudden. If NDSolve, taking a discrete step, happens to pass over the axis far enough, the solution will continue in the same direction. This will lead to a computed solution that is oscillatory. With the OP's coefficients, the regions near the axes where the phase vectors turn are very narrow and the OP's NDSolve code often jumps over it. When it doesn't jump over an axis, the oscillations are interrupted, until it does again. In short it's a step-size issue.

What should happen to a solution that starts in the first quadrant is that it should converge on the equilibrium position.

Solution

Here is the set-up of the problem. I put a variable prec to specify the precision, since it took some experimentation to find a robust setting. It may be that prec needs to be adjusted further to handle a particular use-case.

(*Solve for dynamics*)
d = 80000;
prec = 40;
de := SetPrecision[
{rc'[t] == rate*B/w1[V]^2/2/rc[t] - rate*wc[V]*A*rc[t]/w1[V]/2 + d*rm[t],
rm'[t] == rate*B/w1[V]^2/2/rm[t] + rate*wm[V]*A*rm[t]/w1[V]/2 - d*rc[t],
rm[0] == rm0, rc[0] == rc0} /. {B -> 3.4*^-4, A -> 3*^-3},
prec]


The equilibria:

NSolve[DeleteCases[de, _?(FreeQ[#, t] &)] /. {rc'[t] -> 0., rm'[t] -> 0.},
{rc[t], rm[t]}, Reals]
(*
{{rc[t] -> -2.24157*10^-7, rm[t] -> -1.14809*10^-8},
{rc[t] -> 2.24157*10^-7,  rm[t] -> 1.14809*10^-8}}
*)


There are several strategies and parameters of NDSolve that figure into obtaining an accurate solution. Some of them required experimentation, and I left some of this experimental code in, just in case someone wants to play with it. It is difficult to describe the interdependence of the elements of the code, so I will explain each one separately.

• WhenEvent -- This can be used to detect either rc or rm changing sign at a step taken by NDSolve. If this happens, NDSolve should back up and start over. The action "RestartIntegration" causes this.

• "LocationMethod" -- The setting "LocationMethod" -> "StepBegin" causes NDSolve to restart the integration from the beginning of the current step. If we restarted at the zero crossing, we would get a divide by zero error if using arbitrary precision. If we were using MachinePrecision, we would get, a more or less randomly, either the solution staying on the correct side of the axis or skipping across.

• "StartingStepSize" -- When NDSolve starts, this is irrelevant. Its purpose is to reset the step size after the WhenEvent action of "RestartIntegration". (This is the clever bit, IMO, if there is any cleverness at all here.)

• "IntegrateEvent" -- This makes NDSolve take smaller steps near an event. This seemed to be important sometimes.

• WorkingPrecision -- This problem seems to have trouble with round-off error and numbers rather close to zero. It is important to have a high working precision, so that the differential equation and events are calculated accurately.

• AccuracyGoal -- It is important that this be high. The coordinates of the equilibrium are around 10^-8 and sometimes smaller with other parameter settings. In some use-cases this may need to be raised.

• PrecisionGoal -- By default, PrecisionGoal is set to be half of WorkingPrecision. When a computation is numerically challenging, it is sometimes important to increase the ratio of WorkingPrecision to PrecisionGoal. This seems to be the case here. On the other hand, the higher the PrecisionGoal is, the more work NDSolve will do (i.e., it will take longer and use more memory). While I used round numbers, AccuracyGoal exceeds PrecisionalGoal by about how many leading zeros there are in smallest coordinate of the equilibrium position.

• MaxStepSize -- This particular problem seems hard to control with MaxStepSize but I left it in, commented out, in case it might be useful in a particular use-case.

• Reap/Sow -- This is also part of the code for experimenting. StepMonitor can be used to see what steps are taken in computing the solution. One can inspect the values of t, rc[t] and rm[t] actually used by NDSolve. I used it to monitor the step size. Uncomment the GridLines options in the Plot to see the steps.

Solution:

{mysolve, tvals} = Reap[NDSolve[{
de,
WhenEvent[rc[t] < 0, "RestartIntegration",
"LocationMethod" -> "StepBegin", "IntegrateEvent" -> True],
WhenEvent[rm[t] < 0, "RestartIntegration",
"LocationMethod" -> "StepBegin", "IntegrateEvent" -> True]},
{rc, rm}, {t, 0., 0.04},
MaxSteps -> Infinity,
StartingStepSize -> 1*^-11, (*MaxStepSize->1*^-6,*)
AccuracyGoal -> 20, PrecisionGoal -> 10, WorkingPrecision -> prec,
StepMonitor :> Sow[t]]];


Plot:

(* to determine the plot range *)
{cPR, mPR} = {Max[Abs@rc[t] /. #], Max[Abs@rm[t] /. #]} &@
NSolve[DeleteCases[de, _?(FreeQ[#, t] &)] /. {rc'[t] -> 0., rm'[t] -> 0.},
{rc[t], rm[t]}, Reals];

(*Plot dynamics*)
Plot[Evaluate[rc[t] /. mysolve], {t, 0.0, 0.04},
PlotRange -> {-0.1 cPR, 1.5 cPR}, AspectRatio -> 0.25,
ImageSize -> Full (*,GridLines->{First@tvals,None}*)]
Plot[Evaluate[rm[t] /. mysolve], {t, 0.0, 0.04},
PlotRange -> {-0.1 mPR, 1.5 mPR}, AspectRatio -> 0.25,
ImageSize -> Full]


• Thanks for your help - that definitely improves things (fewer glitches) but I am still finding lots of combinations give me the sort of output I posted above. For example, any value of d<97000 with the above parameters gives me a bad output. – Joe Goodwin Sep 5 '14 at 18:46
• Very nice answer! – user21 Sep 8 '14 at 8:13
• @JoeGoodwin Is the new answer any better? – Michael E2 Sep 13 '14 at 19:06