trapezoidalrule.cs

Labview s trapezoital algorithm

/////////////////////////////////////////////////////////////////////////////
// Copyright (c) 2003 CenterSpace Software, LLC                            //
//                                                                         //
// This code is free software under the Artistic license.                  //
//                                                                         //
// CenterSpace Software                                                    //
// 2098 NW Myrtlewood Way                                                  //
// Corvallis, Oregon, 97330                                                //
// USA                                                                     //
// http://www.centerspace.net                                              //
/////////////////////////////////////////////////////////////////////////////

using System;

namespace Trapezoid
{
  /// <summary>
  /// Simple implementation of the trapezoid rule. Maximum number of iterations 
  /// and/or tolerance can be changed.  Not thread safe. Adapted from source 
  /// code in "Numerical Recipes in C.", 2nd Edition, by Press, Teulosky, 
  /// Vetterling and Flannery.
  /// </summary>
  public class TrapezoidalRule
  {
    private static int maximumIterations_ = 20;
    private static double tolerance_ = 1e-6;
    private static double trap_s;

    /// <summary>
    /// Gets or sets the maximum number of iterations.
    /// </summary>
    public static int MaximumIterations
    {
      get
      {
        return maximumIterations_;
      }
      set
      {
        maximumIterations_ = value;
      }
    }

    /// <summary>
    /// Gets or sets the desired tolerance.
    /// </summary>
    public static double Tolerance
    {
      get
      {
        return tolerance_;
      }
      set
      {
        tolerance_ = value;
      }
    }
    
    /// <summary>
    /// Repeatedly calls the trapezoidal rule on finer and finer intervals until
    /// either the maximum number of iterations is reached or the desired tolerance
    /// is reached. If the maximum number of iterations does not produce a result
    /// within the tolerance, an exception is thrown.
    /// </summary>
    /// <param name="function">A delegate to a function that takes a double and
    /// returns a double.</param>
    /// <param name="lower">Lower bound.</param>
    /// <param name="upper">Upper bound.</param>
    /// <exception>Throws an exception when the rule does not converge.</exception>
    /// <returns></returns>
    public static double Integrate( DoubleFunction function, double lower, double upper )
    {
      double olds = -1.0e30;
      double s;

      for ( int i = 1; i <= maximumIterations_; i++ ) 
      {
        s = trapzd( function, lower, upper, i );
        if ( Math.Abs( s - olds ) < tolerance_ * Math.Abs( olds )) 
        {
          return s;
        }
        olds = s;
      }
      throw new Exception( "Did not converge" );
    }

    private static double trapzd( DoubleFunction function, double lower, double upper, int iterations )
    {
      double x, tnm, sum, del;
      int it, j;

      if ( iterations == 1 ) 
      {
        trap_s = 0.5 * ( upper - lower ) * ( function(lower) + function(upper));
        return trap_s;
      }
      else
      {
        it = 1;
        for ( int i = 1; i < iterations - 1; i++ )
        {
          it <<= 1;
        }
        tnm = it;
        del = ( upper - lower ) / tnm;
        x = lower + 0.5 * del;
        sum = 0.0;
        for ( j = 1; j <= it; j++ )
        {
          sum += function( x );
          x += del;
        }
        trap_s = 0.5 * ( trap_s + ( upper - lower ) * sum / tnm );
        return trap_s;
      }
    }
  }

  /// <summary>
  /// A function that takes double and returns a double.
  /// </summary>
  public delegate double DoubleFunction( double d );
}