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.
Fastest numeric method for ODE
-
0Also, FYI this question is probably better suited for either [SciComp](http://scicomp.stackexchange.com) or [Stack Overflow](http://stackoverflow.com). 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. – 2012-07-21
2 Answers
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.
- 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$.)
- 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.
- 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).
-
0Have you profiled the code? I have had many surprises with c++. (Try inlining *set_torus*, for example.) – 2012-07-22
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.