001/* 002 * Licensed to the Apache Software Foundation (ASF) under one or more 003 * contributor license agreements. See the NOTICE file distributed with 004 * this work for additional information regarding copyright ownership. 005 * The ASF licenses this file to You under the Apache License, Version 2.0 006 * (the "License"); you may not use this file except in compliance with 007 * the License. You may obtain a copy of the License at 008 * 009 * http://www.apache.org/licenses/LICENSE-2.0 010 * 011 * Unless required by applicable law or agreed to in writing, software 012 * distributed under the License is distributed on an "AS IS" BASIS, 013 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. 014 * See the License for the specific language governing permissions and 015 * limitations under the License. 016 */ 017package org.apache.commons.math3.analysis.integration.gauss; 018 019import org.apache.commons.math3.exception.DimensionMismatchException; 020import org.apache.commons.math3.util.Pair; 021import org.apache.commons.math3.util.FastMath; 022 023/** 024 * Factory that creates a 025 * <a href="http://en.wikipedia.org/wiki/Gauss-Hermite_quadrature"> 026 * Gauss-type quadrature rule using Hermite polynomials</a> 027 * of the first kind. 028 * Such a quadrature rule allows the calculation of improper integrals 029 * of a function 030 * <p> 031 * \(f(x) e^{-x^2}\) 032 * </p><p> 033 * Recurrence relation and weights computation follow 034 * <a href="http://en.wikipedia.org/wiki/Abramowitz_and_Stegun"> 035 * Abramowitz and Stegun, 1964</a>. 036 * </p><p> 037 * The coefficients of the standard Hermite polynomials grow very rapidly. 038 * In order to avoid overflows, each Hermite polynomial is normalized with 039 * respect to the underlying scalar product. 040 * The initial interval for the application of the bisection method is 041 * based on the roots of the previous Hermite polynomial (interlacing). 042 * Upper and lower bounds of these roots are provided by </p> 043 * <blockquote> 044 * I. Krasikov, 045 * <em>Nonnegative quadratic forms and bounds on orthogonal polynomials</em>, 046 * Journal of Approximation theory <b>111</b>, 31-49 047 * </blockquote> 048 * 049 * @since 3.3 050 */ 051public class HermiteRuleFactory extends BaseRuleFactory<Double> { 052 /** π<sup>1/2</sup> */ 053 private static final double SQRT_PI = 1.77245385090551602729; 054 /** π<sup>-1/4</sup> */ 055 private static final double H0 = 7.5112554446494248286e-1; 056 /** π<sup>-1/4</sup> √2 */ 057 private static final double H1 = 1.0622519320271969145; 058 059 /** {@inheritDoc} */ 060 @Override 061 protected Pair<Double[], Double[]> computeRule(int numberOfPoints) 062 throws DimensionMismatchException { 063 064 if (numberOfPoints == 1) { 065 // Break recursion. 066 return new Pair<Double[], Double[]>(new Double[] { 0d }, 067 new Double[] { SQRT_PI }); 068 } 069 070 // Get previous rule. 071 // If it has not been computed yet it will trigger a recursive call 072 // to this method. 073 final int lastNumPoints = numberOfPoints - 1; 074 final Double[] previousPoints = getRuleInternal(lastNumPoints).getFirst(); 075 076 // Compute next rule. 077 final Double[] points = new Double[numberOfPoints]; 078 final Double[] weights = new Double[numberOfPoints]; 079 080 final double sqrtTwoTimesLastNumPoints = FastMath.sqrt(2 * lastNumPoints); 081 final double sqrtTwoTimesNumPoints = FastMath.sqrt(2 * numberOfPoints); 082 083 // Find i-th root of H[n+1] by bracketing. 084 final int iMax = numberOfPoints / 2; 085 for (int i = 0; i < iMax; i++) { 086 // Lower-bound of the interval. 087 double a = (i == 0) ? -sqrtTwoTimesLastNumPoints : previousPoints[i - 1].doubleValue(); 088 // Upper-bound of the interval. 089 double b = (iMax == 1) ? -0.5 : previousPoints[i].doubleValue(); 090 091 // H[j-1](a) 092 double hma = H0; 093 // H[j](a) 094 double ha = H1 * a; 095 // H[j-1](b) 096 double hmb = H0; 097 // H[j](b) 098 double hb = H1 * b; 099 for (int j = 1; j < numberOfPoints; j++) { 100 // Compute H[j+1](a) and H[j+1](b) 101 final double jp1 = j + 1; 102 final double s = FastMath.sqrt(2 / jp1); 103 final double sm = FastMath.sqrt(j / jp1); 104 final double hpa = s * a * ha - sm * hma; 105 final double hpb = s * b * hb - sm * hmb; 106 hma = ha; 107 ha = hpa; 108 hmb = hb; 109 hb = hpb; 110 } 111 112 // Now ha = H[n+1](a), and hma = H[n](a) (same holds for b). 113 // Middle of the interval. 114 double c = 0.5 * (a + b); 115 // P[j-1](c) 116 double hmc = H0; 117 // P[j](c) 118 double hc = H1 * c; 119 boolean done = false; 120 while (!done) { 121 done = b - a <= Math.ulp(c); 122 hmc = H0; 123 hc = H1 * c; 124 for (int j = 1; j < numberOfPoints; j++) { 125 // Compute H[j+1](c) 126 final double jp1 = j + 1; 127 final double s = FastMath.sqrt(2 / jp1); 128 final double sm = FastMath.sqrt(j / jp1); 129 final double hpc = s * c * hc - sm * hmc; 130 hmc = hc; 131 hc = hpc; 132 } 133 // Now h = H[n+1](c) and hm = H[n](c). 134 if (!done) { 135 if (ha * hc < 0) { 136 b = c; 137 hmb = hmc; 138 hb = hc; 139 } else { 140 a = c; 141 hma = hmc; 142 ha = hc; 143 } 144 c = 0.5 * (a + b); 145 } 146 } 147 final double d = sqrtTwoTimesNumPoints * hmc; 148 final double w = 2 / (d * d); 149 150 points[i] = c; 151 weights[i] = w; 152 153 final int idx = lastNumPoints - i; 154 points[idx] = -c; 155 weights[idx] = w; 156 } 157 158 // If "numberOfPoints" is odd, 0 is a root. 159 // Note: as written, the test for oddness will work for negative 160 // integers too (although it is not necessary here), preventing 161 // a FindBugs warning. 162 if (numberOfPoints % 2 != 0) { 163 double hm = H0; 164 for (int j = 1; j < numberOfPoints; j += 2) { 165 final double jp1 = j + 1; 166 hm = -FastMath.sqrt(j / jp1) * hm; 167 } 168 final double d = sqrtTwoTimesNumPoints * hm; 169 final double w = 2 / (d * d); 170 171 points[iMax] = 0d; 172 weights[iMax] = w; 173 } 174 175 return new Pair<Double[], Double[]>(points, weights); 176 } 177}