# Explicit Euler method too slow for reaction-diffusion problem

I am solving Turing's reaction-diffusion system with following C++ code. It is too slow: for 128x128 pixel texture, acceptable number of iterations is 200 – which results in 2.5 seconds of delay. I need 400 iterations to obtain interesting image – but 5 seconds of waiting is too much. Also, size of the texture should be in fact 512x512 – but this results in huge waiting time. The devices are iPad, iPod.

Is there any chance to do this faster? Euler method converges slowly (wikipedia) – having quicker method would allow to drop number of iterations?

EDIT: As Thomas Klimpel pointed out, the lines: "if( m_An[i][j] < 0.0 ) { ... }", "if( m_Bn[i][j] < 0.0 ) { ... }" are delaying convergence: after removing, meaningful image appears after 75 iterations. I have commented out the lines in code below.

void TuringSystem::solve( int iterations, double CA, double CB ) {
m_iterations = iterations;
m_CA = CA;
m_CB = CB;

solveProcess();
}

void set_torus( int & x_plus1, int & x_minus1, int x, int size ) {
// Wrap "edges"
x_plus1 = x+1;
x_minus1 = x-1;
if( x == size - 1 ) { x_plus1 = 0; }
if( x == 0 ) { x_minus1 = size - 1; }
}

void TuringSystem::solveProcess() {
double DiA, ReA, DiB, ReB;

// uses Euler's method to solve the diff eqns
for( n=0; n < m_iterations; ++n ) {
for( i=0; i < m_height; ++i ) {

for( j=0; j < m_width; ++j ) {

// Component A
DiA = m_CA * ( m_Ao[i_add1][j] - 2.0 * m_Ao[i][j] + m_Ao[i_sub1][j]   +   m_Ao[i][j_add1] - 2.0 * m_Ao[i][j] + m_Ao[i][j_sub1] );
ReA = m_Ao[i][j] * m_Bo[i][j] - m_Ao[i][j] - 12.0;
m_An[i][j] = m_Ao[i][j] + 0.01 * (ReA + DiA);
// if( m_An[i][j] < 0.0 ) { m_An[i][j] = 0.0; }

// Component B
DiB = m_CB * ( m_Bo[i_add1][j] - 2.0 * m_Bo[i][j] + m_Bo[i_sub1][j]   +   m_Bo[i][j_add1] - 2.0 * m_Bo[i][j] + m_Bo[i][j_sub1] );
ReB = 16.0 - m_Ao[i][j] * m_Bo[i][j];
m_Bn[i][j] = m_Bo[i][j] + 0.01 * (ReB + DiB);
// if( m_Bn[i][j] < 0.0 ) { m_Bn[i][j]=0.0; }
}
}

// Swap Ao for An, Bo for Bn
swapBuffers();
}
}

-
Also, I want to mention that it is preferred that you don't cross-post questions, since it appears you have asked very similar questions both here and here. – Godric Seer Jul 22 '12 at 0:16
Have you already seen Greg Turk's work on this, by any chance? – J. M. Jul 22 '12 at 3:44
@J.M.: Not yet. I just tried run his code: it requires X server with PseudoColor, i.e. 8 bit color depth. I think I cannot provide this on OSX. I tried various VNC servers but no luck. – AllCoder Jul 22 '12 at 15:30
I think you should still be able to adapt Turk's approach to the matter at hand; reaction-diffusion patterns seem to be used a fair bit in computer graphics nowadays. – J. M. Jul 22 '12 at 15:32
I might be wrong, but the part with m_An[i][j] = 0.0; might actually add an element to this system that cannot be modeled by a differential equation with a continuous right hand side. This makes it a bit difficult to come up with a faster solver. – Thomas Klimpel Jul 22 '12 at 21:56

You seem to be limited by stability, which is expected since diffusion is stiff as you refine the grid. Good methods for stiff systems are at least partly implicit. It will take some effort, but you can implement a simple multigrid algorithm (or use a library) to solve this system with a cost of less than ten "work units" (essentially the cost of one of your time steps). When you refine the grid, the number of iterations will not increase.

-
If only diffusion would be stiff here, he could use an ADI method like Douglas-Gunn and everything would be fine. However, in my own experience, the reaction part is often much worse with respect to stiffness in addition to being badly nonlinear. – Thomas Klimpel Jul 22 '12 at 14:27
ADI unfortunately has terrible memory locality. Also note that reaction can be treated implicitly regardless of whether diffusion is. Under grid refinement, diffusion will eventually become dominant, but we can't tell where the threshold is without knowing the constants. – Jed Brown Jul 22 '12 at 15:05
Example code implementing backward Euler for this (in Python) is here: scicomp.stackexchange.com/a/2247/123 – David Ketcheson Jul 22 '12 at 20:08
@DavidKetcheson: Using implicit methods requires solving an equation? This is why there is linalg.spsolve() in the code? – AllCoder Jul 22 '12 at 22:00
@AllCoder Yes, it requires a solve, but the solve can be done much faster than all time steps required for an explicit method to be stable. – Jed Brown Jul 22 '12 at 22:05