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 */ 017package org.apache.commons.math3.analysis.integration.gauss; 018 019import org.apache.commons.math3.exception.DimensionMismatchException; 020import org.apache.commons.math3.util.Pair; 021 022/** 023 * Factory that creates Gauss-type quadrature rule using Legendre polynomials. 024 * In this implementation, the lower and upper bounds of the natural interval 025 * of integration are -1 and 1, respectively. 026 * The Legendre polynomials are evaluated using the recurrence relation 027 * presented in <a href="http://en.wikipedia.org/wiki/Abramowitz_and_Stegun"> 028 * Abramowitz and Stegun, 1964</a>. 029 * 030 * @since 3.1 031 */ 032public class LegendreRuleFactory extends BaseRuleFactory<Double> { 033 /** {@inheritDoc} */ 034 @Override 035 protected Pair<Double[], Double[]> computeRule(int numberOfPoints) 036 throws DimensionMismatchException { 037 038 if (numberOfPoints == 1) { 039 // Break recursion. 040 return new Pair<Double[], Double[]>(new Double[] { 0d }, 041 new Double[] { 2d }); 042 } 043 044 // Get previous rule. 045 // If it has not been computed yet it will trigger a recursive call 046 // to this method. 047 final Double[] previousPoints = getRuleInternal(numberOfPoints - 1).getFirst(); 048 049 // Compute next rule. 050 final Double[] points = new Double[numberOfPoints]; 051 final Double[] weights = new Double[numberOfPoints]; 052 053 // Find i-th root of P[n+1] by bracketing. 054 final int iMax = numberOfPoints / 2; 055 for (int i = 0; i < iMax; i++) { 056 // Lower-bound of the interval. 057 double a = (i == 0) ? -1 : previousPoints[i - 1].doubleValue(); 058 // Upper-bound of the interval. 059 double b = (iMax == 1) ? 1 : previousPoints[i].doubleValue(); 060 // P[j-1](a) 061 double pma = 1; 062 // P[j](a) 063 double pa = a; 064 // P[j-1](b) 065 double pmb = 1; 066 // P[j](b) 067 double pb = b; 068 for (int j = 1; j < numberOfPoints; j++) { 069 final int two_j_p_1 = 2 * j + 1; 070 final int j_p_1 = j + 1; 071 // P[j+1](a) 072 final double ppa = (two_j_p_1 * a * pa - j * pma) / j_p_1; 073 // P[j+1](b) 074 final double ppb = (two_j_p_1 * b * pb - j * pmb) / j_p_1; 075 pma = pa; 076 pa = ppa; 077 pmb = pb; 078 pb = ppb; 079 } 080 // Now pa = P[n+1](a), and pma = P[n](a) (same holds for b). 081 // Middle of the interval. 082 double c = 0.5 * (a + b); 083 // P[j-1](c) 084 double pmc = 1; 085 // P[j](c) 086 double pc = c; 087 boolean done = false; 088 while (!done) { 089 done = b - a <= Math.ulp(c); 090 pmc = 1; 091 pc = c; 092 for (int j = 1; j < numberOfPoints; j++) { 093 // P[j+1](c) 094 final double ppc = ((2 * j + 1) * c * pc - j * pmc) / (j + 1); 095 pmc = pc; 096 pc = ppc; 097 } 098 // Now pc = P[n+1](c) and pmc = P[n](c). 099 if (!done) { 100 if (pa * pc <= 0) { 101 b = c; 102 pmb = pmc; 103 pb = pc; 104 } else { 105 a = c; 106 pma = pmc; 107 pa = pc; 108 } 109 c = 0.5 * (a + b); 110 } 111 } 112 final double d = numberOfPoints * (pmc - c * pc); 113 final double w = 2 * (1 - c * c) / (d * d); 114 115 points[i] = c; 116 weights[i] = w; 117 118 final int idx = numberOfPoints - i - 1; 119 points[idx] = -c; 120 weights[idx] = w; 121 } 122 // If "numberOfPoints" is odd, 0 is a root. 123 // Note: as written, the test for oddness will work for negative 124 // integers too (although it is not necessary here), preventing 125 // a FindBugs warning. 126 if (numberOfPoints % 2 != 0) { 127 double pmc = 1; 128 for (int j = 1; j < numberOfPoints; j += 2) { 129 pmc = -j * pmc / (j + 1); 130 } 131 final double d = numberOfPoints * pmc; 132 final double w = 2 / (d * d); 133 134 points[iMax] = 0d; 135 weights[iMax] = w; 136 } 137 138 return new Pair<Double[], Double[]>(points, weights); 139 } 140}