LegendreRuleFactory.java

  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.math3.analysis.integration.gauss;

  18. import org.apache.commons.math3.exception.DimensionMismatchException;
  19. import org.apache.commons.math3.util.Pair;

  20. /**
  21.  * Factory that creates Gauss-type quadrature rule using Legendre polynomials.
  22.  * In this implementation, the lower and upper bounds of the natural interval
  23.  * of integration are -1 and 1, respectively.
  24.  * The Legendre polynomials are evaluated using the recurrence relation
  25.  * presented in <a href="http://en.wikipedia.org/wiki/Abramowitz_and_Stegun"
  26.  * Abramowitz and Stegun, 1964</a>.
  27.  *
  28.  * @since 3.1
  29.  */
  30. public class LegendreRuleFactory extends BaseRuleFactory<Double> {
  31.     /** {@inheritDoc} */
  32.     @Override
  33.     protected Pair<Double[], Double[]> computeRule(int numberOfPoints)
  34.         throws DimensionMismatchException {

  35.         if (numberOfPoints == 1) {
  36.             // Break recursion.
  37.             return new Pair<Double[], Double[]>(new Double[] { 0d },
  38.                                                 new Double[] { 2d });
  39.         }

  40.         // Get previous rule.
  41.         // If it has not been computed yet it will trigger a recursive call
  42.         // to this method.
  43.         final Double[] previousPoints = getRuleInternal(numberOfPoints - 1).getFirst();

  44.         // Compute next rule.
  45.         final Double[] points = new Double[numberOfPoints];
  46.         final Double[] weights = new Double[numberOfPoints];

  47.         // Find i-th root of P[n+1] by bracketing.
  48.         final int iMax = numberOfPoints / 2;
  49.         for (int i = 0; i < iMax; i++) {
  50.             // Lower-bound of the interval.
  51.             double a = (i == 0) ? -1 : previousPoints[i - 1].doubleValue();
  52.             // Upper-bound of the interval.
  53.             double b = (iMax == 1) ? 1 : previousPoints[i].doubleValue();
  54.             // P[j-1](a)
  55.             double pma = 1;
  56.             // P[j](a)
  57.             double pa = a;
  58.             // P[j-1](b)
  59.             double pmb = 1;
  60.             // P[j](b)
  61.             double pb = b;
  62.             for (int j = 1; j < numberOfPoints; j++) {
  63.                 final int two_j_p_1 = 2 * j + 1;
  64.                 final int j_p_1 = j + 1;
  65.                 // P[j+1](a)
  66.                 final double ppa = (two_j_p_1 * a * pa - j * pma) / j_p_1;
  67.                 // P[j+1](b)
  68.                 final double ppb = (two_j_p_1 * b * pb - j * pmb) / j_p_1;
  69.                 pma = pa;
  70.                 pa = ppa;
  71.                 pmb = pb;
  72.                 pb = ppb;
  73.             }
  74.             // Now pa = P[n+1](a), and pma = P[n](a) (same holds for b).
  75.             // Middle of the interval.
  76.             double c = 0.5 * (a + b);
  77.             // P[j-1](c)
  78.             double pmc = 1;
  79.             // P[j](c)
  80.             double pc = c;
  81.             boolean done = false;
  82.             while (!done) {
  83.                 done = b - a <= Math.ulp(c);
  84.                 pmc = 1;
  85.                 pc = c;
  86.                 for (int j = 1; j < numberOfPoints; j++) {
  87.                     // P[j+1](c)
  88.                     final double ppc = ((2 * j + 1) * c * pc - j * pmc) / (j + 1);
  89.                     pmc = pc;
  90.                     pc = ppc;
  91.                 }
  92.                 // Now pc = P[n+1](c) and pmc = P[n](c).
  93.                 if (!done) {
  94.                     if (pa * pc <= 0) {
  95.                         b = c;
  96.                         pmb = pmc;
  97.                         pb = pc;
  98.                     } else {
  99.                         a = c;
  100.                         pma = pmc;
  101.                         pa = pc;
  102.                     }
  103.                     c = 0.5 * (a + b);
  104.                 }
  105.             }
  106.             final double d = numberOfPoints * (pmc - c * pc);
  107.             final double w = 2 * (1 - c * c) / (d * d);

  108.             points[i] = c;
  109.             weights[i] = w;

  110.             final int idx = numberOfPoints - i - 1;
  111.             points[idx] = -c;
  112.             weights[idx] = w;
  113.         }
  114.         // If "numberOfPoints" is odd, 0 is a root.
  115.         // Note: as written, the test for oddness will work for negative
  116.         // integers too (although it is not necessary here), preventing
  117.         // a FindBugs warning.
  118.         if (numberOfPoints % 2 != 0) {
  119.             double pmc = 1;
  120.             for (int j = 1; j < numberOfPoints; j += 2) {
  121.                 pmc = -j * pmc / (j + 1);
  122.             }
  123.             final double d = numberOfPoints * pmc;
  124.             final double w = 2 / (d * d);

  125.             points[iMax] = 0d;
  126.             weights[iMax] = w;
  127.         }

  128.         return new Pair<Double[], Double[]>(points, weights);
  129.     }
  130. }