Numerical integration

From Rosetta Code

(Redirected from Numerical Integration)
Jump to: navigation, search
Task
Numerical integration
You are encouraged to solve this task according to the task description, using any language you may know.

Write functions to calculate the definite integral of a function (f(x)) using rectangular (left, right, and midpoint), trapezium, and Simpson's methods. Your functions should take in the upper and lower bounds (a and b) and the number of approximations to make in that range (n). Assume that your example already has a function that gives values for f(x).

Simpson's method is defined by the following pseudocode:

h := (b - a) / n
sum1 := f(a + h/2)
sum2 := 0

loop on i from 1 to (n - 1)
    sum1 := sum1 + f(a + h * i + h/2)
    sum2 := sum2 + f(a + h * i)

answer := (h / 6) * (f(a) + f(b) + 4*sum1 + 2*sum2)

What it will be doing is computing integrals for multiple quadratics (like the one shown in Wikipedia) and summing them.

See also


Contents

[edit] ActionScript

Integration functions:

function leftRect(f:Function, a:Number, b:Number, n:uint):Number 
{
var sum:Number = 0;
var dx:Number = (b-a)/n;
for (var x:Number = a; n > 0; n--, x += dx)
sum += f(x);
return sum * dx;
}
 
function rightRect(f:Function, a:Number, b:Number, n:uint):Number
{
var sum:Number = 0;
var dx:Number = (b-a)/n;
for (var x:Number = a + dx; n > 0; n--, x += dx)
sum += f(x);
return sum * dx;
}
 
function midRect(f:Function, a:Number, b:Number, n:uint):Number
{
var sum:Number = 0;
var dx:Number = (b-a)/n;
for (var x:Number = a + (dx / 2); n > 0; n--, x += dx)
sum += f(x);
return sum * dx;
}
function trapezium(f:Function, a:Number, b:Number, n:uint):Number
{
var dx:Number = (b-a)/n;
var x:Number = a;
var sum:Number = f(a);
for(var i:uint = 1; i < n; i++)
{
a += dx;
sum += f(a)*2;
}
sum += f(b);
return 0.5 * dx * sum;
}
function simpson(f:Function, a:Number, b:Number, n:uint):Number
{
var dx:Number = (b-a)/n;
var sum1:Number = f(a + dx/2);
var sum2:Number = 0;
for(var i:uint = 1; i < n; i++)
{
sum1 += f(a + dx*i + dx/2);
sum2 += f(a + dx*i);
}
return (dx/6) * (f(a) + f(b) + 4*sum1 + 2*sum2);
}

Usage:

function f1(n:Number):Number {
return (2/(1+ 4*(n*n)));
}
trace(leftRect(f1, -1, 2, 4));
trace(rightRect(f1, -1, 2, 4));
trace(midRect(f1, -1, 2, 4));
trace(trapezium(f1, -1, 2 ,4 ));
trace(simpson(f1, -1, 2 ,4 ));
 

[edit] Ada

This solution creates a generic package into which the function F(X) is passed during generic instantiation. The first part is the package specification. The second part is the package body.

generic
with function F(X : Long_Float) return Long_Float;
package Integrate is
function Left_Rect(A, B, N : Long_Float) return Long_Float;
function Right_Rect(A, B, N : Long_Float) return Long_Float;
function Mid_Rect(A, B, N : Long_Float) return Long_Float;
function Trapezium(A, B, N : Long_Float) return Long_Float;
function Simpson(A, B, N : Long_Float) return Long_Float;
end Integrate;
package body Integrate is
 
---------------
-- Left_Rect --
---------------
 
function Left_Rect (A, B, N : Long_Float) return Long_Float is
H : constant Long_Float := (B - A) / N;
Sum : Long_Float := 0.0;
X : Long_Float := A;
begin
while X <= B - H loop
Sum := Sum + (H * F(X));
X := X + H;
end loop;
return Sum;
end Left_Rect;
 
----------------
-- Right_Rect --
----------------
 
function Right_Rect (A, B, N : Long_Float) return Long_Float is
H : constant Long_Float := (B - A) / N;
Sum : Long_Float := 0.0;
X : Long_Float := A + H;
begin
while X <= B - H loop
Sum := Sum + (H * F(X));
X := X + H;
end loop;
return Sum;
end Right_Rect;
 
--------------
-- Mid_Rect --
--------------
 
function Mid_Rect (A, B, N : Long_Float) return Long_Float is
H : constant Long_Float := (B - A) / N;
Sum : Long_Float := 0.0;
X : Long_Float := A;
begin
while X <= B - H loop
Sum := Sum + (H / 2.0) * (F(X) + F(X + H));
X := X + H;
end loop;
return Sum;
end Mid_Rect;
 
---------------
-- Trapezium --
---------------
 
function Trapezium (A, B, N : Long_Float) return Long_Float is
H : constant Long_Float := (B - A) / N;
Sum : Long_Float := F(A) + F(B);
X : Long_Float := 1.0;
begin
while X <= N - 1.0 loop
Sum := Sum + 2.0 * F(A + X * (B - A) / N);
X := X + 1.0;
end loop;
return (B - A) / (2.0 * N) * Sum;
end Trapezium;
 
-------------
-- Simpson --
-------------
 
function Simpson (A, B, N : Long_Float) return Long_Float is
H : constant Long_Float := (B - A) / N;
Sum1 : Long_Float := 0.0;
Sum2 : Long_Float := 0.0;
Limit : Integer := Integer(N) - 1;
begin
for I in 0..Limit loop
Sum1 := Sum1 + F(A + H * Long_Float(I) + H / 2.0);
end loop;
for I in 1..Limit loop
Sum2 := Sum2 + F(A + H * Long_Float(I));
end loop;
return H / 6.0 * (F(A) + F(B) + 4.0 * Sum1 + 2.0 * Sum2);
end Simpson;
 
end Integrate;

[edit] ALGOL 68

MODE F = PROC(LONG REAL)LONG REAL;
 
###############
#
# left rect ##
#
##############
 
PROC left rect = (F f, LONG REAL a, b, INT n) LONG REAL:
BEGIN
LONG REAL h= (b - a) / n;
LONG REAL sum:= 0;
LONG REAL x:= a;
WHILE x <= b - h DO
sum := sum + (h * f(x));
x +:= h
OD;
sum
END # left rect #;
 
#################
#
# right rect ##
#
################
 
PROC right rect = (F f, LONG REAL a, b, INT n) LONG REAL:
BEGIN
LONG REAL h= (b - a) / n;
LONG REAL sum:= 0;
LONG REAL x:= a + h;
WHILE x <= b - h DO
sum := sum + (h * f(x));
x +:= h
OD;
sum
END # right rect #;
 
###############
#
# mid rect ##
#
##############
 
PROC mid rect = (F f, LONG REAL a, b, INT n) LONG REAL:
BEGIN
LONG REAL h= (b - a) / n;
LONG REAL sum:= 0;
LONG REAL x:= a;
WHILE x <= b - h DO
sum := sum + (h / 2) * (f(x) + f(x + h));
x +:= h
OD;
sum
END # mid rect #;
 
###############
#
# trapezium ##
#
##############
 
PROC trapezium = (F f, LONG REAL a, b, INT n) LONG REAL:
BEGIN
LONG REAL h= (b - a) / n;
LONG REAL sum:= f(a) + f(b);
LONG REAL x:= 1;
WHILE x <= n - 1 DO
sum := sum + 2 * f(a + x * h );
x +:= 1
OD;
(b - a) / (2 * n) * sum
END # trapezium #;
 
#############
#
# simpson ##
#
############
 
PROC simpson = (F f, LONG REAL a, b, INT n) LONG REAL:
BEGIN
LONG REAL h= (b - a) / n;
LONG REAL sum1:= 0;
LONG REAL sum2:= 0;
INT limit:= n - 1;
FOR i FROM 0 TO limit DO
sum1 := sum1 + f(a + h * LONG REAL(i) + h / 2)
OD;
FOR i FROM 1 TO limit DO
sum2 +:= f(a + h * LONG REAL(i))
OD;
h / 6 * (f(a) + f(b) + 4 * sum1 + 2 * sum2)
END # simpson #;
SKIP

[edit] AutoHotkey

ahk discussion

MsgBox % Rect("fun", 0, 1, 10,-1) ; 0.45 left
MsgBox % Rect("fun", 0, 1, 10) ; 0.50 mid
MsgBox % Rect("fun", 0, 1, 10, 1) ; 0.55 right
MsgBox % Trapez("fun", 0, 1, 10) ; 0.50
MsgBox % Simpson("fun", 0, 1, 10) ; 0.50
 
Rect(f,a,b,n,side=0) { ; side: -1=left, 0=midpoint, 1=right
h := (b - a) / n
sum := 0, a += (side-1)*h/2
Loop %n%
sum += %f%(a + h*A_Index)
Return h*sum
}
 
Trapez(f,a,b,n) {
h := (b - a) / n
sum := 0
Loop % n-1
sum += %f%(a + h*A_Index)
Return h/2 * (%f%(a) + %f%(b) + 2*sum)
}
 
Simpson(f,a,b,n) {
h := (b - a) / n
sum1 := sum2 := 0, ah := a - h/2
Loop %n%
sum1 += %f%(ah + h*A_Index)
Loop % n-1
sum2 += %f%(a + h*A_Index)
Return h/6 * (%f%(a) + %f%(b) + 4*sum1 + 2*sum2)
}
 
fun(x) { ; linear test function
Return x
}

[edit] BASIC

Works with: QuickBasic version 4.5 Translation of: Java

FUNCTION leftRect(a, b, n)
h = (b - a) / n
sum = 0
FOR x = a TO b - h STEP h
sum = sum + h * (f(x))
NEXT x
leftRect = sum
END FUNCTION
 
FUNCTION rightRect(a, b, n)
h = (b - a) / n
sum = 0
FOR x = a + h TO b - h STEP h
sum = sum + h * (f(x))
NEXT x
rightRect = sum
END FUNCTION
 
FUNCTION midRect(a, b, n)
h = (b - a) / n
sum = 0
FOR x = a TO b - h STEP h
sum = sum + (h / 2) * (f(x) + f(x + h))
NEXT x
midRect = sum
END FUNCTION
 
FUNCTION trap(a, b, n)
h = (b - a) / n
sum = f(a) + f(b)
FOR i = 1 TO n-1
sum = sum + 2 * f((a + i * h))
NEXT i
trap = h / 2 * sum
END FUNCTION
 
FUNCTION simpson(a, b, n)
h = (b - a) / n
sum1 = 0
sum2 = 0
 
FOR i = 0 TO n-1
sum1 = sum + f(a + h * i + h / 2)
NEXT i
 
FOR i = 1 TO n - 1
sum2 = sum2 + f(a + h * i)
NEXT i
 
simpson = h / 6 * (f(a) + f(b) + 4 * sum1 + 2 * sum2)
END FUNCTION

[edit] C

Translation of: Java

#include <math.h>
 
double int_leftrect(double from, double to, double n, double (*func)())
{
double h = (to-from)/n;
double sum = 0, x;
for(x=from; x <= (to-h); x += h)
sum += func(x);
return h*sum;
}
 
double int_rightrect(double from, double to, double n, double (*func)())
{
double h = (to-from)/n;
double sum = 0, x;
for(x=from+h; x <= (to-h); x += h)
sum += func(x);
return h*sum;
}
 
double int_midrect(double from, double to, double n, double (*func)())
{
double h = (to-from)/n;
double sum = 0.0, x;
for(x=from; x <= (to-h); x += h)
sum += (func(x) + func(x + h));
return h*sum/2.0;
}
 
double int_trapezium(double from, double to, double n, double (*func)())
{
double h = (to - from) / n;
double sum = func(from) + func(to);
int i;
for(i = 1;i < n;i++)
sum += func(from + i * h);
return h * sum;
}
 
double int_simpson(double from, double to, double n, double (*func)())
{
double h = (to - from) / n;
double sum1 = 0.0;
double sum2 = 0.0;
int i;
 
for(i = 0;i < n;i++)
sum1 += func(from + h * i + h / 2.0);
 
for(i = 1;i < n;i++)
sum2 += func(from + h * i);
 
return h / 6.0 * (func(from) + func(to) + 4.0 * sum1 + 2.0 * sum2);
}

Usage example:

#include <stdio.h>
#include <stdlib.h>
 
double f1(double x)
{
return 2.0*x*log(x*x+1.0);
}
 
double f1a(double x)
{
double y = x*x+1;
return y*(log(y)-1.0);
}
 
double f2(double x)
{
return cos(x);
}
 
double f2a(double x)
{
return sin(x);
}
 
typedef double (*pfunc)(double, double, double, double (*)());
typedef double (*rfunc)(double);
 
#define INTG(F,A,B) (F((B))-F((A)))
 
int main()
{
int i,j;
double ic;
 
pfunc f[5] = { &int_leftrect, &int_rightrect,
&int_midrect, &int_trapezium,
&int_simpson };
const char *names[5] = {"leftrect", "rightrect", "midrect",
"trapezium", "simpson" };
rfunc rf[2] = { &f1, &f2 };
rfunc If[2] = { &f1a, &f2a };
 
for(j=0; j<2; j++)
{
for(i=0; i < 5 ; i++)
{
ic = (*f[i])(0.0, 1.0, 100.0, rf[j]);
printf("%10s [ 0,1] num: %+lf, an: %lf\n",
names[i], ic, INTG((*If[j]), 0.0, 1.0));
}
printf("\n");
}
}

Output of the test (with a little bit of enhancing by hand):

  leftrect [ 0,1] num: +0.365865, an: 0.386294
 rightrect [ 0,1] num: +0.365865
   midrect [ 0,1] num: +0.372628
 trapezium [ 0,1] num: +0.393254
   simpson [ 0,1] num: +0.386294

  leftrect [ 0,1] num: +0.838276, an: 0.841471
 rightrect [ 0,1] num: +0.828276
   midrect [ 0,1] num: +0.836019
 trapezium [ 0,1] num: +0.849165
   simpson [ 0,1] num: +0.841471

[edit] C++

Due to their similarity, it makes sense to make the integration method a policy.

// the integration routine
template<typename Method, typename F, typename Float>
double integrate(F f, Float a, Float b, int steps, Method m)
{
double s = 0;
double h = (b-a)/steps;
for (int i = 0; i < steps; ++i)
s += m(f, a + h*i, h);
return h*s;
}
 
// methods
class rectangular
{
public:
enum position_type { left, middle, right };
rectangular(position_type pos): position(pos) {}
template<typename F, typename Float>
double operator()(F f, Float x, Float h)
{
switch(position)
{
case left:
return f(x);
case middle:
return f(x+h/2);
case right:
return f(x+h);
}
}
private:
position_type position;
};
 
class trapezium
{
public:
template<typename F, typename Float>
double operator()(F f, Float x, Float h)
{
return (f(x) + f(x+h))/2;
}
};
 
class simpson
{
public:
template<typename F, typename Float>
double operator()(F f, Float x, Float h)
{
return (f(x) + 4*f(x+h/2) + f(x+h))/6;
}
};
 
// sample usage
double f(double x) { return x*x; )
 
// inside a function somewhere:
double rl = integrate(f, 0.0, 1.0, 10, rectangular(rectangular::left));
double rm = integrate(f, 0.0, 1.0, 10, rectangular(rectangular::middle));
double rr = integrate(f, 0.0, 1.0, 10, rectangular(rectangular::right));
double t = integrate(f, 0.0, 1.0, 10, trapezium());
double s = integrate(f, 0.0, 1.0, 10, simpson());

[edit] Common Lisp

(defun left-rectangle (f a b n &aux (d (/ (- b a) n)))
(* d (loop for x from a below b by d summing (funcall f x))))
 
(defun right-rectangle (f a b n &aux (d (/ (- b a) n)))
(* d (loop for x from b above a by d summing (funcall f x))))
 
(defun midpoint-rectangle (f a b n &aux (d (/ (- b a) n)))
(* d (loop for x from (+ a (/ d 2)) below b by d summing (funcall f x))))
 
(defun trapezium (f a b n &aux (d (/ (- b a) n)))
(* (/ d 2)
(+ (funcall f a)
(* 2 (loop for x from (+ a d) below b by d summing (funcall f x)))
(funcall f b))))
 
(defun simpson (f a b n)
(loop with h = (/ (- b a) n)
with sum1 = (funcall f (+ a (/ h 2)))
with sum2 = 0
for i from 1 below n
do (incf sum1 (funcall f (+ a (* h i) (/ h 2))))
do (incf sum2 (funcall f (+ a (* h i))))
finally (return (* (/ h 6)
(+ (funcall f a)
(funcall f b)
(* 4 sum1)
(* 2 sum2))))))

[edit] D

The function f(x) is Func1/Func2 in this program.

module integrate ;
import std.stdio ;
import std.math ;
 
alias real function(real) realFn ;
 
class Sigma{
int n ;
real a, b , h ;
realFn fn ;
string desc ;
enum Method {LEFT = 0, RGHT, MIDD, TRAP, SIMP} ;
static string[5] methodName =
["LeftRect ", "RightRect", "MidRect ", "Trapezium", "Simpson "] ;
static Sigma opCall() {
Sigma s = new Sigma() ;
return s ;
}
static Sigma opCall(realFn f, int n, real a, real b) {
Sigma s = new Sigma(f, n, a, b) ;
return s ;
}
static real opCall(realFn f, int n, real a, real b, Method m) {
return Sigma(f,n,a,b).getSum(m) ;
}
private this() {} ;
this(realFn f, int n, real a, real b) {
setFunction(f) ;
setStep(n) ;
setRange(a,b) ;
}
Sigma opCall(Method m) {
return doSum(m) ;
}
Sigma setFunction(realFn f) {
this.fn = f ;
return this ;
}
Sigma setRange(real a, real b) {
this.a = a ; this.b = b ;
setInterval() ;
return this ;
}
Sigma setStep(int n) {
this.n = n ;
setInterval() ;
return this ;
}
Sigma setDesc(string d) {
this.desc = d ;
return this ;
}
private void setInterval() {
this.h = (b - a) / n ;
}
private real partSum(int i, Method m) {
real x = a + h * i ;
switch(m) {
case Method.LEFT:
return fn(x) ;
case Method.RGHT:
return fn(x + h) ;
case Method.MIDD:
return fn(x + h/2) ;
case Method.TRAP:
return (fn(x) + fn(x + h))/2 ;
default:
}
//case SIMPSON:
return (fn(x) + 4 * fn(x + h/2) + fn(x + h))/6 ;
}
real getSum(Method m) {
real sum = 0 ;
for(int i = 0; i < n ; i++)
sum += partSum(i, m) ;
return sum * h ;
}
Sigma doSum(Method m) {
writefln("%10s = %9.6f", methodName[m], getSum(m)) ;
return this ;
}
Sigma showSetting() {
writefln("\n%s\nA = %9.6f, B = %9.6f, N = %s, h = %s", desc, a, b, n, h) ;
return this ;
}
Sigma doLeft() { return doSum(Method.LEFT) ; }
Sigma doRight() { return doSum(Method.RGHT) ; }
Sigma doMid() { return doSum(Method.MIDD) ; }
Sigma doTrap() { return doSum(Method.TRAP) ; }
Sigma doSimp() { return doSum(Method.SIMP) ; }
Sigma doAll() {
showSetting() ;
doLeft() ; doRight() ; doMid() ; doTrap() ; doSimp() ;
return this ;
}
}
real Func1(real x) {
return cos(x) + sin(x) ;
}
real Func2(real x) {
return 2.0L/(1 + 4*x*x) ;
}
void main() {
// use as a re-usable and manageable object
Sigma s = Sigma(&Func1, 10, -PI/2, PI/2).showSetting()
.doLeft().doRight().doMid().doTrap()(Sigma.Method.SIMP) ;
s.setFunction(&Func2).setStep(4).setRange(-1.0L,2.0L) ;
s.setDesc("Function(x) = 2 / (1 + 4x^2)").doAll() ;
// use as a single function call
writefln("\nLeftRect Integral of FUNC2 =") ;
writefln("%12.9f (%3dstep)\n%12.9f (%3dstep)\n%12.9f (%3dstep).",
Sigma(&Func2, 8, -1.0L, 2.0L, Sigma.Method.LEFT), 8,
Sigma(&Func2, 64, -1.0L, 2.0L, Sigma.Method.LEFT), 64,
Sigma(&Func2, 512, -1.0L, 2.0L, Sigma.Method.LEFT),512) ;
writefln("\nSimpson Integral of FUNC2 =") ;
writefln("%12.9f (%3dstep).",
Sigma(&Func2, 512, -1.0L, 2.0L, Sigma.Method.SIMP),512) ;
}

Parts of the output:

Function(x) = 2 / (1 + 4x^2)
A = -1.000000, B =  2.000000, N = 4, h = 0.75
 LeftRect  =  2.456897
 RightRect =  2.245132
 MidRect   =  2.496091
 Trapezium =  2.351014
 Simpson   =  2.447732

[edit] E

Translation of: Python

pragma.enable("accumulator")
 
def leftRect(f, x, h) {
return f(x)
}
 
def midRect(f, x, h) {
return f(x + h/2)
}
 
def rightRect(f, x, h) {
return f(x + h)
}
 
def trapezium(f, x, h) {
return (f(x) + f(x+h)) / 2
}
 
def simpson(f, x, h) {
return (f(x) + 4 * f(x + h / 2) + f(x+h)) / 6
}
 
def integrate(f, a, b, steps, meth) {
def h := (b-a) / steps
return h * accum 0 for i in 0..!steps { _ + meth(f, a+i*h, h) }
}
? integrate(fn x { x ** 2 }, 3.0, 7.0, 30, simpson) 
# value: 105.33333333333334
 
? integrate(fn x { x ** 9 }, 0, 1, 300, simpson)
# value: 0.10000000002160479

[edit] Forth

fvariable step
 
defer method ( fn F: x -- fn[x] )
 
: left execute ;
: right step f@ f+ execute ;
: mid step f@ 2e f/ f+ execute ;
: trap
dup fdup left
fswap right f+ 2e f/ ;
: simpson
dup fdup left
dup fover mid 4e f* f+
fswap right f+ 6e f/ ;
 
: set-step ( n F: a b -- n F: a )
fover f- dup 0 d>f f/ step f! ;
 
: integrate ( xt n F: a b -- F: sigma )
set-step
0e
0 do
dup fover method f+
fswap step f@ f+ fswap
loop
drop fnip
step f@ f* ;
\ testing similar to the D example
: test
' is method ' 4 -1e 2e integrate f. ;
 
: fn1 fsincos f+ ;
: fn2 fdup f* 4e f* 1e f+ 2e fswap f/ ;
 
7 set-precision
test left fn2 \ 2.456897
test right fn2 \ 2.245132
test mid fn2 \ 2.496091
test trap fn2 \ 2.351014
test simpson fn2 \ 2.447732

[edit] Fortran

In ISO Fortran 95 and later if function f() is not already defined to be "elemental", define an elemental wrapper function around it to allow for array-based initialization:

elemental function elemf(x)
real :: elemf, x
elemf = f(x)
end function elemf

Use Array Initializers, Pointers, Array invocation of Elemental functions, Elemental array-array and array-scalar arithmetic, and the SUM intrinsic function. Methods are collected into a single function in a module.

module Integration
implicit none
 
contains
 
! function, lower limit, upper limit, steps, method
function integrate(f, a, b, in, method)
real :: integrate
real, intent(in) :: a, b
integer, optional, intent(in) :: in
character(len=*), intent(in), optional :: method
interface
elemental function f(ra)
real :: f
real, intent(in) :: ra
end function f
end interface
 
integer :: n, i, m
real :: h
real, dimension(:), allocatable :: xpoints
real, dimension(:), target, allocatable :: fpoints
real, dimension(:), pointer :: fleft, fmid, fright
 
if ( present(in) ) then
n = in
else
n = 20
end if
 
if ( present(method) ) then
select case (method)
case ('leftrect')
m = 1
case ('midrect')
m = 2
case ('rightrect')
m = 3
case ( 'trapezoid' )
m = 4
case default
m = 0
end select
else
m = 0
end if
 
h = (b - a) / n
 
allocate(xpoints(0:2*n), fpoints(0:2*n))
 
xpoints = (/ (a + h*i/2, i = 0,2*n) /)
 
fpoints = f(xpoints)
fleft => fpoints(0 : 2*n-2 : 2)
fmid => fpoints(1 : 2*n-1 : 2)
fright => fpoints(2 : 2*n : 2)
 
select case (m)
case (0) ! simpson
integrate = h / 6.0 * sum(fleft + fright + 4.0*fmid)
case (1) ! leftrect
integrate = h * sum(fleft)
case (2) ! midrect
integrate = h * sum(fmid)
case (3) ! rightrect
integrate = h * sum(fright)
case (4) ! trapezoid
integrate = h * sum(fleft + fright) / 2
end select
 
deallocate(xpoints, fpoints)
end function integrate
 
end module Integration

Usage example:

program IntegrationTest
use Integration
use FunctionHolder
implicit none
 
print *, integrate(afun, 0., 3**(1/3.), method='simpson')
print *, integrate(afun, 0., 3**(1/3.), method='leftrect')
print *, integrate(afun, 0., 3**(1/3.), method='midrect')
print *, integrate(afun, 0., 3**(1/3.), method='rightrect')
print *, integrate(afun, 0., 3**(1/3.), method='trapezoid')
 
end program IntegrationTest

The FunctionHolder module:

module FunctionHolder
implicit none
 
contains
 
pure function afun(x)
real :: afun
real, intent(in) :: x
 
afun = x**2
end function afun
 
end module FunctionHolder

[edit] Haskell

Different approach from most of the other examples: First, the function f might be expensive to calculate, and so it should not be evaluated several times. So, ideally, we want to have positions x and weights w for each method and then just calculate the approximation of the integral by

approx f xs ws = sum [w * f x | (x,w) <- zip xs ws]

Second, let's to generalize all integration methods into one scheme. The methods can all be characterized by the coefficients vs they use in a particular interval. These will be fractions, and for terseness, we extract the denominator as an extra argument v.

Now there are the closed formulas (which include the endpoints) and the open formulas (which exclude them). Let's do the open formulas first, because then the coefficients don't overlap:

integrateOpen :: Fractional a => a -> [a] -> (a -> a) -> a -> a -> Int -> a
integrateOpen v vs f a b n = approx f xs ws * h / v where
m = fromIntegral (length vs) * n
h = (b-a) / fromIntegral m
ws = concat $ replicate n vs
c = a + h/2
xs = [c + h * fromIntegral i | i <- [0..m-1]]

Similarly for the closed formulas, but we need an additional function overlap which sums the coefficients overlapping at the interior interval boundaries:

integrateClosed :: Fractional a => a -> [a] -> (a -> a) -> a -> a -> Int -> a
integrateClosed v vs f a b n = approx f xs ws * h / v where
m = fromIntegral (length vs - 1) * n
h = (b-a) / fromIntegral m
ws = overlap n vs
xs = [a + h * fromIntegral i | i <- [0..m]]
 
overlap :: Num a => Int -> [a] -> [a]
overlap n [] = []
overlap n (x:xs) = x : inter n xs where
inter 1 ys = ys
inter n [] = x : inter (n-1) xs
inter n [y] = (x+y) : inter (n-1) xs
inter n (y:ys) = y : inter n ys

And now we can just define

intLeftRect  = integrateClosed  1 [1,0]     
intRightRect = integrateClosed 1 [0,1]
intMidRect = integrateOpen 1 [1]
intTrapezium = integrateClosed 2 [1,1]
intSimpson = integrateClosed 3 [1,4,1]

or, as easily, some additional schemes:

intMilne     = integrateClosed 45 [14,64,24,64,14]
intOpen1 = integrateOpen 2 [3,3]
intOpen2 = integrateOpen 3 [8,-4,8]

Some examples:

*Main> intLeftRect  (\x -> x*x) 0 1 10  
0.2850000000000001
*Main> intRightRect (\x -> x*x) 0 1 10  
0.38500000000000006
*Main> intMidRect   (\x -> x*x) 0 1 10 
0.3325
*Main> intTrapezium (\x -> x*x) 0 1 10 
0.3350000000000001
*Main> intSimpson   (\x -> x*x) 0 1 10 
0.3333333333333334

[edit] J

Solution:

integrate=: adverb define
'a b steps'=. 3{.y,128
size=. (b - a)%steps
size * +/ 2 u\ a + size * i.>:steps
)
 
rectangle=: adverb def 'u -: +/ y'
 
trapezium=: adverb def '-: +/ u y'
 
simpson =: adverb def '6 %~ +/ 1 1 4 * u y, -:+/y'

Example usage
Integrate square (*:) from 0 to π in 10 steps using various methods.

   *: rectangle integrate 0 1p1 10        
10.3095869962
*: trapezium integrate 0 1p1 10
10.3871026879
*: simpson integrate 0 1p1 10
10.3354255601

Integrate sin from 0 to π in 10 steps using various methods.

   sin=: 1&o.
sin rectangle integrate 0 1p1 10
2.00824840791
sin trapezium integrate 0 1p1 10
1.98352353751
sin simpson integrate 0 1p1 10
2.00000678444

Note that J has a primitive verb p.. for integrating polynomials. For example the integral of x2 (which can be described in terms of its coefficients as 0 0 1) is:

   0 p.. 0 0 1
0 0 0 0.333333333333
0 p.. 0 0 1x NB. or using rationals
0 0 0 1r3

That is: 0x^0 + 0x^1 + 0x^2 + \tfrac{1}{3}x^3
So to integrate x2 from 0 to π :

   0 0 1 (0&p..@[ -~/@:p. ]) 0 1p1 
10.3354255601

[edit] Java

The function in this example is assumed to be f(double x).

public class Integrate{
public static double leftRect(double a, double b, double n){
double h = (b - a) / n;
double sum = 0;
for(double x = a;x <= b - h;x += h)
sum += f(x);
return h * sum;
}
 
public static double rightRect(double a, double b, double n){
double h = (b - a) / n;
double sum = 0;
for(double x = a + h;x <= b - h;x += h)
sum += f(x);
return h * sum;
}
 
public static double midRect(double a, double b, double n){
double h = (b - a) / n;
double sum = 0;
for(double x = a;x <= b - h;x += h)
sum += (f(x) + f(x + h));
return (h / 2) * sum;
}
 
public static double trap(double a, double b, double n){
double h = (b - a) / n;
double sum = f(a) + f(b);
for(int i = 1;i < n;i++)
sum += f(a + i * h);
return h * sum;
}
 
public static double simpson(double a, double b, double n){
double h = (b - a) / n;
double sum1 = 0;
double sum2 = 0;
 
for(int i = 0;i < n;i++)
sum1 += f(a + h * i + h / 2);
 
for(int i = 1;i < n;i++)
sum2 += f(a + h * i);
 
return h / 6 * (f(a) + f(b) + 4 * sum1 + 2 * sum2);
}
//assume f(double x) elsewhere in the class
}

[edit] Logo

to i.left :fn :x :step
output invoke :fn :x
end
to i.right :fn :x :step
output invoke :fn :x + :step
end
to i.mid :fn :x :step
output invoke :fn :x + :step/2
end
to i.trapezium :fn :x :step
output ((i.left :fn :x :step) + (i.right :fn :x :step)) / 2
end
to i.simpsons :fn :x :step
output ( (i.left :fn :x :step)
+ (i.mid :fn :x :step) * 4
+ (i.right :fn :x :step) ) / 6
end
 
to integrate :method :fn :steps :a :b
localmake "step (:b - :a) / :steps
localmake "sigma 0
 ; for [x :a :b-:step :step] [make "sigma :sigma + apply :method (list :fn :x :step)]
repeat :steps [
make "sigma :sigma + (invoke :method :fn :a :step)
make "a :a + :step ]
output :sigma * :step
end
 
to fn2 :x
output 2 / (1 + 4 * :x * :x)
end
print integrate "i.left "fn2 4 -1 2  ; 2.456897
print integrate "i.right "fn2 4 -1 2  ; 2.245132
print integrate "i.mid "fn2 4 -1 2  ; 2.496091
print integrate "i.trapezium "fn2 4 -1 2  ; 2.351014
print integrate "i.simpsons "fn2 4 -1 2  ; 2.447732

[edit] MATLAB

For all of the examples given, the function that is passed to the method as parameter f is a function handle.

Function for performing left rectangular integration: leftRectIntegration.m

function integral = leftRectIntegration(f,a,b,n)
 
format long;
width = (b-a)/n; %calculate the width of each devision
x = linspace(a,b,n); %define x-axis
integral = width * sum( f(x(1:n-1)) );
 
end

Function for performing right rectangular integration: rightRectIntegration.m

function integral = rightRectIntegration(f,a,b,n)
 
format long;
width = (b-a)/n; %calculate the width of each devision
x = linspace(a,b,n); %define x-axis
integral = width * sum( f(x(2:n)) );
 
end

Function for performing mid-point rectangular integration: midPointRectIntegration.m

function integral = midPointRectIntegration(f,a,b,n)
 
format long;
width = (b-a)/n; %calculate the width of each devision
x = linspace(a,b,n); %define x-axis
integral = width * sum( f( (x(1:n-1)+x(2:n))/2 ) );
 
end

Function for performing trapezoidal integration: trapezoidalIntegration.m

function integral = trapezoidalIntegration(f,a,b,n)
 
format long;
x = linspace(a,b,n); %define x-axis
integral = trapz( x,f(x) );
 
end

Simpson's rule for numerical integration is already included in MATLAB as "quad()". It is not the same as the above examples, instead of specifying the amount of points to divide the x-axis into, the programmer passes the acceptable error tolerance for the calculation (parameter "tol").

integral = quad(f,a,b,tol)

Sample inputs

Using anonymous functions

trapezoidalIntegration(@(x)( exp(-(x.^2)) ),0,10,100000)
 
ans =
 
0.886226925452753

Using predefined functions

Built-in MATLAB function sin(x):

quad(@sin,0,pi,1/1000000000000)
 
ans =
 
2.000000000000000

User defined scripts and functions: fermiDirac.m

function answer = fermiDirac(x)
k = 8.617343e-5; %Boltazmann's Constant in eV/K
answer = 1./( 1+exp( (x)/(k*2000) ) ); %Fermi-Dirac distribution with mu = 0 and T = 2000K
end
 rightRectIntegration(@fermiDirac,-1,1,1000000)
 
ans =
 
0.999998006023282

[edit] OCaml

let integrate f a b steps meth =
let h = (b -. a) /. float_of_int steps in
let rec helper i s =
if i >= steps then s
else helper (succ i) (s +. meth f (a +. h *. float_of_int i) h)
in
h *. helper 0 0.
 
let left_rect f x _ =
f x
 
let mid_rect f x h =
f (x +. h /. 2.)
 
let right_rect f x h =
f (x +. h)
 
let trapezium f x h =
(f x +. f (x +. h)) /. 2.
 
let simpson f x h =
(f x +. 4. *. f (x +. h /. 2.) +. f (x +. h)) /. 6.
 
let square x = x *. x
 
 
let rl = integrate square 0. 1. 10 left_rect
let rm = integrate square 0. 1. 10 mid_rect
let rr = integrate square 0. 1. 10 right_rect
let t = integrate square 0. 1. 10 trapezium
let s = integrate square 0. 1. 10 simpson

[edit] Pascal

function RectLeft(function f(x: real): real; xl, xr: real): real;
begin
RectLeft := f(xl)
end;
 
function RectMid(function f(x: real): real; xl, xr: real) : real;
begin
RectMid := f((xl+xr)/2)
end;
 
function RectRight(function f(x: real): real; xl, xr: real): real;
begin
RectRight := f(xr)
end;
 
function Trapezium(function f(x: real): real; xl, xr: real): real;
begin
Trapezium := (f(xl) + f(xr))/2
end;
 
function Simpson(function f(x: real): real; xl, xr: real): real;
begin
Simpson := (f(xl) + 4*f((xl+xr)/2) + f(xr))/6
end;
 
function integrate(function method(function f(x: real): real; xl, xr: real): real;
function f(x: real): real;
a, b: real;
n: integer);
var
integral, h: real;
k: integer;
begin
integral := 0;
h := (b-a)/n;
for k := 0 to n-1 do
begin
integral := integral + method(f, a + k*h, a + (k+1)*h)
end;
integrate := integral
end;

[edit] PL/I

 
integrals: procedure options (main);
 
/* The function to be integrated */
f: procedure (x) returns (float);
declare x float;
return (3*x**2 + 2*x);
end f;
 
declare (a, b) float;
declare (rect_area, trap_area, Simpson) float;
declare (d, dx) fixed decimal (10,2);
declare (l, r) float;
declare (S1, S2) float;
 
l = 0; r = 5;
a = 0; b = 5; /* bounds of integration */
dx = 0.05;
 
/* Rectangle method */
rect_area = 0;
do d = a to b by dx;
rect_area = rect_area + dx*f(d);
end;
put skip data (rect_area);
 
/* trapezoid method */
trap_area = 0;
do d = a to b by dx;
trap_area = trap_area + dx*(f(d) + f(d+dx))/2;
end;
put skip data (trap_area);
 
/* Simpson's */
S1 = f(a+dx/2);
S2 = 0;
do d = a to b by dx;
S1 = S1 + f(d+dx+dx/2);
S2 = S2 + f(d+dx);
end;
Simpson = dx * (f(a) + f(b) + 4*S1 + 2*S2) / 6;
put skip data (Simpson);
 
end integrals;
 

[edit] PicoLisp

(scl 6)
 
(de leftRect (Fun X)
(Fun X) )
 
(de rightRect (Fun X H)
(Fun (+ X H)) )
 
(de midRect (Fun X H)
(Fun (+ X (/ H 2))) )
 
(de trapezium (Fun X H)
(/ (+ (Fun X) (Fun (+ X H))) 2) )
 
(de simpson (Fun X H)
(*/
(+
(Fun X)
(* 4 (Fun (+ X (/ H 2))))
(Fun (+ X H)) )
6 ) )
 
(de square (X)
(*/ X X 1.0) )
 
(de integrate (Fun From To Steps Meth)
(let (H (/ (- To From) Steps) Sum 0)
(for (X From (>= (- To H) X) (+ X H))
(inc 'Sum (Meth Fun X H)) )
(*/ H Sum 1.0) ) )
 
(prinl (round (integrate square 3.0 7.0 30 simpson) 3))

Output:

105.333

[edit] PureBasic

Translation of: Ada

Procedure.d leftRect(a.d, b.d, n.d)
Protected.d h, sum, x
 
h = (b - a) / n
x = a
While x <= b - h
sum + h * (f(x))
x + h
Wend
 
ProcedureReturn sum
EndProcedure
 
Procedure.d rightRect(a.d, b.d, n.d)
Protected.d h, sum, x
 
h = (b - a) / n
x = a + h
While x <= b - h
sum + h * (f(x))
x + h
Wend
 
ProcedureReturn sum
EndProcedure
 
Procedure.d midRect(a.d, b.d, n.d)
Protected.d h, sum, x
 
h = (b - a) / n
x = a
While x <= b - h
sum + (h / 2) * (f(x) + f(x + h))
x + h
Wend
 
ProcedureReturn sum
EndProcedure
 
Procedure.d trapezium(a.d, b.d, n.d)
Protected.d h, sum
Protected i
 
h = (b - a) / n
sum = f(a) + f(b)
For i = 1 To n - 1
sum + 2 * f((a + i * h))
Next
 
ProcedureReturn h / 2 * sum
EndProcedure
 
Procedure.d simpson(a.d, b.d, n.d)
Protected.d h, sum1, sum2
Protected i
 
h = (b - a) / n
 
For i = 0 To n - 1
sum1 + f(a + h * i + h / 2)
Next
 
For i = 1 To n - 1
sum2 + f(a + h * i)
Next
 
ProcedureReturn h / 6 * (f(a) + f(b) + 4 * sum1 + 2 * sum2)
EndProcedure

[edit] Python

def left_rect(f,x,h):
return f(x)
 
def mid_rect(f,x,h):
return f(x + h/2)
 
def right_rect(f,x,h):
return f(x+h)
 
def trapezium(f,x,h):
return (f(x) + f(x+h))/2.0
 
def simpson(f,x,h):
return (f(x) + 4*f(x + h/2) + f(x+h))/6.0
 
def square(x):
return x*x
 
def integrate( f, a, b, steps, meth):
h = (b-a)/steps
ival = h * sum(meth(f, a+i*h, h) for i in xrange(steps))
return ival
 
t = integrate( square, 3.0, 7.0, 30, simpson )

A faster Simpson's rule integrator is

def faster_simpson(f, a, b, steps):
h = (b-a)/steps
a1 = a+h/2
s1 = sum( f(a1+i*h) for i in range(0,steps))
s2 = sum( f(a+i*h) for i in range(1,steps))
return (h/6.0)*(f(a)+f(b)+4.0*s1+2.0*s2)

[edit] Ruby

Translation of: Tcl

def leftrect(f, left, right)
f.call(left)
end
 
def midrect(f, left, right)
f.call((left+right)/2.0)
end
 
def rightrect(f, left, right)
f.call(right)
end
 
def trapezium(f, left, right)
(f.call(left) + f.call(right)) / 2.0
end
 
def simpson(f, left, right)
(f.call(left) + 4*f.call((left+right)/2.0) + f.call(right)) / 6.0
end
 
def integrate(f, a, b, steps, method)
delta = 1.0 * (b - a) / steps
total = 0.0
steps.times do |i|
left = a + i*delta
right = left + delta
total += delta * send(method, f, left, right)
end
total
end
 
def square(x)
x**2
end
 
def def_int(f, a, b)
l = case f.to_s
when /sin>/
lambda {|x| -Math.cos(x)}
when /square>/
lambda {|x| (x**3)/3.0}
end
l.call(b) - l.call(a)
end
 
a = 0
b = Math::PI
steps = 10
 
for func in [method(:square), Math.method(:sin)]
puts "integral of #{func} from #{a} to #{b} in #{steps} steps"
actual = def_int(func, a, b)
for method in [:leftrect, :midrect, :rightrect, :trapezium, :simpson]
int = integrate(func, a, b, steps, method)
diff = (int - actual) * 100.0 / actual
printf "  %-10s  %s\t(%.1f%%)\n", method, int, diff
end
end

outputs

integral of #<Method: Object#square> from 0 to 3.14159265358979 in 10 steps
   leftrect    8.83678885388545	(-14.5%)
   midrect     10.3095869961997	(-0.2%)
   rightrect   11.9374165219154	(15.5%)
   trapezium   10.3871026879004	(0.5%)
   simpson     10.3354255600999	(0.0%)
integral of #<Method: Math.sin> from 0 to 3.14159265358979 in 10 steps
   leftrect    1.98352353750945	(-0.8%)
   midrect     2.00824840790797	(0.4%)
   rightrect   1.98352353750945	(-0.8%)
   trapezium   1.98352353750945	(-0.8%)
   simpson     2.0000067844418	(0.0%)

[edit] Scheme

(define (integrate f a b steps meth)
(define h (/ (- b a) steps))
(* h
(let loop ((i 0) (s 0))
(if (>= i steps)
s
(loop (+ i 1) (+ s (meth f (+ a (* h i)) h)))))))
 
(define (left-rect f x h) (f x))
(define (mid-rect f x h) (f (+ x (/ h 2))))
(define (right-rect f x h) (f (+ x h)))
(define (trapezium f x h) (/ (+ (f x) (f (+ x h))) 2))
(define (simpson f x h) (/ (+ (f x) (* 4 (f (+ x (/ h 2)))) (f (+ x h))) 6))
 
(define (square x) (* x x))
 
(define rl (integrate square 0 1 10 left-rect))
(define rm (integrate square 0 1 10 mid-rect))
(define rr (integrate square 0 1 10 right-rect))
(define t (integrate square 0 1 10 trapezium))
(define s (integrate square 0 1 10 simpson))

[edit] Standard ML

fun integrate (f, a, b, steps, meth) = let
val h = (b - a) / real steps
fun helper (i, s) =
if i >= steps then s
else helper (i+1, s + meth (f, a + h * real i, h))
in
h * helper (0, 0.0)
end
 
fun leftRect (f, x, _) = f x
fun midRect (f, x, h) = f (x + h / 2.0)
fun rightRect (f, x, h) = f (x + h)
fun trapezium (f, x, h) = (f x + f (x + h)) / 2.0
fun simpson (f, x, h) = (f x + 4.0 * f (x + h / 2.0) + f (x + h)) / 6.0
 
fun square x = x * x
 
 
val rl = integrate (square, 0.0, 1.0, 10, left_rect )
val rm = integrate (square, 0.0, 1.0, 10, mid_rect )
val rr = integrate (square, 0.0, 1.0, 10, right_rect)
val t = integrate (square, 0.0, 1.0, 10, trapezium )
val s = integrate (square, 0.0, 1.0, 10, simpson )

[edit] Tcl

package require Tcl 8.5
 
proc leftrect {f left right} {
$f $left
}
proc midrect {f left right} {
set mid [expr {($left + $right) / 2.0}]
$f $mid
}
proc rightrect {f left right} {
$f $right
}
proc trapezium {f left right} {
expr {([$f $left] + [$f $right]) / 2.0}
}
proc simpson {f left right} {
set mid [expr {($left + $right) / 2.0}]
expr {([$f $left] + 4*[$f $mid] + [$f $right]) / 6.0}
}
 
proc integrate {f a b steps method} {
set delta [expr {1.0 * ($b - $a) / $steps}]
set total 0.0
for {set i 0} {$i < $steps} {incr i} {
set left [expr {$a + $i * $delta}]
set right [expr {$left + $delta}]
set total [expr {$total + $delta * [$method $f $left $right]}]
}
return $total
}
 
interp alias {} sin {} ::tcl::mathfunc::sin
proc square x {expr {$x*$x}}
proc def_int {f a b} {
switch -- $f {
sin {set lambda {x {expr {-cos($x)}}}}
square {set lambda {x {expr {$x**3/3.0}}}}
}
return [expr {[apply $lambda $b] - [apply $lambda $a]}]
}
 
set a 0
set b [expr {4*atan(1)}]
set steps 10
 
foreach func {square sin} {
puts "integral of ${func}(x) from $a to $b in $steps steps"
set actual [def_int $func $a $b]
foreach method {leftrect midrect rightrect trapezium simpson} {
set int [integrate $func $a $b $steps $method]
set diff [expr {($int - $actual) * 100.0 / $actual}]
puts [format "  %-10s  %s\t(%.1f%%)" $method $int $diff]
}
}
integral of square(x) from 0 to 3.141592653589793 in 10 steps
   leftrect    8.836788853885448	(-14.5%)
   midrect     10.30958699619969	(-0.2%)
   rightrect   11.93741652191543	(15.5%)
   trapezium   10.387102687900438	(0.5%)
   simpson     10.335425560099939	(0.0%)
integral of sin(x) from 0 to 3.141592653589793 in 10 steps
   leftrect    1.9835235375094544	(-0.8%)
   midrect     2.0082484079079745	(0.4%)
   rightrect   1.9835235375094544	(-0.8%)
   trapezium   1.9835235375094546	(-0.8%)
   simpson     2.0000067844418012	(0.0%)

[edit] TI-89 BASIC

TI-89 BASIC has built-in numerical integration with the ∫ operator, but no control over the method used is available so it doesn't really correspond to this task.

\overbrace{{\textstyle\int}(\mathrm{expr},x,a,b)}^{\mathrm{TI-89}} = \overbrace{\int_{a}^b \mathrm{expr}\ dx}^{\text{Math notation}}

An explicit numerical integration program should be written here.

[edit] Ursala

A higher order function parameterized by a method m returns a function that integrates by that method. The method m is meant to specify whether it's rectangular, trapezoidal, etc.. The integrating function constructed from a given method takes a quadruple (f,a,b,n) containing the integrand f, the bounds (a,b), and the number of intervals n.

#import std
#import nat
#import flo
 
(integral_by "m") ("f","a","b","n") =
 
iprod ^(* ! div\float"n" minus/"b" "a",~&) ("m" "f")*ytp (ari successor "n")/"a" "b"

An alternative way of defining this function shown below prevents redundant evaluations of the integrand at the cost of building a table-driven finite map in advance.

(integral_by "m") ("f","a","b","n") =
 
iprod ^(* ! div\float"n" minus/"b" "a",~&) ^H(*+ "m"+ -:"f"+ * ^/~& "f",~&ytp) (ari successor "n")/"a" "b"

As mentioned in the Haskell solution, the latter choice is preferable if evaluating the integrand is expensive. An integrating function is defined for each method as follows.

left       = integral_by "f". ("l","r"). "f" "l"
right = integral_by "f". ("l","r"). "f" "r"
midpoint = integral_by "f". ("l","r"). "f" div\2. plus/"l" "r"
trapezium = integral_by "f". ("l","r"). div\2. plus "f"~~/"l" "r"
simpson = integral_by "f". ("l","r"). div\6. plus:-0. <"f" "l",times/4. "f" div\2. plus/"l" "r","f" "r">

As shown above, the method passed to the integral_by function is itself a higher order function taking an integrand f as an argument and returning a function that operates on the pair of left and right interval enpoints. Here is a test program showing the results of integrating the square from zero to π in ten intervals by all five methods.

#cast %eL
 
examples = <.left,midpoint,rignt,trapezium,simpson> (sqr,0.,pi,10)

output:

<
   8.836789e+00,
   1.030959e+01,
   1.193742e+01,
   1.038710e+01,
   1.033543e+01>

(The GNU Scientific Library integration routines are also callable in Ursala, and are faster and more accurate.)

Personal tools