How do I apply BDF2 in a STRANG splitting I have a 3D diffusion equation that I want to solve using a time splitting (2D+1D). Assume that $A$ is the 2D discrete diffusion operator and $B$ is the 1D discrete diffusion operator.
I want to use a STRANG splitting to ensure a global 2nd order of precision :
\begin{equation}
\frac{\partial f}{\partial t}=\frac{1}{2}Af\\
\frac{\partial f}{\partial t}=Bf\\
\frac{\partial f}{\partial t}=\frac{1}{2}Af\\
\end{equation}
And for each time integration, I want to apply the BDF2 implicit formula to ensure stability and a local 2nd order precision:
\begin{equation}
\frac{\partial f}{\partial t}\approx \frac{f^{n+1}-\frac{4}{3}f^n+\frac{1}{3}f^{n-1}}{\frac{2}{3}\Delta t}
\end{equation}
My question is : for the time iteration $n\Delta t$, what is the $f^{n-1}$ of each stage of the splitting ? Is it the state of $f$ at the previous time iteration or the previous splitting stage ?
To be more literal, do I follow this algorithm ?
\begin{equation}
f0,f1\rightarrow(\frac{\Delta t}{2} A) f21\\
f1,f21 \rightarrow(\Delta t B) f22\\
f21,f22\rightarrow(\frac{\Delta t}{2} A) f2
\end{equation}
and for the following iteration we initialize $f0=f22$ and $f1=f2$.
Or this one ?
\begin{equation}
f0,f1\rightarrow(\frac{\Delta t}{2}A) f21\\
f0,f21 \rightarrow(\Delta t B) f22\\
f0,f22\rightarrow(\frac{\Delta t}{2} A) f2\\
\end{equation}
and for the following iteration we initialize $f0=f1$ and $f1=f2$.
Is there another way to ensure an unconditionnaly stable time integration splitting with a global and local 2nd order ?
 A: The idea behind operator splitting is to separate the problem $\partial_t f = (A+B)f$ into $\partial_t = Af$ and $\partial_t = Bf$, which are then solved successively (discrete operators $\mathcal{H}_A^{\Delta t}$ and $\mathcal{H}_B^{\Delta t}$). The Strang splitting scheme
$$f^{n+1} = \mathcal{H}_A^{\Delta t/2} \mathcal{H}_B^{\Delta t} \mathcal{H}_A^{\Delta t/2} f^n $$
has a second-order splitting error. If integration was performed by using the backward Euler method BDF1 in each step, we would write
$$
\begin{aligned}
f^{(1)} &= f^n + \tfrac12\Delta t A f^{(1)} \\
f^{(2)} &= f^{(1)} + \Delta t B f^{(2)} \\
f^{n+1} &= f^{(2)} + \tfrac12\Delta t A f^{n+1}
\end{aligned}
$$
I don't know how linear multistep methods with more steps (such as BDF2) could be implemented in the splitting framework. Nevertheless, one-step methods with sub-steps such as RK schemes can be used.
You may be interested in using the symmetric splitting scheme by Strang
$$
f^{n+1} = \tfrac12 \left( \mathcal{H}_A^{\Delta t}\mathcal{H}_B^{\Delta t} + \mathcal{H}_B^{\Delta t}\mathcal{H}_A^{\Delta t}\right) f^{n}
$$
along with Heun's RK2 method (aka trapezoidal rule). The scheme reads
$$
\begin{aligned}
f^{(1a)} &= f^n + \tfrac12\Delta t A (f^n + f^{(1a)}) &  f^{(1b)} &= f^n + \tfrac12\Delta t B (f^n + f^{(1b)})\\
f^{(2a)} &= f^{(1a)} + \tfrac12\Delta t B (f^{(1a)} + f^{(2a)}) & f^{(2b)} &= f^{(1b)} + \tfrac12\Delta t A (f^{(1b)} + f^{(2b)}) \\
f^{n+1} &= \tfrac12 \big(f^{(2a)} + f^{(2b)}\big)
\end{aligned}
$$
