@Test(expectedExceptions = IllegalArgumentException.class) public void testBadN() { ORTHONORMAL.getPolynomials(-3); }
@Override public DoubleFunction1D[] getPolynomials(int n) { ArgChecker.isTrue(n >= 0); DoubleFunction1D[] polynomials = new DoubleFunction1D[n + 1]; for (int i = 0; i <= n; i++) { if (i == 0) { polynomials[i] = F0; } else if (i == 1) { polynomials[i] = polynomials[0].multiply(Math.sqrt(2)).multiply(getX()); } else { polynomials[i] = polynomials[i - 1].multiply(getX()).multiply(Math.sqrt(2. / i)) .subtract(polynomials[i - 2].multiply(Math.sqrt((i - 1.) / i))); } } return polynomials; }
@Override public Pair<DoubleFunction1D, DoubleFunction1D>[] getPolynomialsAndFirstDerivative(int n) { ArgChecker.isTrue(n >= 0); @SuppressWarnings("unchecked") Pair<DoubleFunction1D, DoubleFunction1D>[] polynomials = new Pair[n + 1]; DoubleFunction1D p, dp, p1, p2; double sqrt2 = Math.sqrt(2); DoubleFunction1D x = getX(); for (int i = 0; i <= n; i++) { if (i == 0) { polynomials[i] = Pair.of((DoubleFunction1D) F0, getZero()); } else if (i == 1) { polynomials[i] = Pair.of(polynomials[0].getFirst().multiply(sqrt2).multiply(x), (DoubleFunction1D) DF1); } else { p1 = polynomials[i - 1].getFirst(); p2 = polynomials[i - 2].getFirst(); p = p1.multiply(x).multiply(Math.sqrt(2. / i)).subtract(p2.multiply(Math.sqrt((i - 1.) / i))); dp = p1.multiply(Math.sqrt(2 * i)); polynomials[i] = Pair.of(p, dp); } } return polynomials; }
@Override public GaussianQuadratureData generate(int n) { ArgChecker.isTrue(n > 0); double[] x = new double[n]; double[] w = new double[n]; boolean odd = n % 2 != 0; int m = (n + 1) / 2 - (odd ? 1 : 0); Pair<DoubleFunction1D, DoubleFunction1D>[] polynomials = HERMITE.getPolynomialsAndFirstDerivative(n); Pair<DoubleFunction1D, DoubleFunction1D> pair = polynomials[n]; DoubleFunction1D function = pair.getFirst(); DoubleFunction1D derivative = pair.getSecond(); double root = 0; for (int i = 0; i < m; i++) { root = getInitialRootGuess(root, i, n, x); root = ROOT_FINDER.getRoot(function, derivative, root); double dp = derivative.applyAsDouble(root); x[i] = -root; x[n - 1 - i] = root; w[i] = 2. / (dp * dp); w[n - 1 - i] = w[i]; } if (odd) { double dp = derivative.applyAsDouble(0.0); w[m] = 2. / dp / dp; } return new GaussianQuadratureData(x, w); }
@Test public void test() { final int n = 15; final DoubleFunction1D[] f1 = HERMITE.getPolynomials(n); final DoubleFunction1D[] f2 = ORTHONORMAL.getPolynomials(n); final double x = 3.4; for (int i = 0; i < f1.length; i++) { assertEquals( f1[i].applyAsDouble(x) / Math.sqrt(CombinatoricsUtils.factorialDouble(i) * Math.pow(2, i) * Math.sqrt(Math.PI)), f2[i].applyAsDouble(x), EPS); } }