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

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?

• Yes, there are many pitfalls. So many, in fact, that I'd say your question is so broad that it's impossible to answer it in a way that would be helpful to you. I would suggest you take a course on numerical PDE solvers, or at least read a book or two on the issue. A short answer such as the ones we can give you here in this forum will not do the issue sufficient justice. – Wolfgang Bangerth Feb 11 '16 at 6:19
• @WolfgangBangerth I am reading Crank's book called "Mathematics of Diffusion" but I am not fully aware of different solvers. There are no well documented and flexible PDE solvers in MATLAB too. I did look at this post and it seems to be a bit helpful. Unfortunately I don't have much time for taking courses at this moment. Even if you can outline some fundamental issues with my problem and provide me some good references, the answer would be helpful. Or else if you can simply suggest the best scheme for solving these systems then nothing like it. – WYSIWYG Feb 11 '16 at 6:39
• @WolfgangBangerth I also found this document which talks about something similar to my approach and calls it "Method of Lines" but does not explain when it can fail. – WYSIWYG Feb 11 '16 at 6:55
• I'd recommend Randy LeVeque's book on finite difference methods. I also agree with Wolfgang; your current question is hard to answer in a canonical way beyond "yes, this is a correct approach" and "yes, there can be issues". Just listing all the possible issues would not make a good answer for this site (look at the help center), so I'd suggest narrowing down your question to your specific situation -- in particular, if you're actually seeing stability issues. – Christian Clason Feb 11 '16 at 8:09
• He does discuss the method of lines on page 184ff; the difficulty with the method of lines is that you need to look at the system of ODEs you actually get after spatial discretization, so it's not quite as easy as "this ODE solver is good for this kind of PDE systems"; you'd need to look at the eigenvalues of the matrix resulting from the spatial discretization. If you have problems with that, you're more than welcome to come back with specific questions! – Christian Clason Feb 11 '16 at 8:36

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.

• I don't need a stiff solver because the rates are more or less similar for all variables. Like you said I can represent the system as you showed above. However, the equation is more like this $u_i''(x,t)=Au_i+f(u)$ i.e. the nonlinear reaction function requires all $u_i$. So it is difficult to use the matrix form. I'm thinking of better ways which I think I'll optimise as I keep working. My current issue is to understand if the approach is good enough. – WYSIWYG Feb 12 '16 at 4:32
• $f$ is applied to all reactants at all points in time? That makes no sense in just about any model. I edited for clarity that $u_i$ is the vector of all reactants at the time point $i \Delta t$. In MATLAB you write $f$ as a vector valued function that takes in all of the reactants at time $i$ and spits out $f$ evaluated at each point in space (you may need to reshape the vector inside of $f$ and use a bunch of .* and ./). Then the code is simply has u(i+1)=A*u(i) + f(u(i)) as the update function with no loop. Let me know if you need more details. – Chris Rackauckas Feb 12 '16 at 4:45
• This vectorized form in MATLAB should be by far the fastest. However, there are some details in the implementation as to why it's not as fast as possible in more efficient langauges I have a blog post about it if you're interested, so if the vectorized version is too slow, then take look at writing your function as a MEX (C) call or in Julia in the devectorized form. Or you can go all the way out to Fortran, but that's a little inconvenient for graphing. – Chris Rackauckas Feb 12 '16 at 14:58
• First try the diffusion equation (no reaction). There is a known solution via Fourier transforms that you can test against. You should check that your order of accuracy is 2 (evaluate by halving/doubling dx a few times and graph it). Then set diffusion to zero and test a reaction equation. Usually I would use some linear reaction equation since there is an actual solution, and you should get the order that matches your method (halving and doubling dt this time). When you have all of that in order, I usually make sure it can solve a Turing patterning problem. – Chris Rackauckas Mar 2 '16 at 15:00
• Something like this or this. If you have the right linear behavior and can reproduce non-trivial nonlinear behavior that's a good sign it's working. – Chris Rackauckas Mar 2 '16 at 15:02

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.

• don't forget the disclaimer. – nicoguaro Feb 12 '16 at 14:37

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..