# NDSolve gets the wrong answer! How get NDSolve to correctly solve this "stiff" equation?

Below, I am trying to solve for x[t], in a sort of complicated differential equation. In this equation, x[t] is a 3-D vector with x, y, and z components. Here is my code with runTime=100 seconds.

runTime = 100;
f[t_] := 0.5907 + 0.56*Exp[-12.2  (t - 15)^2]
func = {Derivative[1][x][t] == -0.85*f[t]* x[t]\[Cross](x[t]\[Cross]{0, 0, 1}) - 176*x[t]\[Cross]{0.003, 0, -1*x[t] . {0, 0, 1}} + 0.1*x[t]\[Cross]Derivative[1][x][t], x[0] == {1, 0, 0}};
sol = NDSolve[func, {x}, {t, 0, runTime}, MaxSteps -> \[Infinity]];

temp = Table[x[t] /. sol[[1]], {t, 0, runTime, 0.001}];
phi = Arg[temp[[All, 1]] + I*temp[[All, 2]] ]; (* Finds phase *)
phi = phi - Accumulate[2*Pi*Unitize[Threshold[Prepend[Differences[phi],0.],0.5]]]; (* Unwraps phase *)
ListLinePlot[-phi, AspectRatio -> 0.2, PlotRange -> All, AxesLabel -> {None, "Phase, radians"}] (*Plots phase*)


Here is the graph it outputs, which is incorrect:

It also kicks back an error: "NDSolve:: Cannot solve to find an explicit formula for the derivatives. NDSolve will try solving the system as differential-algebraic equations."

If I instead use runtime=25;' seconds, it returns what I feel is the correct solution:

See the 2\pi$jump in phase? That's supposed to happen above. I tried Method -> "StiffnessSwitching", which worked for a similar system, but that kicks back a bunch of errors, including: NDSolve:: The method NDSolveStiffnessSwitching is not currently implemented to solve differential-algebraic equations. Use Method -> Automatic instead. Question: How can I get NDSolve to provide the correct solution for a longer run time? Would it help if we somehow informed Mathematica that x[t] has 3 dimensions? ## 2 Answers Just set a MaxStepSize -> 1 in NDSolve: sol = NDSolve[func, {x}, {t, 0, runTime}, MaxSteps -> Infinity, MaxStepSize -> 1];  Remaining part of the code is the same as in the question so omitted here. Alternatively, setting a single StartingStepSize -> 1 in NDSolve (yeah we don't need MaxStepSize in this case) can also fix the code: sol = NDSolve[func, {x}, {t, 0, runTime}, MaxSteps -> Infinity, StartingStepSize -> 1];  Remaining part of the code is the same as in the question so omitted here. But based on my limited experience I fear that this solution may not be that stable, so personally I don't quite recommend it. • Could you extend this answer with picture and code? Commented Aug 1 at 8:33 • @AlexTrounev Well, I just find it not that necessary. Anyway, updated. Commented Aug 1 at 11:16 • Could you add option StartingStepSize -> 10^-5? :) Commented Aug 1 at 13:23 • @AlexTrounev Actually for this problem, a single StartingStepSize->1 (yeah we don't need MaxStepSize) can also fix the code, but based on my limited experience I fear that this solution may not be that stable. But I agree that I should at least mention this, updated. Commented Aug 1 at 14:16 ## Analysis Analyzing the ODE, I'd say the issue is numerics. The coefficient f[t_] := 0.5907 + 0.56*Exp[-12.2 (t-15)^2] is constant at machine precision except in an interval $$[13.26\mskip-3mu\dots, 16.73\mskip-3mu\dots]$$. This is probably an important point that can be observed in other problems that integrate functions with small regions of support. Solve[SetPrecision[0.56, Infinity] E^( SetPrecision[-12.2, Infinity] (-15 + t)^2) == 2^(-54), t, Reals] /. t_?NumericQ :> N@If[t < 15, Floor, Ceiling][t, 2^(3 - 53)] f[10] - f[t /. %] (* {{t -> 13.262040737546037}, {t -> 16.737959262453963}} {0., 0.} *)  And f[t] varies relatively by less than $$\sqrt{\varepsilon}$$, where $$\varepsilon$$ is the machine epsilon $$2^{-52}$$ for t outside the interval $$[13.78, 16.22]$$. The relative variation $$\pm\sqrt{\varepsilon}$$ is at the boundary of the tolerance prescribed by the default PrecisionGoal. This is not how PrecisionGoal is used by NDSolve[], and I doubt this interval exactly determines where the issue arises (see next para.); however, it gives an idea of where f[t] might be perceived as constant by the error estimator used by NDSolve[]. Namely, f[t] is definitely constant up to t == 13.26, and it is probably treated as constant a little farther, probably not past t == 13.78. And similarly on the other side of t == 15. The issue I referred to (sorry for the order) is that if the integrator steps over the interval, the blip in the variation of f[t] is missed. This is what happens in the OP's code. It is an accident. The accidental nature probably is what leads @xzczd to remark about the StartingStepSize solution, "I fear that this solution may not be that stable." Below are the steps taken by the OP's code. There are 373 in all. The integrator steps from t == 13.3445, which lies in limbo between constant and non-constant as described above, and t == 17.009, where f[t] is constant. However the step-error was calculated, it is clear the step was accepted and the phase shift the OP expected was missed. x["Grid"] /. sol // Flatten // Short[#, 10] & (* {0., 6.25*10^-6, 0.00001875, 0.00003125, 0.00004375, 0.00006875, 0.00009375, 0.00011875, 0.00014375, 0.00019375, 0.00029375, 0.00039375, 0.00059375, 0.00079375, 0.00099375, <<343>>, 9.68007, 10.1381, 10.5962, 11.5123, 13.3445, 17.009, 24.3379, 34.3379, 44.3379, 54.3379, 64.3379, 74.3379, 84.3379, 94.3379, 100.} *)  ## Workaround An easy fix is to restart the integration near the beginning of the blip. It will start with small steps, and it won't jump over the blip: sol = NDSolve[ {func, WhenEvent[t == 13.5, "RestartIntegration"]}, {x}, {t, 0, runTime}, MaxSteps -> \[Infinity]]  Since the system is solved as a DAE and hence must use the IDA method, which is available only in machine precision, workarounds are limited. The OP asked about 3D, which opens other possibilities. ## Formulating as an explicit 3D system We can change to an explicit system rather easily with the replacement rule: xTo3d = x -> Function[t, {x1[t], x2[t], x3[t]}];  The basic call is sol3d = NDSolve[func /. xTo3d, {x1, x2, x3}, {t, 0, runTime}]  It still fails at machine precision, more or less just how it did in the OP's code. However, we can use arbitrary precision. Either of these work: sol3d = NDSolve[Rationalize[Rationalize@func /. xTo3d, 0], {x1, x2, x3}, {t, 0, runTime} , PrecisionGoal -> 8, AccuracyGoal -> 8, WorkingPrecision ->$MachinePrecision];

sol3d = NDSolve[Rationalize[Rationalize@func /. xTo3d, 0],
{x1, x2, x3}, {t, 0, runTime}
, WorkingPrecision -> 16];


Note that the precision and accuracy need to be at least 8, which is interesting. This give me the feeling, as @xzczd expressed it, "I fear that this solution may not be that stable." A higher precision might be more reliable, such as WorkingPrecision -> 24: The higher the precision, the greater the range over which f[t] is non-constant. Note the higher precision method in this case is a bit faster than the machine-precision IDA. This is because IDA uses a nonlinear implicit solver and the 3D version can use an explicit method.

plotTime = runTime/5;
temp = Table[x[t] /. xTo3d /. sol3d[[1]], {t, 0, plotTime, 0.001}];
phi = Arg[temp[[All, 1]] + I*temp[[All, 2]]]; (*Finds phase*)
phi = phi -
Accumulate[
2*Pi*Unitize[
Threshold[Prepend[Differences[phi], 0.],
0.5]]]; (*Unwraps phase*)
ListLinePlot[-phi, AspectRatio -> 0.2,
PlotRange -> {{0, plotTime}, All},