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: 825919 $ $Date: 2009-10-16 10:51:55 -0400 (Fri, 16 Oct 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 clearResult(); 070 verifyInterval(min, max); 071 verifyIterationCount(); 072 073 TrapezoidIntegrator qtrap = new TrapezoidIntegrator(); 074 if (minimalIterationCount == 1) { 075 final double s = (4 * qtrap.stage(f, min, max, 1) - qtrap.stage(f, min, max, 0)) / 3.0; 076 setResult(s, 1); 077 return result; 078 } 079 // Simpson's rule requires at least two trapezoid stages. 080 double olds = 0; 081 double oldt = qtrap.stage(f, min, max, 0); 082 for (int i = 1; i <= maximalIterationCount; ++i) { 083 final double t = qtrap.stage(f, min, max, i); 084 final double s = (4 * t - oldt) / 3.0; 085 if (i >= minimalIterationCount) { 086 final double delta = Math.abs(s - olds); 087 final double rLimit = 088 relativeAccuracy * (Math.abs(olds) + Math.abs(s)) * 0.5; 089 if ((delta <= rLimit) || (delta <= absoluteAccuracy)) { 090 setResult(s, i); 091 return result; 092 } 093 } 094 olds = s; 095 oldt = t; 096 } 097 throw new MaxIterationsExceededException(maximalIterationCount); 098 } 099 100 /** {@inheritDoc} */ 101 @Override 102 protected void verifyIterationCount() throws IllegalArgumentException { 103 super.verifyIterationCount(); 104 // at most 64 bisection refinements 105 if (maximalIterationCount > 64) { 106 throw MathRuntimeException.createIllegalArgumentException( 107 "invalid iteration limits: min={0}, max={1}", 108 0, 64); 109 } 110 } 111 }