I'm studying some new numerical integration algorithms using NDSolve Method Framework to specify them. When I try to set WorkingPrecision option to other value than MachinePrecision I get the following error message:

NDSolve::bdstep: The Step function for Method -> CRK4 returned {0.03529465924388010,0.07025569022824054}, which is not an acceptable form. >>

Here I apply the first CRK4 definition from the official tutorial for the Van der Pol equation:

system = 
   {{Y1'[T] == Y2'[T], Y2'[T] == -Y1[T] + 2 (1 - Y1[T]^2) Y2[T]}, 
    {Y1[0] == 0, Y2[0] == 2}, {Y1[T], Y2[T]}, {T, 0, 30}, {}, {}, {}}];
vars = system["DependentVariables"]; 
time = system["TimeData"];
sol1 = 
  Flatten[vars /. NDSolve[system, Method -> CRK4, WorkingPrecision -> 16]] 

I have little experience working with Wolfram Mathematica and I will be very grateful to get a clue how to fix this issue.


Change the definitions so that the step output is a list:

CRK4[]["Step"[rhs_, h_, t_, x_, xp_]] := Module[{k0, k1, k2, k3},
  k0 = h xp;
  k1 = h rhs[t + h/2, x + k0/2];
  k2 = h rhs[t + h/2, x + k1/2];
  k3 = h rhs[t + h, x + k2];
  {(k0 + 2  k1 + 2  k2 + k3)/6}]        (* <-- List *)

CRK4[___]["StepInput"] = {"F"["T", "X"], "H", "T", "X", "XP"};
CRK4[___]["StepOutput"] = {"XI"};       (* <-- List *)

Then both the MachinePrecision and WorkingPrecision -> 16 code return comparable solutions:

sol1mp = Flatten[vars /. NDSolve[system, Method -> CRK4, StartingStepSize -> 0.01]];
sol1 = Flatten[vars /. NDSolve[system, Method -> CRK4, StartingStepSize -> 0.01, 
     WorkingPrecision -> 16]];

Plot[sol1mp - sol1 // Evaluate, Evaluate@system["TimeData"], PlotRange -> All]

Mathematica graphics

Now, why does it work? Well, the check takes place in internal code that cannot be traced. So I don't know. I don't know if it's a bug or an error in the example. There is quite a bit of documentation on NDSolve[], but there seems to be a lot that is not explained, too. There may be a reason for wanting the two return types.

I can share this tip. If you use the form

CRK4[___]["StepOutput"] = "XI";

it seems the output of the "Step" method needs to be a packed array or you'll get the error. That means machine precision only. If the output form is {"XI"}, then NDSolve seems to tolerate any numeric array of the appropriate dimensions.

  • $\begingroup$ Thanks a lot! That solves the problem. $\endgroup$ – Sheldon Nov 26 '15 at 14:03
  • $\begingroup$ @Sheldon You're welcome. $\endgroup$ – Michael E2 Nov 26 '15 at 14:05

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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