# Numerical integration
# 
# The algorithms are described in:
# Numerical Analysis, 5th edition
# by Richard L. Burden and J. Douglas Faires
# PWS Publishing Company, Boston, 1993.
# Library of congress: QA 297 B84 1993

# In the below, f indicates the function whose integral we wish to determine,
# a,b indicate the left and right endpoints of the interval over which
# we wish to integrate, and n is the number of intervals into which we
# divide [a,b]

# These methods all return one value, the value of the integral

# Currently only works for real functions of a real variable

# Composite Simpson's Rule, Section 4.4, Algorithm 4.1, p. 186
# Note that this has error term = max(f'''')*h^4*(b-a)/180,
# where h=(b-a)/n
# If we can get maximums and derivatives, this would allow us to determine
# automatically what n should be.

# Composite simpson's rule is implemented as a built in function

SetHelp ("CompositeSimpsonsRuleTolerance", "calculus", "Integration of f by Composite Simpson's Rule on the interval [a,b] with the number of steps calculated by the fourth derivative bound and the desired tolerance")
function CompositeSimpsonsRuleTolerance(f,a,b,FourthDerivativeBound,Tolerance) =
(
	# Error term = max(f'''')*h^4*(b-a)/180,
	# where h=(b-a)/n
	n = ceil(|FourthDerivativeBound*(b-a)^5 / (180*Tolerance)|^(1/4));

	# Note that this is done automatically by CompositeSimpsonsRule
	# function:
	#if IsOdd(n) then n = n+1;

	CompositeSimpsonsRule (f, a, b, n)
)
protect ("CompositeSimpsonsRuleTolerance")

# FIXME: reference, though this is really stupid
SetHelp ("MidpointRule", "calculus", "Integration by midpoint rule")
function MidpointRule(f,a,b,n) =
(
  if(not IsFunction(f)) then
	(error("MidpointRule: argument 1 must be a function");bailout)
  else if(not IsReal(a) or not IsReal(b)) then
	(error("MidpointRule: arguments 2, 3 must be real values");bailout)
  else if(not IsInteger(n)) then
	(error("MidpointRule: argument 4 must be an integer");bailout);
  ## check bounds
  if(a>b) then (error("MidpointRule: argument 2 must be less than or equal to argument 3");bailout)
  else if(n<= 0) then (error("MidpointRule: argument 4 must be positive");bailout);

  len = float(b-a);
  s = sum i=1 to n do float(f(a+(len*(i-0.5))/n));
  (s*len)/n
)
protect ("MidpointRule")

SetHelp ("NumericalIntegralSteps", "parameters", "Steps to perform in NumericalIntegral")
parameter NumericalIntegralSteps = 1000

SetHelp ("NumericalIntegralFunction", "parameters", "The function used for numerical integration in NumericalIntegral")
parameter NumericalIntegralFunction = `CompositeSimpsonsRule

SetHelp ("NumericalIntegral", "calculus", "Integration by rule set in NumericalIntegralFunction of f from a to b using NumericalIntegralSteps steps")
function NumericalIntegral(f,a,b) = (NumericalIntegralFunction call (f,a,b,NumericalIntegralSteps))
protect ("NumericalIntegral")