A weighted five-stage strang splitting scheme for the viscous burgers equation

Introduction

Burgers' equation, uₜ + u uₓ = ν uₓₓ, x ∈ ℝ, t > 0, (1) is the simplest nonlinear partial differential equation that combines advective steepening with diffusive smoothing, and it has served for decades as a testbed for numerical methods aimed at the Navier-Stokes equations and other advection-diffusion systems. A natural way to solve it numerically is operator splitting: write the right-hand side as the sum of a nonlinear advection operator A(u) = -u uₓ and a linear diffusion operator D(u) = ν uₓₓ, and approximate the flow of A + D by composing the flows of A and D separately, each of which is much easier to solve on its own than the full equation. The simplest composition, A(h)D(h), is only first order accurate because A and D do not commute. Strang's well-known remedy (Strang, 1968) [4] is the symmetric composition

S(h) = A(h/2) D(h) A(h/2),

which is second order accurate because it is palindromic: reading the stages in reverse order and reversing the sign of h reproduces the same map. It is a classical fact (see e.g. Hairer, Lubich & Wanner, 2006, Thm. II.3.3) [1] that any symmetric one-step method has even order, and this is the only reason the two outer half-steps of A are given equal weight 1/2: equal weighting is the simplest way to make the composition symmetric, not a requirement of symmetry itself.

In this paper we point out that this freedom is normally left on the table. We use a five-stage palindromic composition

S a (h) = A(ah) D(h/2) A((1-2a)h) D(h/2) A(ah),

which is symmetric -and hence at least second order -for every choice of the weight a ∈ (0, 1/2), because reversing the five stages and sending h to -h reproduces Sa(h) regardless of the value of a. The classical scheme is recovered at a = 1/2, which collapses the middle advection stage to zero width. For any other a, the scheme costs exactly the same two diffusion solves and the same total amount of advection integration -the work is simply distributed differently across the three advection sub-steps. The question we answer is: how should a be chosen to make the most of this freedom?

We answer this with a short, explicit calculation. A Baker-Campbell-Hausdorff (BCH) expansion of Sa(h) shows that the leading local error, of size h³, is a linear combination of two independent double commutators, with coefficients c1(a) and c2(a) that we compute explicitly. Because one commutator is weighted more heavily by advection and the other by diffusion, the two coefficients respond differently to changes in a, and balancing them against each other gives a simple closed-form rule for a nearoptimal weight, expressed in terms of the relative strength of diffusion and advection in the problem -in effect, a Reynolds-or Péclet-number-dependent splitting weight. The resulting scheme is no more expensive than classical Strang splitting; it simply spends the same work differently. This idea is related to, but distinct from, the balanced splitting of Speth, Green, MacNamara & Strang (2013) [3] , who improve the accuracy of Strang splitting near a steady state by shifting mass between the two operators, A → A + cn, D → D -cn, with an additive constant recomputed at every step. Our scheme instead leaves the operators untouched and rebalances only the time fractions assigned to a fixed five-stage composition, with a single weight chosen once from the problem data; the two ideas address different mechanisms and could in principle be combined. Section 2 sets up the splitting and the finite-difference discretization used in our experiments. Section 3 derives the order-2 property and the closed-form leading error, and proposes the weight-balancing rule. Section 4 reports numerical experiments confirming the predicted order of convergence and the practical benefit of the asymmetric weight. Section 5 discusses the scope.

We write A(τ) for the exact (or numerically accurate) flow of uₜ = A(u) over a time interval of length τ, and likewise D(τ) for the flow of uₜ = D(u). Because the inviscid Burgers flow A(τ) and the heat-equation flow D(τ) do not commute, any composition of them that approximates the flow of A + D carries a splitting error in addition to whatever error is introduced by the numerical scheme used inside each sub-step.

A weighted five-stage composition

Instead of the classical three-stage Strang composition (2), we use the five-stage composition (3), repeated here for convenience:

The three advection weights a, 1-2a, a sum to 1, as required for consistency, and the composition is symmetric under reversal for every a: reversing the stage order and replacing h by -h returns Sa(h). This is the only property of (2) that Strang's order-2 argument actually uses, so (5) inherits the same order-2 guarantee -a fact we re-derive directly in Section 3, because we also need the size of the leading error term, which the abstract symmetry argument alone does not provide. Setting a = 1/2 collapses the middle advection stage to width zero and reduces (5) to two consecutive copies of the classical Strang scheme (2) run over half a step each, which is algebraically equivalent to (2) over a full step (a standard fact about Strang-type compositions; see McLachlan & Quispel, 2002, §2.3); we verify this numerically in Section 4 and use it as our baseline. The point of this paper is that a need not be 1/2: any other value still gives a valid, second-order scheme at the same cost, and the freedom in a can be spent to reduce the size of the leading error term.

Discretization used in the experiments

To isolate the splitting error from other sources of error, we solve each sub-step as accurately as practical.

Advection sub-step. The inviscid equation uₜ + u uₓ = 0 is discretized in space with the standard conservative central difference for uₜ + (u²/2)ₓ = 0 on a periodic grid, and integrated in time with classical fourth-order Runge-Kutta using many small internal sub-steps, so that the temporal error of this stage is negligible compared with the splitting error under study. Diffusion sub-step. The linear equation uₜ = ν uₓₓ, discretized in space with the standard second-order central Laplacian, is a constant-coefficient linear ODE system that is diagonalized exactly by the discrete Fourier transform on a periodic grid; we advance it by multiplying each Fourier mode by its exact exponential factor. This removes the time-discretization error of the diffusion sub-step entirely, so that any error remaining in our experiments is due only to the splitting itself and to the (negligible) temporal error of the advection sub-step. Neither choice is essential to the method: on a non-periodic domain the diffusion sub-step can be solved with any accurate implicit scheme (e.g. Crank-Nicolson), and the advection sub-step with any accurate scheme for the inviscid Burgers equation (e.g. a high-resolution Godunov-type or WENO scheme near shocks). The splitting analysis of Section 3 is independent of these choices and concerns only the time-splitting error between the two operators.

Error analysis 3.1 Order of accuracy for every weight

We treat A and D as abstract, generally non-commuting operators and expand both the exact flow E(h) = exp(h(A + D)) and the split flow Sa(h) of ( 5) as formal power series in h, using the standard Baker-Campbell-Hausdorff calculus for products of operator exponentials. Collecting terms order by order in h (we verified this expansion to h³ symbolically), one finds that the h⁰ and h¹ terms of Sa(h) -E(h) vanish identically, as required for consistency, and that the h² term is

where the commutator is [A, D] = AD DA. This term disappears for each and every option of a, not only a = 1/2 because the symmetry argument and the composition (5) are symmetric for each a the direct BCH computation agree, as they must. This demonstrates that Sa(h) is correct at the second order. for the entire weight family, with no upper limit on a beyond a

The leading error term

The h³ term of Sa(h) -E(h) is a sum of triple products of A and D, and it can be written exactly as a combination of two independent double commutators,

We verified ( 6)-( 8) by symbolic expansion, multiplying out the operator exponentials of ( 5) and ( 1) as formal series in noncommuting symbols and matching coefficients order by order; we omit the routine but lengthy intermediate algebra. As a check, setting a = 1/2 gives c1(1/2) = -1/24 and c2(1/2) = -1/12, which is exactly the classical leading-order error constant of the threestage Strang scheme (2) (Hairer, Lubich & Wanner, 2006, Ch. II.3) [1] ; this is a useful sanity check on the longer five-stage calculation.

The two commutators in (7) are not interchangeable in size. Roughly speaking, [A, [A, D]] involves two applications of the advection operator and one of diffusion, so its norm scales like (advection strength)² × (diffusion strength); symmetrically, [D, [A, D]] scales like (diffusion strength)² × (advection strength). When advection dominates diffusion -the high-Reynolds-number regime that makes Burgers' equation interesting -the first commutator is generally the more dangerous one, and there is no reason that the value of a minimizing its coefficient c1(a) should coincide with a = 1/2.

A balancing rule for the weight

To turn (8) into a practical recommendation we introduce a single dimensionless ratio r that compares the strength of diffusion to advection. We use the straightforward, easy-to-calculate proxy throughout our experiments.

the opposite of the cell Péclet number, which is the natural local measure of how much diffusion occurs in a given area. on a particular grid, competes with advection. We then choose a so that the two terms in (7) are equal in magnitude and opposite in sign -a simple balance condition,

which, using (8), is a quadratic equation in a with explicit solution a*(r) = 1/2 + r/4 -√(9r² + 24r + 12) / 12.

(The other root of the quadratic lies above 1/2 and is discarded, since we restrict a to (0, 1/2) by symmetry of labeling.) The formula (11) is admittedly only a heuristic: it balances two terms of a multi-term error expression using a single scalar proxy for the relative operator strengths, rather than minimizing a true operator norm of the full local error. Because it is closed-form, its value is nothing to evaluate and close to the empirically best weight, as Section 4 demonstrates. It is important to note two limitations of (11) a*(0) = 1/2 as r 0 (advection dominates completely). a*(r) = 1/6 + 1/(18r) = 1/6 = 0.167 as r (where diffusion is predominant). In both areas, Across the entire range, a* stays well below the classical value 1/2 and only slightly changes with r-The balancing rule's primary message is not so much the exact value of a as it is the fact that In virtually every regime of interest for Burgers, the ideal weight is significantly less than half.

Numerical results

Set-up

With 256 grid points and initial data of u0(x), we use the periodic domain [0, 2] to solve Burgers' equation. sin x + 0.5 sin 2x, viscosity = 0.05, and the final time T = 0.5, when the profile is clearly visible steepened (Figure 3) without yet establishing a true discontinuity, assuming smoothness Section 3 in the BCH expansion is still a reasonable approximation. An example solution iscomputed by integrating the unsplit equation directly with fourth-order Runge-Kutta at a very small time step (Δt = 10⁻⁴), well below the step sizes used in the splitting experiments.

Convergence in time

The classical scheme (a = 1/2) and the asymmetric scheme (a = 0.16), which is the over a sequence of halving time steps, the empirically best weight was found in 4.3 below. Both strategies converge at the predicted second order, with an empirical order that is close to 2 as t 0 confirming the 3.1 analysis. The asymmetric scheme is significantly more symmetric at each step size tested. more precise, by a ratio that increases from approximately 16 at the largest step size to approximately 38 at the smallest. at small t, they are parallel to the reference line, indicating second-order convergence; the asymmetric scheme consistently falls below the classical scheme (a = 1/2) with a value of 0.16.

Fig 1:

Log-log plot of the data in Table 1, together with a reference slope-2 line. Both curves are parallel to the reference line at small Δt, confirming second-order convergence; the asymmetric scheme (a = 0.16) sits consistently below the classical scheme (a = 1/2). The error has a pronounced, a well-defined minimum close to a 0.16, which is far from the traditional value of a = 1/2, and rises on either side of the bare minimum, which is in line with the quadratic dependence on a that was predicted by The balancing formula (11) for this problem uses the grid and U 1.5 to estimate r. spacing x = 2/256) yields a* 0.19, which is close to but not exactly the same as the empirical optimum. Given that (11) uses only a single scalar proxy for the relative strength of two operators that act on a whole function, this level of agreement is about as much as one should expect from a closed-form heuristic, and it is sufficient to identify the right neighborhood for a without any tuning.

Locating the best weight

Sensitivity to viscosity

Table 2 repeats the weight scan of §4.3 for several values of ν, spanning a twenty-fold range, and reports the empirically best weight found at each ν, alongside the value predicted by (11). Both the empirical optimum and the theoretical prediction stay confined to a narrow band, roughly a ∈ [0.10, 0.20], across the whole range of ν tested. This is itself a useful practical finding: rather than re-tuning a for every new value of the viscosity, a single fixed choice in this range -we used a = 0.16 throughout §4.2-4.3 -already captures most of the available benefit.

Solution behavior

The initial data and the solution at T = 0.5 used in this section are depicted in Figure 3. the characteristic, partially smoothed advection-induced steepening of the leading edge of each lobe by having a viscosity of 0.05. This is the regime, which is well short of triggering a shock but has significant local steepening was the setting for our experiments. The BCH expansion and the commutators that are assumed in this analysis to have adequate smoothness appear in it-to have significance, which is reasonable away from shock formation but becomes doubtful once ux gets big and the advection sub-step needs to use shock-capturing. equipment, such as limiters and WENO reconstructions, whose own truncation error is not smooth in u. Coupling the weight a to a local measure of or extending the analysis through shock formation solution steepness rather than a single global ratio r, is a natural next step that we have no followed here. In a second sense, the balancing rule (11) is also just a heuristic because it equates the magnitudes of two. scalar coefficients are the sum of two different operators multiplied by a single number, r, acting as a substitute for the ratio of their actual solution's operator norms. A more rigorous (but more expensive) alternative would estimate the norms of the two commutators directly from the discrete solution at each time step and minimize their combined contribution numerically; our experiments suggest that the extra effort buys little here, since both the theoretical and the empirical optimum stay in a narrow range across a wide span of viscosities (Table 2), but this may not hold for every problem in this family. Finally, the same five-stage idea is not specific to Burgers' equation. Any advection-diffusion or, more generally, any twooperator splitting that currently uses the classical equal-weight Strang composition could be re-examined in the same way: the order-2 property survives for any palindromic weighting, and the leading-error balancing argument of §3.3 carries over essentially unchanged, with A and D re-interpreted as whatever two operators are being split. The reaction-diffusion and convection-reaction systems studied via balanced splitting by Speth, Green, MacNamara & Strang (2013) [3] would be natural candidates for this comparison.

Conclusion

The classical Strang splitting of Burgers' equation gives equal weight to its two advection half-steps because equal weight is the simplest way to make the composition symmetric, not because it minimizes the splitting error. Writing the same idea as a fivestage composition with a free weight a preserves second-order accuracy for every a, at no extra cost, and exposes a genuine free parameter that can be set by a short, explicit calculation. The resulting closed-form rule (11) is simple enough to compute once per problem from the viscosity, the velocity scale, and the grid spacing, and our experiments show that using it -or simply moving the weight to a ≈ 0.15-0.20 instead of 1/2 -reduces the error by more than an order of magnitude at the same computational cost as the classical scheme.

References

1. Hairer E, Lubich C, Wanner G. Geometric Numerical Integration: Structure-Preserving Algorithms for Ordinary Differential Equations. 2nd ed. Berlin: Springer; 2006.
2. McLachlan RI, Quispel GRW. Splitting methods. Acta Numerica. 2002;11:341-434.
3. Speth RL, Green WH, MacNamara S, Strang G. Balanced splitting and rebalanced splitting. SIAM Journal on Numerical Analysis. 2013;51(6):3084-3105.
4. Strang G. On the construction and comparison of difference schemes. SIAM Journal on Numerical Analysis. 1968;5(3):506517.
5. Trotter HF. On the product of semi-groups of operators. Proceedings of the American Mathematical Society. 1959;10(4):545551.