Setting the DifferenceOrder Option

Question

I've been playing around with Method in NDSolve[...] and can't quite seem to figure out how to force NDSolve to use a certain difference order with a certain method.

For instance, I am trying to validate my work against results of solving an ODE with the Backward Euler method which essentially is BDF of order 1.

Now, Mathematica does have the BDF method for solving stiff equations but when I set the following:

Method->{"BDF", "DifferenceOrder"->1}

I get the following error:

NDSolve::moptx: Method option DifferenceOrder in {NDSolve`BDF,{DifferenceOrder->Automatic}} is not one of {ImplicitSolver,MaxDifferenceOrder,VariableStepCoefficients}. >>

NDSolve::initf: The initialization of the method NDSolve`BDF failed. >>

But when I run NDSolve with just Method->{"BDF"}, lo and behold, I get my solution.

Just to make sure that BDF does take DifferenceOrder options, I checked with

Options[NDSolve`BDF] as prompted in a previous question.

Is this not how I should set the order for BDF? What am I missing?

EDIT 1:

I realized I should be using the IDA method as that apparently uses the BDF method to solve DAEs. Still trying to figure out how I might go about doing this though.

EDIT 2:

I think I might have figured it out. This is what I changed in Method

Method -> {"IDA", "ImplicitSolver" -> "Newton", "MaxDifferenceOrder" -> 5},

From IDA's documentation provided in Advanced Numerical Differential Equation Solving in Mathematica, I think with IDA with Newton Implicit Solver of Order 5, I am actually using a BDF5 method.

EDIT 3:

Per rcollyer's comment, it should be "MaxDifferenceOrder" and NOT "DifferenceOrder" in case BDF is used as the Method

Example: Method->{"BDF", "MaxDifferenceOrder->5"} executes BDF5.

Any thoughts? Comments? Anything else that I might look for, for say other methods such as the LSODA etc?

EDIT 4:

As per rcollyer's comment, BDF with "MaxDifferenceOrder"->5, NDSolve[..] changes the order as required UP TO 5 and doesn't just use order 5.

Answer 1

This is not an answer, but it is really too long for a comment.

I don't know if IDA is the same as BDF, I would have to look closer. However, as you wish to have explicit control over the order of the method and that is not available, it is likely that you can extract what order it does use out of the InterpolationFunction returned.

If you set InterpolationOrder -> All on NDSolve, it's supposed to use the underlying method's order as the interpolation order. Experimenting with

Interpolation[..., InterpolationOrder -> ...]

reveals you may be able to extract the order of interpolation used. For instance,

Interpolation[..., InterpolationOrder -> 2][[2]]

gives something like

{4, 7, 0, {11}, {3}, 0, 0, 0, 0, Automatic}

The fifth term, {3} in this case, is always one more than the order. Using that technique on

sol = NDSolve[ {y''[t] + y'[t] + y[t] == 0, y'[0] == y[0] == 1}, 
               y, 
               {t, 0, 2 \[Pi]}, 
               Method -> {"BDF"}, 
               InterpolationOrder -> All
      ]

via sol[[1, 1, 2, 2]] gives

{4, 17, 5, {119}, {4}, {InterpolatingFunction[{{0., 6.28319}}, "<>"]}, 
   0, 0, 0, Automatic}

which implies that it uses a 3$^\text{rd}$ order method.

---

Edit:

As I intimated, the above method seems workable, but if the format of InterpolatingFunction changes, it will fail. However, there is a package specifically designed to extract this type of information from an InterpolatingFunction: DifferentialEquations`InterpolatingFunctionAnatomy`. Using sol from above,

Needs["DifferentialEquations`InterpolatingFunctionAnatomy`"];
InterpolatingFunctionInterpolationOrder @ sol[[1, 1, 2]]

{3}

which confirms what the other method came up with.

Edit 2:

In addition to the package, above, some opaque objects like SparseArray and InterpolatingFunction can have their properties queried. Thanks to Oleksandr, for InterpolatingFunction these are {"Domain", "Coordinates", "DerivativeOrder", "InterpolationOrder", "Grid", "ValuesOnGrid"}. So, then you can determine the order used, via

sol[[1, 1, 2]]["InterpolationOrder"]