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.math4.legacy.analysis.integration.gauss; 018 019import org.apache.commons.math4.core.jdkmath.JdkMath; 020import org.apache.commons.math4.legacy.core.Pair; 021 022/** 023 * Factory that creates a 024 * <a href="http://en.wikipedia.org/wiki/Gauss-Hermite_quadrature"> 025 * Gauss-type quadrature rule using Hermite polynomials</a> 026 * of the first kind. 027 * Such a quadrature rule allows the calculation of improper integrals 028 * of a function 029 * <p> 030 * \(f(x) e^{-x^2}\) 031 * </p><p> 032 * Recurrence relation and weights computation follow 033 * <a href="http://en.wikipedia.org/wiki/Abramowitz_and_Stegun"> 034 * Abramowitz and Stegun, 1964</a>. 035 * </p><p> 036 * The coefficients of the standard Hermite polynomials grow very rapidly. 037 * In order to avoid overflows, each Hermite polynomial is normalized with 038 * respect to the underlying scalar product. 039 * The initial interval for the application of the bisection method is 040 * based on the roots of the previous Hermite polynomial (interlacing). 041 * Upper and lower bounds of these roots are provided by </p> 042 * <blockquote> 043 * I. Krasikov, 044 * <em>Nonnegative quadratic forms and bounds on orthogonal polynomials</em>, 045 * Journal of Approximation theory <b>111</b>, 31-49 046 * </blockquote> 047 * 048 * @since 3.3 049 */ 050public class HermiteRuleFactory extends BaseRuleFactory<Double> { 051 /** π<sup>1/2</sup>. */ 052 private static final double SQRT_PI = 1.77245385090551602729; 053 /** π<sup>-1/4</sup>. */ 054 private static final double H0 = 7.5112554446494248286e-1; 055 /** π<sup>-1/4</sup> √2. */ 056 private static final double H1 = 1.0622519320271969145; 057 058 /** {@inheritDoc} */ 059 @Override 060 protected Pair<Double[], Double[]> computeRule(int numberOfPoints) { 061 if (numberOfPoints == 1) { 062 // Break recursion. 063 return new Pair<>(new Double[] { 0d }, 064 new Double[] { SQRT_PI }); 065 } 066 067 // Get previous rule. 068 // If it has not been computed yet it will trigger a recursive call 069 // to this method. 070 final int lastNumPoints = numberOfPoints - 1; 071 final Double[] previousPoints = getRuleInternal(lastNumPoints).getFirst(); 072 073 // Compute next rule. 074 final Double[] points = new Double[numberOfPoints]; 075 final Double[] weights = new Double[numberOfPoints]; 076 077 final double sqrtTwoTimesLastNumPoints = JdkMath.sqrt(2.0 * lastNumPoints); 078 final double sqrtTwoTimesNumPoints = JdkMath.sqrt(2.0 * numberOfPoints); 079 080 // Find i-th root of H[n+1] by bracketing. 081 final int iMax = numberOfPoints / 2; 082 for (int i = 0; i < iMax; i++) { 083 // Lower-bound of the interval. 084 double a = (i == 0) ? -sqrtTwoTimesLastNumPoints : previousPoints[i - 1].doubleValue(); 085 // Upper-bound of the interval. 086 double b = (iMax == 1) ? -0.5 : previousPoints[i].doubleValue(); 087 088 // H[j-1](a) 089 double hma = H0; 090 // H[j](a) 091 double ha = H1 * a; 092 // H[j-1](b) 093 double hmb = H0; 094 // H[j](b) 095 double hb = H1 * b; 096 for (int j = 1; j < numberOfPoints; j++) { 097 // Compute H[j+1](a) and H[j+1](b) 098 final double jp1 = j + 1.0; 099 final double s = JdkMath.sqrt(2 / jp1); 100 final double sm = JdkMath.sqrt(j / jp1); 101 final double hpa = s * a * ha - sm * hma; 102 final double hpb = s * b * hb - sm * hmb; 103 hma = ha; 104 ha = hpa; 105 hmb = hb; 106 hb = hpb; 107 } 108 109 // Now ha = H[n+1](a), and hma = H[n](a) (same holds for b). 110 // Middle of the interval. 111 double c = 0.5 * (a + b); 112 // P[j-1](c) 113 double hmc = H0; 114 // P[j](c) 115 double hc = H1 * c; 116 boolean done = false; 117 while (!done) { 118 done = b - a <= Math.ulp(c); 119 hmc = H0; 120 hc = H1 * c; 121 for (int j = 1; j < numberOfPoints; j++) { 122 // Compute H[j+1](c) 123 final double jp1 = j + 1.0; 124 final double s = JdkMath.sqrt(2 / jp1); 125 final double sm = JdkMath.sqrt(j / jp1); 126 final double hpc = s * c * hc - sm * hmc; 127 hmc = hc; 128 hc = hpc; 129 } 130 // Now h = H[n+1](c) and hm = H[n](c). 131 if (!done) { 132 if (ha * hc < 0) { 133 b = c; 134 hmb = hmc; 135 hb = hc; 136 } else { 137 a = c; 138 hma = hmc; 139 ha = hc; 140 } 141 c = 0.5 * (a + b); 142 } 143 } 144 final double d = sqrtTwoTimesNumPoints * hmc; 145 final double w = 2 / (d * d); 146 147 points[i] = c; 148 weights[i] = w; 149 150 final int idx = lastNumPoints - i; 151 points[idx] = -c; 152 weights[idx] = w; 153 } 154 155 // If "numberOfPoints" is odd, 0 is a root. 156 // Note: as written, the test for oddness will work for negative 157 // integers too (although it is not necessary here), preventing 158 // a FindBugs warning. 159 if ((numberOfPoints & 1) != 0) { 160 double hm = H0; 161 for (int j = 1; j < numberOfPoints; j += 2) { 162 final double jp1 = j + 1.0; 163 hm = -JdkMath.sqrt(j / jp1) * hm; 164 } 165 final double d = sqrtTwoTimesNumPoints * hm; 166 final double w = 2 / (d * d); 167 168 points[iMax] = 0d; 169 weights[iMax] = w; 170 } 171 172 return new Pair<>(points, weights); 173 } 174}