# Mathematical foundations of the fractional step method

The well-known method for solving multi-physics time-evolution problems called "the method of fractional steps" or "the method of splitting" usually applies to PDEs which have a sum of terms on the right-hand side, and basically it consists of using partial updates, one term at a time. So for a PDE

$$\partial_{t} f = F(f,t) + G(f,t)$$

the fractional-step approach would result in a time-stepping scheme like this

$$\tilde{f} = f(t) + F(f(t),t) dt$$

$$f(t+dt) = \tilde{f} + G(\tilde{f},t) dt,$$

and there are some variations to improve the accuracy (Strang splitting etc.)

Of course there is a great advantage in using dedicated methods (codes) tuned up for specific parts of physics in the full problem, e.g., for dealing with non-ideal hydrodynamics one could combine a method specific for ideal hydrodynamics and a diffusion equation solver. But would this be always guaranteed to work?

If we have two (or more) complicated models $$F$$ and $$G$$ available as “black boxes” (possibly implemented as separate complicated codes), and we want to study self-consistent time evolution with both $$F$$ and $$G$$ in the right-hand side, how can we convince ourselves that the solution is correct? Some confidence can be perhaps based on convergence studies varying the time-step $$dt$$. But even if the results do converge in the time step $$dt$$, does this really guarantee the solution is correct? I guess this question more broadly applies to any time-integration method for nonlinear equations. There is the Lax equivalence theorem that guarantees convergence to the exact solution for a consistent finite-difference approximation, if the time-integration method is stable, but it is formulated only for linear time-evolution problems.

For what classes of problems can it be proven that the fractional-step time-integration techniques solution converges to the exact solution of the equation? For what classes of problems can it be demonstrated that this techniques would not lead to correct solution?

## 1 Answer

The short answer to your question is that it works because of the Zassenhaus formula. Briefly, if $$A$$ and $$B$$ are linear operators, then

$$e^{t(A + B)} = e^{tA}e^{tB}e^{-\frac{t^2}{2}[A, B]}\cdot\ldots$$

where the ellipses denote terms of higher order in $$t$$. If both parts of your dynamics are linear, then for small enough timesteps you can split the integration and the splitting error will be of order $$\delta t^2$$ by this formula. Those "variations to improve the accuracy" as you put it can be derived by symmetrizing the scheme and then looking at the higher-order terms in the Zassenhaus formula. For nonlinear problems, you can linearize around the current value and assume the timestep is sufficiently small.

This approach can work well for many problems but there are lots of caveats. Suppose you were solving an advection-diffusion equation by splitting the two parts. If you solve the unsplit form with an implicit method, the timestepping scheme and the diffusive dynamics will provide all the stability you need to integrate the advective dynamics right. Now say you solved the split form using an implicit method for the diffusion and an explicit method for the advection, but you weren't very careful about adding artifical diffusion or other stabilization to the advective solver. If the timestep is too long, you could end up introducing instabilities via the advection solver that the diffusion solver can't adequately suppress.

That's more of an easy mistake to make, but there are problems that are just fundamentally more difficult to solve via splitting. A good example is problems involving combustion, where the timescale for chemical reactions is orders of magnitude shorter than the timescales for fluid transport. It would be really convenient if you could use a stiff ODE solver for the chemistry and a shock-capturing conservation law solver for the fluid dynamics. But even though you might have a stable integrator for the chemistry, the splitting error can become gigantic. This paper has a nice analysis of why that's such a hard problem.

If you want a reference, there's a wonderful description of this in Hairer, Lubich, and Wanner's book Geometric Numerical Integration in chapter 3. They present things the fancier way using Lie derivatives, which is more rigorous than my slapdash argument using linearization.