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/SimpsonsRule.html"> 026 * Simpson's Rule</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 * This implementation employs basic trapezoid rule as building blocks to 031 * calculate the Simpson's rule of alternating 2/3 and 4/3.</p> 032 * 033 * @version $Revision: 799857 $ $Date: 2009-08-01 09:07:12 -0400 (Sat, 01 Aug 2009) $ 034 * @since 1.2 035 */ 036 public class SimpsonIntegrator extends UnivariateRealIntegratorImpl { 037 038 /** 039 * Construct an integrator for the given function. 040 * 041 * @param f function to integrate 042 * @deprecated as of 2.0 the integrand function is passed as an argument 043 * to the {@link #integrate(UnivariateRealFunction, double, double)}method. 044 */ 045 @Deprecated 046 public SimpsonIntegrator(UnivariateRealFunction f) { 047 super(f, 64); 048 } 049 050 /** 051 * Construct an integrator. 052 */ 053 public SimpsonIntegrator() { 054 super(64); 055 } 056 057 /** {@inheritDoc} */ 058 @Deprecated 059 public double integrate(final double min, final double max) 060 throws MaxIterationsExceededException, FunctionEvaluationException, IllegalArgumentException { 061 return integrate(f, min, max); 062 } 063 064 /** {@inheritDoc} */ 065 public double integrate(final UnivariateRealFunction f, 066 final double min, final double max) 067 throws MaxIterationsExceededException, FunctionEvaluationException, IllegalArgumentException { 068 069 int i = 1; 070 double s, olds, t, oldt; 071 072 clearResult(); 073 verifyInterval(min, max); 074 verifyIterationCount(); 075 076 TrapezoidIntegrator qtrap = new TrapezoidIntegrator(); 077 if (minimalIterationCount == 1) { 078 s = (4 * qtrap.stage(f, min, max, 1) - qtrap.stage(f, min, max, 0)) / 3.0; 079 setResult(s, 1); 080 return result; 081 } 082 // Simpson's rule requires at least two trapezoid stages. 083 olds = 0; 084 oldt = qtrap.stage(f, min, max, 0); 085 while (i <= maximalIterationCount) { 086 t = qtrap.stage(f, min, max, i); 087 s = (4 * t - oldt) / 3.0; 088 if (i >= minimalIterationCount) { 089 final double delta = Math.abs(s - olds); 090 final double rLimit = 091 relativeAccuracy * (Math.abs(olds) + Math.abs(s)) * 0.5; 092 if ((delta <= rLimit) || (delta <= absoluteAccuracy)) { 093 setResult(s, i); 094 return result; 095 } 096 } 097 olds = s; 098 oldt = t; 099 i++; 100 } 101 throw new MaxIterationsExceededException(maximalIterationCount); 102 } 103 104 /** {@inheritDoc} */ 105 @Override 106 protected void verifyIterationCount() throws IllegalArgumentException { 107 super.verifyIterationCount(); 108 // at most 64 bisection refinements 109 if (maximalIterationCount > 64) { 110 throw MathRuntimeException.createIllegalArgumentException( 111 "invalid iteration limits: min={0}, max={1}", 112 0, 64); 113 } 114 } 115 }