001 /*
002 * Licensed to the Apache Software Foundation (ASF) under one or more
003 * contributor license agreements. See the NOTICE file distributed with
004 * this work for additional information regarding copyright ownership.
005 * The ASF licenses this file to You under the Apache License, Version 2.0
006 * (the "License"); you may not use this file except in compliance with
007 * the License. You may obtain a copy of the License at
008 *
009 * http://www.apache.org/licenses/LICENSE-2.0
010 *
011 * Unless required by applicable law or agreed to in writing, software
012 * distributed under the License is distributed on an "AS IS" BASIS,
013 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
014 * See the License for the specific language governing permissions and
015 * limitations under the License.
016 */
017 package org.apache.commons.math3.analysis.integration.gauss;
018
019 import org.apache.commons.math3.util.Pair;
020 import org.apache.commons.math3.exception.NotStrictlyPositiveException;
021 import org.apache.commons.math3.exception.util.LocalizedFormats;
022
023 /**
024 * Factory that creates Gauss-type quadrature rule using Legendre polynomials.
025 * In this implementation, the lower and upper bounds of the natural interval
026 * of integration are -1 and 1, respectively.
027 * The Legendre polynomials are evaluated using the recurrence relation
028 * presented in <a href="http://en.wikipedia.org/wiki/Abramowitz_and_Stegun"
029 * Abramowitz and Stegun, 1964</a>.
030 *
031 * @since 3.1
032 * @version $Id: LegendreRuleFactory.java 1382197 2012-09-07 22:35:01Z erans $
033 */
034 public class LegendreRuleFactory extends BaseRuleFactory<Double> {
035 /**
036 * {@inheritDoc}
037 */
038 @Override
039 protected Pair<Double[], Double[]> computeRule(int numberOfPoints)
040 throws NotStrictlyPositiveException {
041 if (numberOfPoints <= 0) {
042 throw new NotStrictlyPositiveException(LocalizedFormats.NUMBER_OF_POINTS,
043 numberOfPoints);
044 }
045
046 if (numberOfPoints == 1) {
047 // Break recursion.
048 return new Pair<Double[], Double[]>(new Double[] { 0d },
049 new Double[] { 2d });
050 }
051
052 // Get previous rule.
053 // If it has not been computed yet it will trigger a recursive call
054 // to this method.
055 final Double[] previousPoints = getRuleInternal(numberOfPoints - 1).getFirst();
056
057 // Compute next rule.
058 final Double[] points = new Double[numberOfPoints];
059 final Double[] weights = new Double[numberOfPoints];
060
061 // Find i-th root of P[n+1] by bracketing.
062 final int iMax = numberOfPoints / 2;
063 for (int i = 0; i < iMax; i++) {
064 // Lower-bound of the interval.
065 double a = (i == 0) ? -1 : previousPoints[i - 1].doubleValue();
066 // Upper-bound of the interval.
067 double b = (iMax == 1) ? 1 : previousPoints[i].doubleValue();
068 // P[j-1](a)
069 double pma = 1;
070 // P[j](a)
071 double pa = a;
072 // P[j-1](b)
073 double pmb = 1;
074 // P[j](b)
075 double pb = b;
076 for (int j = 1; j < numberOfPoints; j++) {
077 final int two_j_p_1 = 2 * j + 1;
078 final int j_p_1 = j + 1;
079 // P[j+1](a)
080 final double ppa = (two_j_p_1 * a * pa - j * pma) / j_p_1;
081 // P[j+1](b)
082 final double ppb = (two_j_p_1 * b * pb - j * pmb) / j_p_1;
083 pma = pa;
084 pa = ppa;
085 pmb = pb;
086 pb = ppb;
087 }
088 // Now pa = P[n+1](a), and pma = P[n](a) (same holds for b).
089 // Middle of the interval.
090 double c = 0.5 * (a + b);
091 // P[j-1](c)
092 double pmc = 1;
093 // P[j](c)
094 double pc = c;
095 boolean done = false;
096 while (!done) {
097 done = b - a <= Math.ulp(c);
098 pmc = 1;
099 pc = c;
100 for (int j = 1; j < numberOfPoints; j++) {
101 // P[j+1](c)
102 final double ppc = ((2 * j + 1) * c * pc - j * pmc) / (j + 1);
103 pmc = pc;
104 pc = ppc;
105 }
106 // Now pc = P[n+1](c) and pmc = P[n](c).
107 if (!done) {
108 if (pa * pc <= 0) {
109 b = c;
110 pmb = pmc;
111 pb = pc;
112 } else {
113 a = c;
114 pma = pmc;
115 pa = pc;
116 }
117 c = 0.5 * (a + b);
118 }
119 }
120 final double d = numberOfPoints * (pmc - c * pc);
121 final double w = 2 * (1 - c * c) / (d * d);
122
123 points[i] = c;
124 weights[i] = w;
125
126 final int idx = numberOfPoints - i - 1;
127 points[idx] = -c;
128 weights[idx] = w;
129 }
130 // If "numberOfPoints" is odd, 0 is a root.
131 // Note: as written, the test for oddness will work for negative
132 // integers too (although it is not necessary here), preventing
133 // a FindBugs warning.
134 if (numberOfPoints % 2 != 0) {
135 double pmc = 1;
136 for (int j = 1; j < numberOfPoints; j += 2) {
137 pmc = -j * pmc / (j + 1);
138 }
139 final double d = numberOfPoints * pmc;
140 final double w = 2 / (d * d);
141
142 points[iMax] = 0d;
143 weights[iMax] = w;
144 }
145
146 return new Pair<Double[], Double[]>(points, weights);
147 }
148 }