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.math3.analysis.integration; 018 019 import org.apache.commons.math3.exception.MaxCountExceededException; 020 import org.apache.commons.math3.exception.NotStrictlyPositiveException; 021 import org.apache.commons.math3.exception.NumberIsTooLargeException; 022 import org.apache.commons.math3.exception.NumberIsTooSmallException; 023 import org.apache.commons.math3.exception.TooManyEvaluationsException; 024 import org.apache.commons.math3.util.FastMath; 025 026 /** 027 * Implements the <a href="http://mathworld.wolfram.com/RombergIntegration.html"> 028 * Romberg Algorithm</a> for integration of real univariate functions. For 029 * reference, see <b>Introduction to Numerical Analysis</b>, ISBN 038795452X, 030 * chapter 3. 031 * <p> 032 * Romberg integration employs k successive refinements of the trapezoid 033 * rule to remove error terms less than order O(N^(-2k)). Simpson's rule 034 * is a special case of k = 2.</p> 035 * 036 * @version $Id: RombergIntegrator.java 1364387 2012-07-22 18:14:11Z tn $ 037 * @since 1.2 038 */ 039 public class RombergIntegrator extends BaseAbstractUnivariateIntegrator { 040 041 /** Maximal number of iterations for Romberg. */ 042 public static final int ROMBERG_MAX_ITERATIONS_COUNT = 32; 043 044 /** 045 * Build a Romberg integrator with given accuracies and iterations counts. 046 * @param relativeAccuracy relative accuracy of the result 047 * @param absoluteAccuracy absolute accuracy of the result 048 * @param minimalIterationCount minimum number of iterations 049 * @param maximalIterationCount maximum number of iterations 050 * (must be less than or equal to {@link #ROMBERG_MAX_ITERATIONS_COUNT}) 051 * @exception NotStrictlyPositiveException if minimal number of iterations 052 * is not strictly positive 053 * @exception NumberIsTooSmallException if maximal number of iterations 054 * is lesser than or equal to the minimal number of iterations 055 * @exception NumberIsTooLargeException if maximal number of iterations 056 * is greater than {@link #ROMBERG_MAX_ITERATIONS_COUNT} 057 */ 058 public RombergIntegrator(final double relativeAccuracy, 059 final double absoluteAccuracy, 060 final int minimalIterationCount, 061 final int maximalIterationCount) 062 throws NotStrictlyPositiveException, NumberIsTooSmallException, NumberIsTooLargeException { 063 super(relativeAccuracy, absoluteAccuracy, minimalIterationCount, maximalIterationCount); 064 if (maximalIterationCount > ROMBERG_MAX_ITERATIONS_COUNT) { 065 throw new NumberIsTooLargeException(maximalIterationCount, 066 ROMBERG_MAX_ITERATIONS_COUNT, false); 067 } 068 } 069 070 /** 071 * Build a Romberg integrator with given iteration counts. 072 * @param minimalIterationCount minimum number of iterations 073 * @param maximalIterationCount maximum number of iterations 074 * (must be less than or equal to {@link #ROMBERG_MAX_ITERATIONS_COUNT}) 075 * @exception NotStrictlyPositiveException if minimal number of iterations 076 * is not strictly positive 077 * @exception NumberIsTooSmallException if maximal number of iterations 078 * is lesser than or equal to the minimal number of iterations 079 * @exception NumberIsTooLargeException if maximal number of iterations 080 * is greater than {@link #ROMBERG_MAX_ITERATIONS_COUNT} 081 */ 082 public RombergIntegrator(final int minimalIterationCount, 083 final int maximalIterationCount) 084 throws NotStrictlyPositiveException, NumberIsTooSmallException, NumberIsTooLargeException { 085 super(minimalIterationCount, maximalIterationCount); 086 if (maximalIterationCount > ROMBERG_MAX_ITERATIONS_COUNT) { 087 throw new NumberIsTooLargeException(maximalIterationCount, 088 ROMBERG_MAX_ITERATIONS_COUNT, false); 089 } 090 } 091 092 /** 093 * Construct a Romberg integrator with default settings 094 * (max iteration count set to {@link #ROMBERG_MAX_ITERATIONS_COUNT}) 095 */ 096 public RombergIntegrator() { 097 super(DEFAULT_MIN_ITERATIONS_COUNT, ROMBERG_MAX_ITERATIONS_COUNT); 098 } 099 100 /** {@inheritDoc} */ 101 @Override 102 protected double doIntegrate() 103 throws TooManyEvaluationsException, MaxCountExceededException { 104 105 final int m = iterations.getMaximalCount() + 1; 106 double previousRow[] = new double[m]; 107 double currentRow[] = new double[m]; 108 109 TrapezoidIntegrator qtrap = new TrapezoidIntegrator(); 110 currentRow[0] = qtrap.stage(this, 0); 111 iterations.incrementCount(); 112 double olds = currentRow[0]; 113 while (true) { 114 115 final int i = iterations.getCount(); 116 117 // switch rows 118 final double[] tmpRow = previousRow; 119 previousRow = currentRow; 120 currentRow = tmpRow; 121 122 currentRow[0] = qtrap.stage(this, i); 123 iterations.incrementCount(); 124 for (int j = 1; j <= i; j++) { 125 // Richardson extrapolation coefficient 126 final double r = (1L << (2 * j)) - 1; 127 final double tIJm1 = currentRow[j - 1]; 128 currentRow[j] = tIJm1 + (tIJm1 - previousRow[j - 1]) / r; 129 } 130 final double s = currentRow[i]; 131 if (i >= getMinimalIterationCount()) { 132 final double delta = FastMath.abs(s - olds); 133 final double rLimit = getRelativeAccuracy() * (FastMath.abs(olds) + FastMath.abs(s)) * 0.5; 134 if ((delta <= rLimit) || (delta <= getAbsoluteAccuracy())) { 135 return s; 136 } 137 } 138 olds = s; 139 } 140 141 } 142 143 }