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 }