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    }