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?


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.


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.


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?


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.

  • 1
    $\begingroup$ Options[NDSolve`BDF] gives {"ImplicitSolver" -> NDSolve`Newton, "MaxDifferenceOrder" -> 5, "VariableStepCoefficients" -> Automatic} for me on v8.0.4 which says that it doesn't take "DifferenceOrder" as an option. $\endgroup$
    – rcollyer
    Commented Mar 28, 2012 at 14:24
  • $\begingroup$ @rcollyer You are right. However, is running IDA with Newtons method and "MaxDifferenceOrder"->5 the same as running BDF with "MaxDifferenceOrder"-> 5? They both run without errors. $\endgroup$
    – dearN
    Commented Mar 28, 2012 at 14:28
  • $\begingroup$ MaxDifferenceOrder->X mean that the maximum difference order is X. NDSolve may, however switch between orders up until X to get to a solution. $\endgroup$
    – user21
    Commented May 15, 2012 at 9:30

1 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}, 
               {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.


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,

InterpolatingFunctionInterpolationOrder @ sol[[1, 1, 2]]

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"]
  • $\begingroup$ This is a really interesting comment. The reason I am interested in fixing the order to a user defined value is that it makes comparison with published data so much more easier. $\endgroup$
    – dearN
    Commented Mar 28, 2012 at 14:59
  • 1
    $\begingroup$ @DNA no joy on the alternative method, but I'll keep poking at it and see if anything comes of it. $\endgroup$
    – rcollyer
    Commented Mar 28, 2012 at 15:58
  • 1
    $\begingroup$ In the name of completion and all that, I was able to reproduce results in this paper by increasing the "MaxStepFraction" to 150. That generated a grid of 151x151 which seemed well and good for the nonlinear partial differential equation as described in the paper.. $\endgroup$
    – dearN
    Commented Mar 29, 2012 at 20:16
  • 1
    $\begingroup$ @DNA already added to the answer. $\endgroup$
    – rcollyer
    Commented Apr 5, 2012 at 15:01
  • 1
    $\begingroup$ @DNA grumble. Fixed now. The issue was ["InterpolationOrder"] was not being applied to the InterpolatingFunction. $\endgroup$
    – rcollyer
    Commented Apr 5, 2012 at 19:05

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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