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/TrapezoidalRule.html"> 026 * Trapezoidal 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 * The function should be integrable.</p> 031 * 032 * @version $Revision: 799857 $ $Date: 2009-08-01 09:07:12 -0400 (Sat, 01 Aug 2009) $ 033 * @since 1.2 034 */ 035 public class TrapezoidIntegrator extends UnivariateRealIntegratorImpl { 036 037 /** Intermediate result. */ 038 private double s; 039 040 /** 041 * Construct an integrator for the given function. 042 * 043 * @param f function to integrate 044 * @deprecated as of 2.0 the integrand function is passed as an argument 045 * to the {@link #integrate(UnivariateRealFunction, double, double)}method. 046 */ 047 @Deprecated 048 public TrapezoidIntegrator(UnivariateRealFunction f) { 049 super(f, 64); 050 } 051 052 /** 053 * Construct an integrator. 054 */ 055 public TrapezoidIntegrator() { 056 super(64); 057 } 058 059 /** 060 * Compute the n-th stage integral of trapezoid rule. This function 061 * should only be called by API <code>integrate()</code> in the package. 062 * To save time it does not verify arguments - caller does. 063 * <p> 064 * The interval is divided equally into 2^n sections rather than an 065 * arbitrary m sections because this configuration can best utilize the 066 * alrealy computed values.</p> 067 * 068 * @param f the integrand function 069 * @param min the lower bound for the interval 070 * @param max the upper bound for the interval 071 * @param n the stage of 1/2 refinement, n = 0 is no refinement 072 * @return the value of n-th stage integral 073 * @throws FunctionEvaluationException if an error occurs evaluating the 074 * function 075 */ 076 double stage(final UnivariateRealFunction f, 077 final double min, final double max, final int n) 078 throws FunctionEvaluationException { 079 080 long i, np; 081 double x, spacing, sum = 0; 082 083 if (n == 0) { 084 s = 0.5 * (max - min) * (f.value(min) + f.value(max)); 085 return s; 086 } else { 087 np = 1L << (n-1); // number of new points in this stage 088 spacing = (max - min) / np; // spacing between adjacent new points 089 x = min + 0.5 * spacing; // the first new point 090 for (i = 0; i < np; i++) { 091 sum += f.value(x); 092 x += spacing; 093 } 094 // add the new sum to previously calculated result 095 s = 0.5 * (s + sum * spacing); 096 return s; 097 } 098 } 099 100 /** {@inheritDoc} */ 101 @Deprecated 102 public double integrate(final double min, final double max) 103 throws MaxIterationsExceededException, FunctionEvaluationException, IllegalArgumentException { 104 return integrate(f, min, max); 105 } 106 107 /** {@inheritDoc} */ 108 public double integrate(final UnivariateRealFunction f, 109 final double min, final double max) 110 throws MaxIterationsExceededException, FunctionEvaluationException, IllegalArgumentException { 111 112 int i = 1; 113 double t, oldt; 114 115 clearResult(); 116 verifyInterval(min, max); 117 verifyIterationCount(); 118 119 oldt = stage(f, min, max, 0); 120 while (i <= maximalIterationCount) { 121 t = stage(f, min, max, i); 122 if (i >= minimalIterationCount) { 123 final double delta = Math.abs(t - oldt); 124 final double rLimit = 125 relativeAccuracy * (Math.abs(oldt) + Math.abs(t)) * 0.5; 126 if ((delta <= rLimit) || (delta <= absoluteAccuracy)) { 127 setResult(t, i); 128 return result; 129 } 130 } 131 oldt = t; 132 i++; 133 } 134 throw new MaxIterationsExceededException(maximalIterationCount); 135 } 136 137 /** {@inheritDoc} */ 138 @Override 139 protected void verifyIterationCount() throws IllegalArgumentException { 140 super.verifyIterationCount(); 141 // at most 64 bisection refinements 142 if (maximalIterationCount > 64) { 143 throw MathRuntimeException.createIllegalArgumentException( 144 "invalid iteration limits: min={0}, max={1}", 145 0, 64); 146 } 147 } 148 }