1 /*
2 * Licensed to the Hipparchus project 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 Hipparchus project 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 package org.hipparchus.analysis.integration.gauss;
18
19 import org.hipparchus.CalculusFieldElement;
20 import org.hipparchus.Field;
21 import org.hipparchus.exception.MathIllegalArgumentException;
22 import org.hipparchus.util.MathArrays;
23 import org.hipparchus.util.Pair;
24
25 /**
26 * Factory that creates Gauss-type quadrature rule using Legendre polynomials.
27 * In this implementation, the lower and upper bounds of the natural interval
28 * of integration are -1 and 1, respectively.
29 * The Legendre polynomials are evaluated using the recurrence relation
30 * presented in <a href="http://en.wikipedia.org/wiki/Abramowitz_and_Stegun">
31 * Abramowitz and Stegun, 1964</a>.
32 *
33 * @param <T> Type of the number used to represent the points and weights of
34 * the quadrature rules.
35 * @since 2.0
36 */
37 public class FieldLegendreRuleFactory<T extends CalculusFieldElement<T>> extends FieldAbstractRuleFactory<T> {
38
39 /** Simple constructor
40 * @param field field to which rule coefficients belong
41 */
42 public FieldLegendreRuleFactory(final Field<T> field) {
43 super(field);
44 }
45
46 /** {@inheritDoc} */
47 @Override
48 public Pair<T[], T[]> computeRule(int numberOfPoints)
49 throws MathIllegalArgumentException {
50
51 final Field<T> field = getField();
52
53 if (numberOfPoints == 1) {
54 // Break recursion.
55 final T[] points = MathArrays.buildArray(field, numberOfPoints);
56 final T[] weights = MathArrays.buildArray(field, numberOfPoints);
57 points[0] = field.getZero();
58 weights[0] = field.getZero().newInstance(2);
59 return new Pair<>(points, weights);
60 }
61
62 // find nodes as roots of Legendre polynomial
63 final Legendre<T> p = new Legendre<>(numberOfPoints);
64 final T[] points = findRoots(numberOfPoints, p::ratio);
65 enforceSymmetry(points);
66
67 // compute weights
68 final T[] weights = MathArrays.buildArray(field, numberOfPoints);
69 for (int i = 0; i <= numberOfPoints / 2; i++) {
70 final T c = points[i];
71 final T[] pKpKm1 = p.pNpNm1(c);
72 final T d = pKpKm1[1].subtract(c.multiply(pKpKm1[0])).multiply(numberOfPoints);
73 weights[i] = c.square().subtract(1).multiply(-2).divide(d.multiply(d));
74
75 // symmetrical point
76 final int idx = numberOfPoints - i - 1;
77 weights[idx] = weights[i];
78
79 }
80
81 return new Pair<>(points, weights);
82
83 }
84
85 /** Legendre polynomial.
86 * @param <T> Type of the field elements.
87 */
88 private static class Legendre<T extends CalculusFieldElement<T>> {
89
90 /** Degree. */
91 private int degree;
92
93 /** Simple constructor.
94 * @param degree polynomial degree
95 */
96 Legendre(int degree) {
97 this.degree = degree;
98 }
99
100 /** Compute ratio P(x)/P'(x).
101 * @param x point at which ratio must be computed
102 * @return ratio P(x)/P'(x)
103 */
104 public T ratio(T x) {
105 T pm = x.getField().getOne();
106 T p = x;
107 T d = x.getField().getOne();
108 for (int n = 1; n < degree; n++) {
109 // apply recurrence relations (n+1) Pₙ₊₁(x) = (2n+1) x Pₙ(x) - n Pₙ₋₁(x)
110 // and P'ₙ₊₁(x) = (n+1) Pₙ(x) + x P'ₙ(x)
111 final T pp = p.multiply(x.multiply(2 * n + 1)).subtract(pm.multiply(n)).divide(n + 1);
112 d = p.multiply(n + 1).add(d.multiply(x));
113 pm = p;
114 p = pp;
115 }
116 return p.divide(d);
117 }
118
119 /** Compute Pₙ(x) and Pₙ₋₁(x).
120 * @param x point at which polynomials are evaluated
121 * @return array containing Pₙ(x) at index 0 and Pₙ₋₁(x) at index 1
122 */
123 private T[] pNpNm1(final T x) {
124 T[] p = MathArrays.buildArray(x.getField(), 2);
125 p[0] = x;
126 p[1] = x.getField().getOne();
127 for (int n = 1; n < degree; n++) {
128 // apply recurrence relation (n+1) Pₙ₊₁(x) = (2n+1) x Pₙ(x) - n Pₙ₋₁(x)
129 final T pp = p[0].multiply(x.multiply(2 * n + 1)).subtract(p[1].multiply(n)).divide(n + 1);
130 p[1] = p[0];
131 p[0] = pp;
132 }
133 return p;
134 }
135
136 }
137
138 }