Utility Packages for Numerical Differential Equation Solving

InterpolatingFunctionAnatomy

(RefLink[NDSolve, paclet:ref/NDSolve]) returns solutions as (RefLink[InterpolatingFunction, paclet:ref/InterpolatingFunction]) objects. Most of the time, simply using these as functions does what is needed, but occasionally it is useful to access the data inside, which includes the actual values and points (RefLink[NDSolve, paclet:ref/NDSolve]) computed when taking steps. The exact structure of an (RefLink[InterpolatingFunction, paclet:ref/InterpolatingFunction]) object is arranged to make the data storage efficient and evaluation at a given point fast. This structure may change between Wolfram Language versions, so code that is written in terms of accessing parts of (RefLink[InterpolatingFunction, paclet:ref/InterpolatingFunction]) objects may not work with new versions of the Wolfram Language. The DifferentialEquations`InterpolatingFunctionAnatomy` package provides an interface to the data in an (RefLink[InterpolatingFunction, paclet:ref/InterpolatingFunction]) object that will be maintained for future Wolfram Language versions.

Anatomy of (RefLink[InterpolatingFunction, paclet:ref/InterpolatingFunction]) objects.

This loads the package:

Needs["DifferentialEquations`InterpolatingFunctionAnatomy`"];

One common situation where the InterpolatingFunctionAnatomy package is useful is when (RefLink[NDSolve, paclet:ref/NDSolve]) cannot compute a solution over the full range of values that you specified, and you want to plot all of the solution that was computed to try to understand better what might have gone wrong.

Here is an example of a differential equation which cannot be computed up to the specified endpoint:

ifun = First[x /. NDSolve[{x'[t] == Exp[x[t]] - x[t], x[0] == 1},x,{t,0,10}]]

This gets the domain:

domain = InterpolatingFunctionDomain[ifun]

Once the domain has been returned in a list, it is easy to use (RefLink[Part, paclet:ref/Part]) to get the desired endpoints and make the plot:

{begin, end} = domain[[1]];
Plot[ifun[t],{t,begin, end}]

From the plot, it is quite apparent that a singularity has formed and it will not be possible to integrate the system any further.

Sometimes it is useful to see where (RefLink[NDSolve, paclet:ref/NDSolve]) took steps. Getting the coordinates is useful for doing this.

This shows the values that (RefLink[NDSolve, paclet:ref/NDSolve]) computed at each step it took. It is quite apparent from this that nearly all of the steps were used to try to resolve the singularity:

coords = First[InterpolatingFunctionCoordinates[ifun]];
ListPlot[Transpose[{coords, ifun[coords]}]]

The package is particularly useful for analyzing the computed solutions of PDEs.

With this initial condition, Burgers' equation forms a steep front:

mdfun = First[u /.NDSolve[{D[u[x,t],t] == 0.01 D[u[x,t],x,x] - u[x,t] D[u[x,t],x], u[0, t] == u[1,t], u[x,0] == Sin[2 Pi x]},u,{x,0,1},{t,0,0.5}]]

This shows the number of points used in each dimension:

Map[Length, InterpolatingFunctionCoordinates[mdfun]]

This shows the interpolation order used in each dimension:

InterpolatingFunctionInterpolationOrder[mdfun]

This shows that the inability to resolve the front has manifested itself as numerical instability:

Max[Abs[InterpolatingFunctionValuesOnGrid[mdfun]]]

This shows the values computed at the spatial grid points at the endpoint of the temporal integration:

end = InterpolatingFunctionDomain[mdfun][[2,-1]];
x = InterpolatingFunctionCoordinates[mdfun][[1]];
ListPlot[Transpose[{x, mdfun[x, end]}], PlotStyle->PointSize[.025], PlotRange->{-1,1}]

It is easily seen from the point plot that the front has not been resolved.

This makes a 3D plot showing the time evolution for each of the spatial grid points. The initial condition is shown in red:

Show[Graphics3D[{Map[Line, MapThread[Append, {InterpolatingFunctionGrid[mdfun], InterpolatingFunctionValuesOnGrid[mdfun]}, 2]], {RGBColor[1,0,0],Line[Transpose[{x, 0. x, mdfun[x,0.]}]]}}], BoxRatios->{1,1,1}, PlotRange->{All, All, {-1,1}}]

When a derivative of an (RefLink[InterpolatingFunction, paclet:ref/InterpolatingFunction]) object is taken, a new (RefLink[InterpolatingFunction, paclet:ref/InterpolatingFunction]) object is returned that gives the requested derivative when evaluated at a point. The InterpolatingFunctionDerivativeOrder is a way of determining what derivative will be evaluated.

The derivative returns a new (RefLink[InterpolatingFunction, paclet:ref/InterpolatingFunction]) object:

dmdfun = Derivative[0,1][mdfun]

This shows what derivative will be evaluated:

InterpolatingFunctionDerivativeOrder[dmdfun]

NDSolveUtilities

A number of utility routines have been written to facilitate the investigation and comparison of various (RefLink[NDSolve, paclet:ref/NDSolve]) methods. These functions have been collected in the package DifferentialEquations`NDSolveUtilities`.

Functions provided in the NDSolveUtilities package.

This loads the package:

Needs["DifferentialEquations`NDSolveUtilities`"]

A useful means of analyzing Runge–Kutta methods is to study how they behave when applied to a scalar linear test problem (see the package FunctionApproximations.m).

This assigns the (exact or infinitely precise) coefficients for the 2-stage implicit Runge–Kutta Gauss method of order 4:

{amat, bvec, cvec} =NDSolve`ImplicitRungeKuttaGaussCoefficients[4, Infinity]

This computes the linear stability function, which corresponds to the (2,2) Padé approximation to the exponential at the origin:

RungeKuttaLinearStabilityFunction[amat, bvec, z]

Examples of the functions CompareMethods, FinalSolutions, RungeKuttaLinearStabilityFunction, and StepDataPlot can be found within "ExplicitRungeKutta Method for NDSolve".

Examples of the function InvariantErrorPlot can be found within "Projection Method for NDSolve".

InvariantErrorPlot Options

The function InvariantErrorPlot has a number of options that can be used to control the form of the result.

Options of the function InvariantErrorPlot.

The default value for InvariantDimensions is to determine the dimensions from the structure of the input, (RefLink[Dimensions, paclet:ref/Dimensions])[invts].

The default value for InvariantErrorFunction is a function to compute the absolute error.

The default value for InvariantErrorSampleRate is to sample all points if there are less than 1000 steps taken. Above this threshold a logarithmic sample rate is used.