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.math3.analysis.integration; 018 019import org.apache.commons.math3.exception.MaxCountExceededException; 020import org.apache.commons.math3.exception.NotStrictlyPositiveException; 021import org.apache.commons.math3.exception.NumberIsTooLargeException; 022import org.apache.commons.math3.exception.NumberIsTooSmallException; 023import org.apache.commons.math3.exception.TooManyEvaluationsException; 024import 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 * @since 1.2 037 */ 038public class RombergIntegrator extends BaseAbstractUnivariateIntegrator { 039 040 /** Maximal number of iterations for Romberg. */ 041 public static final int ROMBERG_MAX_ITERATIONS_COUNT = 32; 042 043 /** 044 * Build a Romberg integrator with given accuracies and iterations counts. 045 * @param relativeAccuracy relative accuracy of the result 046 * @param absoluteAccuracy absolute accuracy of the result 047 * @param minimalIterationCount minimum number of iterations 048 * @param maximalIterationCount maximum number of iterations 049 * (must be less than or equal to {@link #ROMBERG_MAX_ITERATIONS_COUNT}) 050 * @exception NotStrictlyPositiveException if minimal number of iterations 051 * is not strictly positive 052 * @exception NumberIsTooSmallException if maximal number of iterations 053 * is lesser than or equal to the minimal number of iterations 054 * @exception NumberIsTooLargeException if maximal number of iterations 055 * is greater than {@link #ROMBERG_MAX_ITERATIONS_COUNT} 056 */ 057 public RombergIntegrator(final double relativeAccuracy, 058 final double absoluteAccuracy, 059 final int minimalIterationCount, 060 final int maximalIterationCount) 061 throws NotStrictlyPositiveException, NumberIsTooSmallException, NumberIsTooLargeException { 062 super(relativeAccuracy, absoluteAccuracy, minimalIterationCount, maximalIterationCount); 063 if (maximalIterationCount > ROMBERG_MAX_ITERATIONS_COUNT) { 064 throw new NumberIsTooLargeException(maximalIterationCount, 065 ROMBERG_MAX_ITERATIONS_COUNT, false); 066 } 067 } 068 069 /** 070 * Build a Romberg integrator with given iteration counts. 071 * @param minimalIterationCount minimum number of iterations 072 * @param maximalIterationCount maximum number of iterations 073 * (must be less than or equal to {@link #ROMBERG_MAX_ITERATIONS_COUNT}) 074 * @exception NotStrictlyPositiveException if minimal number of iterations 075 * is not strictly positive 076 * @exception NumberIsTooSmallException if maximal number of iterations 077 * is lesser than or equal to the minimal number of iterations 078 * @exception NumberIsTooLargeException if maximal number of iterations 079 * is greater than {@link #ROMBERG_MAX_ITERATIONS_COUNT} 080 */ 081 public RombergIntegrator(final int minimalIterationCount, 082 final int maximalIterationCount) 083 throws NotStrictlyPositiveException, NumberIsTooSmallException, NumberIsTooLargeException { 084 super(minimalIterationCount, maximalIterationCount); 085 if (maximalIterationCount > ROMBERG_MAX_ITERATIONS_COUNT) { 086 throw new NumberIsTooLargeException(maximalIterationCount, 087 ROMBERG_MAX_ITERATIONS_COUNT, false); 088 } 089 } 090 091 /** 092 * Construct a Romberg integrator with default settings 093 * (max iteration count set to {@link #ROMBERG_MAX_ITERATIONS_COUNT}) 094 */ 095 public RombergIntegrator() { 096 super(DEFAULT_MIN_ITERATIONS_COUNT, ROMBERG_MAX_ITERATIONS_COUNT); 097 } 098 099 /** {@inheritDoc} */ 100 @Override 101 protected double doIntegrate() 102 throws TooManyEvaluationsException, MaxCountExceededException { 103 104 final int m = getMaximalIterationCount() + 1; 105 double previousRow[] = new double[m]; 106 double currentRow[] = new double[m]; 107 108 TrapezoidIntegrator qtrap = new TrapezoidIntegrator(); 109 currentRow[0] = qtrap.stage(this, 0); 110 incrementCount(); 111 double olds = currentRow[0]; 112 while (true) { 113 114 final int i = getIterations(); 115 116 // switch rows 117 final double[] tmpRow = previousRow; 118 previousRow = currentRow; 119 currentRow = tmpRow; 120 121 currentRow[0] = qtrap.stage(this, i); 122 incrementCount(); 123 for (int j = 1; j <= i; j++) { 124 // Richardson extrapolation coefficient 125 final double r = (1L << (2 * j)) - 1; 126 final double tIJm1 = currentRow[j - 1]; 127 currentRow[j] = tIJm1 + (tIJm1 - previousRow[j - 1]) / r; 128 } 129 final double s = currentRow[i]; 130 if (i >= getMinimalIterationCount()) { 131 final double delta = FastMath.abs(s - olds); 132 final double rLimit = getRelativeAccuracy() * (FastMath.abs(olds) + FastMath.abs(s)) * 0.5; 133 if ((delta <= rLimit) || (delta <= getAbsoluteAccuracy())) { 134 return s; 135 } 136 } 137 olds = s; 138 } 139 140 } 141 142}