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.
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).
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.