Numerical integration: Difference between revisions
(add E example) |
(added standard ml) |
||
Line 18: | Line 18: | ||
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. |
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> |
<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 |
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> |
|||
=={{header|ALGOL 68}}== |
=={{header|ALGOL 68}}== |
||
<lang algol68> |
<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 - 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</lang> |
|||
=={{header|BASIC}}== |
=={{header|BASIC}}== |
||
Line 205: | Line 204: | ||
{{trans|Java}} |
{{trans|Java}} |
||
<lang qbasic> |
<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> |
|||
=={{header|C}}== |
=={{header|C}}== |
||
Line 391: | Line 390: | ||
Due to their similarity, it makes sense to make the integration method a policy. |
Due to their similarity, it makes sense to make the integration method a policy. |
||
<lang cpp> |
<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> |
|||
=={{header|D}}== |
=={{header|D}}== |
||
Line 623: | Line 622: | ||
=={{header|Forth}}== |
=={{header|Forth}}== |
||
<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* ; |
|||
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 |
\ 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> |
|||
=={{header|Fortran}}== |
=={{header|Fortran}}== |
||
Line 713: | Line 711: | ||
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: |
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> |
<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: |
Similarly for the closed formulas, but we need an additional function ''overlap'' which sums the coefficients overlapping at the interior interval boundaries: |
||
<lang haskell> |
<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 |
And now we can just define |
||
<lang haskell> |
<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: |
or, as easily, some additional schemes: |
||
<lang haskell> |
<lang haskell>intMilne = integrateClosed 45 [14,64,24,64,14] |
||
intOpen1 = integrateOpen 2 [3,3] |
|||
intOpen2 = integrateOpen 3 [8,-4,8]</lang> |
|||
Some examples: |
Some examples: |
||
Line 767: | Line 765: | ||
=={{header|Java}}== |
=={{header|Java}}== |
||
The function in this example is assumed to be <tt>f(double x)</tt>. |
The function in this example is assumed to be <tt>f(double x)</tt>. |
||
<lang java5> |
<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> |
|||
=={{header|Logo}}== |
=={{header|Logo}}== |
||
<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> |
|||
=={{header|OCaml}}== |
=={{header|OCaml}}== |
||
<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> |
|||
=={{header|Pascal}}== |
=={{header|Pascal}}== |
||
Line 970: | Line 968: | ||
=={{header|Scheme}}== |
=={{header|Scheme}}== |
||
<lang 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> |
|||
=={{header|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> |
Revision as of 01:11, 26 March 2009
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 = sum2 = 0 loop on i from 0 to (n - 1) sum1 = sum1 + f(a + h * i + h / 2) loop on i from 1 to (n - 1) 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.
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 - 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</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>
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 <lang fortran> program integrate
integer, parameter :: n = 20 ! or whatever real, parameter :: a = 0.0, b = 15.0 ! or whatever real, parameter :: h = (b - a) / n real, parameter, dimension(0:2*n) :: xpoints = (/ (a + h*i/(n*2), i = 0, 2*n) /) real, dimension(0:2*n), target :: fpoints ! gather up all the f(x) values needed for all methods once and distribute via pointers real, dimension(:), pointer :: fleft, fmid, fright real :: leftrect, midrect, rightrect, trapezoid, simpson fpoints = elemf(xpoints) ! elemental function wrapper for F runs once for each element of XPOINTS fleft => fpoints(0 : 2*n-2 : 2) fmid => fpoints(1 : 2*n-1 : 2) fright => fpoints(2 : 2*n : 2) ! Left rectangular rule integral leftrect = h * sum(fleft) ! Middle rectangular rule integral midrect = h * sum(fmid) ! Right rectangular rule integral rightrect = h * sum(fright) ! Trapezoid rule integral trapezoid = h / 2 * sum(fleft + fright) ! Simpson's rule integral simpson = h / 6 * sum(fleft + fright + 4*fmid) end program integrate</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
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>
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>
Python
<lang 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 ) </lang> 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>
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>