"Projection" Method for NDSolve

Introduction

When a differential system has a certain structure, it is advantageous if a numerical integration method preserves the structure. In certain situations it is useful to solve differential equations in which solutions are constrained. Projection methods work by taking a time step with a numerical integration method and then projecting the approximate solution onto the manifold on which the true solution evolves.

(RefLink[NDSolve, paclet:ref/NDSolve]) includes a differential algebraic solver which may be appropriate and is described in more detail within "Numerical Solution of Differential-Algebraic Equations".

Sometimes the form of the equations may not be reduced to the form required by a DAE solver. Furthermore so-called index reduction techniques can destroy certain structural properties, such as symplecticity, that the differential system may possess (see [HW96] and [HLW02]). An example that illustrates this can be found in the documentation for DAEs.

In such cases it is often possible to solve a differential system and then use a projective procedure to ensure that the constraints are conserved. This is the idea behind the method "Projection".

If the differential system is ρ-reversible then a symmetric projection process can be advantageous (see [H00]). Symmetric projection is generally more costly than projection and has not yet been implemented in (RefLink[NDSolve, paclet:ref/NDSolve]).

Invariants

Consider a differential equation

(OverDot[y])==f(y),y((Subscript[t, 0]))==(Subscript[y, 0]),

where y may be a vector or a matrix.

Definition: A nonconstant function I(y) is called an invariant of () if I'(y)f(y) == 0 for all y.

This implies that every solution y(t) of () satisfies I(y(t))==I((Subscript[y, 0]))==Constant.

Synonymous with invariant, the terms first integral, conserved quantity, or constant of the motion are also common.

Manifolds

Given an (n-m)-dimensional submanifold of (ℝ^n) with g:(ℝ^n)↦(ℝ^m):

ℳ=={y;g(y)==0}.

Given a differential equation () then (Subscript[y, 0])∈ℳ implies y(t)∈ℳ for all t. This is a weaker assumption than invariance and g(y) is called a weak invariant (see [HLW02]).

Projection Algorithm

Let (Subscript[(Overscript[y, ∼]), n+1]) denote the solution from a one-step numerical integrator. Considering a constrained minimization problem leads to the following system (see [AP91], [HW96] and [HLW02]):

({{(Subscript[y, n+1]), ==, (Subscript[(Overscript[y, ∼]), n+1])+g'((((Subscript[y, n+1])))^T)λ}, {0, ==, g((Subscript[y, n+1])).}})

To save work g((Subscript[y, n+1])) is approximated as g((Subscript[(Overscript[y, ∼]), n+1])). Substituting the first relation into the second relation in () leads to the following simplified Newton scheme for λ:

({{(Subscript[Δλ, i]) , ==, -((( g'((Subscript[(Overscript[y, ∼]), n+1])) g'((((Subscript[(Overscript[y, ∼]), n+1])))^T)))^(-1))g((Subscript[(Overscript[y, ∼]), n+1]) + g'((((Subscript[(Overscript[y, ∼]), n+1])))^T) (Subscript[λ, i])) ,}, {(Subscript[λ, i+1]), ==, (Subscript[λ, i])+ (Subscript[Δλ, i]) }})

with (Subscript[λ, 0])==0.

The first increment (Subscript[Δλ, 0]) is of size O((Subscript[h, n]^(p+1))) so that () usually converges quickly.

The added expense of using a higher-order integration method can be offset by fewer Newton iterations in the projective step.

For the termination criterion in the method "Projection", the option "IterationSafetyFactor" is combined with one Unit in the Last Place in the working precision used by (RefLink[NDSolve, paclet:ref/NDSolve]).

Examples

Load some utility packages:

Needs["DifferentialEquations`NDSolveProblems`"];
Needs["DifferentialEquations`NDSolveUtilities`"];

Linear Invariants

Define a stiff system modeling a chemical reaction:

system = GetNDSolveProblem["Robertson"];
vars = system["DependentVariables"];

This system has a linear invariant:

invariant = system["Invariants"]

Linear invariants are generally conserved by numerical integrators (see [S86]), including the default (RefLink[NDSolve, paclet:ref/NDSolve]) method, as can be observed in a plot of the error in the invariant:

sol = NDSolve[system];

InvariantErrorPlot[invariant,vars,T,sol]

Therefore in this example there is no need to use the method "Projection".

Certain numerical methods preserve quadratic invariants exactly (see for example [C87]). The implicit midpoint rule, or one-stage Gauss implicit Runge–Kutta method, is one such method.

Harmonic Oscillator

Define the harmonic oscillator:

system = GetNDSolveProblem["HarmonicOscillator"];
vars = system["DependentVariables"];

The harmonic oscillator has the following invariant:

invariant = system["Invariants"]

Solve the system using the method "ExplicitRungeKutta". The error in the invariant grows roughly linearly, which is typical behavior for a dissipative method applied to a Hamiltonian system:

erksol = NDSolve[system, Method->"ExplicitRungeKutta"];

InvariantErrorPlot[invariant,vars,T,erksol]

This also solves the system using the method "ExplicitRungeKutta" but it projects the solution at the end of each step. A plot of the error in the invariant shows that it is conserved up to roundoff:

projerksol = NDSolve[system,Method->{"Projection",Method->"ExplicitRungeKutta","Invariants"->invariant}];

InvariantErrorPlot[invariant,vars,T,projerksol]

Since the system is Hamiltonian (the invariant is the Hamiltonian), a symplectic integrator performs well on this problem, giving a small bounded error:

projerksol = NDSolve[system,Method->{"SymplecticPartitionedRungeKutta","DifferenceOrder"->8,"PositionVariables"->{(Subscript[Y, 1])[T]}},StartingStepSize->1/5];

InvariantErrorPlot[invariant,vars,T,projerksol]

Perturbed Kepler Problem

This loads a Hamiltonian system known as the perturbed Kepler problem, sets the integration interval and the step size to take, as well as defining the position variables in the Hamiltonian formalism:

system = GetNDSolveProblem["PerturbedKepler"];
time = system["TimeData"];
step = 3/100;
pvars = Take[system["DependentVariables"], 2]

The system has two invariants, which are defined as H and L:

{H,L} = system["Invariants"]

An experiment now illustrates the importance of using all the available invariants in the projective process (see [HLW02]). Consider the solutions obtained using:

sol = NDSolve[system, Method->"ExplicitEuler", StartingStepSize->step];

ParametricPlot[Evaluate[pvars /. First[sol]], Evaluate[time]]

sol = NDSolve[system, Method->{"Projection", Method->"ExplicitEuler",
"Invariants"->{H}}, StartingStepSize->step];

ParametricPlot[Evaluate[pvars /. First[sol]], Evaluate[time]]

sol = NDSolve[system, Method->{"Projection", Method->"ExplicitEuler",
"Invariants"->{L}}, StartingStepSize->step];

ParametricPlot[Evaluate[pvars /. First[sol]], Evaluate[time]]

sol = NDSolve[system, Method->{"Projection", Method->"ExplicitEuler",
"Invariants"->{H,L}}, StartingStepSize->step];

ParametricPlot[Evaluate[pvars /. First[sol]], Evaluate[time]]

It can be observed that only the solution with projection onto both invariants gives the correct qualitative behavior—for comparison, results using an efficient symplectic solver can be found in "SymplecticPartitionedRungeKutta Method for NDSolve".

Lotka Volterra

An example of constraint projection for the Lotka–Volterra system is given within "Numerical Methods for Solving the Lotka–Volterra Equations".

Euler's Equations

An example of constraint projection for Euler's equations is given within "Rigid Body Solvers".

Option Summary

Options of the method "Projection".