Transformer to Nordsieck vectors for Adams integrators.

This class is 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>0 (j+1) (-i)j sj+1(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+1) (-i)^j terms with i being the row number starting from 1 and j being the column number starting from 1:

       [  -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.

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

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

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