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 <a href="http://mathworld.wolfram.com/SimpsonsRule.html"> 028 * Simpson's Rule</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 * This implementation employs the basic trapezoid rule to calculate Simpson's 033 * rule.</p> 034 * 035 * @version $Id: SimpsonIntegrator.java 1364387 2012-07-22 18:14:11Z tn $ 036 * @since 1.2 037 */ 038 public class SimpsonIntegrator extends BaseAbstractUnivariateIntegrator { 039 040 /** Maximal number of iterations for Simpson. */ 041 public static final int SIMPSON_MAX_ITERATIONS_COUNT = 64; 042 043 /** 044 * Build a Simpson 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 #SIMPSON_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 #SIMPSON_MAX_ITERATIONS_COUNT} 056 */ 057 public SimpsonIntegrator(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 > SIMPSON_MAX_ITERATIONS_COUNT) { 064 throw new NumberIsTooLargeException(maximalIterationCount, 065 SIMPSON_MAX_ITERATIONS_COUNT, false); 066 } 067 } 068 069 /** 070 * Build a Simpson 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 #SIMPSON_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 #SIMPSON_MAX_ITERATIONS_COUNT} 080 */ 081 public SimpsonIntegrator(final int minimalIterationCount, 082 final int maximalIterationCount) 083 throws NotStrictlyPositiveException, NumberIsTooSmallException, NumberIsTooLargeException { 084 super(minimalIterationCount, maximalIterationCount); 085 if (maximalIterationCount > SIMPSON_MAX_ITERATIONS_COUNT) { 086 throw new NumberIsTooLargeException(maximalIterationCount, 087 SIMPSON_MAX_ITERATIONS_COUNT, false); 088 } 089 } 090 091 /** 092 * Construct an integrator with default settings. 093 * (max iteration count set to {@link #SIMPSON_MAX_ITERATIONS_COUNT}) 094 */ 095 public SimpsonIntegrator() { 096 super(DEFAULT_MIN_ITERATIONS_COUNT, SIMPSON_MAX_ITERATIONS_COUNT); 097 } 098 099 /** {@inheritDoc} */ 100 @Override 101 protected double doIntegrate() 102 throws TooManyEvaluationsException, MaxCountExceededException { 103 104 TrapezoidIntegrator qtrap = new TrapezoidIntegrator(); 105 if (getMinimalIterationCount() == 1) { 106 return (4 * qtrap.stage(this, 1) - qtrap.stage(this, 0)) / 3.0; 107 } 108 109 // Simpson's rule requires at least two trapezoid stages. 110 double olds = 0; 111 double oldt = qtrap.stage(this, 0); 112 while (true) { 113 final double t = qtrap.stage(this, iterations.getCount()); 114 iterations.incrementCount(); 115 final double s = (4 * t - oldt) / 3.0; 116 if (iterations.getCount() >= getMinimalIterationCount()) { 117 final double delta = FastMath.abs(s - olds); 118 final double rLimit = 119 getRelativeAccuracy() * (FastMath.abs(olds) + FastMath.abs(s)) * 0.5; 120 if ((delta <= rLimit) || (delta <= getAbsoluteAccuracy())) { 121 return s; 122 } 123 } 124 olds = s; 125 oldt = t; 126 } 127 128 } 129 130 }