LegendreHighPrecisionRuleFactory.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 java.math.BigDecimal;
  19. import java.math.MathContext;

  20. import org.apache.commons.math3.exception.DimensionMismatchException;
  21. import org.apache.commons.math3.util.Pair;

  22. /**
  23.  * Factory that creates Gauss-type quadrature rule using Legendre polynomials.
  24.  * In this implementation, the lower and upper bounds of the natural interval
  25.  * of integration are -1 and 1, respectively.
  26.  * The Legendre polynomials are evaluated using the recurrence relation
  27.  * presented in <a href="http://en.wikipedia.org/wiki/Abramowitz_and_Stegun"
  28.  * Abramowitz and Stegun, 1964</a>.
  29.  *
  30.  * @since 3.1
  31.  */
  32. public class LegendreHighPrecisionRuleFactory extends BaseRuleFactory<BigDecimal> {
  33.     /** Settings for enhanced precision computations. */
  34.     private final MathContext mContext;
  35.     /** The number {@code 2}. */
  36.     private final BigDecimal two;
  37.     /** The number {@code -1}. */
  38.     private final BigDecimal minusOne;
  39.     /** The number {@code 0.5}. */
  40.     private final BigDecimal oneHalf;

  41.     /**
  42.      * Default precision is {@link MathContext#DECIMAL128 DECIMAL128}.
  43.      */
  44.     public LegendreHighPrecisionRuleFactory() {
  45.         this(MathContext.DECIMAL128);
  46.     }

  47.     /**
  48.      * @param mContext Precision setting for computing the quadrature rules.
  49.      */
  50.     public LegendreHighPrecisionRuleFactory(MathContext mContext) {
  51.         this.mContext = mContext;
  52.         two = new BigDecimal("2", mContext);
  53.         minusOne = new BigDecimal("-1", mContext);
  54.         oneHalf = new BigDecimal("0.5", mContext);
  55.     }

  56.     /** {@inheritDoc} */
  57.     @Override
  58.     protected Pair<BigDecimal[], BigDecimal[]> computeRule(int numberOfPoints)
  59.         throws DimensionMismatchException {

  60.         if (numberOfPoints == 1) {
  61.             // Break recursion.
  62.             return new Pair<BigDecimal[], BigDecimal[]>(new BigDecimal[] { BigDecimal.ZERO },
  63.                                                         new BigDecimal[] { two });
  64.         }

  65.         // Get previous rule.
  66.         // If it has not been computed yet it will trigger a recursive call
  67.         // to this method.
  68.         final BigDecimal[] previousPoints = getRuleInternal(numberOfPoints - 1).getFirst();

  69.         // Compute next rule.
  70.         final BigDecimal[] points = new BigDecimal[numberOfPoints];
  71.         final BigDecimal[] weights = new BigDecimal[numberOfPoints];

  72.         // Find i-th root of P[n+1] by bracketing.
  73.         final int iMax = numberOfPoints / 2;
  74.         for (int i = 0; i < iMax; i++) {
  75.             // Lower-bound of the interval.
  76.             BigDecimal a = (i == 0) ? minusOne : previousPoints[i - 1];
  77.             // Upper-bound of the interval.
  78.             BigDecimal b = (iMax == 1) ? BigDecimal.ONE : previousPoints[i];
  79.             // P[j-1](a)
  80.             BigDecimal pma = BigDecimal.ONE;
  81.             // P[j](a)
  82.             BigDecimal pa = a;
  83.             // P[j-1](b)
  84.             BigDecimal pmb = BigDecimal.ONE;
  85.             // P[j](b)
  86.             BigDecimal pb = b;
  87.             for (int j = 1; j < numberOfPoints; j++) {
  88.                 final BigDecimal b_two_j_p_1 = new BigDecimal(2 * j + 1, mContext);
  89.                 final BigDecimal b_j = new BigDecimal(j, mContext);
  90.                 final BigDecimal b_j_p_1 = new BigDecimal(j + 1, mContext);

  91.                 // Compute P[j+1](a)
  92.                 // ppa = ((2 * j + 1) * a * pa - j * pma) / (j + 1);

  93.                 BigDecimal tmp1 = a.multiply(b_two_j_p_1, mContext);
  94.                 tmp1 = pa.multiply(tmp1, mContext);
  95.                 BigDecimal tmp2 = pma.multiply(b_j, mContext);
  96.                 // P[j+1](a)
  97.                 BigDecimal ppa = tmp1.subtract(tmp2, mContext);
  98.                 ppa = ppa.divide(b_j_p_1, mContext);

  99.                 // Compute P[j+1](b)
  100.                 // ppb = ((2 * j + 1) * b * pb - j * pmb) / (j + 1);

  101.                 tmp1 = b.multiply(b_two_j_p_1, mContext);
  102.                 tmp1 = pb.multiply(tmp1, mContext);
  103.                 tmp2 = pmb.multiply(b_j, mContext);
  104.                 // P[j+1](b)
  105.                 BigDecimal ppb = tmp1.subtract(tmp2, mContext);
  106.                 ppb = ppb.divide(b_j_p_1, mContext);

  107.                 pma = pa;
  108.                 pa = ppa;
  109.                 pmb = pb;
  110.                 pb = ppb;
  111.             }
  112.             // Now pa = P[n+1](a), and pma = P[n](a). Same holds for b.
  113.             // Middle of the interval.
  114.             BigDecimal c = a.add(b, mContext).multiply(oneHalf, mContext);
  115.             // P[j-1](c)
  116.             BigDecimal pmc = BigDecimal.ONE;
  117.             // P[j](c)
  118.             BigDecimal pc = c;
  119.             boolean done = false;
  120.             while (!done) {
  121.                 BigDecimal tmp1 = b.subtract(a, mContext);
  122.                 BigDecimal tmp2 = c.ulp().multiply(BigDecimal.TEN, mContext);
  123.                 done = tmp1.compareTo(tmp2) <= 0;
  124.                 pmc = BigDecimal.ONE;
  125.                 pc = c;
  126.                 for (int j = 1; j < numberOfPoints; j++) {
  127.                     final BigDecimal b_two_j_p_1 = new BigDecimal(2 * j + 1, mContext);
  128.                     final BigDecimal b_j = new BigDecimal(j, mContext);
  129.                     final BigDecimal b_j_p_1 = new BigDecimal(j + 1, mContext);

  130.                     // Compute P[j+1](c)
  131.                     tmp1 = c.multiply(b_two_j_p_1, mContext);
  132.                     tmp1 = pc.multiply(tmp1, mContext);
  133.                     tmp2 = pmc.multiply(b_j, mContext);
  134.                     // P[j+1](c)
  135.                     BigDecimal ppc = tmp1.subtract(tmp2, mContext);
  136.                     ppc = ppc.divide(b_j_p_1, mContext);

  137.                     pmc = pc;
  138.                     pc = ppc;
  139.                 }
  140.                 // Now pc = P[n+1](c) and pmc = P[n](c).
  141.                 if (!done) {
  142.                     if (pa.signum() * pc.signum() <= 0) {
  143.                         b = c;
  144.                         pmb = pmc;
  145.                         pb = pc;
  146.                     } else {
  147.                         a = c;
  148.                         pma = pmc;
  149.                         pa = pc;
  150.                     }
  151.                     c = a.add(b, mContext).multiply(oneHalf, mContext);
  152.                 }
  153.             }
  154.             final BigDecimal nP = new BigDecimal(numberOfPoints, mContext);
  155.             BigDecimal tmp1 = pmc.subtract(c.multiply(pc, mContext), mContext);
  156.             tmp1 = tmp1.multiply(nP);
  157.             tmp1 = tmp1.pow(2, mContext);
  158.             BigDecimal tmp2 = c.pow(2, mContext);
  159.             tmp2 = BigDecimal.ONE.subtract(tmp2, mContext);
  160.             tmp2 = tmp2.multiply(two, mContext);
  161.             tmp2 = tmp2.divide(tmp1, mContext);

  162.             points[i] = c;
  163.             weights[i] = tmp2;

  164.             final int idx = numberOfPoints - i - 1;
  165.             points[idx] = c.negate(mContext);
  166.             weights[idx] = tmp2;
  167.         }
  168.         // If "numberOfPoints" is odd, 0 is a root.
  169.         // Note: as written, the test for oddness will work for negative
  170.         // integers too (although it is not necessary here), preventing
  171.         // a FindBugs warning.
  172.         if (numberOfPoints % 2 != 0) {
  173.             BigDecimal pmc = BigDecimal.ONE;
  174.             for (int j = 1; j < numberOfPoints; j += 2) {
  175.                 final BigDecimal b_j = new BigDecimal(j, mContext);
  176.                 final BigDecimal b_j_p_1 = new BigDecimal(j + 1, mContext);

  177.                 // pmc = -j * pmc / (j + 1);
  178.                 pmc = pmc.multiply(b_j, mContext);
  179.                 pmc = pmc.divide(b_j_p_1, mContext);
  180.                 pmc = pmc.negate(mContext);
  181.             }

  182.             // 2 / pow(numberOfPoints * pmc, 2);
  183.             final BigDecimal nP = new BigDecimal(numberOfPoints, mContext);
  184.             BigDecimal tmp1 = pmc.multiply(nP, mContext);
  185.             tmp1 = tmp1.pow(2, mContext);
  186.             BigDecimal tmp2 = two.divide(tmp1, mContext);

  187.             points[iMax] = BigDecimal.ZERO;
  188.             weights[iMax] = tmp2;
  189.         }

  190.         return new Pair<BigDecimal[], BigDecimal[]>(points, weights);
  191.     }
  192. }