Components and Data Structures

Introduction

(RefLink[NDSolve, paclet:ref/NDSolve]) is broken up into several basic steps. For advanced usage, it can sometimes be advantageous to access components to carry out each of these steps separately.

(RefLink[NDSolve, paclet:ref/NDSolve]) performs each of these steps internally, hiding the details from a casual user.

Here are the low-level functions that are used to break up these steps.

NDSolve`ProcessEquations classifies the differential system into initial value problem, boundary value problem, differential-algebraic problem, partial differential problem, etc. It also chooses appropriate default integration methods and constructs the main NDSolve`StateData data structure.

NDSolve`Iterate advances the numerical solution. The first invocation (there can be several) initializes the numerical integration methods.

NDSolve`ProcessSolutions converts numerical data into an (RefLink[InterpolatingFunction, paclet:ref/InterpolatingFunction]) to represent each solution.

Note that NDSolve`ProcessEquations can take a significant portion of the overall time to solve a differential system. In such cases, it can be useful to perform this step only once and use NDSolve`Reinitialize to repeatedly solve for different options or initial conditions.

Example

Process equations and set up data structures for solving the differential system:

ndssdata = 
First[NDSolve`ProcessEquations[{y''[t] + y[t] == 0, y[0]==1, y'[0]==0}, {y, y'}, t, Method->"ExplicitRungeKutta"]]

Initialize the method "ExplicitRungeKutta" and integrate the system up to time 10. The return value of NDSolve`Iterate is (RefLink[Null, paclet:ref/Null]) in order to avoid extra references, which would lead to undesirable copying:

NDSolve`Iterate[ndssdata, 10]

Convert each set of solution data into an (RefLink[InterpolatingFunction, paclet:ref/InterpolatingFunction]):

ndsol = NDSolve`ProcessSolutions[ndssdata]

Representing the solution as an (RefLink[InterpolatingFunction, paclet:ref/InterpolatingFunction]) allows continuous output even for points that are not part of the numerical solution grid:

ParametricPlot[{y[t], y'[t]} /. ndsol, {t, 0, 10}]

ProcessEquations

The first stage of any solution using (RefLink[NDSolve, paclet:ref/NDSolve]) is processing the equations specified into a form that can be efficiently accessed by the actual integration algorithms. This stage minimally involves determining the differential order of each variable, making substitutions needed to get a first-order system, solving for the time derivatives of the functions in terms of the functions, and forming the result into a "NumericalFunction" object. If you want to save the time of repeating this process for the same set of equations, or if you want more control over the numerical integration process, the processing stage can be executed separately with NDSolve`ProcessEquations.

Processing equations for (RefLink[NDSolve, paclet:ref/NDSolve]).

This creates a list of two NDSolve`StateData objects because there are two possible solutions for the y' in terms of y:

NDSolve`ProcessEquations[{y'[x]^2 == y[x] + x, y[0] == 1},y,x]

Reinitialize

It is not uncommon that the solution to a more sophisticated problem involves solving the same differential equation repeatedly, but with different initial conditions. In some cases, processing equations may be as time consuming as numerically integrating the differential equations. In these situations, it is a significant advantage to be able to simply give new initial values.

Reusing processed equations.

This creates an NDSolve`StateData object for the harmonic oscillator:

state = First[NDSolve`ProcessEquations[{x''[t] + x[t] == 0, x[0] == 0, x'[0] == 1},x,t]]

This creates three new NDSolve`StateData objects, each with a different initial condition:

newstate = NDSolve`Reinitialize[state, {x[1]^3 == 1, x'[1] == 0}]

Using NDSolve`Reinitialize may save computation time when you need to solve the same differential equation for many different initial conditions, as you might in a shooting method for boundary value problems.

A subset of (RefLink[NDSolve, paclet:ref/NDSolve]) options can be specified as options to NDSolve`Reinitialize.

This creates a new NDSolve`StateData object, specifying a starting step size:

newstate = NDSolve`Reinitialize[state,{x[0] == 0, x'[0] == 1}, StartingStepSize->1/10]

Iterating Solutions

One important use of NDSolve`StateData objects is to have more control of the integration. For some problems, it is appropriate to check the solution and start over or change parameters, depending on certain conditions.

Iterating solutions to differential equations.

This creates an NDSolve`StateData object that contains the information needed to solve the equation for an oscillator with a varying coefficient using an explicit Runge–Kutta method:

state = First[NDSolve`ProcessEquations[{x''[t] + (1 + 4 UnitStep[Sin[t]])x[t] == 0, x[0] == 1, x'[0] == 0},x,t, Method->"ExplicitRungeKutta"]]

Note that when you use NDSolve`ProcessEquations, you do not need to give the range of the t variable explicitly because that information is not needed to set up the equations in a form ready to solve. (For PDEs, you do have to give the ranges of all spatial variables, however, since that information is essential for determining an appropriate discretization.)

This computes the solution out to time t==1:

NDSolve`Iterate[state, 1]

NDSolve`Iterate does not return a value because it modifies the NDSolve`StateData object assigned to the variable state. Thus, the command affects the value of the variable in a manner similar to setting parts of a list, as described in "Manipulating Lists by Their Indices". You can see that the value of state has changed, since it now displays the current time to which it is integrated.

The output form of state shows the range of times over which the solution has been integrated:

state

If you want to integrate further, you can call NDSolve`Iterate again, but with a larger value for time.

This computes the solution out to time t==3:

NDSolve`Iterate[state, 3]

You can specify a time that is earlier than the first current time, in which case the integration proceeds backward with respect to time.

This computes the solution from the initial condition backward to t==-π/2:

NDSolve`Iterate[state, -Pi/2]

NDSolve`Iterate allows you to specify intermediate times at which to stop. This can be useful, for example, to avoid discontinuities. Typically, this strategy is more effective with so-called one-step methods, such as the explicit Runge–Kutta method used in this example. However, it generally works with the default (RefLink[NDSolve, paclet:ref/NDSolve]) method as well.

This computes the solution out to t==10 π, making sure that the solution does not have problems with the points of discontinuity in the coefficients at t==π, 2 π, …:

NDSolve`Iterate[state,π Range[10]]

Getting Solution Functions

Once you have integrated a system up to a certain time, typically you want to be able to look at the current solution values and generate an approximate function representing the solution computed so far. The command NDSolve`ProcessSolutions allows you to do both.

Getting solutions as (RefLink[InterpolatingFunction, paclet:ref/InterpolatingFunction]) objects.

This extracts the solution computed in the previous section as an (RefLink[InterpolatingFunction, paclet:ref/InterpolatingFunction]) object:

sol = NDSolve`ProcessSolutions[state]

This plots the solution:

Plot[Evaluate[x[t] /. sol],{t,0, 10 Pi}]

Just as when using (RefLink[NDSolve, paclet:ref/NDSolve]) directly, there will be a rule for each function you specified in the second argument to NDSolve`ProcessEquations. Only the specified components of the solutions are saved in such a way that an (RefLink[InterpolatingFunction, paclet:ref/InterpolatingFunction]) object can be created.

Obtaining the current solution values.

This gives the current solution values and derivatives in the forward direction:

sol = NDSolve`ProcessSolutions[state,"Forward"]

The choices you can give for the direction dir are "Forward" and "Backward", which refer to the integration forward and backward from the initial condition.

Integration direction specifications.

The output given by NDSolve`ProcessSolution is always given in terms of the dependent variables, either at a specific value of the independent variable, or interpolated over all of the saved values. This means that when a partial differential equation is being integrated, you will get results representing the dependent variables over the spatial variables.

This computes the solution to the heat equation from time t==-1/8 to t==2:

state = First[NDSolve`ProcessEquations[{D[u[t,x],t] == D[u[t,x],x,x],u[0,x] == Cos[π/2 x], u[t,0] == 1 , u[t, 1] == 0},u,t,{x,0,1}]];
NDSolve`Iterate[state, {-1/8,2}]

This gives the solution at t==2:

NDSolve`ProcessSolutions[state, "Forward"]

The solution is given as an (RefLink[InterpolatingFunction, paclet:ref/InterpolatingFunction]) object that interpolates over the spatial variable x.

This gives the solution at t==-1/8:

NDSolve`ProcessSolutions[state, "Backward"]

When you process the current solution for partial differential equations, the spatial error estimate is checked. (It is not generally checked except when solutions are produced because doing so would be quite time consuming.) Since it is excessive, the (RefLink[NDSolve, paclet:ref/NDSolve])::eerr message is issued. The typical association of the word "backward" with the heat equation as implying instability gives a clue to what is wrong in this example.

Here is a plot of the solution at t==-1/8:

Plot[Evaluate[u[-0.125, x] /. %],{x,0,1}]

The plot of the solution shows that instability is indeed the problem.

Even though the heat equation example is simple enough to know that the solution backward in time is problematic, using NDSolve`Iterate and NDSolve`ProcessSolutions to monitor the solution of a PDE can be done to save computing a solution that turns out not to be as accurate as desired. Another simple form of monitoring follows.

Entering the following commands generates a sequence of plots showing the solution of a generalization of the sine-Gordon equation as it is being computed:

L =-10;
state = First[NDSolve`ProcessEquations[{D[u[t,x,y],t,t] == D[u[t,x,y],x,x] + D[u[t,x,y],y,y] - Sin[u[t,x,y]],
u[0,x,y] == Exp[-((x^2)+(y^2))], Derivative[1,0,0][u][0,x,y] == 0,
u[t,-L,y] == u[t,L, y], u[t,x,-L] == u[t,x,L]},u,t,{x,-L,L}, {y,-L,L}, Method->{"MethodOfLines", "SpatialDiscretization"->{"TensorProductGrid", "DifferenceOrder"->"Pseudospectral"}}]];
GraphicsGrid[Partition[Table[
NDSolve`Iterate[state, τ];
Plot3D[Evaluate[u[τ, x, y] /. NDSolve`ProcessSolutions[state, "Forward"]],{x,-L,L},{y,-L,L}, PlotRange->{-1/4,1/4}],
{τ, 0., 20.,5.}],2]]

When you monitor a solution in this way, it is usually possible to interrupt the computation if you see that the solution found is sufficient. You can still use the NDSolve`StateData object to get the solutions that have been computed.

NDSolve`StateData Properties

An NDSolve`StateData object contains a lot of information, but it is arranged in a manner that makes it easy to iterate solutions, and not in a manner that makes it easy to understand where the information is kept. However, sometimes you will want to get information from the state data object: for this reason several properties have been defined to make accessing the information easy.

One of the most important uses of the information from an NDSolve`StateData object is to set up integration methods. Examples are shown in "The NDSolve Method Plugin Framework".

Properties for an NDSolve`StateData object state.

Here is an NDSolve`StateData object for the pendulum equation:

state = First[NDSolve`ProcessEquations[{x''[t] + Sin[x[t]]==0, x[0]==3, x'[0]==0},x,t]]

Get the system size:

nf = state@"SystemSize"

The system size is 2 because there are two scalar ODEs.

Much of the available information depends on the current solution values. (RefLink[NDSolve, paclet:ref/NDSolve]) saves different data for the forward and backward (from the initial condition time) integration directions. These can be accessed by the following directional properties.

Directional properties for an NDSolve`StateData object state.

If the direction argument is omitted, the functions will return a list with the data for both directions (a list with a single element at the initial condition). Otherwise, the direction can be "Forward", "Backward", or "Active" as specified in the previous subsection.

The time step is not initialized until integration has been started unless it was explicitly specified with the (RefLink[StartingStepSize, paclet:ref/StartingStepSize]) option:

state@"TimeStep"

Integrate to t==10 and get the time step in the forward direction:

NDSolve`Iterate[state, 10];
state@"TimeStep"["Forward"]

The number of time steps used in the backward and forward directions:

state@"TimeStepsUsed"

Get the solution data for the pendulum example:

state@"SolutionData"

Solution Data

The solution data consists of a list of the lists of values for each type of solution component. The solution data components correspond to different types of variables, such as time, dependent variables, and discrete variables. Any component that is not being used for a particular problem will be given as an empty list ({}) or (RefLink[None, paclet:ref/None]).

Getting and setting solution data components.

The components of solution data are given in the following table.

Solution data components.

Here is an NDSolve`StateData object for a solution of the pendulum equation with a discrete variable that has been computed up to t==10:

state = First[NDSolve`ProcessEquations[{x''[t] + a[t] Sin[x[t]]==0, x[0]==3, x'[0]==0, a[0]==1, WhenEvent[x[t]==0,a[t]->a[t]+1]},{a,x},t, DiscreteVariables->a]];
NDSolve`Iterate[state, 10];
state

This shows the variables that are used in the problem:

vars = state@"Variables"

This shows the variables that are used in the problem:

vars = state@"VariableData"

To extract components from both variable data and solution data, NDSolve`SolutionDataComponent can be used.

Show the dependent variables:

NDSolve`SolutionDataComponent[vars, "DependentVariables"]

This gets the current solution data in the backward and forward directions:

{sdb,sdf} = state@"SolutionData"

This gets the values of the dependent variables component for the forward direction:

NDSolve`SolutionDataComponent[sdf, "DependentVariables"]

This gets the discrete variables for the backward direction:

NDSolve`SolutionDataComponent[sdb, "DiscreteVariables"]

This gets the indexed discrete variables for the backward direction. The result is an empty list because no indexed discrete variables are used in this computation:

NDSolve`SolutionDataComponent[sdb, "IndexedDiscreteVariables"]

Note that the solution data is effectively the raw data that (RefLink[NDSolve, paclet:ref/NDSolve]) uses to compute the solution. In the example above, the dependent variables component is a list of length 2 because the second-order equation is automatically reduced to a first-order system. If you want to have a clear identification of the parts of the current solution, you can use NDSolve`ProcessSolutions.

Get solution data associated with the variables via rules:

NDSolve`ProcessSolutions[state, "Forward"]

NumericalFunction

When doing time integration, the methods (RefLink[NDSolve, paclet:ref/NDSolve]) uses work with a right-hand side function for ODEs such that x'==f(t,x) or a residual function for DAEs such that f(t,x,x')==0, where t is real and x, x', and other solution data components used are vectors. The actual variables represented by the vectors x and x' may have more complicated structure, so (RefLink[NDSolve, paclet:ref/NDSolve]) uses an intermediate (RefLinkPlain[NumericalFunction, paclet:tutorial/NDSolveStateData#1128081347]) object that takes vector arguments and returns vectors based on any underlying function.

Get the NDSolve`StateData object for a system with vector and scalar variables:

state = First[NDSolve`ProcessEquations[{y'[t]== {{.1,.2},{.3,.4}}.y[t] -  10 Sign[z[t]], z'[t]==y[t].y[t],y[0]=={2,1},z[0]==-1},{y,z},t]]

Get the (RefLinkPlain[NumericalFunction, paclet:tutorial/NDSolveStateData#1128081347]) that evaluates the derivatives:

nf = state@"NumericalFunction"

The (RefLinkPlain[NumericalFunction, paclet:tutorial/NDSolveStateData#1128081347]) that (RefLink[NDSolve, paclet:ref/NDSolve]) constructs for evaluation only uses as arguments the solution data components that are actually used in the problem. You can find what components it actually uses by using the "SolutionDataComponents" property.

Find which solution data components are required for evaluation

comp = nf@"SolutionDataComponents"

Get the initial solution data and the evaluation arguments

sd = state@"SolutionData"[0]

args = NDSolve`SolutionDataComponent[sd, comp]

Note that the indexed discrete variable was introduced as a "DiscontinuitySignature" variable by the automatic discontinuity handling for the (RefLink[Sign, paclet:ref/Sign]) function. (See the solution shown below.)

Evaluating the function gives the time derivatives, since the system is an ODE:

Apply[nf, args]

% == NDSolve`SolutionDataComponent[sd,"X'"]

To make it easier to evaluate with the solution data for a given problem, two convenience functions have been included that evaluate with the solution data.

Evaluating with solution data.

Evaluate the function using the solution data list:

NDSolve`EvaluateWithSolutionData[nf,sd]

One of the things that (RefLinkPlain[NumericalFunction, paclet:tutorial/NDSolveStateData#1128081347]) does is to handle derivatives automatically. By default, symbolic derivatives are used if they can be found, and otherwise finite differences are used.

Evaluate the Jacobian with respect to the dependent variables:

NDSolve`EvaluateJacobianWithSolutionData[nf,sd,"X"]

In this example, finite differences are used because of the vector arguments and (RefLink[Dot, paclet:ref/Dot]) products.

Integrate and plot the solution for z and z':

NDSolve`Iterate[state, 1];
Plot[Evaluate[{z[t],z'[t]} /. NDSolve`ProcessSolutions[state]],{t,0,.2}, PlotRange->All]

It is sometimes useful to set up (RefLinkPlain[NumericalFunction, paclet:tutorial/NDSolveStateData#1128081347]) to evaluate some function of the (RefLink[NDSolve, paclet:ref/NDSolve]) variables, which can be done using "ProcessExpression".

Many of the built-in methods do this during method initialization so that the function can be evaluated efficiently repeated times. For example, the "projection" method sets up a function to evaluate the deviation of the invariants from their initial values and uses the derivatives to implement the projection after each step.

Set up an NDSolve`StateData object for an oscillator:

state = First[NDSolve`ProcessEquations[{x''[t] +a[t]^2 x[t]==0,WhenEvent[x[t]>0,a[t]->a[t]+1],x[0]==1,x'[0]==0,a[0]==1},x,t,DiscreteVariables->a, MaxSteps->Infinity]];

Make a (RefLinkPlain[NumericalFunction, paclet:tutorial/NDSolveStateData#1128081347]) for the energy integral:

efun = state@"ProcessExpression"[Automatic, (1/2)(((x'[t])^2) +(((a[t] x[t]))^2)) ,{}]

Evaluate at the initial condition:

einit = NDSolve`EvaluateWithSolutionData[efun, state@"SolutionData"["Forward"]]

Make a plot of the energy monitored at regular intervals during the integration:

ListPlot[edata = Table[NDSolve`Iterate[state,τ];{τ, NDSolve`EvaluateWithSolutionData[efun, state@"SolutionData"["Forward"]] - einit},{τ,0,50}]]

The energy deviation is actually small considering how quickly the solution is oscillating by t==50:

Plot[Evaluate[x[t] /. NDSolve`ProcessSolutions[state]],{t,49.99,50}]