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.legacy.core.Pair;
20
21 /**
22 * Factory that creates Gauss-type quadrature rule using Legendre polynomials.
23 * In this implementation, the lower and upper bounds of the natural interval
24 * of integration are -1 and 1, respectively.
25 * The Legendre polynomials are evaluated using the recurrence relation
26 * presented in <a href="http://en.wikipedia.org/wiki/Abramowitz_and_Stegun">
27 * Abramowitz and Stegun, 1964</a>.
28 *
29 * @since 3.1
30 */
31 public class LegendreRuleFactory extends BaseRuleFactory<Double> {
32 /** {@inheritDoc} */
33 @Override
34 protected Pair<Double[], Double[]> computeRule(int numberOfPoints) {
35 if (numberOfPoints == 1) {
36 // Break recursion.
37 return new Pair<>(new Double[] { 0d },
38 new Double[] { 2d });
39 }
40
41 // Get previous rule.
42 // If it has not been computed yet it will trigger a recursive call
43 // to this method.
44 final Double[] previousPoints = getRuleInternal(numberOfPoints - 1).getFirst();
45
46 // Compute next rule.
47 final Double[] points = new Double[numberOfPoints];
48 final Double[] weights = new Double[numberOfPoints];
49
50 // Find i-th root of P[n+1] by bracketing.
51 final int iMax = numberOfPoints / 2;
52 for (int i = 0; i < iMax; i++) {
53 // Lower-bound of the interval.
54 double a = (i == 0) ? -1 : previousPoints[i - 1].doubleValue();
55 // Upper-bound of the interval.
56 double b = (iMax == 1) ? 1 : previousPoints[i].doubleValue();
57 // P[j-1](a)
58 double pma = 1;
59 // P[j](a)
60 double pa = a;
61 // P[j-1](b)
62 double pmb = 1;
63 // P[j](b)
64 double pb = b;
65 for (int j = 1; j < numberOfPoints; j++) {
66 final int two_j_p_1 = 2 * j + 1;
67 final int j_p_1 = j + 1;
68 // P[j+1](a)
69 final double ppa = (two_j_p_1 * a * pa - j * pma) / j_p_1;
70 // P[j+1](b)
71 final double ppb = (two_j_p_1 * b * pb - j * pmb) / j_p_1;
72 pma = pa;
73 pa = ppa;
74 pmb = pb;
75 pb = ppb;
76 }
77 // Now pa = P[n+1](a), and pma = P[n](a) (same holds for b).
78 // Middle of the interval.
79 double c = 0.5 * (a + b);
80 // P[j-1](c)
81 double pmc = 1;
82 // P[j](c)
83 double pc = c;
84 boolean done = false;
85 while (!done) {
86 done = b - a <= Math.ulp(c);
87 pmc = 1;
88 pc = c;
89 for (int j = 1; j < numberOfPoints; j++) {
90 // P[j+1](c)
91 final double ppc = ((2 * j + 1) * c * pc - j * pmc) / (j + 1);
92 pmc = pc;
93 pc = ppc;
94 }
95 // Now pc = P[n+1](c) and pmc = P[n](c).
96 if (!done) {
97 if (pa * pc <= 0) {
98 b = c;
99 pmb = pmc;
100 pb = pc;
101 } else {
102 a = c;
103 pma = pmc;
104 pa = pc;
105 }
106 c = 0.5 * (a + b);
107 }
108 }
109 final double d = numberOfPoints * (pmc - c * pc);
110 final double w = 2 * (1 - c * c) / (d * d);
111
112 points[i] = c;
113 weights[i] = w;
114
115 final int idx = numberOfPoints - i - 1;
116 points[idx] = -c;
117 weights[idx] = w;
118 }
119 // If "numberOfPoints" is odd, 0 is a root.
120 // Note: as written, the test for oddness will work for negative
121 // integers too (although it is not necessary here), preventing
122 // a FindBugs warning.
123 if ((numberOfPoints & 1) != 0) {
124 double pmc = 1;
125 for (int j = 1; j < numberOfPoints; j += 2) {
126 pmc = -j * pmc / (j + 1);
127 }
128 final double d = numberOfPoints * pmc;
129 final double w = 2 / (d * d);
130
131 points[iMax] = 0d;
132 weights[iMax] = w;
133 }
134
135 return new Pair<>(points, weights);
136 }
137 }