# Solve steady state reaction-diffusion/Helmholtz equation numerically

I am solving a problem of the form:

$\dfrac{\partial u(x,y,t)}{\partial t} = \nabla^2 u(x,y,t) - f(x,y,t)u(x,y,t) - \kappa(x,y,t)$

At the moment, I am solving this at each time step by assuming a quasi-steady-state:

$\nabla^2 u(x,y,t) = f(x,y,t)u(x,y,t) + \kappa(x,y,t)$

To do this, I use finite differences, which means solving the following (forgetting boundary conditions for now, but for reference they are Neumann) at each time step:

$\dfrac{u^t_{i, j+1} + u^t_{i+1, j} + u^t_{i, j-1} + u^t_{i-1, j} - 4u^t_{i,j}}{h^2} = f^t_{i,j} u^t_{i,j} + \kappa^t_{i,j}$

This means I have to solve a linear system of the form:

$Ax = b$

Where $A$ is a matrix made of the laplacian, as well as the $f^t_{i,j} u^t_{i,j}$ term. For my 100x100 system (resulting in a laplacian of size 10000x10000), this takes about 0.4 seconds for each time step using Eigen's C++ library solver:SimplicialLLT. However, I would like to speed up solving by any means possible. The issue is that I am limited in the pre-conditioning that I can do since the matrix $A$ changes each time step.

Does anyone have any ideas? Finite element, spectral methods, non-steady state approximations all welcome.

Before anyone asks, the terms $f^t_{i,j}$ and $\kappa^t_{i,j}$ cannot be forecast ahead of time, as they depend on $u^{t-1}_{i,j}$. Hence parallelisation is not going to be possible I think.

Best,

Ben

• This seems way too slow. Octave, on my MBP laptop, can solve such a system in $\sim0.0004$ seconds (try n=100; A=rand(n,n); b=rand(n,1); tic; x=A\b; toc). Are you timing just the solve or the whole code? Your matrix formation routine might be the slow part. – Bill Barth Dec 8 '14 at 19:54
• @BillBarth Thanks for your comment. However, the system I solve is much larger than that. It is 100x100 spatial dimensions, which results in a finite difference matrix formulation which is 10000x10000. This takes a time of about a minute on my computer to solve in Matlab. – ben18785 Dec 9 '14 at 0:13
• Sorry, your post said that you were trying to solve a 100x100 system and it was taking 0.4 seconds, you can see how I might have assumed that was the dimensions of your matrix. You should update the post. – Bill Barth Dec 9 '14 at 0:16
• @BillBarth I have updated the post to remove any ambiguity. – ben18785 Dec 9 '14 at 0:19

One approach to acceleration would go as follows: Assume that $A^t=A(u^t)$ is the matrix you try to solve with, i.e., you are looking to solve the linear systems $$A^t x^t = b^t.$$ Let me assume for a moment that $f^t \ge 0$, then $A^t$ is a symmetric and positive definite matrix. (If my assumption should be wrong, then it is still symmetric but may no longer be positive definite.) Given that, you can solve the linear system with the Conjugate Gradients (CG) method. The question of course is how to precondition it.
My approach would be to compute a Cholesky decomposition of $A^0$ in the first time step and use this as the preconditioner for the next time steps. You will be able to solve time step zero in exactly one CG iteration since the preconditioner is perfect. In subsequent time steps, $A^t$ will start to deviate from the preconditioner and you will need more iterations. Once the number of iterations becomes too large, you re-compute the Cholesky decomposition in the next time step. This works because the expensive step is to compute the Cholesky decomposition, but you only do that every few time steps.