Which ode solver to use

Because it is a one-step solver, it may be more efficient than ode15s at crude tolerances. It can solve some kinds of stiff problems for which ode15s is not effective. Use this solver if the problem is only moderately stiff and you need a solution without numerical damping.

By construction, the same iteration matrix is used in evaluating both stages. Like ode23s , this solver may be more efficient than ode15s at crude tolerances. See Also deval , odeset , odeget , function handle. References [1] Bank, R. Coughran, Jr.

Fichtner, E. Grosse, D. Rose, and R. CAD , 4 , pp Shampine, "A 3 2 pair of Runge-Kutta formulas," Appl. Letters , Vol. Prince, "A family of embedded Runge-Kutta formulae," J.

Malcolm, and C. Moler, and S. Freeman, San Francisco, Reichelt, and J. A function that evaluates the right-hand side of the differential equations.

A vector specifying the interval of integration, [t0,tf]. Optional integration argument created using the odeset function. For the i th event function in events : value i is the value of the function. The structure sol always includes these fields: sol.

Each column sol. If you specify the Events option and events are detected, sol also includes these fields: sol. Points at which events, if any, occurred. Indices into the vector returned by the function specified in the Events option. Problems with state-dependent mass matrices are more difficult: If the mass matrix does not depend on the state variable and the function MASS is to be called with one input argument, t , set the MStateDependence property to ' none '.

If the mass matrix depends weakly on , set MStateDependence to ' weak ' the default and otherwise, to ' strong '. In either case, the function MASS is called with the two arguments t , y. If there are many differential equations, it is important to exploit sparsity: Return a sparse. Supply the sparsity pattern of using the JPattern property or a sparse using the Jacobian property.

If some components of y ' are missing, then the equations are called differential algebraic equations , or DAEs, and the system of DAEs contains some algebraic variables. Algebraic variables are dependent variables whose derivatives do not appear in the equations. A system of DAEs can be rewritten as an equivalent system of first-order ODEs by taking derivatives of the equations to eliminate the algebraic variables. The ode15s and ode23t solvers can solve index-1 DAEs. Fully implicit ODEs cannot be rewritten in an explicit form, and might also contain some algebraic variables.

The ode15i solver is designed for fully implicit problems, including index-1 DAEs. You can supply additional information to the solver for some types of problems by using the odeset function to create an options structure. You can specify any number of coupled ODE equations to solve, and in principle the number of equations is only limited by available computer memory.

If the system of equations has n equations,. For example, consider the system of two equations. You must rewrite higher-order ODEs as an equivalent system of first-order equations using the generic substitutions. The result of these substitutions is a system of n first-order equations. To solve it, separate the real and imaginary parts into different solution components, then recombine the results at the end.

Conceptually, this looks like. When you run a solver to obtain the solution, the initial condition y0 is also separated into real and imaginary parts to provide an initial condition for each solution component.

Once you obtain the solution, combine the real and imaginary components together to obtain the final result. However, ode23 , ode78 , ode89 and ode can be more efficient than ode45 for problems with looser or tighter accuracy requirements. Some ODE problems exhibit stiffness , or difficulty in evaluation.

Stiffness is a term that defies a precise definition, but in general, stiffness occurs when there is a difference in scaling somewhere in the problem.

For example, if an ODE has two solution components that vary on drastically different time scales, then the equation might be stiff. You can identify a problem as stiff if nonstiff solvers such as ode45 are unable to solve the problem or are extremely slow.

If you observe that a nonstiff solver is very slow, try using a stiff solver such as ode15s instead. When using a stiff solver, you can improve reliability and efficiency by supplying the Jacobian matrix or its sparsity pattern. This table provides general guidelines on when to use each of the different solvers. It also has an implicit RKC method as well. Altogether, there's some interesting stuff here, but the interfaces are not too uniform. Like the other Fortran suite, it has "event handling" in terms of "you can control it step-by-step and have a rootfinder, so you can write your own event handling".

They do hold your hand a little bit more by providing driver functions which do some cool stuff , but this is still definitely the old-school way of having a feature. While this is one solver I haven't used myself, I am sure that its efficient like the other Fortran methods, but it hasn't updated with implicit Runge-Kutta or Rosenbrock methods which are known to be more efficient in many cases in their modern form.

Yes, this has a long history. The interface for these is very similar to the ODEPACK interface, which means you can control it step by step and use the rootfinding capabilities to write your own event handling interface. Many different scientific computing software internally are making use of Sundials because it can handle just about anything. Well, it does have many limitations. There's no stochasticity or delays allowed. But more importantly, Hairer and DifferentialEquations.

Multistep methods are also not very stable at their higher orders, and so at higher tolerances lower accuracy these methods may fail to have their steps converge on standard test problems see this note in the ROBER tests , meaning that you may have to increase the accuracy and thus computational cost due to stiffness issues. But, since multistep methods only require a single function evaluation per step or are implicit in only single steps , they are the most efficient for asymptotically hard problems i.

For this reason, these methods excel at solving large discretizations of PDEs. This means that it covers the methods which I mentioned earlier are more efficient in many of the cases though it is a bit lacking on the explicit tableaus and thus could be more efficient, but that's just details. Thus it should be considered separate from the rest of the "Sundials brand", giving a form of completeness but not with the same level of performance as you would expect.

And if you have an asymtopically large problem or very expensive function evaluations, this will be the most efficient as well. Plus the parallel versions are cool! You do have to live with the limitations of low-level software forcing specific number and array types, along with the fact that you need to write your own event handling, but if you're "hardcore" and writing in a compiled language this suite is a good bet.

Now we come to SciPy's suite. SciPy 1. However, its basic Runge-Kutta integrator is written directly in Python with loops, utilizes an old school timestepping method instead of newer more efficient ones this makes it less stable and more likely to diverge than most implementations , it doesn't have very many options, etc.

So this is a definite step backwards in terms of "hardcore efficiency" and the features for optimizing an RK method to a given problem. However, this is a big step forward because it allows the methods to solve on matrices and tensors instead of just arrays, allows for a few of the methods to have dense output, and more flexibility.

These same changes were made to its BDF method as well. The developer chats seem to say it's around a 10x performance loss for problems with around ODEs, smaller problems getting more overhead and larger ones getting less due to the use of vectorization when possible.

Using a Python interface with classes and the like, it does offer a step-by-step interface. It now has some basic event handling, but for example it only checks the endpoints so if your function isn't monotonic it can easily miss it along with other things like floating point issues which come from not doing a pullback , so it's still at the point where it's a decent undergraduate homework problem to ask them to find a way to break it there are several quick ways.

This in't just theoretical, users have posted more than issue with easy ways to break its event handling , so even for simple cases use caution.

It has dop and dopri5 hidden in an old interface. In the developer's tests these are more efficient, but they don't have these nice extra features like post-solution interpolation. It does have a big options list, so if you were planning on swapping out the linear solver for PETSc you're out of luck, but it does expose some of the internal options of LSODA for banded Jacobians.

It only has ODE solvers, no differential-algebraic, delay, or stochastic solvers. It does include a single BVP solver though. The tableaus are all stored at double precision so even if higher precision numbers are accepted it wouldn't give a higher precision result.

As with the methods analysis, its only higher order Runge-Kutta method for efficient solving of nonstiff equations is dop which is now berried in the legacy interface without extra features and it's missing Rosenbrock and SDIRK methods entirely, opting to only provide the multistep methods.

Since the derivative function is where the ODE solver spends most of its time for sufficiently difficult problems , this means that even though you are calling Fortran code, you will lose out on a lot of efficiency. Still, if efficiency isn't a big deal and you don't need bells and whistles, this suite will do the basics. It handles a much smaller domain and is missing a bunch of methods, but if your problem is simple enough it'll do it and can handle some basic events as well.

It's nice to teach a class with and make some plots, and then you can go grab a Sundials wrapper when you need it. There's not much to say other than that deSolve is very much like SciPy's suite.

It wraps almost the same solvers, has pretty much the same limitations, and has the same efficiency problem since in this case it's calling the user-provided R function most of the time. One advantage is that it does have event handling.

Less vanilla with a little bit more features, but generally the same as SciPy. I do want to make a note about its delay differential equation solvers dede. The documentation specifically says:.

They end with saying it's only for "simple" delay equations without mentioning these as the reasons. I think that's going too far: there is no good guarantee that the error is low or the method converges well. I will mark it as "Poor" in the table, but I would advise against using this method and would recommend the developers to remove this from the docs because this is definitely an easy place for users to be mislead by a numerical result.

PyDSTool is an odd little beast. Part of the software is for analytic continuation i. But another part of it is for ODE solvers. It contains one ODE solver which is written in Python itself and it recommends against actually using this for efficiency reasons. Instead, it wraps some of the Hairer methods, specifically dopri5 and radau, and recommends these. By doing so, it's much more efficient. We still note that its array of available methods is small and it offers radau which is great for high accuracy ODEs and DAEs, but is not the most efficient at lower accuracy so it would've been nice to see Rosenbrock and ESDIRK methods.

It has some basic event handling and methods for DDEs again as wrappers to a Hairer method. This is a pretty good suite if you can get it working, though I do caution that getting the extra non-Python methods setup and compiled is nontrivial. One big point to note is that I find the documentation spectacularly difficult to parse. Together, it's pretty efficient and has a good set of methods which will do basic event handling and solve problems at double precision. I haven't tried it out myself but I'll assume this will get you as efficient as though you used it from Fortran.

However, it is lacking in the features department, not offering advanced things like arbitrary number types, event handling, etc. But if you have a vanilla ODE to solve and you want to easily do it efficiently in Python, this is a good option to look at.

It uses the high order strong order 1. Thus it can be pretty efficient for solving the most standard stiff and nonstiff ODEs. However, its implementations do not make use of advanced timestepping techniques with its explicit RK method PI-controllers which makes it require more steps to achieve the same accuracy as some of the more advanced software, making it not really any more efficient in practice though it does do Gustafsson acceleration in its Rosenbrock scheme to prevent "the Hump" behavior.

It and DifferentialEquations. Thus if you are familiar with templates and really want to make use of them, this might be the library to look at, otherwise you're probably better off looking elsewhere like Sundials.

Actually, they are just kind of weird. When comparing their choices against what is known to be efficient according to the ODE research and benchmarks, the methods that they choose to implement are pretty bizarre like extrapolation methods which have repeatedly been shown to not be very efficient, while not included other more important methods.

But they do have some of the basics like multistep Adams and BDF methods. This benchmark , while not measuring runtime and only uses function evaluations which can be very very different according to more sophisticated benchmarks like the Hairer and DifferentialEquations. Mathematica's ODE solvers are very sophisticated.

These methods can be almost 5x faster than the older high order explicit RK methods which themselves are the most efficient class of methods for many nonstiff ODEs, and thus these do quite well.

It includes interpolations to make their solutions continuous functions that plot nicely. Its native methods are able to make full use of Mathematica and its arbitrary precision, but sadly most of the methods it uses are wrappers for the classic solvers.

It uses a method of steps implementation over its explicit Runge-Kutta methods for solving nonstiff DDEs efficiently, and includes high order Runge-Kutta methods for stochastic differential equations though it doesn't do adaptive timestepping in this case.

One nice feature is that all solutions come with an interpolation to make them continuous. Additionally, it can use the symbolic form of the user-provided equation in order to create a function for the Jacobian, and use that instead of a numerical differentiation fallback like all of the previous methods in the stiff solvers for more accuracy and efficiency.

It also has symplectic methods for solving second order ODEs, and its event handling is very expressive. It's very impressive, but since it's using a lot of wrapped methods you cannot always make use of Mathematica's arbitrary precision inside of these numerical methods.

Additionally, its interface doesn't give you control over all of the fine details of the solvers that it wraps. Maple's ODE solvers has the basics.

It has a classic higher order explicit RK method due to Verner not the newer more efficient ones though , and for stiff problems it uses a high order Rosenbrock method. Its native methods can make use of arbitrary arithmetic and complex numbers. It has a BVP solver, but it explicitly says that it's not for stiff equations or ones with singularities it relies on deferred correction or extrapolation. It can solve state-dependent delay differential equations with its explicit RK and Rosenbrock methods.

It has event handling which doesn't seem to be well documented. As another symbolic programming language, it computes Jacobians analytically to pass to the stiff solvers like Mathematica, helping it out in that respect. Thus while some of its basics are a little odd why a Fehlberg method? The Dormand-Prince paper is one of the most famous papers about increasing efficiency in explicit RK tableaus and not as fleshed out, it does have a good variety with quite a bit of capabilities.

I would actually go as far as to say that, if I had to create a suite with a small set of solvers I'd probably make something similar to what Maple has. It matches what seems to benchmark best "for most problems in normal tolerances", which is what most people are probably solving. Thus it has something that's pretty efficient for pretty much every case. What makes it special is that it includes the ability to do sensitivity analysis calcuations.

It can't do anything but double-precision numbers and doesn't have event handling, but the sensitivity calculations makes it quite special. Okay, time for DifferentialEquations. I left it for last because it is by far the most complex of the solver suites, and pulls ideas from many of them. While most of the other suite offer no more than about 15 methods on the high end with most offering about 8 or less , DifferentialEquations.

Thus all of the standard methods mentioned before are available in this suite. But then there are the native Julia methods. In each of these categories it has a huge amount of variety, offering pretty much every method from the other suites along with some unique methods. Some unique methods to point out are that it has the only 5th order Rosenbrock method, it has the efficient Verner methods discussed in the Mathematica part, it has newer 5th order methods i.

It has methods specialized to reduce interpolation error the OwrenZen methods , and methods which are strong-stability preserving for hyperbolic PDEs. It by default the solutions as continuous functions via a high order interpolation though this can be turned off to make the saving more efficient.

Each of these implementations employ a lot of extra tricks for efficiency. For example, the interpolation is "lazy", meaning that if the method requires extra function evaluations for the continuous form, it will only do those extra calculations when the continuous function is used so when you directly ask for it or when you plot.

This is just a peek at the special things the library does to gain an edge. The native Julia methods benchmark very well as well, and all of the benchmarks are openly available. They even allow you to tweak a lot of the internals and swap out the linear algebra routines to use parallel solvers like PETSc. Its event handling is the most advanced out of all of the offerings. You can change just about anything. You can make your ODE do things like change its size during the solving if you want, and you can make the event handling adjust and adapt internal solver parameters.

It's not a hyperbole to say that the user is given "full control" since the differential equation solvers themselves are written as essentially a method on the event handling interface, meaning that anything it can do internally you can do. The variety of methods is also much more diverse than the other offerings. It has symplectic integrators like Harier's suite, but has more high and low order methods.

It has DDE solvers for constant-lag and state-dependent delays, and it has stiff solvers for each of these cases. The stiff solvers also all allow for solving DAEs in mass matrix form though fully implicit ODEs are possible, but can only be solved using a few methods like a wrapper to Sundials' IDA and doesn't include event handling here quite yet. It allows arbitrary numbers and arrays like Boost. This means you can use arbitrary precision numbers in the native Julia methods, or you can use "array-like" objects to model multiscale biological organisms instead of always having to use simple contiguous arrays.

It has addons for things like sensitivity analysis and parameter estimation. It hits the strong points of each of the previously mentioned suites because it was designed from the get-go to do so.

And it benchmarks extremely well. The only downside is that, because it is so feature and performance focused, its documentation is heavy. The beginning tutorial will give you the basics making it hopefully as easy as MATLAB to pick up , but pretty much every control knob is available, making the extended portion of the documentation a long read. Let's take a step back and summarize this information. Since it also features solvers for the non-ordinary differential equations and its unique Julia methods also benchmarks well, I think it's clear that DifferentialEquations.

As for the other choices from scripting languages, MATLAB wasn't designed to have all of the most efficient methods, but it'll handle basic equations with delays and events and output good plots.

Mathematica and Maple will do symbolic pre-calculations to speed things up and can JiT compile functions, along with offering pretty good event handling, and thus their wrappers are more like DifferentialEquations. So in a pinch when not all of the bells and whistles are necessary, each of these scripting language suites will get you by. Behind DifferentialEquations. Still, you're going to have to write a lot of stuff yourself to make the rootfinding into an event handling interface, but if you put the work in this will do you well.

Any of these will do you well if you want to really get down-and-dirty in a compiled language and write a lot of the interfaces yourself, but they will be a sacrifice in productivity with no clear performance gain over the scripting language methods which also include some form of JIT compilation.

Shortly after being posted this has received quite a bit of buzz. I am thankful so many people are interested! A few things were pointed out to me that I have since corrected:. There are probably a few other corrections needed. For one, I may have put PyDSTool too high in some areas like state-dependent delay solvers I had notes that say it can't do it, but I can't find it in the docs anymore More edits today.

This time for Maple. Some more refinements.