Operator Splitting, Strang Splitting, and IMEX Schemes for PDE Solvers

Time-dependent partial differential equations often combine several physical processes. A transport model may include advection, diffusion, chemical reactions, external sources, and nonlinear feedback. Each part can have different mathematical properties and may require a different numerical treatment.

Advection is commonly handled with methods designed for wave-like transport. Diffusion often creates stiffness and benefits from implicit integration. Reaction terms may be inexpensive and non-stiff, or they may contain very fast chemical processes that require a specialized ordinary differential equation solver.

A single monolithic method can solve all terms together, but the resulting system may be large and difficult to implement. Operator splitting provides another option. It separates the full equation into smaller subproblems, solves them independently, and combines their results over each time step.

Lie-Trotter splitting offers a simple first-order method. Strang splitting improves temporal accuracy through a symmetric sequence of substeps. IMEX schemes pursue a related goal by treating selected terms implicitly and others explicitly within one additive time-integration method.

What Is Operator Splitting?

After spatial discretization, a time-dependent PDE often becomes a large system of ordinary differential equations:

du/dt = A(u) + B(u)

The operator A may represent advection, while B represents diffusion or reaction. More complicated systems may contain three or more operators.

Operator splitting replaces the combined problem with a sequence of simpler subproblems. Instead of integrating A + B simultaneously, the method advances the solution under A and then under B.

For a linear autonomous problem, the exact evolution over a time step h can be written formally as:

u(t + h) = exp(h(A + B))u(t)

If the operators commute, meaning:

[A, B] = AB - BA = 0

then the exponential separates exactly:

exp(h(A + B)) = exp(hA) exp(hB)

In that special case, sequential integration introduces no splitting error. In most practical PDEs, however, the operators do not commute. Their order then matters, and the separated solution only approximates the combined evolution.

Why Split a PDE Into Separate Operators?

The main benefit is modularity. Each physical process can use the numerical method best suited to it.

An advection operator may use an explicit finite volume method with a flux limiter. A diffusion operator may use an implicit linear solver. A reaction operator may use a local stiff ODE integrator. These components can be developed, tested, and improved separately.

Splitting can also reduce memory requirements. A monolithic implicit method may require one large matrix containing every coupled term. A split method can solve smaller systems or reuse existing operator-specific solvers.

The approach is especially attractive in multiphysics applications where mature solvers already exist for each process. Instead of rewriting them as one system, developers can connect them through a controlled time-stepping sequence.

Lie-Trotter Splitting

The simplest sequential method is commonly called Lie-Trotter splitting. For one time step of length h, it applies one operator followed by the other:

u*      = SolveA(uⁿ, h)
uⁿ⁺¹    = SolveB(u*, h)

The order can also be reversed:

u*      = SolveB(uⁿ, h)
uⁿ⁺¹    = SolveA(u*, h)

When the operators do not commute, the two sequences generally produce different results. Both are first-order accurate in time under standard assumptions. The local splitting error is typically proportional to h², while the accumulated global error over a fixed interval is proportional to h.

Lie-Trotter splitting is easy to implement and requires only one solve for each operator per time step. It is useful for prototypes, low-accuracy calculations, and applications where the time step is already restricted by another stability or resolution condition.

Its main weakness is that achieving a small temporal error may require many short steps.

Strang Splitting

Strang splitting uses a symmetric half-step, full-step, half-step sequence:

u*      = SolveA(uⁿ, h / 2)
u**     = SolveB(u*, h)
uⁿ⁺¹    = SolveA(u**, h / 2)

The alternative arrangement places B on the outside:

u*      = SolveB(uⁿ, h / 2)
u**     = SolveA(u*, h)
uⁿ⁺¹    = SolveB(u**, h / 2)

The symmetric composition gives second-order global accuracy when the operators and their solutions are sufficiently regular. Its local error is generally proportional to h³.

This improvement makes Strang splitting a common default for production simulations. It provides substantially better temporal accuracy than first-order sequential splitting without requiring a fully coupled solver.

The method is not automatically accurate for every time step. Each subsolver must also resolve its own process adequately. A formally second-order splitting sequence cannot compensate for an inaccurate first-order subsolver or a time step that fails to capture fast dynamics.

Why Symmetry Improves Accuracy

The error can be studied with the Baker-Campbell-Hausdorff expansion. For two linear operators, a simple sequential product has the form:

exp(hA) exp(hB)
= exp(h(A + B) + h²[A, B] / 2 + higher-order terms)

The commutator term shows why applying the operators independently does not normally reproduce the exact combined solution.

The Strang composition is:

exp(hA / 2) exp(hB) exp(hA / 2)

Because this sequence is symmetric in time, the leading first-order global splitting error cancels. The remaining leading terms involve nested commutators such as:

[A, [A, B]]
[B, [B, A]]

The exact coefficients depend on the chosen arrangement, but the practical conclusion is clear: splitting accuracy depends not only on the time-step size but also on how strongly the operators fail to commute.

Understanding Splitting Error

Splitting error is separate from spatial discretization error and from the error introduced by each time integrator. A simulation may therefore contain several error sources at once.

The splitting contribution tends to be small when the operators interact weakly or vary smoothly. It can become larger when coefficients change sharply, nonlinear feedback is strong, or one process immediately changes the coefficients used by another.

Consider a reaction-diffusion problem in which reaction rates depend strongly on local temperature. If the reaction step changes the temperature or concentration rapidly, performing diffusion before reaction may produce a noticeably different intermediate state than performing reaction first.

Reducing the time step usually reduces this disagreement. Comparing both operator orders can also provide a simple indication that splitting effects are significant, although it is not a complete error estimate.

Operator Ordering Matters

For non-commuting operators, there is no universally best order. The choice should reflect the physics, the relative time scales, and the required output.

In Strang splitting, the operator placed on the outside is evaluated twice per complete step, although the final half-step of one step can sometimes be combined with the first half-step of the next. The more expensive operator may therefore be placed in the middle to reduce repeated setup work.

The outer operator also acts last, which can influence which constraints are satisfied most accurately at the end of a time step. For example, a final reaction half-step may preserve a local chemical equilibrium differently from a final transport half-step.

Developers should test plausible orderings against a reference solution or a much smaller time step rather than assuming that one arrangement is always superior.

What Are IMEX Schemes?

IMEX means implicit-explicit. An IMEX method divides the right-hand side into a non-stiff part and a stiff part:

du/dt = F(u) + G(u)

The term F is evaluated explicitly, while G is treated implicitly. This avoids solving the entire nonlinear system implicitly while retaining better stability for the stiff contribution.

IMEX schemes are often constructed as additive Runge-Kutta or multistep methods. Unlike fractional-step operator splitting, the explicit and implicit terms participate in a shared set of intermediate stages.

This coupling can reduce some errors caused by solving complete physical processes one after another. However, the implementation usually requires a compatible IMEX framework and implicit solves at individual stages.

Operator Splitting vs. IMEX

Operator splitting is often the natural choice when a codebase already contains independent transport, diffusion, and reaction solvers. IMEX is attractive when the software already supports additive time integrators or when simultaneous stage coupling produces better accuracy.

A Reaction-Diffusion Example

A common model combines diffusion with a nonlinear reaction:

∂u/∂t = D∇²u + k u(1 - u)

The diffusion operator is:

A(u) = D∇²u

The reaction operator is:

B(u) = k u(1 - u)

A Strang step can advance reaction for half a time step, diffusion for a full time step, and reaction again for half a step.

The following FiPy example demonstrates this pattern:

from fipy import Grid1D, CellVariable, TransientTerm, DiffusionTerm

# Spatial mesh
nx = 100
length = 1.0
dx = length / nx
mesh = Grid1D(nx=nx, dx=dx)

# Solution variable
phi = CellVariable(
    name="phi",
    mesh=mesh,
    value=0.0
)

# Initial condition
x = mesh.cellCenters[0]
phi.setValue(
    1.0,
    where=(x > 0.4) & (x < 0.6)
)

# Model coefficients
diffusion_coefficient = 1.0
reaction_rate = 5.0

# Implicit diffusion equation
diffusion_equation = (
    TransientTerm(var=phi)
    == DiffusionTerm(
        coeff=diffusion_coefficient,
        var=phi
    )
)

dt = 0.001
number_of_steps = 1000

for step in range(number_of_steps):
    # First reaction half-step
    reaction = (
        reaction_rate
        * phi.value
        * (1.0 - phi.value)
    )
    phi.setValue(
        phi.value + 0.5 * dt * reaction
    )

    # Full implicit diffusion step
    diffusion_equation.solve(
        var=phi,
        dt=dt
    )

    # Second reaction half-step
    reaction = (
        reaction_rate
        * phi.value
        * (1.0 - phi.value)
    )
    phi.setValue(
        phi.value + 0.5 * dt * reaction
    )

The example uses an explicit Euler update for each reaction half-step and an implicit FiPy solve for diffusion. It illustrates the splitting sequence, but the explicit reaction update still has its own stability and accuracy restrictions.

If the reaction is strongly stiff, a local implicit method, an exact reaction solution, or a dedicated stiff ODE solver may be required. Strang splitting determines how the operators are composed; it does not determine which numerical method must be used inside each substep.

Choosing an Appropriate Time Step

A time step must satisfy more than one requirement. It must resolve the physical processes, keep explicit subsolvers stable, and make the splitting error acceptably small.

For explicit advection, the time step may be limited by a Courant condition. Explicit diffusion methods often have an even stronger restriction connected to the square of the spatial cell size. Explicit reactions may require a small step when reaction rates are large.

Implicit treatment removes some stability restrictions, but it does not remove accuracy requirements. A very large implicit step may remain stable while producing a poor approximation of rapid transients.

A practical convergence study should repeat the simulation with smaller time steps and compare the quantities that matter, such as peak concentration, front position, total mass, or reaction yield.

Adaptive Splitting Strategies

Adaptive methods adjust the time step according to an estimated local error. One practical strategy compares a first-order split result with a second-order Strang result over the same interval.

Another option compares one full step with two half steps. If the solutions differ by more than a selected tolerance, the method rejects the step and tries again with a smaller value.

Adaptive control is useful when the model contains quiet periods followed by rapid reactions, fronts, or other short events. Fixed small steps may waste computation during slow phases, while fixed large steps may miss important dynamics.

The error calculation should include appropriate scaling so that both small and large solution components are evaluated fairly.

Subcycling and Multiple Time Scales

Some operators evolve much faster than others. Operator splitting allows the fast process to use several short internal steps while the slower process advances once.

For example, a reaction solver might take ten small substeps during one larger transport interval:

Reaction: 10 × h/10
Transport: 1 × h

This approach is called subcycling or multirate integration. It can reduce cost when applying the small time step to every operator would be unnecessary.

Subcycling introduces additional design questions. Information exchanged between operators may need interpolation, and the slow process may still influence the fast one during the larger interval. The method should be tested carefully when coupling is strong.

When Operator Splitting Works Well

Splitting is particularly effective when the physical processes can be separated cleanly and specialized solvers already exist.

• The operators interact weakly over one time step.
• Their coefficients and solution fields vary smoothly.
• Different processes require very different numerical methods.
• A monolithic matrix would be too large or expensive.
• The model contains clearly separated time scales.
• The codebase benefits from modular physics components.

Reaction-diffusion systems, reactive transport, atmospheric chemistry, combustion, plasma models, and multiphase simulations frequently use some form of splitting.

When Splitting Becomes Difficult

Operator splitting may require very small time steps when the separated terms are strongly coupled.

Sharp spatial coefficient changes, moving interfaces, rapid nonlinear feedback, and nearly instantaneous equilibrium constraints can increase splitting error. The intermediate state produced after one subsolver may also be physically invalid for the next process.

Conservation can become another concern. Although each subsolver may conserve a quantity independently, the complete composition may not preserve every coupled invariant.

In these situations, possible alternatives include a smaller time step, iterative coupling within each step, an IMEX method, or a fully coupled implicit solver.

Common Implementation Mistakes

A frequent mistake is assuming that Strang splitting makes the complete algorithm second order automatically. Every subsolver must have sufficient accuracy, and boundary conditions must be applied consistently during each substep.

Other common problems include:

• Using an explicit subsolver outside its stability limit
• Applying the wrong operator order without testing alternatives
• Failing to recompute coefficients after another operator changes the solution
• Comparing results only at one time-step size
• Ignoring conservation changes between substeps
• Reusing stale boundary or source data
• Confusing splitting error with spatial discretization error
• Calling a reaction or diffusion term stiff without examining its actual time scale

Choosing a Time-Integration Strategy

A Practical Implementation Workflow

Begin by writing the PDE as a sum of physically meaningful operators. Identify which terms are stiff, which can be treated explicitly, and which already have reliable specialized solvers.

Implement the simplest split method first and verify every subproblem independently. Test conservation, boundary conditions, and expected limiting behavior.

Next, implement Strang splitting and compare it with the first-order sequence. Run time-step refinement tests and, when possible, compare the results with a fully coupled reference solution.

Measure both accuracy and cost. A method that requires fewer steps may still be slower if each split operation performs expensive matrix assembly or data transfer.

Document the operator order, subsolver methods, internal tolerances, and coupling assumptions. These choices are part of the scientific model and should be reproducible.

Conclusion

Operator splitting transforms a coupled PDE system into a sequence of smaller subproblems. This allows advection, diffusion, reaction, and other processes to use numerical methods suited to their individual behavior.

Lie-Trotter splitting is simple but first-order accurate. Strang splitting uses a symmetric half-step, full-step, half-step composition to achieve second-order global accuracy under suitable conditions. IMEX schemes separate stiff and non-stiff terms within a shared additive time integrator rather than solving complete physical operators sequentially.

The effectiveness of splitting depends on the time step, operator ordering, subsolver accuracy, and strength of the coupling. Non-commuting operators introduce splitting error, while sharp gradients and rapid feedback can make that error significant.

For many multiphysics PDEs, Strang splitting provides a practical balance between modularity, computational efficiency, and temporal accuracy. When coupling is too strong for a sequential method, IMEX or monolithic approaches may provide more reliable results.