/** Simple constructor. * @param n number of steps of the multistep method * (excluding the one being computed) */ private AdamsNordsieckTransformer(final int n) { final int rows = n - 1; // compute exact coefficients FieldMatrix<BigFraction> bigP = buildP(rows); FieldDecompositionSolver<BigFraction> pSolver = new FieldLUDecomposition<BigFraction>(bigP).getSolver(); BigFraction[] u = new BigFraction[rows]; Arrays.fill(u, BigFraction.ONE); BigFraction[] bigC1 = pSolver.solve(new ArrayFieldVector<BigFraction>(u, false)).toArray(); // update coefficients are computed by combining transform from // Nordsieck to multistep, then shifting rows to represent step advance // then applying inverse transform BigFraction[][] shiftedP = bigP.getData(); for (int i = shiftedP.length - 1; i > 0; --i) { // shift rows shiftedP[i] = shiftedP[i - 1]; } shiftedP[0] = new BigFraction[rows]; Arrays.fill(shiftedP[0], BigFraction.ZERO); FieldMatrix<BigFraction> bigMSupdate = pSolver.solve(new Array2DRowFieldMatrix<BigFraction>(shiftedP, false)); // convert coefficients to double update = MatrixUtils.bigFractionMatrixToRealMatrix(bigMSupdate); c1 = new double[rows]; for (int i = 0; i < rows; ++i) { c1[i] = bigC1[i].doubleValue(); } }
/** Simple constructor. * @param field field to which the time and state vector elements belong * @param n number of steps of the multistep method * (excluding the one being computed) */ private AdamsNordsieckFieldTransformer(final Field<T> field, final int n) { this.field = field; final int rows = n - 1; // compute coefficients FieldMatrix<T> bigP = buildP(rows); FieldDecompositionSolver<T> pSolver = new FieldLUDecomposition<T>(bigP).getSolver(); T[] u = MathArrays.buildArray(field, rows); Arrays.fill(u, field.getOne()); c1 = pSolver.solve(new ArrayFieldVector<T>(u, false)).toArray(); // update coefficients are computed by combining transform from // Nordsieck to multistep, then shifting rows to represent step advance // then applying inverse transform T[][] shiftedP = bigP.getData(); for (int i = shiftedP.length - 1; i > 0; --i) { // shift rows shiftedP[i] = shiftedP[i - 1]; } shiftedP[0] = MathArrays.buildArray(field, rows); Arrays.fill(shiftedP[0], field.getZero()); update = new Array2DRowFieldMatrix<T>(pSolver.solve(new Array2DRowFieldMatrix<T>(shiftedP, false)).getData()); }
/** Simple constructor. * @param nSteps number of steps of the multistep method * (excluding the one being computed) */ private AdamsNordsieckTransformer(final int nSteps) { // compute exact coefficients FieldMatrix<BigFraction> bigP = buildP(nSteps); FieldDecompositionSolver<BigFraction> pSolver = new FieldLUDecomposition<BigFraction>(bigP).getSolver(); BigFraction[] u = new BigFraction[nSteps]; Arrays.fill(u, BigFraction.ONE); BigFraction[] bigC1 = pSolver .solve(new ArrayFieldVector<BigFraction>(u, false)).toArray(); // update coefficients are computed by combining transform from // Nordsieck to multistep, then shifting rows to represent step advance // then applying inverse transform BigFraction[][] shiftedP = bigP.getData(); for (int i = shiftedP.length - 1; i > 0; --i) { // shift rows shiftedP[i] = shiftedP[i - 1]; } shiftedP[0] = new BigFraction[nSteps]; Arrays.fill(shiftedP[0], BigFraction.ZERO); FieldMatrix<BigFraction> bigMSupdate = pSolver.solve(new Array2DRowFieldMatrix<BigFraction>(shiftedP, false)); // convert coefficients to double update = MatrixUtils.bigFractionMatrixToRealMatrix(bigMSupdate); c1 = new double[nSteps]; for (int i = 0; i < nSteps; ++i) { c1[i] = bigC1[i].doubleValue(); } }
/** * Determinant of a matrix. * @param a matrix * @return determinant of a */ public C determinant(GenMatrix<C> a) { FieldMatrix<CMFieldElement<C>> am = CMFieldElementUtil.<C> toCMFieldMatrix(a); final FieldLUDecomposition<CMFieldElement<C>> lu = new FieldLUDecomposition<CMFieldElement<C>>(am); CMFieldElement<C> dm = lu.getDeterminant(); C d = dm.val; return d; }
/** * Inverse of a matrix. * @param a matrix * @return inverse matrix of a */ public GenMatrix<C> inverse(GenMatrix<C> a) { FieldMatrix<CMFieldElement<C>> am = CMFieldElementUtil.<C> toCMFieldMatrix(a); final FieldLUDecomposition<CMFieldElement<C>> lu = new FieldLUDecomposition<CMFieldElement<C>>(am); FieldDecompositionSolver<CMFieldElement<C>> fds = lu.getSolver(); FieldMatrix<CMFieldElement<C>> bm = fds.getInverse(); GenMatrix<C> g = CMFieldElementUtil.<C> matrixFromCMFieldMatrix(a.ring, bm); return g; }
/** * Solve a linear system: a x = b. * @param a matrix * @param b vector of right hand side * @return a solution vector x */ public GenVector<C> solve(GenMatrix<C> a, GenVector<C> b) { FieldMatrix<CMFieldElement<C>> am = CMFieldElementUtil.<C> toCMFieldMatrix(a); FieldVector<CMFieldElement<C>> bv = CMFieldElementUtil.<C> toCMFieldElementVector(b); final FieldLUDecomposition<CMFieldElement<C>> lu = new FieldLUDecomposition<CMFieldElement<C>>(am); FieldDecompositionSolver<CMFieldElement<C>> fds = lu.getSolver(); FieldVector<CMFieldElement<C>> xv = fds.solve(bv); GenVector<C> xa = CMFieldElementUtil.<C> vectorFromCMFieldVector(b.modul, xv); return xa; }
/** Initialize the high order scaled derivatives at step start. * @param h step size to use for scaling * @param t first steps times * @param y first steps states * @param yDot first steps derivatives * @return Nordieck vector at start of first step (h<sup>2</sup>/2 y''<sub>n</sub>, * h<sup>3</sup>/6 y'''<sub>n</sub> ... h<sup>k</sup>/k! y<sup>(k)</sup><sub>n</sub>) */ public Array2DRowFieldMatrix<T> initializeHighOrderDerivatives(final T h, final T[] t, final T[][] y, final T[][] yDot) { // using Taylor series with di = ti - t0, we get: // y(ti) - y(t0) - di y'(t0) = di^2 / h^2 s2 + ... + di^k / h^k sk + O(h^k) // y'(ti) - y'(t0) = 2 di / h^2 s2 + ... + k di^(k-1) / h^k sk + O(h^(k-1)) // we write these relations for i = 1 to i= 1+n/2 as a set of n + 2 linear // equations depending on the Nordsieck vector [s2 ... sk rk], so s2 to sk correspond // to the appropriately truncated Taylor expansion, and rk is the Taylor remainder. // The goal is to have s2 to sk as accurate as possible considering the fact the sum is // truncated and we don't want the error terms to be included in s2 ... sk, so we need // to solve also for the remainder final T[][] a = MathArrays.buildArray(field, c1.length + 1, c1.length + 1); final T[][] b = MathArrays.buildArray(field, c1.length + 1, y[0].length); final T[] y0 = y[0]; final T[] yDot0 = yDot[0]; for (int i = 1; i < y.length; ++i) { final T di = t[i].subtract(t[0]); final T ratio = di.divide(h); T dikM1Ohk = h.reciprocal(); // linear coefficients of equations // y(ti) - y(t0) - di y'(t0) and y'(ti) - y'(t0) final T[] aI = a[2 * i - 2]; final T[] aDotI = (2 * i - 1) < a.length ? a[2 * i - 1] : null; for (int j = 0; j < aI.length; ++j) { dikM1Ohk = dikM1Ohk.multiply(ratio); aI[j] = di.multiply(dikM1Ohk); if (aDotI != null) { aDotI[j] = dikM1Ohk.multiply(j + 2); } } // expected value of the previous equations final T[] yI = y[i]; final T[] yDotI = yDot[i]; final T[] bI = b[2 * i - 2]; final T[] bDotI = (2 * i - 1) < b.length ? b[2 * i - 1] : null; for (int j = 0; j < yI.length; ++j) { bI[j] = yI[j].subtract(y0[j]).subtract(di.multiply(yDot0[j])); if (bDotI != null) { bDotI[j] = yDotI[j].subtract(yDot0[j]); } } } // solve the linear system to get the best estimate of the Nordsieck vector [s2 ... sk], // with the additional terms s(k+1) and c grabbing the parts after the truncated Taylor expansion final FieldLUDecomposition<T> decomposition = new FieldLUDecomposition<T>(new Array2DRowFieldMatrix<T>(a, false)); final FieldMatrix<T> x = decomposition.getSolver().solve(new Array2DRowFieldMatrix<T>(b, false)); // extract just the Nordsieck vector [s2 ... sk] final Array2DRowFieldMatrix<T> truncatedX = new Array2DRowFieldMatrix<T>(field, x.getRowDimension() - 1, x.getColumnDimension()); for (int i = 0; i < truncatedX.getRowDimension(); ++i) { for (int j = 0; j < truncatedX.getColumnDimension(); ++j) { truncatedX.setEntry(i, j, x.getEntry(i, j)); } } return truncatedX; }
