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.legacy.analysis.integration.gauss; 018 019import org.apache.commons.math4.legacy.core.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.legacy.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.legacy.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.legacy.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.legacy.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.legacy.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.legacy.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.legacy.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.legacy.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.legacy.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.legacy.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.legacy.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.legacy.exception.NotStrictlyPositiveException if number of points is not positive 146 * @throws org.apache.commons.math4.legacy.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}