IDA Method for NDSolve

The IDA package is part of the SUNDIALS (SUite of Nonlinear and DIfferential/ALgebraic equation Solvers) developed at the Center for Applied Scientific Computing of Lawrence Livermore National Laboratory. As described in the IDA user guide [HT99], "IDA is a general purpose solver for the initial value problem for systems of differential-algebraic equations (DAEs). The name IDA stands for Implicit Differential-Algebraic solver. IDA is based on DASPK ...". DASPK [BHP94], [BHP98] is a Fortran code for solving large-scale differential-algebraic systems.

In the Wolfram Language, an interface has been provided to the IDA package so that rather than needing to write a function in C for evaluating the residual and compiling the program, the Wolfram Language generates the function automatically from the equations you input to (RefLink[NDSolve, paclet:ref/NDSolve]).

IDA solves the system (1) with Backward Differentiation Formula (BDF) methods of orders 1 through 5, implemented using a variable-step form. The BDF of order k at time (Subscript[t, n])==(Subscript[t, n-1])+(Subscript[h, n]) is given by the formula

The coefficients (Subscript[a, n,i]) depend on the order k and past step sizes. Applying the BDF to the DAE (1) gives a system of nonlinear equations to solve:

The solution of the system is achieved by Newton-type methods, typically using an approximation to the Jacobian

J==((∂F)/(∂x))+(Subscript[c, n])((∂F)/(∂x')), where (Subscript[c, n])==((Subscript[α, n,0])/(Subscript[h, n])).

"Its [IDA's] most notable feature is that, in the solution of the underlying nonlinear system at each time step, it offers a choice of Newton/direct methods or an Inexact Newton/Krylov (iterative) method." [HT99] In the Wolfram Language, you can access these solvers using method options or use the default Wolfram Language (RefLink[LinearSolve, paclet:ref/LinearSolve]), which switches automatically to direct sparse solvers for large problems.

At each step of the solution, IDA computes an estimate (Subscript[E, n]) of the local truncation error, and the step size and order are chosen so that the weighted norm (RefLink[Norm, paclet:ref/Norm])[(Subscript[E, n])/(Subscript[w, n])] is less than 1. The j(^th)component, (Subscript[w, n,j]), of (Subscript[w, n]) is given by

The values prec and acc are taken from the (RefLink[NDSolve, paclet:ref/NDSolve]) settings for the (RefLink[PrecisionGoal, paclet:ref/PrecisionGoal])->prec and (RefLink[AccuracyGoal, paclet:ref/AccuracyGoal])->acc.

Because IDA provides a great deal of flexibility, particularly in the way nonlinear equations are solved, there are a number of method options which allow you to control how this is done. You can use the method options to IDA by giving (RefLink[NDSolve, paclet:ref/NDSolve]) the option (RefLink[Method, paclet:ref/Method])->{IDA,ida method options}.

The options for the IDA method are associated with the symbol IDA in the NDSolve` context:

Options[NDSolve`IDA]

IDA method options.

When strict accuracy of intermediate values computed with the (RefLink[InterpolatingFunction, paclet:ref/InterpolatingFunction]) object returned from (RefLink[NDSolve, paclet:ref/NDSolve]) is important, you will want to use the (RefLink[NDSolve, paclet:ref/NDSolve]) method option setting (RefLink[InterpolationOrder, paclet:ref/InterpolationOrder])->(RefLink[All, paclet:ref/All]) that uses interpolation based on the order of the method, sometimes called dense output, to represent the solution between time steps. By default (RefLink[NDSolve, paclet:ref/NDSolve]) stores a minimal amount of data to represent the solution well enough for graphical purposes. Keeping the amount of data small saves on both memory and time for more complicated solutions.

As an example which highlights the difference between minimal output and full method interpolation order, consider tracking a quantity f(t)==((x(t))^2)+((y(t))^2) derived from the solution of a simple linear equation where the exact solution can be computed using (RefLink[DSolve, paclet:ref/DSolve]).

This defines the function f giving the quantity as a function of time based on solutions x[t] and y[t]:

f[t_]  :=((x[t])^2)+((y[t])^2);

This defines the linear equations along with initial conditions:

eqns = {x'[t]==x[t]-2 y[t],y'[t]==x[t]+y[t]};
ics = {x[0] == 1, y[0] == 1};

The exact value of f as a function of time can be computed symbolically using (RefLink[DSolve, paclet:ref/DSolve]):

(Subscript[f, exact])[t_] = First[f[t]/.DSolve[{eqns, ics},{x,y},t]]

The exact solution will be compared with solutions computed with and without dense output.

A simple way to track the quantity is to create a function which derives it from the numerical solution of the differential equation:

f1[t_] = First[f[t] /. NDSolve[{eqns, ics}, {x,y},{t,0,1}]]

It can also be computed with dense output:

(Subscript[f1, dense])[t_] = First[f[t] /. NDSolve[{eqns, ics}, {x,y},{t,0,1}, InterpolationOrder->All]]

This plot shows the error in the two computed solutions. The computed solution at the time steps are indicated by black dots. The default output error is shown in gray and the dense output error in black:

Needs["DifferentialEquations`InterpolatingFunctionAnatomy`"];t1 =Cases[f1[t], (if_InterpolatingFunction)[t]->InterpolatingFunctionCoordinates[if], Infinity][[1,1]];
pode = Show[Block[{$DisplayFunction = Identity},{ListPlot[Transpose[{t1, (Subscript[f, exact])[t1] - f1[t1]}], PlotStyle->PointSize[.02]],Plot[(Subscript[f, exact])[t]- f1[t],{t,0,1}, PlotStyle->RGBColor[.8,.8,.8]],Plot[(Subscript[f, exact])[t]- (Subscript[f1, dense])[t],{t,0,1}]}], PlotRange->All]

From the plot, it is quite apparent that when the time steps get large, the default solution output has a much larger error between time steps. The dense output solution represents the solution computed by the solver even between time steps. Keep in mind, however, that the dense output solution takes up much more space.

This compares the sizes of the default and dense output solutions:

ByteCount /@ {f1[t], (Subscript[f1, dense])[t]}

When the quantity you want to derive from the solution is complicated, you can ensure that it is locally kept within tolerances by giving it as an algebraic quantity, forcing the solver to keep its error in control.

This adds a dependent variable with an algebraic equation that sets the dependent variable equal to the quantity of interest:

f2[t_] = First[g[t] /. NDSolve[{eqns, ics, g[t] == f[t]}, {x,y,g},{t,0,1}]]

This computes the same solution with dense output:

(Subscript[f2, dense])[t_] = First[g[t] /. NDSolve[{eqns, ics, g[t] == f[t]}, {x,y,g},{t,0,1}, InterpolationOrder->All]]

This makes a plot comparing the error for all four solutions. The time steps for IDA are shown as blue points and the dense output from IDA is shown in blue with the default output shown in light blue:

t2 =InterpolatingFunctionCoordinates[Head[f2[t]]][[1]];
Show[{pode, ListPlot[Transpose[{t2, (Subscript[f, exact])[t2] - f2[t2]}], PlotStyle->{RGBColor[0,0,1],PointSize[0.02]}],Plot[(Subscript[f, exact])[t]- f2[t],{t,0,1}, PlotStyle->RGBColor[.7,.7,1]],Plot[(Subscript[f, exact])[t]- (Subscript[f2, dense])[t],{t,0,1}, PlotStyle->RGBColor[0,0,1]]}, PlotRange->{{0,1}, 1*^-6 {-1,1}}]

You can see from the plot that the error is somewhat smaller when the quantity is computed algebraically along with the solution.

The remainder of this documentation will focus on suboptions of the two possible settings for the "ImplicitSolver" option, which can be "Newton" or "GMRES". With "Newton", the Jacobian or an approximation to it is computed and the linear system is solved to find the Newton step. On the other hand, "GMRES" is an iterative method and, rather than computing the entire Jacobian, a directional derivative is computed for each iterative step.

The "Newton" method has one method option, "LinearSolveMethod", which you can use to tell the Wolfram Language how to solve the linear equations required to find the Newton step. There are several possible values for this option.

Possible settings for the "LinearSolveMethod" option.

The "GMRES" method may be substantially faster, but is typically quite a bit more tricky to use because to really be effective it typically requires a preconditioner, which can be supplied via a method option. There are also some other method options that control the Krylov subspace process. To use these, refer to the IDA user guide [HT99].

"GMRES" method options.

As an example problem, consider a two-dimensional Burgers' equation

This can typically be solved with an ordinary differential equation solver, but in this example two things are achieved by using the DAE solver. First, boundary conditions are enforced as algebraic conditions. Second, (RefLink[NDSolve, paclet:ref/NDSolve]) is forced to use conservative differencing by using an algebraic term. For comparison, a known exact solution will be used for initial and boundary conditions.

This defines a function that satisfies Burger's equation:

Bsol[t_, x_, y_] = 1/(1+Exp[(x+y-t)/(2 ν)]);

This defines initial and boundary conditions for the unit square consistent with the exact solution:

ic = u[0,x,y] == Bsol[0,x,y];
bc = {
u[t,0,y] == Bsol[t,0,y], u[t,1,y] == Bsol[t,1,y],
u[t,x,0] == Bsol[t,x,0], u[t,x,1] == Bsol[t,x,1]};

This defines the differential equation:

de = D[u[t,x,y],t] == ν( D[u[t,x,y],x,x] +D[u[t,x,y],y,y]) -u[t,x,y](D[u[t,x,y],x] + D[u[t,x,y],y]);

This sets the diffusion constant ν to a value for which you can find a solution in a reasonable amount of time and shows a plot of the solution at t==1:

ν = 0.025;
Plot3D[Bsol[1,x,y],{x,0,1},{y,0,1}]

You can see from the plot that with ν==0.025, there is a fairly steep front. This moves with constant speed.

This solves the problem using the default settings for (RefLink[NDSolve, paclet:ref/NDSolve]) and the IDA method with the exception of the "DifferentiateBoundaryConditions" option for "MethodOfLines", which causes (RefLink[NDSolve, paclet:ref/NDSolve]) to treat the boundary conditions as algebraic:

Timing[sol = NDSolve[{de, ic, bc},u,{t,0,1},{x,0,1},{y,0,1}, Method->{"MethodOfLines", "DifferentiateBoundaryConditions"->False}]]

Since there is an exact solution to compare to, the overall solution accuracy can be compared as well.

This defines a function that finds the maximum deviation between the exact and computed solutions at the grid points at all of the time steps:

errfun[sol_] := Module[{ifun = First[u /. sol], grid, exvals, gvals},
grid = InterpolatingFunctionGrid[ifun];
gvals = InterpolatingFunctionValuesOnGrid[ifun];
exvals = Apply[Bsol,Transpose[grid, RotateLeft[Range[ArrayDepth[grid]],1]]];
Max[Abs[exvals - gvals]]]

This computes the maximal error for the computed solution:

errfun[sol]

In the following, a comparison will be made with different settings for the options of the IDA method. To emphasize the option settings, a function will be defined to time the computation of the solution and give the maximal error.

This defines a function for comparing different IDA option settings:

TimeSolution[idaopts___] := Module[{time, err, steps},
time = First[Timing[sol = NDSolve[{de, ic, bc},u,{t,0,1},{x,0,1},{y,0,1}, Method->{"MethodOfLines", "DifferentiateBoundaryConditions"->False, Method->{IDA, idaopts}}]]];
err = errfun[sol];
steps = Length[First[InterpolatingFunctionCoordinates[First[u /. sol]]]] "Steps";
{time, err, steps}]

No options use the previous defaults:

TimeSolution[]

This uses the "Band" method:

TimeSolution["ImplicitSolver"->{"Newton", "LinearSolveMethod"->"Band"}]

The "Band" method is not very effective because the bandwidth of the Jacobian is relatively large, partly because of the fourth-order derivatives used and partly because the one-sided stencils used near the boundary add width at the top and bottom. You can specify the bandwidth explicitly.

This uses the "Band" method with the width set to include the stencil of the differences in only one direction:

TimeSolution["ImplicitSolver"->{"Newton", "LinearSolveMethod"->{"Band", "BandWidth"->3}}]

While the solution time was smaller, notice that the error is slightly greater and the total number of time steps is a lot greater. If the problem had been more stiff, the iterations likely would not have converged because it was missing information from the other direction. Ideally, the bandwidth should not eliminate information from an entire dimension.

This computes the grids used in the X and Y directions and shows the number of points used:

{X,Y} = InterpolatingFunctionCoordinates[First[u /. sol]][[{2,3}]];
{nx,ny} = {Length[X], Length[Y]}

This uses the "Band" method with the width set to include at least part of the stencil in both directions:

TimeSolution["ImplicitSolver"->{"Newton", "LinearSolveMethod"->{"Band", "BandWidth"->51}}]

With the more appropriate setting of the bandwidth, the solution is still slightly slower than in the default case. The "Band" method can sometimes be effective on two-dimensional problems, but is usually most effective on one-dimensional problems.

This computes the solution using the "GMRES" implicit solver without a preconditioner:

TimeSolution["ImplicitSolver"->"GMRES"]

This is incredibly slow! Using the "GMRES" method without a preconditioner is not recommended for this very reason. However, finding a good preconditioner is not usually trivial. For this example, a diagonal preconditioner will be used.

The setting of the "Preconditioner" option should be a function f, which accepts four arguments that will be given to it by (RefLink[NDSolve, paclet:ref/NDSolve]) such that f[t,x,x',c] returns another function that will apply the preconditioner to the residual vector. (See IDA user guide [HT99] for details on how the preconditioner is used.) The arguments t, x, x', c are the current time, solution vector, solution derivative vector, and the constant c in formula (2) above. For example, if you can determine a procedure that would generate an appropriate preconditioner matrix P as a function of these arguments, you could use

"Preconditioner"->(RefLink[Function, paclet:ref/Function])[{t,x,xp,c},(RefLink[LinearSolve, paclet:ref/LinearSolve])[P[t,x,xp,c]]]

to produce a (RefLink[LinearSolveFunction, paclet:ref/LinearSolveFunction]) object which will effectively invert the preconditioner matrix P. Typically, for each time the preconditioner function is set up, it is applied to the residual vector several times, so using some sort of factorization such as is contained in a (RefLink[LinearSolveFunction, paclet:ref/LinearSolveFunction]) is a good idea.

For the diagonal case, the inverse can be affected simply by using the reciprocal. The most difficult part of setting up a diagonal preconditioner is keeping in mind that values on the boundary should not be affected by it.

This finds the diagonal elements of the differentiation matrix for computing the preconditioner:

DM = NDSolve`FiniteDifferenceDerivative[{2,0},{X,Y}]@"DifferentiationMatrix" + NDSolve`FiniteDifferenceDerivative[{0,2},{X,Y}]@"DifferentiationMatrix";
Short[diag = Tr[DM,List]]

This gets the positions where elements at the boundary that satisfy a simple algebraic condition are in the flattened solution vector:

bound = SparseArray[{{i_, 1}->1., {i_, ny}->1.,{1,i_}->1.,{nx,i_}->1.},{nx,ny},0.];
Short[pos = Drop[ArrayRules[Flatten[bound]],-1][[All,1,1]]]

This defines the function that sets up the function called to get the effective inverse of the preconditioner. For the diagonal case, the inverse is done simply by taking the reciprocal:

pfun[t_, x_, xp_, c_] := Module[{d,dd},
d = 1./(c - ν diag);
d[[pos]] = 1.;
Function[# dd] /. dd->d]

This uses the preconditioned "GMRES" method to compute the solution:

TimeSolution["ImplicitSolver"->{"GMRES", "Preconditioner"->pfun}]

Thus, even with a crude preconditioner, the "GMRES" method computes the solution faster than using the direct sparse solvers.

For PDE discretizations with higher-order temporal derivatives or systems of PDEs, you may need to look at the corresponding (RefLinkPlain[NDSolve`StateData, paclet:tutorial/NDSolveStateData]) object to determine how the variables are ordered so that you can get the structural form of the preconditioner correctly.