From b6703d294ca04c5fee554a580693d1bea77893ca Mon Sep 17 00:00:00 2001 From: Luc Maisonobe Date: Thu, 9 Apr 2009 15:23:11 +0000 Subject: [PATCH] added a transformer between Nordsieck form and multistep form for state history in multistep ODE integrators (this will help implementing adaptive stepsize multistep integrators) git-svn-id: https://svn.apache.org/repos/asf/commons/proper/math/trunk@763692 13f79535-47bb-0310-9956-ffa450edef68 --- .../math/ode/NordsieckTransformer.java | 511 ++++++++++++++++++ .../math/ode/NordsieckTransformerTest.java | 268 +++++++++ 2 files changed, 779 insertions(+) create mode 100644 src/java/org/apache/commons/math/ode/NordsieckTransformer.java create mode 100644 src/test/org/apache/commons/math/ode/NordsieckTransformerTest.java diff --git a/src/java/org/apache/commons/math/ode/NordsieckTransformer.java b/src/java/org/apache/commons/math/ode/NordsieckTransformer.java new file mode 100644 index 000000000..450d0a104 --- /dev/null +++ b/src/java/org/apache/commons/math/ode/NordsieckTransformer.java @@ -0,0 +1,511 @@ +/* + * Licensed to the Apache Software Foundation (ASF) under one or more + * contributor license agreements. See the NOTICE file distributed with + * this work for additional information regarding copyright ownership. + * The ASF licenses this file to You under the Apache License, Version 2.0 + * (the "License"); you may not use this file except in compliance with + * the License. You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ + +package org.apache.commons.math.ode; + +import java.io.Serializable; +import java.math.BigInteger; +import java.util.Arrays; + +import org.apache.commons.math.fraction.BigFraction; +import org.apache.commons.math.linear.RealMatrix; +import org.apache.commons.math.linear.RealMatrixImpl; + +/** + * This class transforms state history between multistep (with or without + * derivatives) and Nordsieck forms. + *

+ * {@link MultistepIntegrator multistep integrators} use state history + * from several previous steps to compute the current state. They may also use + * the first derivative of current state. All states are separated by a fixed + * step size h from each other. Since these methods are based on polynomial + * interpolation, the information from the previous state may be represented + * in another equivalent way: using the state higher order derivatives at + * current step rather. This class transforms state history between these three + * equivalent forms. + *

+ *

+ * The supported forms for a dimension n history are: + *

+ * In these expressions, yk is the state at the current step. For each p, + * yk-p is the state at the pth previous step. y'k, + * y''k ... yn-1k are respectively the first, second, ... + * (n-1)th derivatives of the state at current step and h is the fixed + * step size. + *

+ *

+ * The transforms are exact for polynomials. + *

+ *

+ * In Nordsieck form, the state history can be converted from step size h to step + * size h' by rescaling each component by 1, h'/h, (h'/h)2 ... + * (h'/h)n-1. + *

+ *

+ * Instances of this class are guaranteed to be immutable. + *

+ * @see org.apache.commons.math.ode.MultistepIntegrator + * @see org.apache.commons.math.ode.nonstiff.AdamsBashforthIntegrator + * @see org.apache.commons.math.ode.nonstiff.AdamsMoultonIntegrator + * @version $Revision$ $Date$ + * @since 2.0 + */ +public class NordsieckTransformer implements Serializable { + + /** Serializable version identifier. */ + private static final long serialVersionUID = -2707468304560314664L; + + /** Nordsieck to Multistep without derivatives matrix. */ + private final RealMatrix matNtoMWD; + + /** Multistep without derivatives to Nordsieck matrix. */ + private final RealMatrix matMWDtoN; + + /** Nordsieck to Multistep matrix. */ + private final RealMatrix matNtoM; + + /** Multistep to Nordsieck matrix. */ + private final RealMatrix matMtoN; + + /** + * Build a transformer for a specified order. + * @param n dimension of the history + */ + public NordsieckTransformer(final int n) { + + // from Nordsieck to multistep without derivatives + final BigInteger[][] bigNtoMWD = buildNordsieckToMultistepWithoutDerivatives(n); + double[][] dataNtoMWD = new double[n][n]; + for (int i = 0; i < n; ++i) { + double[] dRow = dataNtoMWD[i]; + BigInteger[] bRow = bigNtoMWD[i]; + for (int j = 0; j < n; ++j) { + dRow[j] = bRow[j].doubleValue(); + } + } + matNtoMWD = new RealMatrixImpl(dataNtoMWD, false); + + // from multistep without derivatives to Nordsieck + final BigFraction[][] bigToN = buildMultistepWithoutDerivativesToNordsieck(n); + double[][] dataMWDtoN = new double[n][n]; + for (int i = 0; i < n; ++i) { + double[] dRow = dataMWDtoN[i]; + BigFraction[] bRow = bigToN[i]; + for (int j = 0; j < n; ++j) { + dRow[j] = bRow[j].doubleValue(); + } + } + matMWDtoN = new RealMatrixImpl(dataMWDtoN, false); + + // from Nordsieck to multistep + final BigInteger[][] bigNtoM = buildNordsieckToMultistep(n); + double[][] dataNtoM = new double[n][n]; + for (int i = 0; i < n; ++i) { + double[] dRow = dataNtoM[i]; + BigInteger[] bRow = bigNtoM[i]; + for (int j = 0; j < n; ++j) { + dRow[j] = bRow[j].doubleValue(); + } + } + matNtoM = new RealMatrixImpl(dataNtoM, false); + + // from multistep to Nordsieck + convertMWDtNtoMtN(bigToN); + double[][] dataMtoN = new double[n][n]; + for (int i = 0; i < n; ++i) { + double[] dRow = dataMtoN[i]; + BigFraction[] bRow = bigToN[i]; + for (int j = 0; j < n; ++j) { + dRow[j] = bRow[j].doubleValue(); + } + } + matMtoN = new RealMatrixImpl(dataMtoN, false); + + } + + /** + * Build the transform from Nordsieck to multistep without derivatives. + * @param n dimension of the history + * @return transform from Nordsieck to multistep without derivatives + */ + public static BigInteger[][] buildNordsieckToMultistepWithoutDerivatives(final int n) { + + final BigInteger[][] array = new BigInteger[n][n]; + + // row 0: [1 0 0 0 ... 0 ] + array[0][0] = BigInteger.ONE; + Arrays.fill(array[0], 1, n, BigInteger.ZERO); + + // the following expressions are direct applications of Taylor series + // rows 1 to n-1: aij = (-i)^j + // [ 1 -1 1 -1 1 ...] + // [ 1 -2 4 -8 16 ...] + // [ 1 -3 9 -27 81 ...] + // [ 1 -4 16 -64 256 ...] + for (int i = 1; i < n; ++i) { + final BigInteger[] row = array[i]; + final BigInteger factor = BigInteger.valueOf(-i); + BigInteger aj = BigInteger.ONE; + for (int j = 0; j < n; ++j) { + row[j] = aj; + aj = aj.multiply(factor); + } + } + + return array; + + } + + /** + * Build the transform from multistep without derivatives to Nordsieck. + * @param n dimension of the history + * @return transform from multistep without derivatives to Nordsieck + */ + public static BigFraction[][] buildMultistepWithoutDerivativesToNordsieck(final int n) { + + final BigInteger[][] iArray = new BigInteger[n][n]; + + // row 0: [1 0 0 0 ... 0 ] + iArray[0][0] = BigInteger.ONE; + Arrays.fill(iArray[0], 1, n, BigInteger.ZERO); + + // We use recursive definitions of triangular integer series for each column. + // For example column 0 of matrices of increasing dimensions are: + // 1/0! for dimension 1 + // 1/1!, 1/1! for dimension 2 + // 2/2!, 3/2!, 1/2! for dimension 3 + // 6/3!, 11/3!, 6/3!, 1/3! for dimension 4 + // 24/4!, 50/4!, 35/4!, 10/4!, 1/4! for dimension 5 + // The numerators are the Stirling numbers of the first kind, (A008275 in + // Sloane's encyclopedia http://www.research.att.com/~njas/sequences/A008275) + // with a multiplicative factor of +/-1 (which we will write +/-binomial(n-1, 0)). + // In the same way, column 1 is A049444 with a multiplicative factor of + // +/-binomial(n-1, 1) and column 2 is A123319 with a multiplicative factor of + // +/-binomial(n-1, 2). The next columns are defined by similar definitions but + // are not identified in Sloane's encyclopedia. + // Another interesting observation is that for each dimension k, the last column + // (except the initial 0) is a copy of the first column of the dimension k-1 matrix, + // possibly with an opposite sign (i.e. these columns are also linked to Stirling + // numbers of the first kind). + for (int i = 1; i < n; ++i) { + + final BigInteger bigI = BigInteger.valueOf(i); + + // row i + BigInteger[] rowK = iArray[i]; + BigInteger[] rowKm1 = iArray[i - 1]; + for (int j = 0; j < i; ++j) { + rowK[j] = BigInteger.ONE; + } + rowK[i] = rowKm1[0]; + + // rows i-1 to 1 + for (int k = i - 1; k > 0; --k) { + + // select rows + rowK = rowKm1; + rowKm1 = iArray[k - 1]; + + // apply recursive defining formula + for (int j = 0; j < i; ++j) { + rowK[j] = rowK[j].multiply(bigI).add(rowKm1[j]); + } + + // initialize new last column + rowK[i] = rowKm1[0]; + + } + rowKm1[0] = rowKm1[0].multiply(bigI); + + } + + // apply column specific factors + final BigInteger factorial = iArray[0][0]; + final BigFraction[][] fArray = new BigFraction[n][n]; + for (int i = 0; i < n; ++i) { + final BigFraction[] fRow = fArray[i]; + final BigInteger[] iRow = iArray[i]; + BigInteger binomial = BigInteger.ONE; + for (int j = 0; j < n; ++j) { + fRow[j] = new BigFraction(binomial.multiply(iRow[j]), factorial); + binomial = binomial.negate().multiply(BigInteger.valueOf(n - j - 1)).divide(BigInteger.valueOf(j + 1)); + } + } + + return fArray; + + } + + /** + * Build the transform from Nordsieck to multistep. + * @param n dimension of the history + * @return transform from Nordsieck to multistep + */ + public static BigInteger[][] buildNordsieckToMultistep(final int n) { + + final BigInteger[][] array = new BigInteger[n][n]; + + // row 0: [1 0 0 0 ... 0 ] + array[0][0] = BigInteger.ONE; + Arrays.fill(array[0], 1, n, BigInteger.ZERO); + + if (n > 1) { + + // row 1: [0 1 0 0 ... 0 ] + array[1][0] = BigInteger.ZERO; + array[1][1] = BigInteger.ONE; + Arrays.fill(array[1], 2, n, BigInteger.ZERO); + + // the following expressions are direct applications of Taylor series + // rows 2 to n-1: aij = (1-i)^j + // [ 1 -1 1 -1 1 ...] + // [ 1 -2 4 -8 16 ...] + // [ 1 -3 9 -27 81 ...] + // [ 1 -4 16 -64 256 ...] + for (int i = 2; i < n; ++i) { + final BigInteger[] row = array[i]; + final BigInteger factor = BigInteger.valueOf(1 - i); + BigInteger aj = BigInteger.ONE; + for (int j = 0; j < n; ++j) { + row[j] = aj; + aj = aj.multiply(factor); + } + } + + } + + return array; + + } + + /** + * Build the transform from multistep to Nordsieck. + * @param n dimension of the history + * @return transform from multistep to Nordsieck + */ + public static BigFraction[][] buildMultistepToNordsieck(final int n) { + final BigFraction[][] array = buildMultistepWithoutDerivativesToNordsieck(n); + convertMWDtNtoMtN(array); + return array; + } + + /** + * Convert a transform from multistep without derivatives to Nordsieck to + * multistep to Nordsieck. + * @param work array, contains tansform from multistep without derivatives + * to Nordsieck on input, will be overwritten with tansform from multistep + * to Nordsieck on output + */ + private static void convertMWDtNtoMtN(BigFraction[][] array) { + + final int n = array.length; + if (n == 1) { + return; + } + + // the second row of the matrix without derivatives represents the linear equation: + // hy' = a0 yk + a1 yk-1 + ... + a(n-1) yk-(n-1) + // we solve it with respect to the oldest state yk-(n-1) and get + // yk-(n-1) = -a0/a(n-1) yk + 1/a(n-1) hy' - a1/a(n-1) yk-1 - ... + final BigFraction[] secondRow = array[1]; + final BigFraction[] solved = new BigFraction[n]; + final BigFraction f = secondRow[n - 1].reciprocal().negate(); + solved[0] = secondRow[0].multiply(f); + solved[1] = f.negate(); + for (int j = 2; j < n; ++j) { + solved[j] = secondRow[j - 1].multiply(f); + } + + // update the matrix so it expects hy' in second element + // rather than yk-(n-1) in last elements when post-multiplied + for (int i = 0; i < n; ++i) { + final BigFraction[] rowI = array[i]; + final BigFraction last = rowI[n - 1]; + for (int j = n - 1; j > 1; --j) { + rowI[j] = rowI[j - 1].add(last.multiply(solved[j])); + } + rowI[1] = last.multiply(solved[1]); + rowI[0] = rowI[0].add(last.multiply(solved[0])); + } + + } + + /** + * Transform a scalar state history from multistep form to Nordsieck form. + *

+ * The input state history must be in multistep form with element 0 for + * current state, element 1 for current state scaled first derivative, element + * 2 for previous state ... element n-1 for (n-2)th previous state. + * The output state history will be in Nordsieck form with element 0 for + * current state, element 1 for current state scaled first derivative, element + * 2 for current state scaled second derivative ... element n-1 for current state + * scaled (n-1)th derivative. + *

+ * @param multistepHistory scalar state history in multistep form + * @return scalar state history in Nordsieck form + */ + public double[] multistepToNordsieck(final double[] multistepHistory) { + return matMtoN.operate(multistepHistory); + } + + /** + * Transform a vectorial state history from multistep form to Nordsieck form. + *

+ * The input state history must be in multistep form with row 0 for + * current state, row 1 for current state scaled first derivative, row + * 2 for previous state ... row n-1 for (n-2)th previous state. + * The output state history will be in Nordsieck form with row 0 for + * current state, row 1 for current state scaled first derivative, row + * 2 for current state scaled second derivative ... row n-1 for current state + * scaled (n-1)th derivative. + *

+ * @param multistepHistory vectorial state history in multistep form + * @return vectorial state history in Nordsieck form + */ + public RealMatrix multistepToNordsieck(final RealMatrix multistepHistory) { + return matMtoN.multiply(multistepHistory); + } + + /** + * Transform a scalar state history from Nordsieck form to multistep form. + *

+ * The input state history must be in Nordsieck form with element 0 for + * current state, element 1 for current state scaled first derivative, element + * 2 for current state scaled second derivative ... element n-1 for current state + * scaled (n-1)th derivative. + * The output state history will be in multistep form with element 0 for + * current state, element 1 for current state scaled first derivative, element + * 2 for previous state ... element n-1 for (n-2)th previous state. + *

+ * @param nordsieckHistory scalar state history in Nordsieck form + * @return scalar state history in multistep form + */ + public double[] nordsieckToMultistep(final double[] nordsieckHistory) { + return matNtoM.operate(nordsieckHistory); + } + + /** + * Transform a vectorial state history from Nordsieck form to multistep form. + *

+ * The input state history must be in Nordsieck form with row 0 for + * current state, row 1 for current state scaled first derivative, row + * 2 for current state scaled second derivative ... row n-1 for current state + * scaled (n-1)th derivative. + * The output state history will be in multistep form with row 0 for + * current state, row 1 for current state scaled first derivative, row + * 2 for previous state ... row n-1 for (n-2)th previous state. + *

+ * @param nordsieckHistory vectorial state history in Nordsieck form + * @return vectorial state history in multistep form + */ + public RealMatrix nordsieckToMultistep(final RealMatrix nordsieckHistory) { + return matNtoM.multiply(nordsieckHistory); + } + + /** + * Transform a scalar state history from multistep without derivatives form + * to Nordsieck form. + *

+ * The input state history must be in multistep without derivatives form with + * element 0 for current state, element 1 for previous state ... element n-1 + * for (n-1)th previous state. + * The output state history will be in Nordsieck form with element 0 for + * current state, element 1 for current state scaled first derivative, element + * 2 for current state scaled second derivative ... element n-1 for current state + * scaled (n-1)th derivative. + *

+ * @param mwdHistory scalar state history in multistep without derivatives form + * @return scalar state history in Nordsieck form + */ + public double[] multistepWithoutDerivativesToNordsieck(final double[] mwdHistory) { + return matMWDtoN.operate(mwdHistory); + } + + /** + * Transform a vectorial state history from multistep without derivatives form + * to Nordsieck form. + *

+ * The input state history must be in multistep without derivatives form with + * row 0 for current state, row 1 for previous state ... row n-1 + * for (n-1)th previous state. + * The output state history will be in Nordsieck form with row 0 for + * current state, row 1 for current state scaled first derivative, row + * 2 for current state scaled second derivative ... row n-1 for current state + * scaled (n-1)th derivative. + *

+ * @param mwdHistory vectorial state history in multistep without derivatives form + * @return vectorial state history in Nordsieck form + */ + public RealMatrix multistepWithoutDerivativesToNordsieck(final RealMatrix mwdHistory) { + return matMWDtoN.multiply(mwdHistory); + } + + /** + * Transform a scalar state history from Nordsieck form to multistep without + * derivatives form. + *

+ * The input state history must be in Nordsieck form with element 0 for + * current state, element 1 for current state scaled first derivative, element + * 2 for current state scaled second derivative ... element n-1 for current state + * scaled (n-1)th derivative. + * The output state history will be in multistep without derivatives form with + * element 0 for current state, element 1 for previous state ... element n-1 + * for (n-1)th previous state. + *

+ * @param nordsieckHistory scalar state history in Nordsieck form + * @return scalar state history in multistep without derivatives form + */ + public double[] nordsieckToMultistepWithoutDerivatives(final double[] nordsieckHistory) { + return matNtoMWD.operate(nordsieckHistory); + } + + /** + * Transform a vectorial state history from Nordsieck form to multistep without + * derivatives form. + *

+ * The input state history must be in Nordsieck form with row 0 for + * current state, row 1 for current state scaled first derivative, row + * 2 for current state scaled second derivative ... row n-1 for current state + * scaled (n-1)th derivative. + * The output state history will be in multistep without derivatives form with + * row 0 for current state, row 1 for previous state ... row n-1 + * for (n-1)th previous state. + *

+ * @param nordsieckHistory vectorial state history in Nordsieck form + * @return vectorial state history in multistep without derivatives form + */ + public RealMatrix nordsieckToMultistepWithoutDerivatives(final RealMatrix nordsieckHistory) { + return matNtoMWD.multiply(nordsieckHistory); + } + +} diff --git a/src/test/org/apache/commons/math/ode/NordsieckTransformerTest.java b/src/test/org/apache/commons/math/ode/NordsieckTransformerTest.java new file mode 100644 index 000000000..e261f2547 --- /dev/null +++ b/src/test/org/apache/commons/math/ode/NordsieckTransformerTest.java @@ -0,0 +1,268 @@ +/* + * Licensed to the Apache Software Foundation (ASF) under one or more + * contributor license agreements. See the NOTICE file distributed with + * this work for additional information regarding copyright ownership. + * The ASF licenses this file to You under the Apache License, Version 2.0 + * (the "License"); you may not use this file except in compliance with + * the License. You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ + +package org.apache.commons.math.ode; + +import java.math.BigInteger; +import java.util.Random; + +import junit.framework.Test; +import junit.framework.TestCase; +import junit.framework.TestSuite; + +import org.apache.commons.math.analysis.polynomials.PolynomialFunction; +import org.apache.commons.math.fraction.BigFraction; +import org.apache.commons.math.linear.RealMatrix; +import org.apache.commons.math.linear.RealMatrixImpl; + +public class NordsieckTransformerTest +extends TestCase { + + public NordsieckTransformerTest(String name) { + super(name); + } + + public void testDimension2() { + NordsieckTransformer transformer = new NordsieckTransformer(2); + double[] nordsieckHistory = new double[] { 1.0, 2.0 }; + double[] mwdHistory = new double[] { 1.0, -1.0 }; + double[] multistepHistory = new double[] { 1.0, 2.0 }; + checkVector(nordsieckHistory, transformer.multistepWithoutDerivativesToNordsieck(mwdHistory)); + checkVector(mwdHistory, transformer.nordsieckToMultistepWithoutDerivatives(nordsieckHistory)); + checkVector(nordsieckHistory, transformer.multistepToNordsieck(multistepHistory)); + checkVector(multistepHistory, transformer.nordsieckToMultistep(nordsieckHistory)); + } + + public void testDimension3() { + NordsieckTransformer transformer = new NordsieckTransformer(3); + double[] nordsieckHistory = new double[] { 1.0, 4.0, 18.0 }; + double[] mwdHistory = new double[] { 1.0, 15.0, 65.0 }; + double[] multistepHistory = new double[] { 1.0, 4.0, 15.0 }; + checkVector(nordsieckHistory, transformer.multistepWithoutDerivativesToNordsieck(mwdHistory)); + checkVector(mwdHistory, transformer.nordsieckToMultistepWithoutDerivatives(nordsieckHistory)); + checkVector(nordsieckHistory, transformer.multistepToNordsieck(multistepHistory)); + checkVector(multistepHistory, transformer.nordsieckToMultistep(nordsieckHistory)); + } + + public void testDimension7() { + NordsieckTransformer transformer = new NordsieckTransformer(7); + RealMatrix nordsieckHistory = + new RealMatrixImpl(new double[][] { + { 1, 2, 3 }, + { -2, 1, 0 }, + { 1, 1, 1 }, + { 0, -1, 1 }, + { 1, -1, 2 }, + { 2, 0, 1 }, + { 1, 1, 2 } + }, false); + RealMatrix mwdHistory = + new RealMatrixImpl(new double[][] { + { 1, 2, 3 }, + { 4, 3, 6 }, + { 25, 60, 127 }, + { 340, 683, 1362 }, + { 2329, 3918, 7635 }, + { 10036, 15147, 29278 }, + { 32449, 45608, 87951 } + }, false); + RealMatrix multistepHistory = + new RealMatrixImpl(new double[][] { + { 1, 2, 3 }, + { -2, 1, 0 }, + { 4, 3, 6 }, + { 25, 60, 127 }, + { 340, 683, 1362 }, + { 2329, 3918, 7635 }, + { 10036, 15147, 29278 } + }, false); + + RealMatrix m = transformer.multistepWithoutDerivativesToNordsieck(mwdHistory); + assertEquals(0.0, m.subtract(nordsieckHistory).getNorm(), 1.0e-11); + m = transformer.nordsieckToMultistepWithoutDerivatives(nordsieckHistory); + assertEquals(0.0, m.subtract(mwdHistory).getNorm(), 1.0e-11); + m = transformer.multistepToNordsieck(multistepHistory); + assertEquals(0.0, m.subtract(nordsieckHistory).getNorm(), 1.0e-11); + m = transformer.nordsieckToMultistep(nordsieckHistory); + assertEquals(0.0, m.subtract(multistepHistory).getNorm(), 1.0e-11); + + } + + public void testInverseWithoutDerivatives() { + for (int n = 1; n < 20; ++n) { + BigInteger[][] nTom = + NordsieckTransformer.buildNordsieckToMultistepWithoutDerivatives(n); + BigFraction[][] mTon = + NordsieckTransformer.buildMultistepWithoutDerivativesToNordsieck(n); + for (int i = 0; i < n; ++i) { + for (int j = 0; j < n; ++j) { + BigFraction s = BigFraction.ZERO; + for (int k = 0; k < n; ++k) { + s = s.add(mTon[i][k].multiply(nTom[k][j])); + } + assertEquals((i == j) ? BigFraction.ONE : BigFraction.ZERO, s); + } + } + } + } + + public void testInverse() { + for (int n = 1; n < 20; ++n) { + BigInteger[][] nTom = + NordsieckTransformer.buildNordsieckToMultistep(n); + BigFraction[][] mTon = + NordsieckTransformer.buildMultistepToNordsieck(n); + for (int i = 0; i < n; ++i) { + for (int j = 0; j < n; ++j) { + BigFraction s = BigFraction.ZERO; + for (int k = 0; k < n; ++k) { + s = s.add(mTon[i][k].multiply(nTom[k][j])); + } + assertEquals((i == j) ? BigFraction.ONE : BigFraction.ZERO, s); + } + } + } + } + + public void testMatrices1() { + checkMatrix(1, new int[][] { { 1 } }, + NordsieckTransformer.buildMultistepWithoutDerivativesToNordsieck(1)); + checkMatrix(new int[][] { { 1 } }, + NordsieckTransformer.buildNordsieckToMultistepWithoutDerivatives(1)); + checkMatrix(1, new int[][] { { 1 } }, + NordsieckTransformer.buildMultistepToNordsieck(1)); + checkMatrix(new int[][] { { 1 } }, + NordsieckTransformer.buildNordsieckToMultistep(1)); + } + + public void testMatrices2() { + checkMatrix(1, new int[][] { { 1, 0 }, { 1, -1 } }, + NordsieckTransformer.buildMultistepWithoutDerivativesToNordsieck(2)); + checkMatrix(new int[][] { { 1, 0 }, { 1, -1 } }, + NordsieckTransformer.buildNordsieckToMultistepWithoutDerivatives(2)); + checkMatrix(1, new int[][] { { 1, 0 }, { 0, 1 } }, + NordsieckTransformer.buildMultistepToNordsieck(2)); + checkMatrix(new int[][] { { 1, 0 }, { 0, 1 } }, + NordsieckTransformer.buildNordsieckToMultistep(2)); + } + + public void testMatrices3() { + checkMatrix(2, new int[][] { { 2, 0, 0 }, { 3, -4, 1 }, { 1, -2, 1 } }, + NordsieckTransformer.buildMultistepWithoutDerivativesToNordsieck(3)); + checkMatrix(new int[][] { { 1, 0, 0 }, { 1, -1, 1 }, { 1, -2, 4 } }, + NordsieckTransformer.buildNordsieckToMultistepWithoutDerivatives(3)); + checkMatrix(1, new int[][] { { 1, 0, 0 }, { 0, 1, 0 }, { -1, 1, 1} }, + NordsieckTransformer.buildMultistepToNordsieck(3)); + checkMatrix(new int[][] { { 1, 0, 0 }, { 0, 1, 0 }, { 1, -1, 1 } }, + NordsieckTransformer.buildNordsieckToMultistep(3)); + } + + public void testMatrices4() { + checkMatrix(6, new int[][] { { 6, 0, 0, 0 }, { 11, -18, 9, -2 }, { 6, -15, 12, -3 }, { 1, -3, 3, -1 } }, + NordsieckTransformer.buildMultistepWithoutDerivativesToNordsieck(4)); + checkMatrix(new int[][] { { 1, 0, 0, 0 }, { 1, -1, 1, -1 }, { 1, -2, 4, -8 }, { 1, -3, 9, -27 } }, + NordsieckTransformer.buildNordsieckToMultistepWithoutDerivatives(4)); + checkMatrix(4, new int[][] { { 4, 0, 0, 0 }, { 0, 4, 0, 0 }, { -7, 6, 8, -1 }, { -3, 2, 4, -1 } }, + NordsieckTransformer.buildMultistepToNordsieck(4)); + checkMatrix(new int[][] { { 1, 0, 0, 0 }, { 0, 1, 0, 0 }, { 1, -1, 1, -1 }, { 1, -2, 4, -8 } }, + NordsieckTransformer.buildNordsieckToMultistep(4)); + } + + public void testPolynomial() { + Random r = new Random(1847222905841997856l); + for (int n = 2; n < 9; ++n) { + + // build a polynomial and its derivatives + double[] coeffs = new double[n + 1]; + for (int i = 0; i < n; ++i) { + coeffs[i] = 2 * r.nextDouble() - 1.0; + } + PolynomialFunction[] polynomials = new PolynomialFunction[n]; + polynomials[0] = new PolynomialFunction(coeffs); + for (int k = 1; k < polynomials.length; ++k) { + polynomials[k] = (PolynomialFunction) polynomials[k - 1].derivative(); + } + double h = 0.01; + + // build a state history in multistep form + double[] multistepHistory = new double[n]; + multistepHistory[0] = polynomials[0].value(1.0); + multistepHistory[1] = h * polynomials[1].value(1.0); + for (int i = 2; i < multistepHistory.length; ++i) { + multistepHistory[i] = polynomials[0].value(1.0 - (i - 1) * h); + } + + // build the same state history in multistep without derivatives form + double[] mwdHistory = new double[n]; + for (int i = 0; i < multistepHistory.length; ++i) { + mwdHistory[i] = polynomials[0].value(1.0 - i * h); + } + + // build the same state history in Nordsieck form + double[] nordsieckHistory = new double[n]; + double scale = 1.0; + for (int i = 0; i < nordsieckHistory.length; ++i) { + nordsieckHistory[i] = scale * polynomials[i].value(1.0); + scale *= h / (i + 1); + } + + // check the transform is exact for these polynomials states + NordsieckTransformer transformer = new NordsieckTransformer(n); + checkVector(nordsieckHistory, transformer.multistepWithoutDerivativesToNordsieck(mwdHistory)); + checkVector(mwdHistory, transformer.nordsieckToMultistepWithoutDerivatives(nordsieckHistory)); + checkVector(nordsieckHistory, transformer.multistepToNordsieck(multistepHistory)); + checkVector(multistepHistory, transformer.nordsieckToMultistep(nordsieckHistory)); + + } + } + + private void checkVector(double[] reference, double[] candidate) { + assertEquals(reference.length, candidate.length); + for (int i = 0; i < reference.length; ++i) { + assertEquals(reference[i], candidate[i], 1.0e-14); + } + } + + private void checkMatrix(int[][] reference, BigInteger[][] candidate) { + assertEquals(reference.length, candidate.length); + for (int i = 0; i < reference.length; ++i) { + int[] rRow = reference[i]; + BigInteger[] cRow = candidate[i]; + assertEquals(rRow.length, cRow.length); + for (int j = 0; j < rRow.length; ++j) { + assertEquals(rRow[j], cRow[j].intValue()); + } + } + } + + private void checkMatrix(int denominator, int[][] reference, BigFraction[][] candidate) { + assertEquals(reference.length, candidate.length); + for (int i = 0; i < reference.length; ++i) { + int[] rRow = reference[i]; + BigFraction[] cRow = candidate[i]; + assertEquals(rRow.length, cRow.length); + for (int j = 0; j < rRow.length; ++j) { + assertEquals(new BigFraction(rRow[j], denominator), cRow[j]); + } + } + } + + public static Test suite() { + return new TestSuite(NordsieckTransformerTest.class); + } + +}