Fastest numeric method for ODE

I am currently using Euler method and it is 4-8 times too slow. Which method will be fastest? I need it to compute Turing's reaction-diffusion system.

• Are you sure it is Euler's method that is slowing you down, and not the implementation? What language are you using? Is it a built-in or library function? What is your data set like? What sort of machine are you performing the calculation on, and with how much RAM? When you say it is 4-8 times too slow, what is this relative to and what timing benchmarks have led you to this conclusion? – ely Jul 21 '12 at 22:18
• What do you mean it is '4-8 times too slow'? Compared with what? – copper.hat Jul 21 '12 at 22:30
• @EMS, cooper.hat: I run it on iPad and iPod. It is implemented in C. By the 4-8 times I mean: user can wait 1-2 seconds for computation of the result, but it takes 8-12 seconds with the Euler method. The data set is 256x256 array of doubles, being direct source of final image. If I do 100 iterations, then there are 100 * 256 * 256 loop executions. – AllCoder Jul 21 '12 at 22:37
• Forward Euler, backward Euler? Have you tried a large step size? – copper.hat Jul 21 '12 at 22:41
• Also, FYI this question is probably better suited for either SciComp or Stack Overflow. It's not clear that this is due to your chosen algorithm. It could easily be due to an inefficient C implementation. If you post some code along with the question at either of those other sites, you'll almost surely get better help... and folks on those sites will also be able to suggest other mathematical algorithms to try if indeed that is the right way to solve the problem. – ely Jul 21 '12 at 22:58

Euler's method is the fastest possible single-step, explicit, non-adaptive method for a given fixed step size $h$.

If this is too slow, I would do the following things, in this order.

1. Perform a fundamental analysis on your code. Insert breakpoints. Ensure that you're not doing things like repeatedly re-sizing arrays. Replace loops with better code (a first-order explicit ODE solver should have one "loop": one walk through the code, from $t=0$ to $t=T$.)
2. Increase step size. Compute what your acceptable error tolerance is, and make $h$ be as large as possible to keep your solution within that error bound.
3. Implement an adaptive-step size solver. These are usually faster for a given accuracy, but are more complicated. Dormand-Prince 4(5) pairs is a standard adaptive RK method.

I'm willing to bet that you can solve your speed problem with step (1) (hint: if you're iterating through 256x256 elements, you don't need to necessarily encapsulate the method in a double loop).

Fellow systems biologist studying reaction-diffusion systems here. There is no fastest numerical method for ODEs, but here's some things to try.

Turing's reaction-diffusion system is actually a PDE, but you can represent it as a system of ODEs by converting the LaPlacian into a matrix (the tridiagonal matrix with tridiags (1,-2,1)). This is an order 2 approximation in space, i.e. halving $\Delta x$ takes a quarter out of your error. Once you have it as a system of ODEs like this, you can apply any numerical ODE method.

What numerical ODE method should you use? Well, first try Dormand-Prince 4(5) (in MATLAB this is ode45). However, if your reactions are stiff (i.e. some reactions are very fast) or you have a large diffusion coefficient, then Runge-Kutta methods (like ode45 and Euler's method) will be required to take REALLY small timesteps. This is most likely the problem you're having. In this case, you need to use a solver for stiff equations. In MATLAB, a built-in solver for stiff equations is ode15s. If that doesn't work, then you may need to use an implicit solver like Implicit Euler or the Midpoint Method.

What I described above is known as the Method of Lines, i.e. convert the PDE into a system of ODEs. However, for some very stiff problems this won't work well. You then have to try other numerical PDE methods like Crank-Nicholson type methods.