The summer continues to allow some clearing of the backlog of projects accumulated during the academic year: Helge Holden, Kenneth Karlsen, Nils Risebro, and myself have uploaded to the arXiv the paper “Operator splitting for the KdV equation“, submitted to Math. Comp..   This paper is concerned with accurate numerical schemes for solving initial value problems for the Korteweg-de Vries equation

u_t + u_{xxx} = u u_x; \quad u(0,x) = u_0(x), (1)

though the analysis here would be valid for a wide range of other semilinear dispersive models as well.  In principle, these equations, which are completely integrable, can be solved exactly by the inverse scattering method, but fast and accurate implementations of this method are still not as satisfactory as one would like.  On the other hand, the linear Korteweg-de Vries equation

u_t + u_{xxx} = 0; \quad u(0,x) = u_0(x) (2)

can be solved exactly (with accurate and fast numerics) via the (fast) Fourier transform, while the (inviscid) Burgers equation

u_t = u u_x; \quad u(0,x) = u_0(x) (3)

can also be solved exactly (and quickly) by the method of characteristics.  Since the KdV equation is in some sense a combination of the equations (2) and (3), it is then reasonable to hope that some combination of the solution schemes for (2) and (3) can be used to solve (1), at least in some approximate sense.

One way to do this is by the method of operator splitting.  Observe from the formal approximation e^{A \Delta t} \approx 1 + A \Delta t + O( \Delta t^2) (where \Delta t should be thought of as small, and A is some matrix or linear operator), that one has

e^{(A+B) \Delta t} u_0 = e^{A \Delta t} e^{B \Delta t} u_0 + O( \Delta t^2 ), (4)

[we do not assume A and B to commute here] and thus we formally have

e^{n (A+B) \Delta t} u_0 = (e^{A \Delta t} e^{B \Delta t})^n u_0 + O( \Delta t ) (5)

if n = T/\Delta t for some fixed time T (thus n = O(1/\Delta t)).  As a consequence, if one wants to solve the linear ODE

u_t = (A+B) u; \quad u(0) = u_0 (1′)

for time T = n \Delta t, one can achieve an approximate solution (accurate to order \Delta t) by alternating n times between evolving the ODE

u_t = A u (2′)

for time \Delta t, and evolving the ODE

u_t = B u (3′)

for time \Delta t, starting with the initial data u_0.

It turns out that this scheme can be formalised, and furthermore generalised to nonlinear settings such as those for the KdV equation (1).  More precisely, we show that if u_0 \in H^s({\Bbb R}) for some s \geq 5, then one can solve (1) to accuracy O(\Delta t) in H^s norm for any fixed time T = n \Delta t by alternating between evolving (2) and (3) for times \Delta t (this scheme is known as Godunov splitting).

Actually, one can obtain faster convergence by modifying the scheme, at the cost of requiring higher regularity on the data; the situation is similar to that of numerical integration (or quadrature), in which the midpoint rule or Simpson’s rule provide more accuracy than the Riemann integral if the integrand is smooth.  For instance, one has the variant

e^{n (A+B) \Delta t} = (e^{A \Delta t/2} e^{B \Delta t/2} e^{B \Delta t/2} e^{A \Delta t/2})^n + O( \Delta t^2 ) (6)

of (5), which can be seen by expansion to second order in \Delta t (or by playing around with the Baker-Campbell-Hausdorff formula).  For KdV, we can rigorously show that the analogous scheme (known as Strang splitting) involving the indicated combination of evolutions of (2) and (3) will also converge to accuracy \Delta t^2 in H^s norm, provided that u_0 \in H^s({\Bbb R}) and s \geq 17.

Our main tools are the energy method (i.e. computations of how H^s({\Bbb R})-norm type energies of the solution evolve in time via differentiation and integration by parts), and the use of two time variables rather than one.  To describe the second method, let us give a PDE-style proof of (4) (assuming that A, B are finite-dimensional matrices and u is vector-valued to avoid the technical complications associated with infinite dimensional vector spaces).  We introduce two time variables t, \tau and consider the PDE problem

\partial_\tau u(t,\tau) = A u(t,\tau)

\partial_t u(t,0) = B u(t,0)

u(0,0) = u_0.

Thus, one first solves the equation (3′) along the t axis with initial data u_0, and then solves the equation (2′) in the \tau direction, thus leading to the explicit solution u(t,\tau) = e^{\tau A} e^{t B} u_0.  In particular, the left-hand side of (4) is u(\Delta t, \Delta t).

Now we look at the evolution of u along the diagonal axis.  Observe from the chain rule that

\frac{d}{dt} u(t,t) = (A+B) u(t,t) + F(t,t) (7)

where F(t,\tau) :=\partial_t u(t,\tau) - B u(t,\tau).  On the other hand, F itself solves an ODE in the \tau direction:

\partial_\tau F = A F + [A,B] u

F(t,0) = 0.

We thus have F(t,t) = O( t ) for small t, and then by Duhamel’s formula and (7) we have u(t,t) = e^{(A+B)t} u_0 + O( t^2 ), which implies (4) as claimed.  Note that this analysis also recovers the classical formula e^{(A+B)t} = e^{At} e^{Bt} in the special case that A and B commute.

About these ads