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 * https://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
18 /*
19 * This is not the original file distributed by the Apache Software Foundation
20 * It has been modified by the Hipparchus project
21 */
22 package org.hipparchus.analysis.integration.gauss;
23
24 import org.hipparchus.exception.MathIllegalArgumentException;
25 import org.hipparchus.util.FastMath;
26 import org.hipparchus.util.Pair;
27
28 /**
29 * Factory that creates a
30 * <a href="http://en.wikipedia.org/wiki/Gauss-Hermite_quadrature">
31 * Gauss-type quadrature rule using Hermite polynomials</a>
32 * of the first kind.
33 * Such a quadrature rule allows the calculation of improper integrals
34 * of a function
35 * <p>
36 * \(f(x) e^{-x^2}\)
37 * </p>
38 * <p>
39 * Recurrence relation and weights computation follow
40 * <a href="http://en.wikipedia.org/wiki/Abramowitz_and_Stegun">
41 * Abramowitz and Stegun, 1964</a>.
42 * </p>
43 *
44 */
45 public class HermiteRuleFactory extends AbstractRuleFactory {
46
47 /** √π. */
48 private static final double SQRT_PI = 1.77245385090551602729;
49
50 /** Empty constructor.
51 * <p>
52 * This constructor is not strictly necessary, but it prevents spurious
53 * javadoc warnings with JDK 18 and later.
54 * </p>
55 * @since 3.0
56 */
57 public HermiteRuleFactory() { // NOPMD - unnecessary constructor added intentionally to make javadoc happy
58 // nothing to do
59 }
60
61 /** {@inheritDoc} */
62 @Override
63 protected Pair<double[], double[]> computeRule(int numberOfPoints)
64 throws MathIllegalArgumentException {
65
66 if (numberOfPoints == 1) {
67 // Break recursion.
68 return new Pair<>(new double[] { 0 } , new double[] { SQRT_PI });
69 }
70
71 // find nodes as roots of Hermite polynomial
72 final double[] points = findRoots(numberOfPoints, new Hermite(numberOfPoints)::ratio);
73 enforceSymmetry(points);
74
75 // compute weights
76 final double[] weights = new double[numberOfPoints];
77 final Hermite hm1 = new Hermite(numberOfPoints - 1);
78 for (int i = 0; i < numberOfPoints; i++) {
79 final double y = hm1.hNhNm1(points[i])[0];
80 weights[i] = SQRT_PI / (numberOfPoints * y * y);
81 }
82
83 return new Pair<>(points, weights);
84
85 }
86
87 /** Hermite polynomial, normalized to avoid overflow.
88 * <p>
89 * The regular Hermite polynomials and associated weights are given by:
90 * <pre>
91 * H₀(x) = 1
92 * H₁(x) = 2 x
93 * Hₙ₊₁(x) = 2x Hₙ(x) - 2n Hₙ₋₁(x), and H'ₙ(x) = 2n Hₙ₋₁(x)
94 * wₙ(xᵢ) = [2ⁿ⁻¹ n! √π]/[n Hₙ₋₁(xᵢ)]²
95 * </pre>
96 * </p>
97 * <p>
98 * In order to avoid overflow with normalize the polynomials hₙ(x) = Hₙ(x) / √[2ⁿ n!]
99 * so the recurrence relations and weights become:
100 * <pre>
101 * h₀(x) = 1
102 * h₁(x) = √2 x
103 * hₙ₊₁(x) = [√2 x hₙ(x) - √n hₙ₋₁(x)]/√(n+1), and h'ₙ(x) = 2n hₙ₋₁(x)
104 * uₙ(xᵢ) = √π/[n Nₙ₋₁(xᵢ)²]
105 * </pre>
106 * </p>
107 */
108 private static class Hermite {
109
110 /** √2. */
111 private static final double SQRT2 = FastMath.sqrt(2);
112
113 /** Degree. */
114 private final int degree;
115
116 /** Simple constructor.
117 * @param degree polynomial degree
118 */
119 Hermite(int degree) {
120 this.degree = degree;
121 }
122
123 /** Compute ratio H(x)/H'(x).
124 * @param x point at which ratio must be computed
125 * @return ratio H(x)/H'(x)
126 */
127 public double ratio(double x) {
128 double[] h = hNhNm1(x);
129 return h[0] / (h[1] * 2 * degree);
130 }
131
132 /** Compute Nₙ(x) and Nₙ₋₁(x).
133 * @param x point at which polynomials are evaluated
134 * @return array containing Nₙ(x) at index 0 and Nₙ₋₁(x) at index 1
135 */
136 private double[] hNhNm1(final double x) {
137 double[] h = { SQRT2 * x, 1 };
138 double sqrtN = 1;
139 for (int n = 1; n < degree; n++) {
140 // apply recurrence relation hₙ₊₁(x) = [√2 x hₙ(x) - √n hₙ₋₁(x)]/√(n+1)
141 final double sqrtNp = FastMath.sqrt(n + 1);
142 final double hp = (h[0] * x * SQRT2 - h[1] * sqrtN) / sqrtNp;
143 h[1] = h[0];
144 h[0] = hp;
145 sqrtN = sqrtNp;
146 }
147 return h;
148 }
149
150 }
151
152 }