GaussIntegratorFactory.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.math4.legacy.analysis.integration.gauss;

  18. import org.apache.commons.math4.legacy.core.Pair;
  19. import java.math.BigDecimal;

  20. /**
  21.  * Class that provides different ways to compute the nodes and weights to be
  22.  * used by the {@link GaussIntegrator Gaussian integration rule}.
  23.  *
  24.  * @since 3.1
  25.  */
  26. public class GaussIntegratorFactory {
  27.     /** Generator of Gauss-Legendre integrators. */
  28.     private final BaseRuleFactory<Double> legendre = new LegendreRuleFactory();
  29.     /** Generator of Gauss-Legendre integrators. */
  30.     private final BaseRuleFactory<BigDecimal> legendreHighPrecision = new LegendreHighPrecisionRuleFactory();
  31.     /** Generator of Gauss-Hermite integrators. */
  32.     private final BaseRuleFactory<Double> hermite = new HermiteRuleFactory();
  33.     /** Generator of Gauss-Laguerre integrators. */
  34.     private final BaseRuleFactory<Double> laguerre = new LaguerreRuleFactory();

  35.     /**
  36.      * Creates a Gauss-Laguerre integrator of the given order.
  37.      * The call to the
  38.      * {@link GaussIntegrator#integrate(org.apache.commons.math4.legacy.analysis.UnivariateFunction)
  39.      * integrate} method will perform an integration on the interval
  40.      * \([0, +\infty)\): the computed value is the improper integral of
  41.      * \(e^{-x} f(x)\)
  42.      * where \(f(x)\) is the function passed to the
  43.      * {@link SymmetricGaussIntegrator#integrate(org.apache.commons.math4.legacy.analysis.UnivariateFunction)
  44.      * integrate} method.
  45.      *
  46.      * @param numberOfPoints Order of the integration rule.
  47.      * @return a Gauss-Legendre integrator.
  48.      * @since 4.0
  49.      */
  50.     public GaussIntegrator laguerre(int numberOfPoints) {
  51.         return new GaussIntegrator(getRule(laguerre, numberOfPoints));
  52.     }

  53.     /**
  54.      * Creates a Gauss-Legendre integrator of the given order.
  55.      * The call to the
  56.      * {@link GaussIntegrator#integrate(org.apache.commons.math4.legacy.analysis.UnivariateFunction)
  57.      * integrate} method will perform an integration on the natural interval
  58.      * {@code [-1 , 1]}.
  59.      *
  60.      * @param numberOfPoints Order of the integration rule.
  61.      * @return a Gauss-Legendre integrator.
  62.      */
  63.     public GaussIntegrator legendre(int numberOfPoints) {
  64.         return new GaussIntegrator(getRule(legendre, numberOfPoints));
  65.     }

  66.     /**
  67.      * Creates a Gauss-Legendre integrator of the given order.
  68.      * The call to the
  69.      * {@link GaussIntegrator#integrate(org.apache.commons.math4.legacy.analysis.UnivariateFunction)
  70.      * integrate} method will perform an integration on the given interval.
  71.      *
  72.      * @param numberOfPoints Order of the integration rule.
  73.      * @param lowerBound Lower bound of the integration interval.
  74.      * @param upperBound Upper bound of the integration interval.
  75.      * @return a Gauss-Legendre integrator.
  76.      * @throws org.apache.commons.math4.legacy.exception.NotStrictlyPositiveException if number of points is not positive
  77.      */
  78.     public GaussIntegrator legendre(int numberOfPoints,
  79.                                     double lowerBound,
  80.                                     double upperBound) {
  81.         return new GaussIntegrator(transform(getRule(legendre, numberOfPoints),
  82.                                              lowerBound, upperBound));
  83.     }

  84.     /**
  85.      * Creates a Gauss-Legendre integrator of the given order.
  86.      * The call to the
  87.      * {@link GaussIntegrator#integrate(org.apache.commons.math4.legacy.analysis.UnivariateFunction)
  88.      * integrate} method will perform an integration on the natural interval
  89.      * {@code [-1 , 1]}.
  90.      *
  91.      * @param numberOfPoints Order of the integration rule.
  92.      * @return a Gauss-Legendre integrator.
  93.      * @throws org.apache.commons.math4.legacy.exception.NotStrictlyPositiveException if number of points is not positive
  94.      */
  95.     public GaussIntegrator legendreHighPrecision(int numberOfPoints) {
  96.         return new GaussIntegrator(getRule(legendreHighPrecision, numberOfPoints));
  97.     }

  98.     /**
  99.      * Creates an integrator of the given order, and whose call to the
  100.      * {@link GaussIntegrator#integrate(org.apache.commons.math4.legacy.analysis.UnivariateFunction)
  101.      * integrate} method will perform an integration on the given interval.
  102.      *
  103.      * @param numberOfPoints Order of the integration rule.
  104.      * @param lowerBound Lower bound of the integration interval.
  105.      * @param upperBound Upper bound of the integration interval.
  106.      * @return a Gauss-Legendre integrator.
  107.      * @throws org.apache.commons.math4.legacy.exception.NotStrictlyPositiveException if number of points is not positive
  108.      */
  109.     public GaussIntegrator legendreHighPrecision(int numberOfPoints,
  110.                                                  double lowerBound,
  111.                                                  double upperBound) {
  112.         return new GaussIntegrator(transform(getRule(legendreHighPrecision, numberOfPoints),
  113.                                              lowerBound, upperBound));
  114.     }

  115.     /**
  116.      * Creates a Gauss-Hermite integrator of the given order.
  117.      * The call to the
  118.      * {@link SymmetricGaussIntegrator#integrate(org.apache.commons.math4.legacy.analysis.UnivariateFunction)
  119.      * integrate} method will perform a weighted integration on the interval
  120.      * \([-\infty, +\infty]\): the computed value is the improper integral of
  121.      * \(e^{-x^2}f(x)\)
  122.      * where \(f(x)\) is the function passed to the
  123.      * {@link SymmetricGaussIntegrator#integrate(org.apache.commons.math4.legacy.analysis.UnivariateFunction)
  124.      * integrate} method.
  125.      *
  126.      * @param numberOfPoints Order of the integration rule.
  127.      * @return a Gauss-Hermite integrator.
  128.      */
  129.     public SymmetricGaussIntegrator hermite(int numberOfPoints) {
  130.         return new SymmetricGaussIntegrator(getRule(hermite, numberOfPoints));
  131.     }

  132.     /**
  133.      * @param factory Integration rule factory.
  134.      * @param numberOfPoints Order of the integration rule.
  135.      * @return the integration nodes and weights.
  136.      * @throws org.apache.commons.math4.legacy.exception.NotStrictlyPositiveException if number of points is not positive
  137.      * @throws org.apache.commons.math4.legacy.exception.DimensionMismatchException if the elements of the rule pair do not
  138.      * have the same length.
  139.      */
  140.     private static Pair<double[], double[]> getRule(BaseRuleFactory<? extends Number> factory,
  141.                                                     int numberOfPoints) {
  142.         return factory.getRule(numberOfPoints);
  143.     }

  144.     /**
  145.      * Performs a change of variable so that the integration can be performed
  146.      * on an arbitrary interval {@code [a, b]}.
  147.      * It is assumed that the natural interval is {@code [-1, 1]}.
  148.      *
  149.      * @param rule Original points and weights.
  150.      * @param a Lower bound of the integration interval.
  151.      * @param b Lower bound of the integration interval.
  152.      * @return the points and weights adapted to the new interval.
  153.      */
  154.     private static Pair<double[], double[]> transform(Pair<double[], double[]> rule,
  155.                                                       double a,
  156.                                                       double b) {
  157.         final double[] points = rule.getFirst();
  158.         final double[] weights = rule.getSecond();

  159.         // Scaling
  160.         final double scale = (b - a) / 2;
  161.         final double shift = a + scale;

  162.         for (int i = 0; i < points.length; i++) {
  163.             points[i] = points[i] * scale + shift;
  164.             weights[i] *= scale;
  165.         }

  166.         return new Pair<>(points, weights);
  167.     }
  168. }