"LocallyExact" Method for NDSolve

Introduction

A differential system can sometimes be solved by analytic means. The function (RefLink[DSolve, paclet:ref/DSolve]) implements many of the known algorithmic techniques.

However, differential systems that can be solved in closed form constitute only a small subset. Despite this fact, when a closed-form solution does not exist for the entire vector field, it is often possible to analytically solve a system of differential equations for part of the vector field. An example of this is the method "Splitting", which breaks up a vector field f into subfields  (Subscript[f, 1]),…,(Subscript[f, n]) such that f==(Subscript[f, 1])+⋯+(Subscript[f, n]).

The idea underlying the method "LocallyExact" is that rather than using a standard numerical integration scheme, when a solution can be found by (RefLink[DSolve, paclet:ref/DSolve]) direct numerical evaluation can be used to locally advance the solution.

Since the method "LocallyExact" makes no attempt to adaptively adjust step sizes, it is primarily intended for use as a submethod between integration steps.

Examples

Load a package with some predefined problems:

Needs["DifferentialEquations`NDSolveProblems`"];

Harmonic Oscillator

Numerically solve the equations of motion for a harmonic oscillator using the method "LocallyExact". The result is two interpolating functions that approximate the solution and the first derivative:

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

sols = vars /. First[NDSolve[system, StartingStepSize->1/10, Method->"LocallyExact"]]

The solution evolves on the unit circle:

ParametricPlot[Evaluate[sols],Evaluate[tdata],AspectRatio->1]

Global versus Local

The method "LocallyExact" is not intended as a substitute for a closed-form (global) solution.

Despite the fact that the method "LocallyExact" uses the analytic solution to advance the solution, it only produces solutions at the grid points in the numerical integration (or even inside grid points if called appropriately). Therefore, there can be errors due to sampling at interpolation points that do not lie exactly on the numerical integration grid.

Plot the error in the first solution component of the harmonic oscillator and compare it with the exact flow:

Plot[Evaluate[First[sols]-Cos[T]],Evaluate[tdata]]

Simplification

The method "LocallyExact" has an option "SimplificationFunction" that can be used to simplify the results of (RefLink[DSolve, paclet:ref/DSolve]).

Here is the linearized component of the differential system that turns up in the splitting of the Lorenz equations using standard values for the parameters:

eqs={(Subscript[Y, 1])'[T]==σ ((Subscript[Y, 2])[T]-(Subscript[Y, 1])[T]),(Subscript[Y, 2])'[T]==r (Subscript[Y, 1])[T]-(Subscript[Y, 2])[T],(Subscript[Y, 3])'[T]==-b (Subscript[Y, 3])[T]}/.{σ->10,r->28,b->8/3};
ics = {(Subscript[Y, 1])[0]==-8,(Subscript[Y, 2])[0]==8,(Subscript[Y, 3])[0]==27};
vars={(Subscript[Y, 1])[T],(Subscript[Y, 2])[T],(Subscript[Y, 3])[T]};

This subsystem is exactly solvable by (RefLink[DSolve, paclet:ref/DSolve]):

DSolve[eqs,vars,T]

Often the results of (RefLink[DSolve, paclet:ref/DSolve]) can be simplified. This defines a function to simplify an expression and also prints out the input and the result:

myfun[x_]:=
Module[{simpx},
Print["Before simplification ",x];
simpx = FullSimplify[ExpToTrig[x]];
Print["After simplification ", simpx];
simpx
];

The function can be passed as an option to the method "LocallyExact":

NDSolve[{eqs, ics}, vars, {T, 0, 1}, StartingStepSize->1/10, Method->{"LocallyExact", "SimplificationFunction"->myfun}]

The simplification is performed only once during the initialization phase that constructs the data object for the numerical integration method.

Option Summary

Option of the method "LocallyExact".