WorkingPrecision and NDSolve Method Plugin Framework

Question

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:

Needs["DifferentialEquations`NDSolveProblems`"];
system = 
 NDSolveProblem[
   {{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.

Answer 1

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]

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.