Numerical integration: Difference between revisions
(Undo revision 90984 by TimToady Integration is valid over the range from zero. (and I've already calculated that :-)) |
|||
| Line 18: | Line 18: | ||
* f(x) = x^3, where x is [0,1], with 100 approximations. The exact result is 1/4, or 0.25. |
* f(x) = x^3, where x is [0,1], with 100 approximations. The exact result is 1/4, or 0.25. |
||
* f(x) = 1/x, where x is [1,100], with 1,000 approximations. The exact result is the natural log of 100, or about 4.605170 |
* f(x) = 1/x, where x is [1,100], with 1,000 approximations. The exact result is the natural log of 100, or about 4.605170 |
||
* f(x) = x, where x is [ |
* f(x) = x, where x is [0,5000], with 5,000,000 approximations. The exact result is 12,500,000. |
||
* f(x) = x, where x is [ |
* f(x) = x, where x is [0,6000], with 6,000,000 approximations. The exact result is 18,000,000. |
||
Revision as of 19:31, 12 September 2010
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)
Demonstrate your function by showing the results for:
- f(x) = x^3, where x is [0,1], with 100 approximations. The exact result is 1/4, or 0.25.
- f(x) = 1/x, where x is [1,100], with 1,000 approximations. The exact result is the natural log of 100, or about 4.605170
- f(x) = x, where x is [0,5000], with 5,000,000 approximations. The exact result is 12,500,000.
- f(x) = x, where x is [0,6000], with 6,000,000 approximations. The exact result is 18,000,000.
See also
- Active object for integrating a function of real time.
ActionScript
Integration functions: <lang ActionScript>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); }</lang> Usage: <lang ActionScript>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 )); </lang>
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.
<lang ada>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;</lang>
ALGOL 68
<lang algol68>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 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 * f(x + h / 2);
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</lang>
AutoHotkey
ahk discussion <lang autohotkey>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
}</lang>
BASIC
<lang qbasic>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</lang>
C
<lang c>#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);
}</lang>
Usage example:
<lang c>#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");
}
}</lang>
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
C++
Due to their similarity, it makes sense to make the integration method a policy. <lang cpp>// 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());</lang>
Common Lisp
<lang 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))))))</lang>
D
The function f(x) is Func1/Func2 in this program. <lang d>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) ;
}</lang> 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
E
<lang e>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) }
}</lang>
<lang e>? 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</lang>
Forth
<lang 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</lang>
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: <lang fortran>elemental function elemf(x)
real :: elemf, x elemf = f(x)
end function elemf</lang>
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. <lang fortran>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</lang>
Usage example: <lang fortran>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</lang>
The FunctionHolder module:
<lang fortran>module FunctionHolder
implicit none
contains
pure function afun(x) real :: afun real, intent(in) :: x
afun = x**2 end function afun
end module FunctionHolder</lang>
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
<lang haskell>approx f xs ws = sum [w * f x | (x,w) <- zip xs ws]</lang>
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:
<lang haskell>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]]</lang>
Similarly for the closed formulas, but we need an additional function overlap which sums the coefficients overlapping at the interior interval boundaries:
<lang haskell>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</lang>
And now we can just define
<lang haskell>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]</lang>
or, as easily, some additional schemes:
<lang haskell>intMilne = integrateClosed 45 [14,64,24,64,14] intOpen1 = integrateOpen 2 [3,3] intOpen2 = integrateOpen 3 [8,-4,8]</lang>
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
J
Solution: <lang j>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'</lang>
Example usage
Integrate square (*:) from 0 to π in 10 steps using various methods.
<lang j> *: rectangle integrate 0 1p1 10
10.3095869962
*: trapezium integrate 0 1p1 10
10.3871026879
*: simpson integrate 0 1p1 10
10.3354255601</lang>
Integrate sin from 0 to π in 10 steps using various methods.
<lang j> 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</lang>
Note that J has a primitive verb p.. for integrating polynomials. For example the integral of (which can be described in terms of its coefficients as 0 0 1) is:
<lang j> 0 p.. 0 0 1
0 0 0 0.333333333333
0 p.. 0 0 1x NB. or using rationals
0 0 0 1r3</lang>
That is:
So to integrate from 0 to π :
<lang j> 0 0 1 (0&p..@[ -~/@:p. ]) 0 1p1
10.3354255601</lang>
That said, J also has d. which can integrate suitable functions.
<lang j> *:d._1]1p1 10.3354</lang>
Java
The function in this example is assumed to be f(double x). <lang java5>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
}</lang>
Logo
<lang 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</lang>
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 <lang MATLAB>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</lang>
Function for performing right rectangular integration: rightRectIntegration.m <lang MATLAB>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</lang>
Function for performing mid-point rectangular integration: midPointRectIntegration.m <lang MATLAB>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</lang>
Function for performing trapezoidal integration: trapezoidalIntegration.m <lang MATLAB>function integral = trapezoidalIntegration(f,a,b,n)
format long; x = linspace(a,b,n); %define x-axis integral = trapz( x,f(x) );
end</lang>
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"). <lang MATLAB>integral = quad(f,a,b,tol)</lang>
Using anonymous functions
<lang MATLAB>trapezoidalIntegration(@(x)( exp(-(x.^2)) ),0,10,100000)
ans =
0.886226925452753</lang>
Using predefined functions
Built-in MATLAB function sin(x): <lang MATLAB>quad(@sin,0,pi,1/1000000000000)
ans =
2.000000000000000</lang>
User defined scripts and functions: fermiDirac.m <lang MATLAB>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</lang>
<lang MATLAB> rightRectIntegration(@fermiDirac,-1,1,1000000)
ans =
0.999998006023282</lang>
OCaml
<lang 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</lang>
Pascal
<lang 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;</lang>
Perl 6
<lang perl6>sub leftrect(&f, $a, $b, $n) {
my $h = ($b - $a) / $n; $h * [+] do f($_) for $a, *+$h ... $b-$h;
}
sub rightrect(&f, $a, $b, $n) {
my $h = ($b - $a) / $n; $h * [+] do f($_) for $a+$h, *+$h ... $b;
}
sub midrect(&f, $a, $b, $n) {
my $h = ($b - $a) / $n; $h * [+] do f($_) for $a+$h/2, *+$h ... $b-$h/2;
}
sub trapez(&f, $a, $b, $n) {
my $h = ($b - $a) / $n; $h / 2 * [+] f($a), f($b), do f($_) * 2 for $a+$h, *+$h ... $b-$h;
}
sub simpsons(&f, $a, $b, $n) {
my $h = ($b - $a) / $n;
my $sum1 = f($a + $h/2);
my $sum2 = 0;
for $a+$h, *+$h ... $b-$h {
$sum1 += f($_ + $h/2);
$sum2 += f($_);
}
($h / 6) * (f($a) + f($b) + 4*$sum1 + 2*$sum2);
}
sub tryem(&f, $a, $b, $n) {
say "rectangle method left: ", leftrect &f, $a, $b, $n; say "rectangle method right: ", rightrect &f, $a, $b, $n; say "rectangle method mid: ", midrect &f, $a, $b, $n; say "composite trapezoidal rule: ", trapez &f, $a, $b, $n; say "quadratic simpsons rule: ", simpsons &f, $a, $b, $n;
}
say '-> $x { 2 / (1 + 4 * $x ** 2) }'; tryem -> $x { 2 / (1 + 4 * $x * $x) }, -1, 2, 4;
say "\n",'{ $_ ** 3 }'; tryem { $_ * $_ * $_ }, 0, 1, 100; # rakudo doesn't do Rat ** Int --> Rat yet
say "\n",'1 / *'; tryem 1 / *, 1, 100, 1000;</lang> Output: <lang>-> $x { 2 / (1 + 4 * $x ** 2) } rectangle method left: 2.45689655172414 rectangle method right: 2.24513184584178 rectangle method mid: 2.49609105850139 composite trapezoidal rule: 2.35101419878296 quadratic simpsons rule: 2.44773210526191
{ $_ ** 3 } rectangle method left: 0.245025 rectangle method right: 0.255025 rectangle method mid: 0.2499875 composite trapezoidal rule: 0.250025 quadratic simpsons rule: 0.25
1 / * rectangle method left: 4.65499105751468 rectangle method right: 4.55698105751468 rectangle method mid: 4.60476254867838 composite trapezoidal rule: 4.60598605751468 quadratic simpsons rule: 4.60517038495714</lang>
Note that these integrations are done with rationals rather than floats, so should be fairly precise (though of course they are not terribly accurate with so few iterations). Some of the sums do overflow into Num (floating point)--currently rakudo allows implements Rat32--but at least all of the interval arithmetic is exact.
PL/I
<lang 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; </lang>
PicoLisp
<lang 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))</lang> Output:
105.333
PureBasic
<lang PureBasic>Prototype.d TestFunction(Arg.d)
Procedure.d LeftIntegral(Start, Stop, Steps, *func.TestFunction)
Protected.d n=(Stop-Start)/Steps, sum, x=Start While x <= Stop-n sum + n * *func(x) x + n Wend ProcedureReturn sum
EndProcedure
Procedure.d MidIntegral(Start, Stop, Steps, *func.TestFunction)
Protected.d n=(Stop-Start)/Steps, sum, x=Start While x <= Stop-n sum + n * *func(x+n/2) x + n Wend ProcedureReturn sum
EndProcedure
Procedure.d RightIntegral(Start, Stop, Steps, *func.TestFunction)
Protected.d n=(Stop-Start)/Steps, sum, x=Start While x < Stop x + n sum + n * *func(x) Wend ProcedureReturn sum
EndProcedure
Procedure.d Trapezium(Start, Stop, Steps, *func.TestFunction)
Protected.d n=(Stop-Start)/Steps, sum, x=Start While x<=Stop sum + n * (*func(x) + *func(x+n))/2 x+n Wend ProcedureReturn sum
EndProcedure
Procedure.d Simpson(Start, Stop, Steps, *func.TestFunction)
Protected.d n=(Stop-Start)/Steps, sum1, sum2, x=Start Protected i For i=0 To steps-1 sum1+ *func(Start+n*i+n/2) Next For i=1 To Steps-1 sum2+ *func(Start+n*i) Next ProcedureReturn n * (*func(Start)+ *func(Stop)+4*sum1+2*sum2) / 6
EndProcedure
- - Set up functions to integrate
Procedure.d Test1(n.d)
ProcedureReturn n*n*n
EndProcedure
Procedure.d Test2(n.d)
ProcedureReturn 1/n
EndProcedure
- This function should be integrated as a integer function, but for
- comparably this will stay as a float.
Procedure.d Test3(n.d)
ProcedureReturn n
EndProcedure
- - Test the code & present the results
CompilerIf #PB_Compiler_Debugger
MessageRequester("Notice!","Running this program in Debug-mode will be slow")
CompilerEndIf
- = 0.25
Define Answer$ Answer$="Left ="+StrD(LeftIntegral (0,1,100,@Test1()))+#CRLF$ Answer$+"Mid ="+StrD(MidIntegral (0,1,100,@Test1()))+#CRLF$ Answer$+"Right ="+StrD(RightIntegral(0,1,100,@Test1()))+#CRLF$ Answer$+"Trapezium="+StrD(Trapezium (0,1,100,@Test1()))+#CRLF$ Answer$+"Simpson ="+StrD(Simpson (0,1,100,@Test1())) MessageRequester("Answer should be 1/4",Answer$)
- = Ln(100) e.g. ~4.60517019...
Answer$="Left ="+StrD(LeftIntegral (1,100,1000,@Test2()))+#CRLF$ Answer$+"Mid ="+StrD(MidIntegral (1,100,1000,@Test2()))+#CRLF$ Answer$+"Right ="+StrD(RightIntegral (1,100,1000,@Test2()))+#CRLF$ Answer$+"Trapezium="+StrD(Trapezium (1,100,1000,@Test2()))+#CRLF$ Answer$+"Simpson ="+StrD(Simpson (1,100,1000,@Test2())) MessageRequester("Answer should be Ln(100), e.g. ~4.60517019",Answer$)
- 12,500,000
Answer$="Left ="+StrD(LeftIntegral (0,5000,5000000,@Test3()))+#CRLF$ Answer$+"Mid ="+StrD(MidIntegral (0,5000,5000000,@Test3()))+#CRLF$ Answer$+"Right ="+StrD(RightIntegral (0,5000,5000000,@Test3()))+#CRLF$ Answer$+"Trapezium="+StrD(Trapezium (0,5000,5000000,@Test3()))+#CRLF$ Answer$+"Simpson ="+StrD(Simpson (0,5000,5000000,@Test3())) MessageRequester("Answer should be 12,500,000",Answer$)
- 18,000,000
Answer$="Left ="+StrD(LeftIntegral (0,6000,6000000,@Test3()))+#CRLF$ Answer$+"Mid ="+StrD(MidIntegral (0,6000,6000000,@Test3()))+#CRLF$ Answer$+"Right ="+StrD(RightIntegral (0,6000,6000000,@Test3()))+#CRLF$ Answer$+"Trapezium="+StrD(Trapezium (0,6000,6000000,@Test3()))+#CRLF$ Answer$+"Simpson ="+StrD(Simpson (0,6000,6000000,@Test3())) MessageRequester("Answer should be 18,000,000",Answer$) </lang>
Left =0.2353220100 Mid =0.2401367513 Right =0.2550250000 Trapezium=0.2500250000 Simpson =0.2500000000 Left =4.6540000764 Mid =4.6037720584 Right =4.5569810575 Trapezium=4.6059860575 Simpson =4.6051703850 Left =12499992.5007297550 Mid =12499995.0007292630 Right =12500002.5007287540 Trapezium=12500000.0007287620 Simpson =12500000.0000000000 Left =17999991.0013914930 Mid =17999994.0013910230 Right =18000003.0013904940 Trapezium=18000000.0013905240 Simpson =17999999.9999999960
Python
Answers are first given using floating point arithmatic, then using fractions, only converted to floating point on output. <lang python>from fractions import Fraction
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 cube(x):
return x*x*x
def reciprocal(x):
return 1/x
def identity(x):
return x
def integrate( f, a, b, steps, meth):
h = (b-a)/steps ival = h * sum(meth(f, a+i*h, h) for i in range(steps)) return ival</lang>
Tests <lang python>for a, b, steps, func in ((0., 1., 100, cube), (1., 100., 1000, reciprocal)):
for rule in (left_rect, mid_rect, right_rect, trapezium, simpson):
print('%s integrated using %s\n from %r to %r (%i steps) = %r' %
(func.__name__, rule.__name__, a, b, steps,
integrate( func, a, b, steps, rule)))
a, b = Fraction.from_float(a), Fraction.from_float(b)
for rule in (left_rect, mid_rect, right_rect, trapezium, simpson):
print('%s integrated using %s\n from %r to %r (%i steps and fractions) = %r' %
(func.__name__, rule.__name__, a, b, steps,
float(integrate( func, a, b, steps, rule))))
- Extra tests (compute intensive)
for a, b, steps, func in ((1., 5000., 5000000, identity),
(1., 6000., 6000000, identity)):
for rule in (left_rect, mid_rect, right_rect, trapezium, simpson):
print('%s integrated using %s\n from %r to %r (%i steps) = %r' %
(func.__name__, rule.__name__, a, b, steps,
integrate( func, a, b, steps, rule)))
a, b = Fraction.from_float(a), Fraction.from_float(b)
for rule in (left_rect, mid_rect, right_rect, trapezium, simpson):
print('%s integrated using %s\n from %r to %r (%i steps and fractions) = %r' %
(func.__name__, rule.__name__, a, b, steps,
float(integrate( func, a, b, steps, rule))))</lang>
Sample test Output
cube integrated using left_rect from 0.0 to 1.0 (100 steps) = 0.24502500000000005 cube integrated using mid_rect from 0.0 to 1.0 (100 steps) = 0.24998750000000006 cube integrated using right_rect from 0.0 to 1.0 (100 steps) = 0.25502500000000006 cube integrated using trapezium from 0.0 to 1.0 (100 steps) = 0.250025 cube integrated using simpson from 0.0 to 1.0 (100 steps) = 0.25 cube integrated using left_rect from Fraction(0, 1) to Fraction(1, 1) (100 steps and fractions) = 0.245025 cube integrated using mid_rect from Fraction(0, 1) to Fraction(1, 1) (100 steps and fractions) = 0.2499875 cube integrated using right_rect from Fraction(0, 1) to Fraction(1, 1) (100 steps and fractions) = 0.255025 cube integrated using trapezium from Fraction(0, 1) to Fraction(1, 1) (100 steps and fractions) = 0.250025 cube integrated using simpson from Fraction(0, 1) to Fraction(1, 1) (100 steps and fractions) = 0.25 reciprocal integrated using left_rect from 1.0 to 100.0 (1000 steps) = 4.65499105751468 reciprocal integrated using mid_rect from 1.0 to 100.0 (1000 steps) = 4.604762548678376 reciprocal integrated using right_rect from 1.0 to 100.0 (1000 steps) = 4.55698105751468 reciprocal integrated using trapezium from 1.0 to 100.0 (1000 steps) = 4.605986057514676 reciprocal integrated using simpson from 1.0 to 100.0 (1000 steps) = 4.605170384957133 reciprocal integrated using left_rect from Fraction(1, 1) to Fraction(100, 1) (1000 steps and fractions) = 4.654991057514676 reciprocal integrated using mid_rect from Fraction(1, 1) to Fraction(100, 1) (1000 steps and fractions) = 4.604762548678376 reciprocal integrated using right_rect from Fraction(1, 1) to Fraction(100, 1) (1000 steps and fractions) = 4.556981057514676 reciprocal integrated using trapezium from Fraction(1, 1) to Fraction(100, 1) (1000 steps and fractions) = 4.605986057514677 reciprocal integrated using simpson from Fraction(1, 1) to Fraction(100, 1) (1000 steps and fractions) = 4.605170384957134 identity integrated using left_rect from 1.0 to 5000.0 (5000000 steps) = 12499997.0009999 identity integrated using mid_rect from 1.0 to 5000.0 (5000000 steps) = 12499999.499999998 identity integrated using right_rect from 1.0 to 5000.0 (5000000 steps) = 12500001.999000099 identity integrated using trapezium from 1.0 to 5000.0 (5000000 steps) = 12499999.499999998 identity integrated using simpson from 1.0 to 5000.0 (5000000 steps) = 12499999.499999998 identity integrated using left_rect from Fraction(1, 1) to Fraction(5000, 1) (5000000 steps and fractions) = 12499997.0009999 identity integrated using mid_rect from Fraction(1, 1) to Fraction(5000, 1) (5000000 steps and fractions) = 12499999.5 identity integrated using right_rect from Fraction(1, 1) to Fraction(5000, 1) (5000000 steps and fractions) = 12500001.9990001 identity integrated using trapezium from Fraction(1, 1) to Fraction(5000, 1) (5000000 steps and fractions) = 12499999.499999998 identity integrated using simpson from Fraction(1, 1) to Fraction(5000, 1) (5000000 steps and fractions) = 12499999.499999998 identity integrated using left_rect from 1.0 to 6000.0 (6000000 steps) = 17999996.50099991 identity integrated using mid_rect from 1.0 to 6000.0 (6000000 steps) = 17999999.500000004 identity integrated using right_rect from 1.0 to 6000.0 (6000000 steps) = 18000002.499000076 identity integrated using trapezium from 1.0 to 6000.0 (6000000 steps) = 17999999.500000004 identity integrated using simpson from 1.0 to 6000.0 (6000000 steps) = 17999999.500000004 identity integrated using left_rect from Fraction(1, 1) to Fraction(6000, 1) (6000000 steps and fractions) = 17999996.500999916 identity integrated using mid_rect from Fraction(1, 1) to Fraction(6000, 1) (6000000 steps and fractions) = 17999999.5 identity integrated using right_rect from Fraction(1, 1) to Fraction(6000, 1) (6000000 steps and fractions) = 18000002.499000084 identity integrated using trapezium from Fraction(1, 1) to Fraction(6000, 1) (6000000 steps and fractions) = 17999999.500000004 identity integrated using simpson from Fraction(1, 1) to Fraction(6000, 1) (6000000 steps and fractions) = 17999999.500000004
A faster Simpson's rule integrator is <lang python>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)</lang>
Ruby
<lang ruby>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</lang> 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%)
Scheme
<lang 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))</lang>
Standard ML
<lang sml>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 )</lang>
Tcl
<lang 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]
}
}</lang>
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%)
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.
An explicit numerical integration program should be written here.
Ursala
A higher order function parameterized by a method returns a function that integrates by that method. The method is meant to specify whether it's rectangular, trapezoidal, etc.. The integrating function constructed from a given method takes a quadruple containing the integrand , the bounds , and the number of intervals .
<lang Ursala>#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"</lang> 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. <lang Ursala>(integral_by "m") ("f","a","b","n") =
iprod ^(* ! div\float"n" minus/"b" "a",~&) ^H(*+ "m"+ -:"f"+ * ^/~& "f",~&ytp) (ari successor "n")/"a" "b"</lang>
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.
<lang Ursala>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"></lang>
As shown above, the method passed to the integral_by function
is itself a higher order function taking an integrand 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.
<lang Ursala>#cast %eL
examples = <.left,midpoint,rignt,trapezium,simpson> (sqr,0.,pi,10)</lang> 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.)
- Programming Tasks
- Arithmetic operations
- ActionScript
- Ada
- Ada examples needing attention
- Examples needing attention
- ALGOL 68
- AutoHotkey
- BASIC
- BASIC examples needing attention
- C
- C examples needing attention
- C++
- Common Lisp
- D
- E
- Forth
- Fortran
- Haskell
- J
- Java
- Java examples needing attention
- Logo
- MATLAB
- Sample inputs
- OCaml
- Pascal
- Perl 6
- PL/I
- PicoLisp
- PureBasic
- Python
- Ruby
- Scheme
- Standard ML
- Tcl
- TI-89 BASIC examples needing attention
- Ursala
- M4/Omit