Class AdamsNordsieckFieldTransformer<T extends CalculusFieldElement<T>>
- Type Parameters:
T
- the type of the field elements
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 si(n) at step n as: \[ \left\{\begin{align} s_1(n) &= h y'_n \text{ for first derivative}\\ s_2(n) &= \frac{h^2}{2} y_n'' \text{ for second derivative}\\ s_3(n) &= \frac{h^3}{6} y_n''' \text{ for third derivative}\\ &\cdots\\ s_k(n) &= \frac{h^k}{k!} y_n^{(k)} \text{ for } k^\mathrm{th} \text{ derivative} \end{align}\right. \]
With the previous definition, the classical representation of multistep methods uses first derivatives only, i.e. it handles yn, s1(n) and qn where qn is defined as: \[ q_n = [ s_1(n-1) s_1(n-2) \ldots s_1(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 yn, s1(n) and rn) where rn is defined as: \[ r_n = [ s_2(n), s_3(n) \ldots s_k(n) ]^T \] (here again we omit the k index in the notation for clarity)
Taylor series formulas show that for any index offset i, s1(n-i) can be computed from s1(n), s2(n) ... sk(n), the formula being exact for degree k polynomials. \[ s_1(n-i) = s_1(n) + \sum_{j\gt 0} (j+1) (-i)^j s_{j+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 rn and qn resulting from the Taylor series formulas above is: \[ q_n = s_1(n) u + P r_n \] 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: \[ P=\begin{bmatrix} -2 & 3 & -4 & 5 & \ldots \\ -4 & 12 & -32 & 80 & \ldots \\ -6 & 27 & -108 & 405 & \ldots \\ -8 & 48 & -256 & 1280 & \ldots \\ & & \ldots\\ \end{bmatrix} \]
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:
- yn+1 = yn + s1(n) + uT rn
- s1(n+1) = h f(tn+1, yn+1)
- rn+1 = (s1(n) - s1(n+1)) P-1 u + P-1 A P rn
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:
- Yn+1 = yn + s1(n) + uT rn
- S1(n+1) = h f(tn+1, Yn+1)
- Rn+1 = (s1(n) - s1(n+1)) P-1 u + P-1 A P rn
- yn+1 = yn + S1(n+1) + [ -1 +1 -1 +1 ... ±1 ] rn+1
- s1(n+1) = h f(tn+1, yn+1)
- rn+1 = Rn+1 + (s1(n+1) - S1(n+1)) P-1 u
where the upper case Yn+1, S1(n+1) and Rn+1 represent the predicted states whereas the lower case yn+1, sn+1 and rn+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.
-
Method Summary
Modifier and TypeMethodDescriptionstatic <T extends CalculusFieldElement<T>>
AdamsNordsieckFieldTransformer<T>getInstance
(Field<T> field, int nSteps) Get the Nordsieck transformer for a given field and number of steps.initializeHighOrderDerivatives
(T h, T[] t, T[][] y, T[][] yDot) Initialize the high order scaled derivatives at step start.updateHighOrderDerivativesPhase1
(Array2DRowFieldMatrix<T> highOrder) Update the high order scaled derivatives for Adams integrators (phase 1).void
updateHighOrderDerivativesPhase2
(T[] start, T[] end, Array2DRowFieldMatrix<T> highOrder) Update the high order scaled derivatives Adams integrators (phase 2).
-
Method Details
-
getInstance
public static <T extends CalculusFieldElement<T>> AdamsNordsieckFieldTransformer<T> getInstance(Field<T> field, int nSteps) Get the Nordsieck transformer for a given field and number of steps.- Type Parameters:
T
- the type of the field elements- Parameters:
field
- field to which the time and state vector elements belongnSteps
- number of steps of the multistep method (excluding the one being computed)- Returns:
- Nordsieck transformer for the specified field and number of steps
-
initializeHighOrderDerivatives
Initialize the high order scaled derivatives at step start.- Parameters:
h
- step size to use for scalingt
- first steps timesy
- first steps statesyDot
- first steps derivatives- Returns:
- Nordieck vector at start of first step (h2/2 y''n, h3/6 y'''n ... hk/k! y(k)n)
-
updateHighOrderDerivativesPhase1
public Array2DRowFieldMatrix<T> updateHighOrderDerivativesPhase1(Array2DRowFieldMatrix<T> highOrder) Update the high order scaled derivatives for Adams integrators (phase 1).The complete update of high order derivatives has a form similar to: \[ r_{n+1} = (s_1(n) - s_1(n+1)) P^{-1} u + P^{-1} A P r_n \] this method computes the P-1 A P rn part.
- Parameters:
highOrder
- high order scaled derivatives (h2/2 y'', ... hk/k! y(k))- Returns:
- updated high order derivatives
- See Also:
-
updateHighOrderDerivativesPhase2
public void updateHighOrderDerivativesPhase2(T[] start, T[] end, Array2DRowFieldMatrix<T> highOrder) Update the high order scaled derivatives Adams integrators (phase 2).The complete update of high order derivatives has a form similar to: \[ r_{n+1} = (s_1(n) - s_1(n+1)) P^{-1} u + P^{-1} A P r_n \] this method computes the (s1(n) - s1(n+1)) P-1 u part.
Phase 1 of the update must already have been performed.
- Parameters:
start
- first order scaled derivatives at step startend
- first order scaled derivatives at step endhighOrder
- high order scaled derivatives, will be modified (h2/2 y'', ... hk/k! y(k))- See Also:
-