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.analysis.integration.gauss;
018
019import org.apache.commons.math4.util.Pair;
020import java.math.BigDecimal;
021
022/**
023 * Class that provides different ways to compute the nodes and weights to be
024 * used by the {@link GaussIntegrator Gaussian integration rule}.
025 *
026 * @since 3.1
027 */
028public class GaussIntegratorFactory {
029    /** Generator of Gauss-Legendre integrators. */
030    private final BaseRuleFactory<Double> legendre = new LegendreRuleFactory();
031    /** Generator of Gauss-Legendre integrators. */
032    private final BaseRuleFactory<BigDecimal> legendreHighPrecision = new LegendreHighPrecisionRuleFactory();
033    /** Generator of Gauss-Hermite integrators. */
034    private final BaseRuleFactory<Double> hermite = new HermiteRuleFactory();
035    /** Generator of Gauss-Laguerre integrators. */
036    private final BaseRuleFactory<Double> laguerre = new LaguerreRuleFactory();
037
038    /**
039     * Creates a Gauss-Laguerre integrator of the given order.
040     * The call to the
041     * {@link GaussIntegrator#integrate(org.apache.commons.math4.analysis.UnivariateFunction)
042     * integrate} method will perform an integration on the interval
043     * \([0, +\infty)\): the computed value is the improper integral of
044     * \(e^{-x} f(x)\)
045     * where \(f(x)\) is the function passed to the
046     * {@link SymmetricGaussIntegrator#integrate(org.apache.commons.math4.analysis.UnivariateFunction)
047     * integrate} method.
048     *
049     * @param numberOfPoints Order of the integration rule.
050     * @return a Gauss-Legendre integrator.
051     * @since 4.0
052     */
053    public GaussIntegrator laguerre(int numberOfPoints) {
054        return new GaussIntegrator(getRule(laguerre, numberOfPoints));
055    }
056
057    /**
058     * Creates a Gauss-Legendre integrator of the given order.
059     * The call to the
060     * {@link GaussIntegrator#integrate(org.apache.commons.math4.analysis.UnivariateFunction)
061     * integrate} method will perform an integration on the natural interval
062     * {@code [-1 , 1]}.
063     *
064     * @param numberOfPoints Order of the integration rule.
065     * @return a Gauss-Legendre integrator.
066     */
067    public GaussIntegrator legendre(int numberOfPoints) {
068        return new GaussIntegrator(getRule(legendre, numberOfPoints));
069    }
070
071    /**
072     * Creates a Gauss-Legendre integrator of the given order.
073     * The call to the
074     * {@link GaussIntegrator#integrate(org.apache.commons.math4.analysis.UnivariateFunction)
075     * integrate} method will perform an integration on the given interval.
076     *
077     * @param numberOfPoints Order of the integration rule.
078     * @param lowerBound Lower bound of the integration interval.
079     * @param upperBound Upper bound of the integration interval.
080     * @return a Gauss-Legendre integrator.
081     * @throws org.apache.commons.math4.exception.NotStrictlyPositiveException if number of points is not positive
082     */
083    public GaussIntegrator legendre(int numberOfPoints,
084                                    double lowerBound,
085                                    double upperBound) {
086        return new GaussIntegrator(transform(getRule(legendre, numberOfPoints),
087                                             lowerBound, upperBound));
088    }
089
090    /**
091     * Creates a Gauss-Legendre integrator of the given order.
092     * The call to the
093     * {@link GaussIntegrator#integrate(org.apache.commons.math4.analysis.UnivariateFunction)
094     * integrate} method will perform an integration on the natural interval
095     * {@code [-1 , 1]}.
096     *
097     * @param numberOfPoints Order of the integration rule.
098     * @return a Gauss-Legendre integrator.
099     * @throws org.apache.commons.math4.exception.NotStrictlyPositiveException if number of points is not positive
100     */
101    public GaussIntegrator legendreHighPrecision(int numberOfPoints) {
102        return new GaussIntegrator(getRule(legendreHighPrecision, numberOfPoints));
103    }
104
105    /**
106     * Creates an integrator of the given order, and whose call to the
107     * {@link GaussIntegrator#integrate(org.apache.commons.math4.analysis.UnivariateFunction)
108     * integrate} method will perform an integration on the given interval.
109     *
110     * @param numberOfPoints Order of the integration rule.
111     * @param lowerBound Lower bound of the integration interval.
112     * @param upperBound Upper bound of the integration interval.
113     * @return a Gauss-Legendre integrator.
114     * @throws org.apache.commons.math4.exception.NotStrictlyPositiveException if number of points is not positive
115     */
116    public GaussIntegrator legendreHighPrecision(int numberOfPoints,
117                                                 double lowerBound,
118                                                 double upperBound) {
119        return new GaussIntegrator(transform(getRule(legendreHighPrecision, numberOfPoints),
120                                             lowerBound, upperBound));
121    }
122
123    /**
124     * Creates a Gauss-Hermite integrator of the given order.
125     * The call to the
126     * {@link SymmetricGaussIntegrator#integrate(org.apache.commons.math4.analysis.UnivariateFunction)
127     * integrate} method will perform a weighted integration on the interval
128     * \([-\infty, +\infty]\): the computed value is the improper integral of
129     * \(e^{-x^2}f(x)\)
130     * where \(f(x)\) is the function passed to the
131     * {@link SymmetricGaussIntegrator#integrate(org.apache.commons.math4.analysis.UnivariateFunction)
132     * integrate} method.
133     *
134     * @param numberOfPoints Order of the integration rule.
135     * @return a Gauss-Hermite integrator.
136     */
137    public SymmetricGaussIntegrator hermite(int numberOfPoints) {
138        return new SymmetricGaussIntegrator(getRule(hermite, numberOfPoints));
139    }
140
141    /**
142     * @param factory Integration rule factory.
143     * @param numberOfPoints Order of the integration rule.
144     * @return the integration nodes and weights.
145     * @throws org.apache.commons.math4.exception.NotStrictlyPositiveException if number of points is not positive
146     * @throws org.apache.commons.math4.exception.DimensionMismatchException if the elements of the rule pair do not
147     * have the same length.
148     */
149    private static Pair<double[], double[]> getRule(BaseRuleFactory<? extends Number> factory,
150                                                    int numberOfPoints) {
151        return factory.getRule(numberOfPoints);
152    }
153
154    /**
155     * Performs a change of variable so that the integration can be performed
156     * on an arbitrary interval {@code [a, b]}.
157     * It is assumed that the natural interval is {@code [-1, 1]}.
158     *
159     * @param rule Original points and weights.
160     * @param a Lower bound of the integration interval.
161     * @param b Lower bound of the integration interval.
162     * @return the points and weights adapted to the new interval.
163     */
164    private static Pair<double[], double[]> transform(Pair<double[], double[]> rule,
165                                                      double a,
166                                                      double b) {
167        final double[] points = rule.getFirst();
168        final double[] weights = rule.getSecond();
169
170        // Scaling
171        final double scale = (b - a) / 2;
172        final double shift = a + scale;
173
174        for (int i = 0; i < points.length; i++) {
175            points[i] = points[i] * scale + shift;
176            weights[i] *= scale;
177        }
178
179        return new Pair<>(points, weights);
180    }
181}