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.math4.legacy.analysis.integration.gauss; 018 019import org.apache.commons.math4.legacy.core.Pair; 020 021/** 022 * Factory that creates Gauss-type quadrature rule using Legendre polynomials. 023 * In this implementation, the lower and upper bounds of the natural interval 024 * of integration are -1 and 1, respectively. 025 * The Legendre polynomials are evaluated using the recurrence relation 026 * presented in <a href="http://en.wikipedia.org/wiki/Abramowitz_and_Stegun"> 027 * Abramowitz and Stegun, 1964</a>. 028 * 029 * @since 3.1 030 */ 031public class LegendreRuleFactory extends BaseRuleFactory<Double> { 032 /** {@inheritDoc} */ 033 @Override 034 protected Pair<Double[], Double[]> computeRule(int numberOfPoints) { 035 if (numberOfPoints == 1) { 036 // Break recursion. 037 return new Pair<>(new Double[] { 0d }, 038 new Double[] { 2d }); 039 } 040 041 // Get previous rule. 042 // If it has not been computed yet it will trigger a recursive call 043 // to this method. 044 final Double[] previousPoints = getRuleInternal(numberOfPoints - 1).getFirst(); 045 046 // Compute next rule. 047 final Double[] points = new Double[numberOfPoints]; 048 final Double[] weights = new Double[numberOfPoints]; 049 050 // Find i-th root of P[n+1] by bracketing. 051 final int iMax = numberOfPoints / 2; 052 for (int i = 0; i < iMax; i++) { 053 // Lower-bound of the interval. 054 double a = (i == 0) ? -1 : previousPoints[i - 1].doubleValue(); 055 // Upper-bound of the interval. 056 double b = (iMax == 1) ? 1 : previousPoints[i].doubleValue(); 057 // P[j-1](a) 058 double pma = 1; 059 // P[j](a) 060 double pa = a; 061 // P[j-1](b) 062 double pmb = 1; 063 // P[j](b) 064 double pb = b; 065 for (int j = 1; j < numberOfPoints; j++) { 066 final int two_j_p_1 = 2 * j + 1; 067 final int j_p_1 = j + 1; 068 // P[j+1](a) 069 final double ppa = (two_j_p_1 * a * pa - j * pma) / j_p_1; 070 // P[j+1](b) 071 final double ppb = (two_j_p_1 * b * pb - j * pmb) / j_p_1; 072 pma = pa; 073 pa = ppa; 074 pmb = pb; 075 pb = ppb; 076 } 077 // Now pa = P[n+1](a), and pma = P[n](a) (same holds for b). 078 // Middle of the interval. 079 double c = 0.5 * (a + b); 080 // P[j-1](c) 081 double pmc = 1; 082 // P[j](c) 083 double pc = c; 084 boolean done = false; 085 while (!done) { 086 done = b - a <= Math.ulp(c); 087 pmc = 1; 088 pc = c; 089 for (int j = 1; j < numberOfPoints; j++) { 090 // P[j+1](c) 091 final double ppc = ((2 * j + 1) * c * pc - j * pmc) / (j + 1); 092 pmc = pc; 093 pc = ppc; 094 } 095 // Now pc = P[n+1](c) and pmc = P[n](c). 096 if (!done) { 097 if (pa * pc <= 0) { 098 b = c; 099 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}