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.math.analysis.integration; 018 019 import org.apache.commons.math.FunctionEvaluationException; 020 import org.apache.commons.math.MathRuntimeException; 021 import org.apache.commons.math.MaxIterationsExceededException; 022 import org.apache.commons.math.analysis.UnivariateRealFunction; 023 024 /** 025 * Implements the <a href="http://mathworld.wolfram.com/RombergIntegration.html"> 026 * Romberg Algorithm</a> for integration of real univariate functions. For 027 * reference, see <b>Introduction to Numerical Analysis</b>, ISBN 038795452X, 028 * chapter 3. 029 * <p> 030 * Romberg integration employs k successive refinements of the trapezoid 031 * rule to remove error terms less than order O(N^(-2k)). Simpson's rule 032 * is a special case of k = 2.</p> 033 * 034 * @version $Revision: 824822 $ $Date: 2009-10-13 11:56:51 -0400 (Tue, 13 Oct 2009) $ 035 * @since 1.2 036 */ 037 public class RombergIntegrator extends UnivariateRealIntegratorImpl { 038 039 /** 040 * Construct an integrator for the given function. 041 * 042 * @param f function to integrate 043 * @deprecated as of 2.0 the integrand function is passed as an argument 044 * to the {@link #integrate(UnivariateRealFunction, double, double)}method. 045 */ 046 @Deprecated 047 public RombergIntegrator(UnivariateRealFunction f) { 048 super(f, 32); 049 } 050 051 /** 052 * Construct an integrator. 053 */ 054 public RombergIntegrator() { 055 super(32); 056 } 057 058 /** {@inheritDoc} */ 059 @Deprecated 060 public double integrate(final double min, final double max) 061 throws MaxIterationsExceededException, FunctionEvaluationException, IllegalArgumentException { 062 return integrate(f, min, max); 063 } 064 065 /** {@inheritDoc} */ 066 public double integrate(final UnivariateRealFunction f, 067 final double min, final double max) 068 throws MaxIterationsExceededException, FunctionEvaluationException, IllegalArgumentException { 069 070 final int m = maximalIterationCount + 1; 071 double previousRow[] = new double[m]; 072 double currentRow[] = new double[m]; 073 074 clearResult(); 075 verifyInterval(min, max); 076 verifyIterationCount(); 077 078 TrapezoidIntegrator qtrap = new TrapezoidIntegrator(); 079 currentRow[0] = qtrap.stage(f, min, max, 0); 080 double olds = currentRow[0]; 081 for (int i = 1; i <= maximalIterationCount; ++i) { 082 083 // switch rows 084 final double[] tmpRow = previousRow; 085 previousRow = currentRow; 086 currentRow = tmpRow; 087 088 currentRow[0] = qtrap.stage(f, min, max, i); 089 for (int j = 1; j <= i; j++) { 090 // Richardson extrapolation coefficient 091 final double r = (1L << (2 * j)) - 1; 092 final double tIJm1 = currentRow[j - 1]; 093 currentRow[j] = tIJm1 + (tIJm1 - previousRow[j - 1]) / r; 094 } 095 final double s = currentRow[i]; 096 if (i >= minimalIterationCount) { 097 final double delta = Math.abs(s - olds); 098 final double rLimit = relativeAccuracy * (Math.abs(olds) + Math.abs(s)) * 0.5; 099 if ((delta <= rLimit) || (delta <= absoluteAccuracy)) { 100 setResult(s, i); 101 return result; 102 } 103 } 104 olds = s; 105 } 106 throw new MaxIterationsExceededException(maximalIterationCount); 107 } 108 109 /** {@inheritDoc} */ 110 @Override 111 protected void verifyIterationCount() throws IllegalArgumentException { 112 super.verifyIterationCount(); 113 // at most 32 bisection refinements due to higher order divider 114 if (maximalIterationCount > 32) { 115 throw MathRuntimeException.createIllegalArgumentException( 116 "invalid iteration limits: min={0}, max={1}", 117 0, 32); 118 } 119 } 120 }