I have a reaction-diffusion equation in 2 dimensions of the typical form:
$\frac{\partial u}{\partial t} = D\nabla^2u - \Phi(u(x))$
I want to stress that $\Phi(u(x))$, is not a constant, but depends on the location, where it can either be a source or a sink term for the field $u$. The rate of consumption for the sink term here depends on the level of the field $u(x)$. The source of the field is constant.
At the moment, I have been assuming that the field is in steady-state at each time step, so have been solving the nonhomogeneous Poisson equation of the form:
$D\nabla^2u = \Phi(u(x))$
I have been using a finite difference approximation (where alterations are made to the laplacian matrix for the sinks) to solve this on a square domain, with Neumann boundary conditions at each side. The size of my domain is 100 x 100 points.
However, I wanted to know if there could be a faster way to solve for the steady state? I am solving for the field at each time step, which means that even if solver are relatively quick (a third of a second at the moment in c++ and Matlab), the time taken to solve for the field is a bottle neck for speeding up my code.
I have done some reading, and there appear to be some suggestions that Spectral methods could be used here, but I don't want to dive in only to find it takes longer. Similarly, some have suggested solving the top parabolic equation. Any advice here would be most appreciative.
Best,
Ben
Why do you assume that at each time step the solution is in steady state? You can avoid solving for the steady state by using a time step $\Delta t$ and using the approximation $$ \frac{\partial u}{\partial t}\approx\frac{u(t+\Delta t)-u(t)}{\Delta t}=\text{ finite diference approximation at time $t$ of }D\nabla^2u+\Phi(u). $$ To improve the stability of the method you can use an implicit or Crank-Nicholson method on the linear part (i.e. evaluating the linear part at $t+\Delta t$) and explicit for the nonlinear part.