1

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?

WYSIWYG
  • 143
  • 1
  • 11
  • 1
    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 06: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 06: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 06:55
  • 3
    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 08:09
  • Thanks @ChristianClason and Wolfgang. I did not realize that this question would be broad. You can VTC it. I'll read the book. Right now my simulations "seem" to work. I'll simulate a known example from the book and compare with the analytical solution to see if my code works fine. Does Randy LeVeque explain which methods are good for what kind of systems? – WYSIWYG Feb 11 '16 at 08:25
  • 1
    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 08:36
  • I think I've corrected a typo on the $u''$ discretisation, if you want to check it – Steve Feb 11 '16 at 10:08
  • Since you are using MATLAB, is there any particular reason why you aren't using the standard pdepe function to solve your PDE system? It generally works well for the type of system you show and using it is certainly more straightforward than writing your own code. – Bill Greene Feb 11 '16 at 13:08
  • @BillGreene That's a valid question: 1. I don't have the pde-toolbox 2. I don't understand the pdepe documentation well enough (looks quite messy compared to ode solvers). 3. I also need to solve for steady state. I'll try pdepe too. Thanks for the suggestion. – WYSIWYG Feb 11 '16 at 13:32
  • @BillGreene I am running pdepe now. It is too damn slow. I just realized that even pdepe runs an ode solver- ode15s – WYSIWYG Feb 11 '16 at 15:01
  • pdepe is an implementation of the method of lines approach so an ODE solver is required. ode15s is usually the MATLAB ODE solver of choice for these types of applications. When you say "too slow", how long is that and what kind of time were you expecting? The number of mesh points in the spatial direction affects the computation time. Perhaps you have more than you need for your required accuracy. Similarly, you can reduce the error tolerances used by ode15s. Finally, you can use the MATLAB profiler to see where the time is spent. You would face the same issues writing your own code. – Bill Greene Feb 11 '16 at 17:21
  • @BillGreene Well.. my code runs in 9 seconds. Whereas pdepe with the same time span and lower number of spatial grid points was taking more than 10 min. Moreover, if pdepe uses ode15s then why does it take the time points. That is anyway a different discussion. Perhaps if you can find time we can discuss this in chat. When I use ode15s (with Abs/RelTol=1e-7) on my FD-discretized system, it takes almost same time as ode45 and gives the same result. I also use a parfor loop in my code (not multicore) so that might make it a bit faster but still. I am reading the book now to know more. – WYSIWYG Feb 11 '16 at 17:45
  • 1
    Solving this shouldn't require looping. You should change the diffusion operator to a sparse tridiagonal matrix. Then by using a matrix equation MATLAB will take care of multi-threading the operation for you, and it will run in C and be much quicker. If you have everything as vectors, simply use .* and ./ in your reaction equation $f$ and it will automatically do it at all points in space (and automatically multi-thread it). – Chris Rackauckas Feb 12 '16 at 02:42
  • @ChrisRackauckas I understand that it is possible without looping and I did try that but the reaction terms (26×1 vector) are such that I have to solve them for each spatial point. I have to either loop for the variables or the points. However my reaction function (similar to an ode function) takes a vector. So the only option I could see is to loop through the gridpoints. There may be better ways, I'm sure. I'll implement this in FORTRAN if necessary. – WYSIWYG Feb 12 '16 at 04:21
  • You should just make one big vector with all the points in reactants and in space, and adjust the construction of the diffusion matrix for this. This will get rid of your loop and make a massive speed up. Or, I think the best solution, port your code to Julia. It doesn't take much time at all to port MATLAB code to Julia, but if you have loops you can get a massive speedup. – Chris Rackauckas Feb 12 '16 at 04:38
  • @ChrisRackauckas I did exactly that way. Yes I can adjust the diffusion matrix accordingly as you said. I was thinking about it. – WYSIWYG Feb 12 '16 at 04:39

3 Answers3

3

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.

Chris Rackauckas
  • 12,320
  • 1
  • 42
  • 68
  • 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 04: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 04:45
  • No I didn't mean time. We can perhaps discuss this better in chat if you are free for that. I just meant that *f* takes all reactants at a given space and time point. If *u* is a long vector (m×n,1) of all reactants(n) at all spatial points(m) either listed reactants first or timepoints first. The *f* should be redefined such that it operates on this huge vector; basically reapply the same reaction rules for all the m-spatial points. This will require a loop. – WYSIWYG Feb 12 '16 at 10:13
  • Yes I can reshape the vector to a matrix and apply *f* on all the columns and again reshape back to a vector. This, I guess, would take more time than running a loop. – WYSIWYG Feb 12 '16 at 10:18
  • I would be highly surprised if that's slower than a loop. Another thing you can do is in your vectorized code you can call "sub-vectors of the reactants". For example, if your vector is all the reactants at space j, then all the reactants at j+1, etc., then you can get the vector of all of the reactant 1 via u(1:26:end). So you can add at all points in space for two reactants via u(1:26:end)+u(2:26:end). If you grouped by space, you would just make the middle number the number of space points. This is pretty much what your loop is likely doing. – Chris Rackauckas Feb 12 '16 at 14:51
  • However, given the way de-referencing works, if you have lots of reactions, it could be more efficient to just save the matrix of u(1:26:end),u(2:26,end), etc., and the act of building this matrix in its most efficient form is reshape. – Chris Rackauckas Feb 12 '16 at 14:58
  • 1
    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
  • I have done some modifications in the code such that there are no loops and the computations are done on the vectorized/matrix forms. I also split the diffusion from the reaction terms as suggested in Randy LeVeque's book. The speed has tremendously increased. However, I need to know if this program is actually solving the problem correctly. Can you suggest me some simple published solutions that I can replicate? Do you think it is a good idea to share my program somewhere for a review ([codereview.se]??) – WYSIWYG Mar 02 '16 at 07:30
  • 1
    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 02 '16 at 15:00
  • 1
    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 02 '16 at 15:02
  • Thanks. I was trying it on a Turing problem. I'll do the basic problems too. Diffusion seems to work as I can see the concentration becoming uniform. I have to check the dynamics though. Reaction works like well mixed system (ODE) when diffusion is set to zero. I didn't get your point on halving and doubling dt. I do not set dt. I use ode15s to solve the spatially discretized problem. Do you mean tolerance? I can change dx and see if the solution is similar. – WYSIWYG Mar 02 '16 at 18:50
  • Oh yeah, if you're using an adpative solver don't worry about changing dt to check the order of accuracy. You can just assume your temporal integrator (ode15s) is correct. I forgot about that and so for example if you coded your own RK4, you would want to check you get 4th order accuracy. – Chris Rackauckas Mar 03 '16 at 06:00
1

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.

David Ketcheson
  • 16,522
  • 4
  • 54
  • 105
0

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

FancyPants
  • 133
  • 5