    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     */
    017    package org.apache.commons.math3.analysis.integration.gauss;
    018
    019    import java.math.BigDecimal;
    020
    021    import org.apache.commons.math3.util.Pair;
    022
    023    /**
    024     * Class that provides different ways to compute the nodes and weights to be
    025     * used by the {@link GaussIntegrator Gaussian integration rule}.
    026     *
    027     * @since 3.1
    028     * @version $Id: GaussIntegratorFactory.java 1364420 2012-07-22 20:01:12Z tn $
    029     */
    030    public class GaussIntegratorFactory {
    031        /** Generator of Gauss-Legendre integrators. */
    032        private final BaseRuleFactory<Double> legendre = new LegendreRuleFactory();
    033        /** Generator of Gauss-Legendre integrators. */
    034        private final BaseRuleFactory<BigDecimal> legendreHighPrecision = new LegendreHighPrecisionRuleFactory();
    035
    036        /**
    037         * Creates an integrator of the given order, and whose call to the
    038         * {@link GaussIntegrator#integrate(org.apache.commons.math3.analysis.UnivariateFunction)
    039         * integrate} method will perform an integration on the natural interval
    040         * {@code [-1 , 1]}.
    041         *
    042         * @param numberOfPoints Order of the integration rule.
    043         * @return a Gauss-Legendre integrator.
    044         */
    045        public GaussIntegrator legendre(int numberOfPoints) {
    046            return new GaussIntegrator(getRule(legendre, numberOfPoints));
    047        }
    048
    049        /**
    050         * Creates an integrator of the given order, and whose call to the
    051         * {@link GaussIntegrator#integrate(org.apache.commons.math3.analysis.UnivariateFunction)
    052         * integrate} method will perform an integration on the given interval.
    053         *
    054         * @param numberOfPoints Order of the integration rule.
    055         * @param lowerBound Lower bound of the integration interval.
    056         * @param upperBound Upper bound of the integration interval.
    057         * @return a Gauss-Legendre integrator.
    058         */
    059        public GaussIntegrator legendre(int numberOfPoints,
    060                                        double lowerBound,
    061                                        double upperBound) {
    062            return new GaussIntegrator(transform(getRule(legendre, numberOfPoints),
    063                                                 lowerBound, upperBound));
    064        }
    065
    066        /**
    067         * Creates an integrator of the given order, and whose call to the
    068         * {@link GaussIntegrator#integrate(org.apache.commons.math3.analysis.UnivariateFunction)
    069         * integrate} method will perform an integration on the natural interval
    070         * {@code [-1 , 1]}.
    071         *
    072         * @param numberOfPoints Order of the integration rule.
    073         * @return a Gauss-Legendre integrator.
    074         */
    075        public GaussIntegrator legendreHighPrecision(int numberOfPoints) {
    076            return new GaussIntegrator(getRule(legendreHighPrecision, numberOfPoints));
    077        }
    078
    079        /**
    080         * Creates an integrator of the given order, and whose call to the
    081         * {@link GaussIntegrator#integrate(org.apache.commons.math3.analysis.UnivariateFunction)
    082         * integrate} method will perform an integration on the given interval.
    083         *
    084         * @param numberOfPoints Order of the integration rule.
    085         * @param lowerBound Lower bound of the integration interval.
    086         * @param upperBound Upper bound of the integration interval.
    087         * @return a Gauss-Legendre integrator.
    088         */
    089        public GaussIntegrator legendreHighPrecision(int numberOfPoints,
    090                                                     double lowerBound,
    091                                                     double upperBound) {
    092            return new GaussIntegrator(transform(getRule(legendreHighPrecision, numberOfPoints),
    093                                                 lowerBound, upperBound));
    094        }
    095
    096        /**
    097         * @param factory Integration rule factory.
    098         * @param numberOfPoints Order of the integration rule.
    099         * @return the integration nodes and weights.
    100         */
    101        private static Pair<double[], double[]> getRule(BaseRuleFactory<? extends Number> factory,
    102                                                        int numberOfPoints) {
    103            return factory.getRule(numberOfPoints);
    104        }
    105
    106        /**
    107         * Performs a change of variable so that the integration can be performed
    108         * on an arbitrary interval {@code [a, b]}.
    109         * It is assumed that the natural interval is {@code [-1, 1]}.
    110         *
    111         * @param rule Original points and weights.
    112         * @param a Lower bound of the integration interval.
    113         * @param b Lower bound of the integration interval.
    114         * @return the points and weights adapted to the new interval.
    115         */
    116        private static Pair<double[], double[]> transform(Pair<double[], double[]> rule,
    117                                                          double a,
    118                                                          double b) {
    119            final double[] points = rule.getFirst();
    120            final double[] weights = rule.getSecond();
    121
    122            // Scaling
    123            final double scale = (b - a) / 2;
    124            final double shift = a + scale;
    125
    126            for (int i = 0; i < points.length; i++) {
    127                points[i] = points[i] * scale + shift;
    128                weights[i] *= scale;
    129            }
    130
    131            return new Pair<double[], double[]>(points, weights);
    132        }
    133    }