Numerical simulation of a reaction-diffusion system on MATLAB with finite difference discretization of spatial derivative

Question

I have a model of a system which consists of diffusing reactants and intermediates. For each variable, $u_i$, the final representative equation looks like the standard reaction-diffusion form:

$$\frac{\partial}{\partial t}u_i(x,t)=f_i(\bar{u}(x,t))-D\frac{\partial^2}{\partial x^2}u_i(x,t)$$

where $f_i$ is a non-linear function that denotes the reaction term. I am using the finite difference method (FDM) to discretize the spatial derivative with Neumann boundary conditions: $$\frac{\partial u_i}{\partial x}{\bigg|}_{x=0, x=L} = 0$$

using the central difference formula:

$$u_i''(x)\approx\frac{u_i(x-h) - 2u_i(x) + u_i(x+h)}{h^2}.$$ After reducing the PDE to an ODE, I am using MATLAB's ode15s (RK4 based solver) to numerically integrate the equations in time.

I have not worked on PDEs much and online resources on FDM (such as even wikipedia), depict simultaneous discretization of the time derivative, also.

As far as I understand, the discretization shown in these cases, depict the simple Euler integration (either explicit or implicit). Even the Crank-Nicolson scheme is shown like that.

Is my approach correct (discretize spatial derivative and solve using an optimized ode solver)? Can there be numerical stability issues?

Answer 1

Systems biologist studying reaction-diffusion systems here. The method that you started with is probably the most common. In fact, I would say that it's the first thing to try for this type of problem. If you just let $A$ be the tridiagonal matrix (1,-2,1), then you can write the system as:

$ u_{i+1} = Au_{i} + f(u_{i})$

where $u_{i}$ is the vector of all reactants and $f$ is in its vectorized form. This could then be thrown into any ODE solver. This is the Method of Lines (MOL) approach that FancyPants noted in a succinct form. For many reaction-diffusion problems this is sufficient. If your instability comes from stiff reaction equations (i.e. the non-spatial model is stiff), then using a stiff ODE solver in this form is usually sufficient. Not only that, but because of the highly vectorized form of this equation, you can easily make use of GPGPU/Xeon Phi/multi-node computing to "brute force" the solution (in MATLAB, this takes little/to no extra effort).

However, if you have high diffusion constants (or your diffusion is a highly variable equation), you may not be able to get the system stable enough for the MOL method to work. If that's the case, there are two ways that you can go. One common way is the Crank-Nicholson method. A quick intuition is that it is simply a midpoint method in space and a midpoint method in time. Just about any computational PDEs text will outline the method.

However, this will result in a large implicit system which you will have to solve using a nonlinear solver like Newton's method. The real problem is that your implicit equations are all coupled, meaning they cannot be solved in parallel which can hurt for large systems. Thus one of the state-of-the-art methods in this field are the Implicit Integration Factor Methods. For example, if you look at Equation 26 in the linked paper, if you were to compute your function at previous values as some large constant $C$, then you see that the implicit equation decouples to that it's independent at each point in space. This means you can run a nonlinear solver at each point in space (meaning a much smaller Jacobian for faster solutions (and you can easily for solve it pen/paper instead of numerically approximating it), and you can solve at each point in space in parallel/GPU/etc.) and get massive practical speedups over the Crank Nicholson Method.

In summary: first try the method of lines with ODE45, if that doesn't work try stiff solvers. Else go to Crank-Nicholson (if you have solve of that code lying around), but if you will be solving lots of these types of problems, try an Implicit Integration Factor Method.

Note: If you're looking at periodic problems, you may also want to try Spectral methods. Also, I am assuming you're using a 1D system. If not, look at operator splitting methods (either ADI or for integration factor methods), or finite element methods.

Answer 2

What you're doing is called the method of lines and it is the most common way of discretizing an evolution PDE. An example of how to do it stably for reaction-diffusion problems (written by me) is available here in Python.

Answer 3

As you said yo can discretize only the space derivative, thus obtainng a set of DAEs (differential-algebraic equations), or discretize both space and time in a way to obtain a set of AEs (algebraic-equations).

As mentioned in a comment, the MOL (method of lines) is a discretization approach which assumes to discretize only the spatial domain and leave the time as a countinuous one. The main advantage is that the time step is chosen automatically by the solver according to the enforcement of some error checking rule.

If you would implement a set of AEs, thus discretizing both space and time domains, then is up to you to choose the time step an/or implement an algorithm able to accept/reject a time step size according to some error-checking rule.

Moreover in case that you would implement this latter formulation, in some cases the numerical scheme can be unstable (in literature there are many examples that show this particular behavior while solving the same set of equations with different numerical methods with fixed stepsize). Finally, when working with nonlinear PDEs (like the one you are using), for a AE implementation you need also to linearize the nonlinear terms around a working point in order to implement some kind of Newton-Raphson method to solve iteratively the nonlinearities. In any case there are many ways to do this : an efficient algorithm was developed by Professor Newman called BAND, which is also available in Fortran routines. This aspect would be automatically managed by the solvers of DAEs such as ode15i, ode23, SUNDIALS suite, DASSL etc..