Optimal use of Strang splitting (for reaction diffusion equation)

Question

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

Answer 1

I'm going to write this as an answer although it doesn't directly answer the question.

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.