Step problem in NDSolve with FixedStep method

Question

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?

Answer 1

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.