finite difference : why should we solve linear equation at each step

Question

I'm trying to simulate the growth of a brain tumor using a 3d reaction-diffusion model $ \partial_{t}u = \nabla.{(D\nabla{u}) } + ku.(1-u)$ , knowing the initial distribution of tumor $u^0$, the non-homogeneous diffusion coefficient $D$ and the non-homogeneous net-proliferation map $k$.

So I first tried to proceed to simulation simply by dividing it into small time steps and directly simulate the process by simply replacing the derivatives what is appearing to me now as forward Euler finite difference. I got serious instability problems which, if I understand well, are due to the difference scheme I'm using.

So I'm trying to learn how to correctly code finite difference algorithms. At this point there's something I don't understand and I hope someone will be able to explain me where I'm wrong.

Let forget the proliferation term and do simply a 1D diffusion ( $ \partial_{t}u = \partial_{x}(D\partial_{x}u) $ ) using a Crank-Nicolson scheme for the time difference and a central difference for the spatial derivatives.

So, if $D$ depends on $x$, we have $ u^{n+1}_{i} = u^{n}_{i} + \frac{\Delta{t}}{2} \{ \partial_{x}(D\partial_{x}u^{n}) + \partial_{x}(D\partial_{x}u^{n+1}) \} $

which develop to

$ \Rightarrow ( 1 - \frac{\Delta t}{2}\partial_{x}(D\partial_{x}\cdot) ) u^{n+1}_{i} = ( 1 + \frac{\Delta t}{2}\partial_{x}(D\partial_{x}\cdot) ) u^{n}_{i}$

and using $ \partial_{x}(D\partial_{x}V) = \partial_{x}D.\partial_{x}V + D.\partial_{xx}V $,

$\partial_{x}V = \frac{V_{i+1}-V_{i-1}}{2\Delta x}$

and $\partial_{xx}V = \frac{V_{i+1}+V_{i-1}-2V_{i}}{{\Delta x}^2}$

one get

$ \partial_{x}(D\partial_{x}u_{i}) = \frac{(D_{i+1}-D_{i-1}).(u_{i+1}-u_{i-1})}{4{\Delta x}^2} + \frac{D_{i}.(u_{i+1}+u_{i-1}-2u_{i})}{{\Delta x}^2}$

ordering it according to the spatial indice of $u$ , one get

$ \partial_{x}(D\partial_{x}u_{i}) = \frac{1}{4{\Delta x}^2}.\{ (D_{i+1}-D_{i-1}+4D_{i})u_{i+1} + 8D_{i}u_{i} - (D_{i+1}-D_{i-1}-4D_{i})u_{i-1} \} $

In matrix form, this can be written

$ \partial_{x}(D\partial_{x}u) = \tilde{D}u $

with $ \tilde{D} \equiv \frac{1}{4{\Delta x}^2}\begin{bmatrix} -8D_1 & (D_2+4D_1) & 0 & \cdots \\ -(D_3 - D_1 -4D_2) & -8D_2 & (D_3-D_1+4D_2) & \cdots \\ \vdots & \vdots & \vdots & \vdots \\ 0 & \cdots & -(- D_N -4D_{N-1}) & -8D_N \end{bmatrix} $

Therefore, one have to solve, at each time step, the equation $ (1 - \frac{\Delta t}{2}\tilde{D})u^{n+1} = (1 + \frac{\Delta t}{2}\tilde{D})u^{n} $

my problem is that on every course/paper I read on this subject, the authors are explaining that one have to solve the equation at each time step which is the time consuming step. However, as it is written here at least, the matrix $\tilde{D}$ is completely determined from the (fixed in time) $D$ distribution and it looks to me that one can just simply generate it, then calculate the total matrix $ M = (1 - \frac{\Delta t}{2}\tilde{D})^{-1} \cdot (1 + \frac{\Delta t}{2}\tilde{D})$ and then simply calculate $u^n = M^nu^0 $ which is very fast.

I tried to code it in MATLAB, and indeed, doing it step by step, or just at once gives me the same results ( and which are coherent with a simple Gaussian blurring with $\sigma = \sqrt{2D.T} $ in the case of an homogeneous D, just to check the validity of the simulation).

So i must have missed something but what ? Maybe finite difference is not the tool to use in case of diffusion-only and are actually used because of the reaction term, but still, if I understand well, one can separate the simulation into smaller sub-time steps to perform only diffusion or only proliferation (solved analytically for the later) with Douglas-Gunn or equivalent approaches. In such a case, the diffusion step can still be solved using a simple matrix multiplication at each step which do not require a long processing ? So I would be very happy if someone can tell me if I'm going on the right track or if (as I suspect) I missed any fundamental issue...

Many thanks !

D. Guez

Answer 1

You're absolutely right that the matrix doesn't change during each time step, so there's no need to invert it more than once. In practice, however, you probably shouldn't be inverting the matrix at all!

The fact that the coefficient matrix is constant is exactly why matrix factorization (e.g. LU) can provide significant speedup. Factorize once at the beginning of the simulation, and you just have to solve triagonal matrix problems, which is $O(N^2)$ at each time step, rather than the $O(N^3)$ for a direct solution (though computing the the factorization is $O(N^3)$, you only do it once). This answer over at Math explains it better the I do, I think.

Also worth noting that tridiagonal systems like your simplified example are especially easy, you can solve them in $O(N)$ using the Thomas Algorithm.

If you really only care about the result at time step $n$, and don't care at all about the intervening values, your $u^n = M^nu^0 $ approach is valid, though I'm not sure how much that will help. In fact, it will probably slow you down since matrix multiplication is between $O(N^{2.3})$ and $O(N^3)$ (i.e. not "very fast", and this suggests Matlab uses an $O(N^3)$ implementation), and the initial matrix inverse itself is costly.

There's nothing wrong with finite difference for diffusion only problems, though gotta be careful with Euler forward to keep your $\Delta t$ small enough for stability issues ($\Delta t < \frac{ \Delta x^2}{2D}$, iirc). Even Crank-Nicholson can be oscillatory (though stable) if too long a step is taken. Also, since your proliferation term is non-linear, if you do an implicit solution, you're going to have to solve a non-linear system of equation at each time step.