1 /*
2 * Licensed to the Apache Software Foundation (ASF) under one or more
3 * contributor license agreements. See the NOTICE file distributed with
4 * this work for additional information regarding copyright ownership.
5 * The ASF licenses this file to You under the Apache License, Version 2.0
6 * (the "License"); you may not use this file except in compliance with
7 * the License. You may obtain a copy of the License at
8 *
9 * http://www.apache.org/licenses/LICENSE-2.0
10 *
11 * Unless required by applicable law or agreed to in writing, software
12 * distributed under the License is distributed on an "AS IS" BASIS,
13 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14 * See the License for the specific language governing permissions and
15 * limitations under the License.
16 */
17 package org.apache.commons.math4.legacy.analysis.integration.gauss;
18
19 import org.apache.commons.math4.core.jdkmath.JdkMath;
20 import org.apache.commons.math4.legacy.core.Pair;
21
22 /**
23 * Factory that creates a
24 * <a href="http://en.wikipedia.org/wiki/Gauss-Hermite_quadrature">
25 * Gauss-type quadrature rule using Hermite polynomials</a>
26 * of the first kind.
27 * Such a quadrature rule allows the calculation of improper integrals
28 * of a function
29 * <p>
30 * \(f(x) e^{-x^2}\)
31 * </p><p>
32 * Recurrence relation and weights computation follow
33 * <a href="http://en.wikipedia.org/wiki/Abramowitz_and_Stegun">
34 * Abramowitz and Stegun, 1964</a>.
35 * </p><p>
36 * The coefficients of the standard Hermite polynomials grow very rapidly.
37 * In order to avoid overflows, each Hermite polynomial is normalized with
38 * respect to the underlying scalar product.
39 * The initial interval for the application of the bisection method is
40 * based on the roots of the previous Hermite polynomial (interlacing).
41 * Upper and lower bounds of these roots are provided by </p>
42 * <blockquote>
43 * I. Krasikov,
44 * <em>Nonnegative quadratic forms and bounds on orthogonal polynomials</em>,
45 * Journal of Approximation theory <b>111</b>, 31-49
46 * </blockquote>
47 *
48 * @since 3.3
49 */
50 public class HermiteRuleFactory extends BaseRuleFactory<Double> {
51 /** π<sup>1/2</sup>. */
52 private static final double SQRT_PI = 1.77245385090551602729;
53 /** π<sup>-1/4</sup>. */
54 private static final double H0 = 7.5112554446494248286e-1;
55 /** π<sup>-1/4</sup> √2. */
56 private static final double H1 = 1.0622519320271969145;
57
58 /** {@inheritDoc} */
59 @Override
60 protected Pair<Double[], Double[]> computeRule(int numberOfPoints) {
61 if (numberOfPoints == 1) {
62 // Break recursion.
63 return new Pair<>(new Double[] { 0d },
64 new Double[] { SQRT_PI });
65 }
66
67 // Get previous rule.
68 // If it has not been computed yet it will trigger a recursive call
69 // to this method.
70 final int lastNumPoints = numberOfPoints - 1;
71 final Double[] previousPoints = getRuleInternal(lastNumPoints).getFirst();
72
73 // Compute next rule.
74 final Double[] points = new Double[numberOfPoints];
75 final Double[] weights = new Double[numberOfPoints];
76
77 final double sqrtTwoTimesLastNumPoints = JdkMath.sqrt(2.0 * lastNumPoints);
78 final double sqrtTwoTimesNumPoints = JdkMath.sqrt(2.0 * numberOfPoints);
79
80 // Find i-th root of H[n+1] by bracketing.
81 final int iMax = numberOfPoints / 2;
82 for (int i = 0; i < iMax; i++) {
83 // Lower-bound of the interval.
84 double a = (i == 0) ? -sqrtTwoTimesLastNumPoints : previousPoints[i - 1].doubleValue();
85 // Upper-bound of the interval.
86 double b = (iMax == 1) ? -0.5 : previousPoints[i].doubleValue();
87
88 // H[j-1](a)
89 double hma = H0;
90 // H[j](a)
91 double ha = H1 * a;
92 // H[j-1](b)
93 double hmb = H0;
94 // H[j](b)
95 double hb = H1 * b;
96 for (int j = 1; j < numberOfPoints; j++) {
97 // Compute H[j+1](a) and H[j+1](b)
98 final double jp1 = j + 1.0;
99 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 }