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