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.util.Pair;
25
26 /**
27 * Factory that creates Gauss-type quadrature rule using Laguerre polynomials.
28 *
29 * @see <a href="http://en.wikipedia.org/wiki/Gauss%E2%80%93Laguerre_quadrature">Gauss-Laguerre quadrature (Wikipedia)</a>
30 */
31 public class LaguerreRuleFactory extends AbstractRuleFactory {
32
33 /** Empty constructor.
34 * <p>
35 * This constructor is not strictly necessary, but it prevents spurious
36 * javadoc warnings with JDK 18 and later.
37 * </p>
38 * @since 3.0
39 */
40 public LaguerreRuleFactory() { // NOPMD - unnecessary constructor added intentionally to make javadoc happy
41 // nothing to do
42 }
43
44 /** {@inheritDoc} */
45 @Override
46 protected Pair<double[], double[]> computeRule(int numberOfPoints) {
47
48 // find nodes as roots of Laguerre polynomial
49 final double[] points = findRoots(numberOfPoints, new Laguerre(numberOfPoints)::ratio);
50
51 // compute weights
52 final double[] weights = new double[numberOfPoints];
53 final int n1 = numberOfPoints + 1;
54 final long n1Squared = n1 * (long) n1;
55 final Laguerre laguerreN1 = new Laguerre(n1);
56 for (int i = 0; i < numberOfPoints; i++) {
57 final double val = laguerreN1.value(points[i]);
58 weights[i] = points[i] / (n1Squared * val * val);
59 }
60
61 return new Pair<>(points, weights);
62
63 }
64
65 /** Laguerre polynomial. */
66 private static class Laguerre {
67
68 /** Degree. */
69 private int degree;
70
71 /** Simple constructor.
72 * @param degree polynomial degree
73 */
74 Laguerre(int degree) {
75 this.degree = degree;
76 }
77
78 /** Evaluate polynomial.
79 * @param x point at which polynomial must be evaluated
80 * @return value of the polynomial
81 */
82 public double value(final double x) {
83 return lNlNm1(x)[0];
84 }
85
86 /** Compute ratio L(x)/L'(x).
87 * @param x point at which ratio must be computed
88 * @return ratio L(x)/L'(x)
89 */
90 public double ratio(double x) {
91 double[] l = lNlNm1(x);
92 return x * l[0] / (degree * (l[0] - l[1]));
93 }
94
95 /** Compute Lₙ(x) and Lₙ₋₁(x).
96 * @param x point at which polynomials are evaluated
97 * @return array containing Lₙ(x) at index 0 and Lₙ₋₁(x) at index 1
98 */
99 private double[] lNlNm1(final double x) {
100 double[] l = { 1 - x, 1 };
101 for (int n = 1; n < degree; n++) {
102 // apply recurrence relation (n+1) Lₙ₊₁(x) = (2n + 1 - x) Lₙ(x) - n Lₙ₋₁(x)
103 final double lp = (l[0] * (2 * n + 1 - x) - l[1] * n) / (n + 1);
104 l[1] = l[0];
105 l[0] = lp;
106 }
107 return l;
108 }
109
110 }
111
112 }