9
$\begingroup$

I want to solve a relatively small system of stiff ODEs (< 10 first-order equations) using high precision floating point arithmetic (using MPFR or alike). What would be the easiest algorithm to implement/port? I'm not only interested in the "stepper" algorithm, but in the combination of stepper/error estimation/step size control.

Improve this question
$\endgroup$

6 Answers 6

Reset to default
4
$\begingroup$

I don't know if your problem would be tractable using a high order Taylor series method, but if so, and you're comfortable with implementing a solution to your problem in Python, then you could try a combination of mpmath, SymPy, and possibly pytaylor (which is not as well-established as the first two Python modules, but does implement 4th order RK in addition to Taylor series methods). Based on a quick read of a paper by Nedialkov in BIT, it seems as though a sufficiently high Taylor series (of order 20-30) could be used for highly accurate, efficient numerical integration of mildly to moderately stiff problems. (Nedialkov studies the differential-algebraic case; presumably, these results must also apply to ODEs.)

Whatever solution method you choose, I encourage you open-source your work, since it sounds like it'd be a useful contribution to the ODE software landscape.

Improve this answer
$\endgroup$
2
  • $\begingroup$ I've decided to look into the Taylor methods; probably using mpmath. Thanks to all for the quick and good answers. $\endgroup$
    – GertVdE
    Commented Jun 21, 2012 at 18:49
  • $\begingroup$ The link to the paper by Nedialkov at springerlink.com is broken. I'm also unable to find any copy saved on the Wayback Machine. $\endgroup$
    – user43608
    Commented Jul 22, 2022 at 5:09
7
$\begingroup$

This may be a bit outdated, but Hairer and Wanner's book recommends their own radau5 code. In my own experience, this code is both extremely robust and efficient. It is also a rather straight-forward Runge-Kutta scheme, and thus not too difficult to implement.

There's a Matlab implementation of radau5 by Christian Engstler somewhere out there.

Improve this answer
$\endgroup$
5
  • $\begingroup$ Thanks Pedro. The RADAU code by H&W is my "reference" code now :-) I've been very happy using it but I wasn't sure I was up to porting it... $\endgroup$
    – GertVdE
    Commented Jun 17, 2012 at 19:58
  • $\begingroup$ Pedro, since the ODE's are stiff, wouldn't this require an extremely small time step? Mightn't it be more advisable to use an implicit method instead? $\endgroup$
    – Paul
    Commented Jun 17, 2012 at 20:07
  • $\begingroup$ @Paul: radau5 uses, if I recall correctly, an implicit Runge-Kutta scheme. In any case, we were using it quite successfully for stiff ODEs. $\endgroup$
    – Pedro
    Commented Jun 17, 2012 at 20:11
  • 2
    $\begingroup$ radau5 is a fully implicit Runge-Kutta method using 3-point Radau quadrature to provide fifth order accuracy. Hairer also wrote a variable-order code named radau that can provide much higher order. SDC is also not a bad choice if you need extremely high order accuracy. $\endgroup$
    – Jed Brown
    Commented Jun 18, 2012 at 5:03
  • $\begingroup$ I know people who have had success hacking DASSL and DVODE, and of course, DVODE was ported to C to become CVODE, but I don't know that it makes those codes "simple" to port. Since I work with extremely stiff systems in combustion, I tend to be partial to BDF methods over IRK methods, because that's what that community uses. $\endgroup$
    – Geoff Oxberry
    Commented Jun 18, 2012 at 13:06
6
$\begingroup$

Obviously, using multiple precision data types like those in MPFR only makes sense if numerical roundoff is at least around the same order of magnitude as the discretization error. You can only achieve this with high order integrators, so anything that has less than, say, 4th order will not be a useful candidate, and even higher order is better.

I'm not an expert in ODEs to recommend an integration package, but the considerations above may narrow your choice.

Improve this answer
$\endgroup$
6
$\begingroup$

I prefer to always take a benchmark-informed approach to this. Look at the SciML Benchmarks:

Generally, the Rodas5 in Julia's OrdinaryDiffEq.jl algorithm (mixed with automatic differentiation for the Jacobian) performs best in this territory (or Rodas5P). Note it already allows for MPFR bigfloats for arbitrary precision, so it could be used for this without modification just by changing the initial condition and time span to big. For example:

using OrdinaryDiffEq
function lorenz(du,u,p,t) # just for show
 du[1] = 10.0(u[2]-u[1])
 du[2] = u[1]*(28.0-u[3]) - u[2]
 du[3] = u[1]*u[2] - (8/3)*u[3]
end
u0 = [1.0;0.0;0.0]
tspan = (0.0,100.0)
prob = ODEProblem(lorenz,big.(u0),big.(tspan))
sol = solve(prob,Rodas5())

MPFR bigfloats are slower than things like ArbFloats and 128-bit floats, both of which are compatible as well, so I'd suggest looking into some of these alternative extended float types if you really want to optimize it.

BDF methods like SUNDIALS CVODE, DASSL, etc. have difficulty scaling down to there (for DASSL, see the DAE benchmarks like this one.

Radau is the only other one that comes close in benchmarks, because it has a pretty good convergence rate.

enter image description here

Note that only the adaptive order one is competitive, radau5 is not. There is currently a pure Julia RadauIIA5 which outperforms the Fortran one and is compatible with high precision arithmetic (in the same way as above). You'll see it in the benchmarks, but the benchmarks show that RadauIIA5 > radau5 but scaling needs adaptive radau, which we don't currently have but wouldn't be too hard to implement. If you need it sooner, we could look into it.

The implicit extrapolation methods (like ImplicitEulerBarycentricExtrapolation) have very good scaling to low tolerances because they can make use of multithreading as the problem gets harder and the order increases. You'll see they start to do very well at the edge of the benchmark, and benchmarking at lower tolerances might reveal them to be the most efficient in that category. So they are worth a shot, but note you will need to do things like tweak the initial and upper bound on the orders to really get them right. I think they can also receive a few general code optimizations to boost them still too: they are quite new.

I don't know if your problem would be tractable using a high order Taylor series method

Generally Taylor methods are not stable for highly stiff equations. At least, I haven't seen one that has a large stability region, the simple idea is a simple explicit extrapolation which has a not so fantastic stability region (on par with RK methods).

Improve this answer
$\endgroup$
3
$\begingroup$

There are existing implementations, though I don't know of any that are publicly available. You could ask Glaser and Rohklin for their implementation of the methods described in this paper.

Improve this answer
$\endgroup$
2
$\begingroup$

  1. Is DLSODE an option?

  2. Is the system linear? If so, DGPADM is probably the best.

They are relatively self-contained, so not too torturous to re-code if that's what you need to do.

Improve this answer
$\endgroup$

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.