Raymond | dee0849 | 2015-04-02 10:43:13 -0700 | [diff] [blame] | 1 | /* |
| 2 | * Licensed to the Apache Software Foundation (ASF) under one or more |
| 3 | * contributor license agreements. See the NOTICE file distributed with |
| 4 | * this work for additional information regarding copyright ownership. |
| 5 | * The ASF licenses this file to You under the Apache License, Version 2.0 |
| 6 | * (the "License"); you may not use this file except in compliance with |
| 7 | * the License. You may obtain a copy of the License at |
| 8 | * |
| 9 | * http://www.apache.org/licenses/LICENSE-2.0 |
| 10 | * |
| 11 | * Unless required by applicable law or agreed to in writing, software |
| 12 | * distributed under the License is distributed on an "AS IS" BASIS, |
| 13 | * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. |
| 14 | * See the License for the specific language governing permissions and |
| 15 | * limitations under the License. |
| 16 | */ |
| 17 | package org.apache.commons.math.analysis.integration; |
| 18 | |
| 19 | import org.apache.commons.math.ConvergenceException; |
| 20 | import org.apache.commons.math.FunctionEvaluationException; |
| 21 | import org.apache.commons.math.MathRuntimeException; |
| 22 | import org.apache.commons.math.MaxIterationsExceededException; |
| 23 | import org.apache.commons.math.analysis.UnivariateRealFunction; |
| 24 | import org.apache.commons.math.exception.util.LocalizedFormats; |
| 25 | import org.apache.commons.math.util.FastMath; |
| 26 | |
| 27 | /** |
| 28 | * Implements the <a href="http://mathworld.wolfram.com/Legendre-GaussQuadrature.html"> |
| 29 | * Legendre-Gauss</a> quadrature formula. |
| 30 | * <p> |
| 31 | * Legendre-Gauss integrators are efficient integrators that can |
| 32 | * accurately integrate functions with few functions evaluations. A |
| 33 | * Legendre-Gauss integrator using an n-points quadrature formula can |
| 34 | * integrate exactly 2n-1 degree polynomials. |
| 35 | * </p> |
| 36 | * <p> |
| 37 | * These integrators evaluate the function on n carefully chosen |
| 38 | * abscissas in each step interval (mapped to the canonical [-1 1] interval). |
| 39 | * The evaluation abscissas are not evenly spaced and none of them are |
| 40 | * at the interval endpoints. This implies the function integrated can be |
| 41 | * undefined at integration interval endpoints. |
| 42 | * </p> |
| 43 | * <p> |
| 44 | * The evaluation abscissas x<sub>i</sub> are the roots of the degree n |
| 45 | * Legendre polynomial. The weights a<sub>i</sub> of the quadrature formula |
| 46 | * integrals from -1 to +1 ∫ Li<sup>2</sup> where Li (x) = |
| 47 | * ∏ (x-x<sub>k</sub>)/(x<sub>i</sub>-x<sub>k</sub>) for k != i. |
| 48 | * </p> |
| 49 | * <p> |
| 50 | * @version $Revision: 1070725 $ $Date: 2011-02-15 02:31:12 +0100 (mar. 15 févr. 2011) $ |
| 51 | * @since 1.2 |
| 52 | */ |
| 53 | |
| 54 | public class LegendreGaussIntegrator extends UnivariateRealIntegratorImpl { |
| 55 | |
| 56 | /** Abscissas for the 2 points method. */ |
| 57 | private static final double[] ABSCISSAS_2 = { |
| 58 | -1.0 / FastMath.sqrt(3.0), |
| 59 | 1.0 / FastMath.sqrt(3.0) |
| 60 | }; |
| 61 | |
| 62 | /** Weights for the 2 points method. */ |
| 63 | private static final double[] WEIGHTS_2 = { |
| 64 | 1.0, |
| 65 | 1.0 |
| 66 | }; |
| 67 | |
| 68 | /** Abscissas for the 3 points method. */ |
| 69 | private static final double[] ABSCISSAS_3 = { |
| 70 | -FastMath.sqrt(0.6), |
| 71 | 0.0, |
| 72 | FastMath.sqrt(0.6) |
| 73 | }; |
| 74 | |
| 75 | /** Weights for the 3 points method. */ |
| 76 | private static final double[] WEIGHTS_3 = { |
| 77 | 5.0 / 9.0, |
| 78 | 8.0 / 9.0, |
| 79 | 5.0 / 9.0 |
| 80 | }; |
| 81 | |
| 82 | /** Abscissas for the 4 points method. */ |
| 83 | private static final double[] ABSCISSAS_4 = { |
| 84 | -FastMath.sqrt((15.0 + 2.0 * FastMath.sqrt(30.0)) / 35.0), |
| 85 | -FastMath.sqrt((15.0 - 2.0 * FastMath.sqrt(30.0)) / 35.0), |
| 86 | FastMath.sqrt((15.0 - 2.0 * FastMath.sqrt(30.0)) / 35.0), |
| 87 | FastMath.sqrt((15.0 + 2.0 * FastMath.sqrt(30.0)) / 35.0) |
| 88 | }; |
| 89 | |
| 90 | /** Weights for the 4 points method. */ |
| 91 | private static final double[] WEIGHTS_4 = { |
| 92 | (90.0 - 5.0 * FastMath.sqrt(30.0)) / 180.0, |
| 93 | (90.0 + 5.0 * FastMath.sqrt(30.0)) / 180.0, |
| 94 | (90.0 + 5.0 * FastMath.sqrt(30.0)) / 180.0, |
| 95 | (90.0 - 5.0 * FastMath.sqrt(30.0)) / 180.0 |
| 96 | }; |
| 97 | |
| 98 | /** Abscissas for the 5 points method. */ |
| 99 | private static final double[] ABSCISSAS_5 = { |
| 100 | -FastMath.sqrt((35.0 + 2.0 * FastMath.sqrt(70.0)) / 63.0), |
| 101 | -FastMath.sqrt((35.0 - 2.0 * FastMath.sqrt(70.0)) / 63.0), |
| 102 | 0.0, |
| 103 | FastMath.sqrt((35.0 - 2.0 * FastMath.sqrt(70.0)) / 63.0), |
| 104 | FastMath.sqrt((35.0 + 2.0 * FastMath.sqrt(70.0)) / 63.0) |
| 105 | }; |
| 106 | |
| 107 | /** Weights for the 5 points method. */ |
| 108 | private static final double[] WEIGHTS_5 = { |
| 109 | (322.0 - 13.0 * FastMath.sqrt(70.0)) / 900.0, |
| 110 | (322.0 + 13.0 * FastMath.sqrt(70.0)) / 900.0, |
| 111 | 128.0 / 225.0, |
| 112 | (322.0 + 13.0 * FastMath.sqrt(70.0)) / 900.0, |
| 113 | (322.0 - 13.0 * FastMath.sqrt(70.0)) / 900.0 |
| 114 | }; |
| 115 | |
| 116 | /** Abscissas for the current method. */ |
| 117 | private final double[] abscissas; |
| 118 | |
| 119 | /** Weights for the current method. */ |
| 120 | private final double[] weights; |
| 121 | |
| 122 | /** |
| 123 | * Build a Legendre-Gauss integrator. |
| 124 | * @param n number of points desired (must be between 2 and 5 inclusive) |
| 125 | * @param defaultMaximalIterationCount maximum number of iterations |
| 126 | * @exception IllegalArgumentException if the number of points is not |
| 127 | * in the supported range |
| 128 | */ |
| 129 | public LegendreGaussIntegrator(final int n, final int defaultMaximalIterationCount) |
| 130 | throws IllegalArgumentException { |
| 131 | super(defaultMaximalIterationCount); |
| 132 | switch(n) { |
| 133 | case 2 : |
| 134 | abscissas = ABSCISSAS_2; |
| 135 | weights = WEIGHTS_2; |
| 136 | break; |
| 137 | case 3 : |
| 138 | abscissas = ABSCISSAS_3; |
| 139 | weights = WEIGHTS_3; |
| 140 | break; |
| 141 | case 4 : |
| 142 | abscissas = ABSCISSAS_4; |
| 143 | weights = WEIGHTS_4; |
| 144 | break; |
| 145 | case 5 : |
| 146 | abscissas = ABSCISSAS_5; |
| 147 | weights = WEIGHTS_5; |
| 148 | break; |
| 149 | default : |
| 150 | throw MathRuntimeException.createIllegalArgumentException( |
| 151 | LocalizedFormats.N_POINTS_GAUSS_LEGENDRE_INTEGRATOR_NOT_SUPPORTED, |
| 152 | n, 2, 5); |
| 153 | } |
| 154 | |
| 155 | } |
| 156 | |
| 157 | /** {@inheritDoc} */ |
| 158 | @Deprecated |
| 159 | public double integrate(final double min, final double max) |
| 160 | throws ConvergenceException, FunctionEvaluationException, IllegalArgumentException { |
| 161 | return integrate(f, min, max); |
| 162 | } |
| 163 | |
| 164 | /** {@inheritDoc} */ |
| 165 | public double integrate(final UnivariateRealFunction f, final double min, final double max) |
| 166 | throws ConvergenceException, FunctionEvaluationException, IllegalArgumentException { |
| 167 | |
| 168 | clearResult(); |
| 169 | verifyInterval(min, max); |
| 170 | verifyIterationCount(); |
| 171 | |
| 172 | // compute first estimate with a single step |
| 173 | double oldt = stage(f, min, max, 1); |
| 174 | |
| 175 | int n = 2; |
| 176 | for (int i = 0; i < maximalIterationCount; ++i) { |
| 177 | |
| 178 | // improve integral with a larger number of steps |
| 179 | final double t = stage(f, min, max, n); |
| 180 | |
| 181 | // estimate error |
| 182 | final double delta = FastMath.abs(t - oldt); |
| 183 | final double limit = |
| 184 | FastMath.max(absoluteAccuracy, |
| 185 | relativeAccuracy * (FastMath.abs(oldt) + FastMath.abs(t)) * 0.5); |
| 186 | |
| 187 | // check convergence |
| 188 | if ((i + 1 >= minimalIterationCount) && (delta <= limit)) { |
| 189 | setResult(t, i); |
| 190 | return result; |
| 191 | } |
| 192 | |
| 193 | // prepare next iteration |
| 194 | double ratio = FastMath.min(4, FastMath.pow(delta / limit, 0.5 / abscissas.length)); |
| 195 | n = FastMath.max((int) (ratio * n), n + 1); |
| 196 | oldt = t; |
| 197 | |
| 198 | } |
| 199 | |
| 200 | throw new MaxIterationsExceededException(maximalIterationCount); |
| 201 | |
| 202 | } |
| 203 | |
| 204 | /** |
| 205 | * Compute the n-th stage integral. |
| 206 | * @param f the integrand function |
| 207 | * @param min the lower bound for the interval |
| 208 | * @param max the upper bound for the interval |
| 209 | * @param n number of steps |
| 210 | * @return the value of n-th stage integral |
| 211 | * @throws FunctionEvaluationException if an error occurs evaluating the |
| 212 | * function |
| 213 | */ |
| 214 | private double stage(final UnivariateRealFunction f, |
| 215 | final double min, final double max, final int n) |
| 216 | throws FunctionEvaluationException { |
| 217 | |
| 218 | // set up the step for the current stage |
| 219 | final double step = (max - min) / n; |
| 220 | final double halfStep = step / 2.0; |
| 221 | |
| 222 | // integrate over all elementary steps |
| 223 | double midPoint = min + halfStep; |
| 224 | double sum = 0.0; |
| 225 | for (int i = 0; i < n; ++i) { |
| 226 | for (int j = 0; j < abscissas.length; ++j) { |
| 227 | sum += weights[j] * f.value(midPoint + halfStep * abscissas[j]); |
| 228 | } |
| 229 | midPoint += step; |
| 230 | } |
| 231 | |
| 232 | return halfStep * sum; |
| 233 | |
| 234 | } |
| 235 | |
| 236 | } |