Norms in NDSolve

(RefLink[NDSolve, paclet:ref/NDSolve]) uses norms of error estimates to determine when solutions satisfy error tolerances. In nearly all cases the norm has been weighted, or scaled, such that it is less than 1 if error tolerances have been satisfied and greater than 1 if error tolerances are not satisfied. One significant advantage of such a scaled norm is that a given method can be written without explicit reference to tolerances: the satisfaction of tolerances is found by comparing the scaled norm to 1, thus simplifying the code required for checking error estimates within methods.

Suppose that v is vector and u is a reference vector to compute weights with (typically u is an approximate solution vector). Then the scaled vector w to which the norm is applied has components:

(Subscript[w, i])==((|(Subscript[v, i])|)/((Subscript[t, a])+(Subscript[t, r])| (Subscript[u, i])|))

where absolute and relative tolerances (Subscript[t, a]) and (Subscript[t, r]) are derived respectively from the (RefLink[AccuracyGoal, paclet:ref/AccuracyGoal])->ag and (RefLink[PrecisionGoal, paclet:ref/PrecisionGoal])->pg options by (Subscript[t, a])==(10^(-ag)) and (Subscript[t, r])==(10^(-pg)).

The actual norm used is determined by the setting for the (RefLink[NormFunction, paclet:ref/NormFunction]) option given to (RefLink[NDSolve, paclet:ref/NDSolve]).

(RefLink[NormFunction, paclet:ref/NormFunction]) option to (RefLink[NDSolve, paclet:ref/NDSolve]).

The setting for the (RefLink[NormFunction, paclet:ref/NormFunction]) option can be any function that returns a scalar for a vector argument and satisfies the properties of a norm. If you specify a function that does not satisfy the required properties of a norm, (RefLink[NDSolve, paclet:ref/NDSolve]) will almost surely run into problems and give an answer, if any, which is incorrect.

The default value of (RefLink[Automatic, paclet:ref/Automatic]) means that (RefLink[NDSolve, paclet:ref/NDSolve]) may use different norms for different methods. Most methods use an infinity-norm, but the IDA method for DAEs uses a 2-norm because that helps maintain smoothness in the merit function for finding roots of the residual. It is strongly recommended that you use (RefLink[Norm, paclet:ref/Norm]) with a particular value of p. For this reason, you can also use the shorthand (RefLink[NormFunction, paclet:ref/NormFunction])->p in place of (RefLink[NormFunction, paclet:ref/NormFunction])->((RefLink[Norm, paclet:ref/Norm])[#,p]/(RefLink[Length, paclet:ref/Length])[#]^(1/p)&). The most commonly used implementations for p==1, p==2, and p==∞ have been specially optimized for speed.

This compares the overall error for computing the solution to the simple harmonic oscillator over 100 cycles with different norms specified:

Map[First[(1 - x[100 π]) /. NDSolve[{x''[t] + x[t] == 0, x[0] == 1, x'[0] == 0}, x, {t,0,100 π}, Method->ExplicitRungeKutta, NormFunction->#]]&, {1, 2, ∞}]

The reason that error decreases with increasing p is because the norms are normalized by multiplying with 1/(n^(1/p)), where n is the length of the vector. This is often important in (RefLink[NDSolve, paclet:ref/NDSolve]) because in many cases, an attempt is being made to check the approximation to a function, where more points should give a better approximation, or less error.

Consider a finite difference approximation to the first derivative of a periodic function u given by (Subscript[u, i]^('))==(((Subscript[u, i+1])-(Subscript[u, i]))/( h)) where (Subscript[u, i])==u((Subscript[x, i])) on a grid with uniform spacing h==(Subscript[x, i+1])-(Subscript[x, i]). In the Wolfram Language, this can easily be computed using (RefLink[ListCorrelate, paclet:ref/ListCorrelate]).

This computes the error of the first derivative approximation for the cosine function on a grid with 16 points covering the interval [0,2π]:

h = 2π/16.;
grid = h Range[16];
err16 = Sin[grid] - ListCorrelate[{1,-1}/h, Cos[grid], {1,1}]

This computes the error of the first derivative approximation for the cosine function on a grid with 32 points covering the interval [0,2π]:

h = 2π/32.;
grid = h Range[32];
err32 = Sin[grid] - ListCorrelate[{1,-1}/h, Cos[grid], {1,1}]

It is quite apparent that the pointwise error is significantly less with a larger number of points.

The 2 norms of the vectors are of the same order of magnitude:

{Norm[err16,2], Norm[err32,2]}

The norms of the vectors are comparable because the number of components in the vector has increased, so the usual linear algebra norm does not properly reflect the convergence. Normalizing by multiplying by 1/(n^(1/p)) reflects the convergence in the function space properly.

The normalized 2 norms of the vectors reflect the convergence to the actual function. Since the approximation is first order, doubling the number of grid points should approximately halve the error:

{Norm[err16,2]/Sqrt[16], Norm[err32,2]/Sqrt[32]}

Note that if you specify a function an option value, and you intend to use it for PDE or function approximation solutions, you should be sure to include a proper normalization in the function.

ScaledVectorNorm

Methods that have error control need to determine whether a step satisfies local error tolerances or not. To simplify the process of checking this, utility function ScaledVectorNorm does the scaling (1) and computes the norm. The table includes the formulas for specific values of p for reference.

ScaledVectorNorm.

This sets up a scaled vector norm object with the default machine-precision tolerances used in (RefLink[NDSolve, paclet:ref/NDSolve]):

svn = NDSolve`ScaledVectorNorm[2, {(10.^(-8)), (10.^(-8))}]

This applies the scaled norm object with a sample error and solution reference vector:

svn[{9. (10.^(-9)), (10.^(-8))}, {2.,1.}]

Because of the absolute tolerance term, the value comes out reasonably even if some of the components of the reference solution are zero:

svn[{9. (10.^(-9)), (10.^(-8)), 2 (10^(-8))}, {1.,0., 0.}]

When setting up a method for (RefLink[NDSolve, paclet:ref/NDSolve]), you can get the appropriate ScaledVectorNorm object to use using the "Norm" method function of the (RefLinkPlain[NDSolve`StateData, paclet:tutorial/NDSolveStateData]) object.

Here is an NDSolve`StateData object:

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

This gets the appropriate scaled norm to use from the state data:

svn = state["Norm"]

This applies it to a sample error vector using the initial condition as reference vector:

svn[{(10.^(-9)), (10.^(-8))}, state@"SolutionVector"["Forward"]]