org.apache.commons.math.ode.nonstiff

Class AdamsNordsieckTransformer

java.lang.Object

org.apache.commons.math.ode.nonstiff.AdamsNordsieckTransformer

public class AdamsNordsieckTransformer
extends java.lang.Object

Transformer to Nordsieck vectors for Adams integrators.

This class i used by Adams-Bashforth and Adams-Moulton integrators to convert between classical representation with several previous first derivatives and Nordsieck representation with higher order scaled derivatives.

We define scaled derivatives s_i(n) at step n as:

 s1(n) = h y'n for first derivative
 s2(n) = h2/2 y''n for second derivative
 s3(n) = h3/6 y'''n for third derivative
 ...
 sk(n) = hk/k! y(k)n for kth derivative
 

With the previous definition, the classical representation of multistep methods uses first derivatives only, i.e. it handles y_n, s₁(n) and q_n where q_n is defined as:

   qn = [ s1(n-1) s1(n-2) ... s1(n-(k-1)) ]T
 

(we omit the k index in the notation for clarity).

Another possible representation uses the Nordsieck vector with higher degrees scaled derivatives all taken at the same step, i.e it handles y_n, s₁(n) and r_n) where r_n is defined as:

 rn = [ s2(n), s3(n) ... sk(n) ]T
 

(here again we omit the k index in the notation for clarity)

Taylor series formulas show that for any index offset i, s₁(n-i) can be computed from s₁(n), s₂(n) ... s_k(n), the formula being exact for degree k polynomials.

 s1(n-i) = s1(n) + ∑j j (-i)j-1 sj(n)
 

The previous formula can be used with several values for i to compute the transform between classical representation and Nordsieck vector at step end. The transform between r_n and q_n resulting from the Taylor series formulas above is:

 qn = s1(n) u + P rn
 

where u is the [ 1 1 ... 1 ]^T vector and P is the (k-1)×(k-1) matrix built with the j (-i)^j-1 terms:

        [  -2   3   -4    5  ... ]
        [  -4  12  -32   80  ... ]
   P =  [  -6  27 -108  405  ... ]
        [  -8  48 -256 1280  ... ]
        [          ...           ]
 

Changing -i into +i in the formula above can be used to compute a similar transform between classical representation and Nordsieck vector at step start. The resulting matrix is simply the absolute value of matrix P.

For Adams-Bashforth method, the Nordsieck vector at step n+1 is computed from the Nordsieck vector at step n as follows:

• y_n+1 = y_n + s₁(n) + u^T r_n
• s₁(n+1) = h f(t_n+1, y_n+1)
• r_n+1 = (s₁(n) - s₁(n+1)) P^-1 u + P^-1 A P r_n

where A is a rows shifting matrix (the lower left part is an identity matrix):

        [ 0 0   ...  0 0 | 0 ]
        [ ---------------+---]
        [ 1 0   ...  0 0 | 0 ]
    A = [ 0 1   ...  0 0 | 0 ]
        [       ...      | 0 ]
        [ 0 0   ...  1 0 | 0 ]
        [ 0 0   ...  0 1 | 0 ]
 

For Adams-Moulton method, the predicted Nordsieck vector at step n+1 is computed from the Nordsieck vector at step n as follows:

• Y_n+1 = y_n + s₁(n) + u^T r_n
• S₁(n+1) = h f(t_n+1, Y_n+1)
• R_n+1 = (s₁(n) - s₁(n+1)) P^-1 u + P^-1 A P r_n

From this predicted vector, the corrected vector is computed as follows:

• y_n+1 = y_n + S₁(n+1) + [ -1 +1 -1 +1 ... ±1 ] r_n+1
• s₁(n+1) = h f(t_n+1, y_n+1)
• r_n+1 = R_n+1 + (s₁(n+1) - S₁(n+1)) P^-1 u

where the upper case Y_n+1, S₁(n+1) and R_n+1 represent the predicted states whereas the lower case y_n+1, s_n+1 and r_n+1 represent the corrected states.

We observe that both methods use similar update formulas. In both cases a P^-1u vector and a P^-1 A P matrix are used that do not depend on the state, they only depend on k. This class handles these transformations.

Since:

2.0

Method Summary

Modifier and Type	Method and Description
static AdamsNordsieckTransformer	getInstance(int nSteps) Get the Nordsieck transformer for a given number of steps.
int	getNSteps() Get the number of steps of the method (excluding the one being computed).
Array2DRowRealMatrix	initializeHighOrderDerivatives(double[] first, double[][] multistep) Initialize the high order scaled derivatives at step start.
Array2DRowRealMatrix	updateHighOrderDerivativesPhase1(Array2DRowRealMatrix highOrder) Update the high order scaled derivatives for Adams integrators (phase 1).
void	updateHighOrderDerivativesPhase2(double[] start, double[] end, Array2DRowRealMatrix highOrder) Update the high order scaled derivatives Adams integrators (phase 2).

Methods inherited from class java.lang.Object

equals, getClass, hashCode, notify, notifyAll, toString, wait, wait, wait

Method Detail

getInstance

public static AdamsNordsieckTransformer getInstance(int nSteps)

Get the Nordsieck transformer for a given number of steps.

Parameters:

nSteps - number of steps of the multistep method (excluding the one being computed)

Returns:

Nordsieck transformer for the specified number of steps

getNSteps

public int getNSteps()

Get the number of steps of the method (excluding the one being computed).

Returns:

number of steps of the method (excluding the one being computed)

initializeHighOrderDerivatives

public Array2DRowRealMatrix initializeHighOrderDerivatives(double[] first,
                                                           double[][] multistep)

Initialize the high order scaled derivatives at step start.

Parameters:

first - first scaled derivative at step start

multistep - scaled derivatives after step start (hy'1, ..., hy'k-1) will be modified

Returns:

high order derivatives at step start

updateHighOrderDerivativesPhase1

public Array2DRowRealMatrix updateHighOrderDerivativesPhase1(Array2DRowRealMatrix highOrder)

Update the high order scaled derivatives for Adams integrators (phase 1).

The complete update of high order derivatives has a form similar to:

 rn+1 = (s1(n) - s1(n+1)) P-1 u + P-1 A P rn
 

this method computes the P^-1 A P r_n part.

Parameters:

highOrder - high order scaled derivatives (h²/2 y'', ... h^k/k! y(k))

Returns:

updated high order derivatives

See Also:

updateHighOrderDerivativesPhase2(double[], double[], Array2DRowRealMatrix)

updateHighOrderDerivativesPhase2

public void updateHighOrderDerivativesPhase2(double[] start,
                                             double[] end,
                                             Array2DRowRealMatrix highOrder)

Update the high order scaled derivatives Adams integrators (phase 2).

The complete update of high order derivatives has a form similar to:

 rn+1 = (s1(n) - s1(n+1)) P-1 u + P-1 A P rn
 

this method computes the (s₁(n) - s₁(n+1)) P^-1 u part.

Phase 1 of the update must already have been performed.

Parameters:

start - first order scaled derivatives at step start

end - first order scaled derivatives at step end

highOrder - high order scaled derivatives, will be modified (h²/2 y'', ... h^k/k! y(k))

See Also:

updateHighOrderDerivativesPhase1(Array2DRowRealMatrix)