Tell me more ×
Computational Science Stack Exchange is a question and answer site for scientists using computers to solve scientific problems. It's 100% free, no registration required.
up vote 6 down vote favorite

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.

share|improve this question

5 Answers

up vote 1 down vote accepted

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.

share|improve this answer
I've decided to look into the Taylor methods; probably using mpmath. Thanks to all for the quick and good answers. – GertVdE Jun 21 '12 at 18:49
up vote 6 down vote

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.

share|improve this answer
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... – GertVdE Jun 17 '12 at 19:58
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? – Paul♦ Jun 17 '12 at 20:07
@Paul: radau5 uses, if I recall correctly, an implicit Runge-Kutta scheme. In any case, we were using it quite successfully for stiff ODEs. – Pedro Jun 17 '12 at 20:11
2  
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. – Jed Brown Jun 18 '12 at 5:03
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. – Geoff Oxberry♦ Jun 18 '12 at 13:06
up vote 3 down vote

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.

share|improve this answer
up vote 3 down vote

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.

share|improve this answer
up vote 2 down vote

  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.

share|improve this answer

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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