# Optimal use of Strang splitting (for reaction diffusion equation)

I made a strange observation while computing the solution to a simple 1D reaction diffusion equation:

$\frac{\partial}{\partial t}a=\frac{\partial^2}{\partial x^2}a-ab$

$\frac{\partial}{\partial t}b=-ab$

$\frac{\partial}{\partial t}c = a$

The initial value of $b$ is a constant ($b(0,x)=b_0$), and I'm only interested in the integral over $a$ from $0$ to $1$ ($\int_0^1a(t,x)dt$). The purpose of $c$ and the equation $\frac{\partial}{\partial t}c = a$ is just to evaluate this integral.

I used a Strang splitting scheme for the coupling between diffusion and reaction (a half step reaction, then a full step diffusion, and then again a half step reaction), a Crank Nicholson scheme for the diffusion, and an analytical solution for the reaction (including the equation $\frac{\partial}{\partial t}c = a$).

Because one step of the analytical solution was more than a factor 3 slower than one step of the Crank Nicholson scheme, I tried to make more than one Crank Nicholson step for each reaction step. I was hoping to get by with fewer step of the Strang splitting scheme, so that I would be faster overall.

However, the opposite effect can be observed, namely that much more steps for the Strang splitting scheme are required if more than one Crank Nicholson step is used. (I'm only concerned with the accuracy of the integral over $a$, which seems to converge faster than $a$ itself.) After wondering for some time, I noticed that the same effect also happens for $b(t,x)=b_0=0$, and that I even understand why for this case. The point is that if I make exactly one Crank Nicholson step, then the overall scheme turns into a trapezoid rule (if $b(t)=0$).

So if I would treat $\frac{\partial}{\partial t}c = a$ as part of the diffusion step, increasing the number of Crank Nicholson steps (probably) wouldn't lead to reduced overall accuracy (as observed). But that seems to defeat the purpose of using an analytical solution for the (non-linear and potentially very stiff) reaction part of the system.

So here is my question: Is there a better way to treat $\frac{\partial}{\partial t}c = a$ in the context of a Strang splitting, than to either treat it as part of the reaction step, or to treat it as part of the diffusion step. I want to avoid being "forced" to use exactly one Crank Nicholson step for the diffusion. (For example in 3D, I would prefer to solve the diffusion analytically by an FFT instead of using Crank Nicholson. Of course I can also combine FFT with Crank Nicholson, so it's not such a big deal.)

• In people.maths.ox.ac.uk/dellar/OperatorLB.html, a similar effect seems to be described. The conclusion is that it's crucial to use Crank Nicholson instead of the exact solution. So perhaps the answer to my question is a simple no. – Thomas Klimpel Apr 5 '12 at 15:21
• Something seems wrong with your equations. $c$ doesn't appear in the first two making the coupling one-way and meaning that you could compute $c$ at any $t$ as a post-processing step. – Bill Barth Apr 5 '12 at 21:13
• @BillBarth I changed the question to clarify the role of $c$. So $c$ is just a means to compute $\int_0^1a(t,x)dt$. Please let me know if you have any suggestion how to compute this integral more accurate (than what I get from the combination of Strang splitting and Crank Nicholson described above), potentially using a post-processing step. – Thomas Klimpel Apr 5 '12 at 23:36
• This is now a long time gone, but did you recognize that this system of equations can be written as a parabolic PDE in $c$ with an exponential reaction term? I guess I wonder if you really want to solve this 3-variable system instead of a simplified one. – Bill Barth Jul 22 '13 at 17:22
• @BillBarth I would be interested to learn how this system can be written as a parabolic PDE with an exponential reaction term. The speed of the solution of this model is a limiting factor during model calibration (which can take several hours), even so the used accuracy with respect to time integration is quite far from full convergence. – Thomas Klimpel Jul 22 '13 at 18:06

Plugging the second equation and the third equation into the first, and plugging the third into the second, together give: \begin{align} \frac{\partial^2 c}{\partial t^2} &= \frac{\partial^2}{\partial x^2}\frac{\partial c}{\partial t} + \frac{\partial b}{\partial t} \\ \frac{\partial b}{\partial t} &= - \left(\frac{\partial c}{\partial t}\right) b \end{align} Rearranging these two gives: \begin{align} \frac{\partial}{\partial t}\left(\frac{\partial c}{\partial t} - \frac{\partial^2 c}{\partial x^2} - b \right) &= 0 \\ \frac{1}{b}\left(\frac{\partial b}{\partial t}\right) &= -\frac{\partial c}{\partial t} \end{align} Now, we can integrate both of these once in $t$, leaving for the first equation: $$\frac{\partial c}{\partial t} - \frac{\partial^2 c}{\partial x^2} = b + A(x)$$ Using the third equation, we can express the "constant" of integration as $A(x)=a_0-\frac{\partial^2 c_0}{\partial x^2}-b_0$. The second equation is a bit more tricky. Rewriting a little, we have: $$\int_0^t \frac{1}{b(x,t')}\left(\frac{\partial b(x,t')}{\partial t'}\right) dt' = -\int_0^t \frac{\partial c(x,t')}{\partial t'} dt'$$ This leads to the solution $$\ln b(x,t) - \ln b_0(x) = - c(x) + c_0(x)$$ Or $$\ln \frac{b}{b_0} = -c + c_0$$ Exponentiating gives: $$b = b_0 e^{c_0-c}$$ And, finally, plugging this into the PDE for $c$ gives $$\frac{\partial c}{\partial t} - \frac{\partial^2 c}{\partial x^2} = b_0 e^{c_0-c}+A(x)$$
Replacing $c$ by $c-c_0$, or equivalently using the initial conditions $c_0(x)=0$, this equation simplifies to $$\frac{\partial c}{\partial t} = \frac{\partial^2 c}{\partial x^2} + a_0 - (1-e^{-c}) b_0$$ Now, you should be able to find considerable literature on the best way to solve this equation. I don't know if Crank-Nicholson is a good choice for the exponential term, but it seems plausible. Some care should probably be taken to guarantee that $c > c_0$ everywhere or the solution may blow up quickly.
I've only walked through this derivation twice, so there may be an error or two in it, but it has the right feel to me. If $b_0 = 0$ everywhere, then this is clearly the right solution, and it has a whiff of plausibility otherwise.
Solving this until $t=1$ and evaluating $c(x,1)$ should give you the answer you are looking for.
• Thanks so much for this answer. I find it quite illuminating, at least it makes it easier for me to understand/predict the behavior of the solution. Another advantage is that the time evolution of $c$ is slower than the time evolution of $a$, so I'm quite optimistic that convergence will be better than before. – Thomas Klimpel Jul 23 '13 at 19:07