"Composition" and "Splitting" Methods for NDSolve

Introduction

In some cases it is useful to split the differential system into subsystems and solve each subsystem using appropriate integration methods. Recombining the individual solutions often allows certain dynamical properties, such as volume, to be conserved. More information on splitting and composition can be found in [MQ02, HLW02], and specific aspects related to (RefLink[NDSolve, paclet:ref/NDSolve]) are discussed in [SS05, SS06].

Definitions

Of concern are initial value problems y'(t)==f(y(t)), where y(0)==(Subscript[y, 0])∈(ℝ^n).

"Composition"

Composition is a useful device for raising the order of a numerical integration scheme.

In contrast to the Aitken–Neville algorithm used in extrapolation, composition can conserve geometric properties of the base integration method (e.g. symplecticity).

Let (Subscript[Φ, (Subscript[f, ,])(Subscript[γ, i])h]^((i))) be a basic integration method that takes a step of size (Subscript[γ, i])h with (Subscript[γ, 1]),…,(Subscript[γ, s]) given real numbers.

Then the s-stage composition method (Subscript[Ψ, f,h]) is given by

Often interest is in composition methods (Subscript[Ψ, f,h]) that involve the same base method Φ==(Φ^((i))),i==1,…,s.

An interesting special case is symmetric composition: (Subscript[γ, i])==(Subscript[γ, s-i+1]), i==1,…,⌊s/2⌋.

The most common types of composition are:

"Splitting"

An s-stage splitting method is a generalization of a composition method in which f is broken up in an additive fashion:

The essential point is that there can often be computational advantages in solving problems involving (Subscript[f, i]) instead of f.

An s-stage splitting method is a composition of the form

with (Subscript[f, 1]),…,(Subscript[f, s]) not necessarily distinct.

Each base integration method now only solves part of the problem, but a suitable composition can still give rise to a numerical scheme with advantageous properties.

If the vector field (Subscript[f, i]) is integrable, then the exact solution or flow (Subscript[φ, (Subscript[f, i]),h]) can be used in place of a numerical integration method.

A splitting method may also use a mixture of flows and numerical methods.

An example is Lie–Trotter splitting [T59]:

Split f==(Subscript[f, 1])+(Subscript[f, 2]) with (Subscript[γ, 1])==(Subscript[γ, 2])==1; then (Subscript[Ψ, f,h])==(Subscript[φ, (Subscript[f, 2]),h]^((2)))◦(Subscript[φ, (Subscript[f, 1]),h]^((1))) yields a first-order integration method.

Computationally it can be advantageous to combine flows using the group property

Implementation

Several changes to the new (RefLink[NDSolve, paclet:ref/NDSolve]) framework were needed in order to implement splitting and composition methods.

A simplified input syntax allows omitted vector fields and methods to be filled in cyclically. These must be defined unambiguously:

{(Subscript[f, 1]),(Subscript[f, 2]),(Subscript[f, 1]),(Subscript[f, 2])} can be input as {(Subscript[f, 1]),(Subscript[f, 2])}.

{(Subscript[f, 1]),(Subscript[f, 2]),(Subscript[f, 3]),(Subscript[f, 2]),(Subscript[f, 1])} cannot be input as {(Subscript[f, 1]),(Subscript[f, 2]),(Subscript[f, 3])} since this corresponds to {(Subscript[f, 1]),(Subscript[f, 2]),(Subscript[f, 3]),(Subscript[f, 1]),(Subscript[f, 2])}.

Nested Methods

The following example constructs a high-order splitting method from a low-order splitting using "Composition".

(RefLink[NDSolve, paclet:ref/NDSolve]) ({{⟶, "Composition" ({{↗, "Splitting"  f=(Subscript[f, 1])+(Subscript[f, 2])  ({{↗, "LocallyExact"  (Subscript[f, 1])}, {⟶, ({{ImplicitMidpoint  (Subscript[f, 2]), }})}, {↘, ({{"LocallyExact"  (Subscript[f, 1]), }})}})}, {⋮, ⋮}, {⟶, "Splitting"  f=(Subscript[f, 1])+(Subscript[f, 2])  ({{↗, "LocallyExact"  (Subscript[f, 1])}, {⟶, ({{ImplicitMidpoint  (Subscript[f, 2]), }})}, {↘, ({{"LocallyExact"  (Subscript[f, 1]), }})}})}, {⋮, ⋮}, {↘, "Splitting"  f=(Subscript[f, 1])+(Subscript[f, 2])  ({{↗, "LocallyExact"  (Subscript[f, 1])}, {⟶, ({{ImplicitMidpoint  (Subscript[f, 2]), }})}, {↘, ({{"LocallyExact"  (Subscript[f, 1]), }})}})}})}})

Simplification

A more efficient integrator can be obtained in the previous example using the group property of flows and calling the "Splitting" method directly.

(RefLink[NDSolve, paclet:ref/NDSolve]) ({{⟶, "Splitting"  f=(Subscript[f, 1])+(Subscript[f, 2])({{↗, ({{"LocallyExact"  (Subscript[f, 1]) }, {ImplicitMidpoint  (Subscript[f, 2])}})}, {⋮, ⋮}, {⟶, ({{"LocallyExact"  (Subscript[f, 1]) }, {ImplicitMidpoint  (Subscript[f, 2])}, {"LocallyExact"  (Subscript[f, 1]) }}) }, {⋮, ⋮}, {↘, ({{ImplicitMidpoint  (Subscript[f, 2])}, {"LocallyExact"  (Subscript[f, 1]) }})  }})}})

Examples

The following examples will use a second-order symmetric splitting known as the Strang splitting [S68], [M68]. The splitting coefficients are automatically determined from the structure of the equations.

Load a package with some useful example problems:

Needs["DifferentialEquations`NDSolveProblems`"];

This defines a method known as symplectic leapfrog in terms of the method "SymplecticPartitionedRungeKutta":

SymplecticLeapfrog = {"SymplecticPartitionedRungeKutta", "DifferenceOrder"->2, "PositionVariables":>qvars};

Symplectic Splitting

Symplectic Leapfrog

"SymplecticPartitionedRungeKutta" is an efficient method for solving separable Hamiltonian systems H(p,q)==T(p)+V(q) with favorable long-time dynamics.

"Splitting" is a more general-purpose method, but it can be used to construct partitioned symplectic methods (though it is somewhat less efficient than "SymplecticPartitionedRungeKutta").

Consider the harmonic oscillator that arises from a linear differential system that is governed by the separable Hamiltonian H(p,q)==(p^2)/2+(q^2)/2:

system = GetNDSolveProblem["HarmonicOscillator"]

Split the Hamiltonian vector field into independent components governing momentum and position. This is done by setting the relevant right-hand sides of the equations to zero:

eqs = system["System"];
Y1 = eqs;
Part[Y1, 1,2] = 0;
Y2 = eqs;
Part[Y2, 2,2] = 0;

This composition of weighted (first-order) Euler integration steps corresponds to the symplectic (second-order) leapfrog method:

tfinal = 1;
time = {T, 0, tfinal};
qvars = {Subscript[Y,1][T]};
splittingsol = NDSolve[system, time, StartingStepSize->1/10,
 Method->{"Splitting","DifferenceOrder"->2, "Equations"->{Y1,Y2,Y1},
"Method"->{"ExplicitEuler", "ExplicitEuler", "ExplicitEuler"}}]

The method "ExplicitEuler" could only have been specified once, since the second and third instances would have been filled in cyclically.

This is the result at the end of the integration step:

InputForm[splittingsol /. T->tfinal]

This invokes the built-in integration method corresponding to the symplectic leapfrog integrator:

sprksol = NDSolve[system, time, StartingStepSize->1/10, Method->SymplecticLeapfrog]

The result at the end of the integration step is identical to the result of the splitting method:

InputForm[sprksol /. T->tfinal]

Composition of Symplectic Leapfrog

This takes the symplectic leapfrog scheme as the base integration method and constructs a fourth-order symplectic integrator using a symmetric composition of Ruth–Yoshida [Y90]:

YoshidaCoefficients = RootReduce[{1/(2-2^(1/3)),-2^(1/3)/(2-2^(1/3)),1/(2-2^(1/3))}];

YoshidaCompositionCoefficients[4, p_]:= N[YoshidaCoefficients, p];

splittingsol = NDSolve[system,time, StartingStepSize->1/10,
 Method->{"Composition", "Coefficients"->YoshidaCompositionCoefficients,
"DifferenceOrder"->4, "Method"->{SymplecticLeapfrog}}]

This is the result at the end of the integration step:

InputForm[splittingsol /. T->tfinal]

This invokes the built-in symplectic integration method using coefficients for the fourth-order methods of Ruth and Yoshida:

SPRK4[4,prec_]:=N[{{Root[-1+12*#1-48*#1^2+48*#1^3&,1,0],Root[1-24*#1^2+48*#1^3&,1,0],Root[1-24*#1^2+48*#1^3&,1,0],Root[-1+12*#1-48*#1^2+48*#1^3&,1,0]},{Root[-1+6*#1-12*#1^2+6*#1^3&,1,0],Root[1-3*#1+3*#1^2+3*#1^3&,1,0],Root[-1+6*#1-12*#1^2+6*#1^3&,1,0],0}},prec];

 sprksol = NDSolve[system,time, StartingStepSize->1/10, 
Method->{"SymplecticPartitionedRungeKutta", "Coefficients"->SPRK4,"DifferenceOrder"->4,"PositionVariables"->qvars}]

The result at the end of the integration step is identical to the result of the composition method:

InputForm[sprksol /. T->tfinal]

Hybrid Methods

While a closed-form solution often does not exist for the entire vector field, in some cases it is possible to analytically solve a system of differential equations for part of the vector field.

When a solution can be found by (RefLink[DSolve, paclet:ref/DSolve]), direct numerical evaluation can be used to locally advance the solution.

This idea is implemented in the method (RefLinkPlain["LocallyExact", paclet:tutorial/NDSolveLocallyExact]).

Harmonic Oscillator Test Example

This example checks that the solution for the exact flows of split components of the harmonic oscillator equations is the same as applying Euler's method to each of the split components:

system = GetNDSolveProblem["HarmonicOscillator"];
eqs = system["System"];
Y1 = eqs;
Part[Y1, 1,2] = 0;
Y2 = eqs;
Part[Y2, 2,2] = 0;
tfinal = 1;
time = {T, 0, tfinal};

solexact = NDSolve[system, time, StartingStepSize->1/10,
 Method->{NDSolve`Splitting, "DifferenceOrder"->2, "Equations"->{Y1,Y2,Y1}, "Method"->{"LocallyExact"}}];

InputForm[solexact /. T->1]

soleuler = NDSolve[system, time, StartingStepSize->1/10,
 Method->{NDSolve`Splitting, "DifferenceOrder"->2, "Equations"->{Y1,Y2,Y1}, "Method"->{"ExplicitEuler"}}];

InputForm[soleuler /. T->tfinal]

Hybrid Numeric-Symbolic Splitting Methods (ABC Flow)

Consider the Arnold, Beltrami, and Childress flow, a widely studied model for volume-preserving three-dimensional flows:

system = GetNDSolveProblem["ArnoldBeltramiChildress"]

When applied directly, a volume-preserving integrator would not in general preserve symmetries. A symmetry-preserving integrator, such as the implicit midpoint rule, would not preserve volume.

This defines a splitting of the system by setting some of the right-hand side components to zero:

eqs = system["System"];
Y1 = eqs;
Part[Y1, 2,2] = 0;
Y2 = eqs;
Part[Y2, {1,3},2] = 0;

Y1

Y2

The system for Y2 is solvable exactly by (RefLink[DSolve, paclet:ref/DSolve]) so that you can use the (RefLinkPlain["LocallyExact", paclet:tutorial/NDSolveLocallyExact]) method.

Y1 is not solvable, however, so you need to use a suitable numerical integrator in order to obtain the desired properties in the splitting method.

This defines a method for computing the implicit midpoint rule in terms of the built-in "ImplicitRungeKutta" method:

ImplicitMidpoint ={"FixedStep",Method-> {"ImplicitRungeKutta", "Coefficients"->"ImplicitRungeKuttaGaussCoefficients","DifferenceOrder"->2, ImplicitSolver->{FixedPoint,AccuracyGoal->MachinePrecision, PrecisionGoal->MachinePrecision,"IterationSafetyFactor"->1}}};

This defines a second-order, volume-preserving, reversing symmetry-group integrator [MQ02]:

splittingsol = NDSolve[system,
StartingStepSize->1/10,
 Method->{"Splitting", "DifferenceOrder"->2,
"Equations"->{Y2,Y1,Y2}, "Method"->{"LocallyExact",ImplicitMidpoint,"LocallyExact"}}]

Lotka–Volterra Equations

Various numerical integrators for this system are compared within "Numerical Methods for Solving the Lotka–Volterra Equations".

Euler's Equations

Various numerical integrators for Euler's equations are compared within "Rigid Body Solvers".

Non-Autonomous Vector Fields

Consider the Duffing oscillator, a forced planar non-autonomous differential system:

system = GetNDSolveProblem["DuffingOscillator"]

This defines a splitting of the system:

Y1 = {(Subscript[Y, 1])'[T]==(Subscript[Y, 2])[T],(Subscript[Y, 2])'[T]==(((Subscript[Y, 2])[T])/4)};
Y2 = {(Subscript[Y, 1])'[T]==0,(Subscript[Y, 2])'[T]==((3 Cos[T])/10)+(Subscript[Y, 1])[T]-(((Subscript[Y, 1])[T])^3)};

The splitting of the time component among the vector fields is ambiguous, so the method issues an error message:

splittingsol = NDSolve[system, StartingStepSize->1/10,
 Method->{"Splitting","DifferenceOrder"->2,
"Equations"->{Y2,Y1,Y1},"Method"->{"LocallyExact"}}]

The equations can be extended by introducing a new "dummy" variable Z[T] such that Z[T]==T and specifying how it should be distributed in the split differential systems:

Y1 = {(Subscript[Y, 1])'[T]==(Subscript[Y, 2])[T],(Subscript[Y, 2])'[T]==(((Subscript[Y, 2])[T])/4), Z'[T]==1};
Y2 = {(Subscript[Y, 1])'[T]==0,(Subscript[Y, 2])'[T]==((3 Cos[Z[T]])/10)+(Subscript[Y, 1])[T]-(((Subscript[Y, 1])[T])^3), Z'[T]==0};
eqs = Join[system["System"],{Z'[T]==1}];
ics = Join[system["InitialConditions"], {Z[0]==1}];
vars = Join[system["DependentVariables"], {Z[T]}];
time = system["TimeData"];

This defines a geometric splitting method that satisfies (Subscript[λ, 1])+(Subscript[λ, 2])==-δ for any finite time interval, where (Subscript[λ, 1]) and (Subscript[λ, 2]) are the Lyapunov exponents [MQ02]:

splittingsol = NDSolve[{eqs, ics}, vars,time, StartingStepSize->1/10,
 Method->{NDSolve`Splitting, "DifferenceOrder"->2,
"Equations"->{Y2,Y1,Y2}, "Method"->{"LocallyExact"}}]

Here is a plot of the solution:

ParametricPlot[Evaluate[system["DependentVariables"[]] /. First[splittingsol]],Evaluate[time],AspectRatio->1]

Option Summary

The default coefficient choice in "Composition" tries to automatically select between "SymmetricCompositionCoefficients" and "SymmetricCompositionSymmetricMethodCoefficients" depending on the properties of the methods specified using the (RefLink[Method, paclet:ref/Method]) option.

Options of the method "Composition".

Options of the method "Splitting".