# Step problem in NDSolve with FixedStep method

We want to solve a set of ODEs with fixed step, whose size is 1/1000, the related codes are given in the following.

The parameters are given by

δt = 1/1000;
finalTime = 5000;

γ = -2 Pi 1.18 10^7;
T10 = 1/0.111;
T20 = 1/0.148;
γse = 0.00248;
PRb = 0.8;
a = (γse + 1/T10)/T20;
b = PRb γse;

ωz[B0_] := γ (B0)
d = T20 (γse + 1/T10);

α = 170 a/b;

setδT2 = ((-3)/(T20 γ) 10^9)/2;
dataB0 = (750 + {-setδT2, setδT2}) 10^-9;


Inintial states are given by

P1x0 = 0.001 Sin[Pi/5 ];
P1y0 = 0;
P1z0 = 0.001 Cos[Pi/5];
P2x0 = 0.001 Sin[Pi/5 ];
P2y0 = 0;
P2z0 = 0.001 Cos[Pi/5];


The equations are

dPxii1dt[Pxii1_, Pyii1_, Pzii1_, Pxii2_, Pyii2_, Pzii2_,
B0ii_] := (Pyii1 ωz[B0ii] +
Pzii1  α (Pxii1 + Pxii2)/2 - Pxii1/T20);
dPyii1dt[Pxii1_, Pyii1_, Pzii1_, Pxii2_, Pyii2_, Pzii2_,
B0ii_] := (Pzii1  α (Pyii1 + Pyii2)/2 -
Pxii1 ωz[B0ii] - Pyii1/T20);
dPzii1dt[Pxii1_, Pyii1_, Pzii1_, Pxii2_, Pyii2_, Pzii2_,
B0ii_] := (-Pxii1  α (Pxii1 + Pxii2)/2 -
Pyii1 α (Pyii1 + Pyii2)/2 -
Pzii1/T10 + γse (PRb - Pzii1)) ;
dPxii2dt[Pxii1_, Pyii1_, Pzii1_, Pxii2_, Pyii2_, Pzii2_,
B0ii_] := (Pyii2 ωz[B0ii] +
Pzii2  α (Pxii1 + Pxii2)/2 - Pxii2/T20);
dPyii2dt[Pxii1_, Pyii1_, Pzii1_, Pxii2_, Pyii2_, Pzii2_,
B0ii_] := (Pzii2  α (Pyii1 + Pyii2)/2 -
Pxii2 ωz[B0ii] - Pyii2/T20) ;
dPzii2dt[Pxii1_, Pyii1_, Pzii1_, Pxii2_, Pyii2_, Pzii2_,
B0ii_] := (-Pxii2  α (Pxii1 + Pxii2)/2 -
Pyii2 α (Pyii1 + Pyii2)/2 -
Pzii2/T10 + γse (PRb - Pzii2));
solvEquation = {Pxii1'[t] ==
dPxii1dt[Pxii1[t], Pyii1[t], Pzii1[t], Pxii2[t], Pyii2[t],
Pzii2[t], dataB0[[1]]],
Pyii1'[t] ==
dPyii1dt[Pxii1[t], Pyii1[t], Pzii1[t], Pxii2[t], Pyii2[t],
Pzii2[t], dataB0[[1]]],
Pzii1'[t] ==
dPzii1dt[Pxii1[t], Pyii1[t], Pzii1[t], Pxii2[t], Pyii2[t],
Pzii2[t], dataB0[[1]]],
Pxii2'[t] ==
dPxii2dt[Pxii1[t], Pyii1[t], Pzii1[t], Pxii2[t], Pyii2[t],
Pzii2[t], dataB0[[2]]],
Pyii2'[t] ==
dPyii2dt[Pxii1[t], Pyii1[t], Pzii1[t], Pxii2[t], Pyii2[t],
Pzii2[t], dataB0[[2]]],
Pzii2'[t] ==
dPzii2dt[Pxii1[t], Pyii1[t], Pzii1[t], Pxii2[t], Pyii2[t],
Pzii2[t], dataB0[[2]]], Pxii1[0] == P1x0, Pyii1[0] == P1y0,
Pzii1[0] == P1z0, Pxii2[0] == P2x0, Pyii2[0] == P2y0,
Pzii2[0] == P2z0};


Solving with

DOPRIamat = {{1/5}, {3/40, 9/40}, {44/45, -56/15,
32/9}, {19372/6561, -25360/2187,
64448/6561, -212/729}, {9017/3168, -355/33, 46732/5247,
49/176, -5103/18656}, {35/384, 0, 500/1113, 125/192, -2187/6784,
11/84}};
DOPRIbvec = {35/384, 0, 500/1113, 125/192, -2187/6784, 11/84, 0};
DOPRIcvec = {1/5, 3/10, 4/5, 8/9, 1, 1};
DOPRIevec = {71/57600, 0, -71/16695, 71/1920, -17253/339200,
22/525, -1/40};
DOPRICoefficients[5, p_] :=
N[{DOPRIamat, DOPRIbvec, DOPRIcvec, DOPRIevec}, p];

solvPxyz =
NDSolve[solvEquation, {Pxii1, Pyii1, Pzii1, Pxii2, Pyii2,
Pzii2}, {t, 0, finalTime},
Method -> {"FixedStep",
Method -> {"ExplicitRungeKutta", "DifferenceOrder" -> 5,
"Coefficients" -> DOPRICoefficients,
"StiffnessTest" -> False}}, StartingStepSize -> δt] //
Flatten;


As shown above, δt = 1/1000;finalTime = 5000; , the length of calculated grid should be 5000001. However,Length[(Pxii1 /. solvPxyz)["Grid"]] gives

5000002


And we have checked that

(Pxii1 /. solvPxyz)["Grid"][[-10 ;; -1]]


outputs

 {{4999.992}, {4999.993}, {4999.994}, {4999.995}, {4999.996}, {4999.997},{4999.9980000000005}, {4999.999}, {4999.9995}, {5000.}}


It can be seen that the extra grid is {4999.9995}. What is the reason here? And how to fix this problem?

The problem can be reproduced with a much simpler sample:

δt = 1/10000;
num = 1000;
tst = NDSolveValue[{x'[t] == 1, x[0] == 0}, x, {t, 0, num δt},
StartingStepSize -> δt, Method -> {"FixedStep", Method -> Automatic}];
Length@tst["Grid"]
(* 1002 *)
tst["Grid"][[-5 ;; All]]
(* {{0.0997}, {0.0998}, {0.0999}, {0.09995}, {0.1}} *)


I'm not sure about why FixedStep method behaves like this (better to contact WRI, I think), but setting WorkingPrecision option resolves the issue:

δt = 1/10000;
num = 1000;
tst = NDSolveValue[{x'[t] == 1, x[0] == 0}, x, {t, 0, num δt},
StartingStepSize -> δt, Method -> {"FixedStep", Method -> Automatic},
WorkingPrecision -> 16];
Length@tst["Grid"]
(* 1001 *)
tst["Grid"][[-5 ;; All]]
(* {{0.09960000000000000}, {0.09970000000000000}, {0.09980000000000000},
{0.09990000000000000}, {0.1000000000000000}} *)


For your original code, setting WorkingPrecision -> 16 will result in precw and mxst warning. The first warning is benign, you can avoid it with Rationalize[…, 0] or SetPrecision[…, Infinity] anyway. The second warning can be resolved with MaxSteps option. (It's a bit surprising to me the mxst warning doesn't show up when the WorkingPrecision is the default MachinePrecision. ) Still, I'll use the sample above to show the usage of Rationalize and WorkingPrecision:

δt = 1/10000;
num = 10000 + 3;
tst = NDSolveValue[{x'[t] == 1., x[0] == 0} // Rationalize[#, 0] &,
x, {t, 0, num δt}, StartingStepSize -> δt,
Method -> {"FixedStep", Method -> Automatic}, MaxStepFraction -> 1/( num),
MaxStepSize -> δt, WorkingPrecision -> 16, MaxSteps -> Infinity];
Length@tst["Grid"]
tst["Grid"][[-5 ;; All]]


Finally, I'd like to emphasize again, the default ODE solver of NDSolve is quite robust, and one should not manually turn to FixedStep method unless there is some special reason.

• Thanks for your reply! @xzczd However, WorkingPrecision -> 16 does not work for my case. The errors NDSolve::precw: The precision of the differential equation ({<<1>>}) is less than WorkingPrecision (16.).and NDSolve::mxst: Maximum number of 10000 steps reached at the point t == 10.000999999999999999999999999999999986616.. are outputted. Commented Nov 11, 2023 at 6:52
• @so_sure The first warning isn't a problem, if you don't want to see it, wrap your system with Rationalize[…, 0] or SetPrecision[…, Infinity]. As to the second warning, add MaxSteps->Infinity. (It's a bit surprising to me that your original code doesn't need manual adjustion for MaxSteps. ) Commented Nov 11, 2023 at 7:01
• @so_sure See the example in my updated answer if you still feel confused. Commented Nov 11, 2023 at 7:09
• It woks well, but consumes too much time. ╮（╯＿╰）╭ Commented Nov 11, 2023 at 7:54
• @so_sure Performance tuning is another topic. And as mentioned at the end of my question, the default ODE solver of NDSolve is quite robust, one should not manually turn to FixedStep method unless there is some special reason, related: mathematica.stackexchange.com/a/274162/1871 Commented Nov 11, 2023 at 8:00