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: 799857 $ $Date: 2009-08-01 09:07:12 -0400 (Sat, 01 Aug 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 int i = 1, j, m = maximalIterationCount + 1; 071 // Array structure here can be improved for better space 072 // efficiency because only the lower triangle is used. 073 double r, t[][] = new double[m][m], s, olds; 074 075 clearResult(); 076 verifyInterval(min, max); 077 verifyIterationCount(); 078 079 TrapezoidIntegrator qtrap = new TrapezoidIntegrator(); 080 t[0][0] = qtrap.stage(f, min, max, 0); 081 olds = t[0][0]; 082 while (i <= maximalIterationCount) { 083 t[i][0] = qtrap.stage(f, min, max, i); 084 for (j = 1; j <= i; j++) { 085 // Richardson extrapolation coefficient 086 r = (1L << (2 * j)) -1; 087 t[i][j] = t[i][j-1] + (t[i][j-1] - t[i-1][j-1]) / r; 088 } 089 s = t[i][i]; 090 if (i >= minimalIterationCount) { 091 final double delta = Math.abs(s - olds); 092 final double rLimit = 093 relativeAccuracy * (Math.abs(olds) + Math.abs(s)) * 0.5; 094 if ((delta <= rLimit) || (delta <= absoluteAccuracy)) { 095 setResult(s, i); 096 return result; 097 } 098 } 099 olds = s; 100 i++; 101 } 102 throw new MaxIterationsExceededException(maximalIterationCount); 103 } 104 105 /** {@inheritDoc} */ 106 @Override 107 protected void verifyIterationCount() throws IllegalArgumentException { 108 super.verifyIterationCount(); 109 // at most 32 bisection refinements due to higher order divider 110 if (maximalIterationCount > 32) { 111 throw MathRuntimeException.createIllegalArgumentException( 112 "invalid iteration limits: min={0}, max={1}", 113 0, 32); 114 } 115 } 116 }