Numerical integration: Difference between revisions
Completed Haskell version |
categories added |
||
Line 3,064: | Line 3,064: | ||
{{omit from|M4}} |
{{omit from|M4}} |
||
[[Category:Arithmetic]] |
|||
[[Category:Mathematics]] |
Revision as of 15:10, 26 April 2011
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
Specification of a generic package implementing the five specified kinds of numerical integration: <lang ada>generic
type Scalar is digits <>; with function F (X : Scalar) return Scalar;
package Integrate is
function Left_Rectangular (A, B : Scalar; N : Positive) return Scalar; function Right_Rectangular (A, B : Scalar; N : Positive) return Scalar; function Midpoint_Rectangular (A, B : Scalar; N : Positive) return Scalar; function Trapezium (A, B : Scalar; N : Positive) return Scalar; function Simpsons (A, B : Scalar; N : Positive) return Scalar;
end Integrate;</lang> An alternative solution is to pass a function reference to the integration function. This solution is probably slightly faster, and works even with Ada83. One could also make each integration function generic, instead of making the whole package generic.
Body of the package implementing numerical integration: <lang ada>package body Integrate is
function Left_Rectangular (A, B : Scalar; N : Positive) return Scalar is H : constant Scalar := (B - A) / Scalar (N); Sum : Scalar := 0.0; X : Scalar; begin for I in 0 .. N - 1 loop X := A + Scalar (I) * H; Sum := Sum + H * F (X); end loop; return Sum; end Left_Rectangular;
function Right_Rectangular (A, B : Scalar; N : Positive) return Scalar is H : constant Scalar := (B - A) / Scalar (N); Sum : Scalar := 0.0; X : Scalar; begin for I in 1 .. N loop X := A + Scalar (I) * H; Sum := Sum + H * F (X); end loop; return Sum; end Right_Rectangular;
function Midpoint_Rectangular (A, B : Scalar; N : Positive) return Scalar is H : constant Scalar := (B - A) / Scalar (N); Sum : Scalar := 0.0; X : Scalar; begin for I in 1 .. N loop X := A + Scalar (I) * H - 0.5 * H; Sum := Sum + H * F (X); end loop; return Sum; end Midpoint_Rectangular;
function Trapezium (A, B : Scalar; N : Positive) return Scalar is H : constant Scalar := (B - A) / Scalar (N); Sum : Scalar := F(A) + F(B); X : Scalar := 1.0; begin while X <= Scalar (N) - 1.0 loop Sum := Sum + 2.0 * F (A + X * (B - A) / Scalar (N)); X := X + 1.0; end loop; return (B - A) / (2.0 * Scalar (N)) * Sum; end Trapezium;
function Simpsons (A, B : Scalar; N : Positive) return Scalar is H : constant Scalar := (B - A) / Scalar (N); Sum_1 : Scalar := 0.0; Sum_2 : Scalar := 0.0; begin for I in 0 .. N - 1 loop Sum_1 := Sum_1 + F (A + H * Scalar (I) + 0.5 * H); Sum_2 := Sum_2 + F (A + H * Scalar (I)); end loop; return H / 6.0 * (F (A) + F (B) + 4.0 * Sum_1 + 2.0 * Sum_2); end Simpsons;
end Integrate;</lang>
Test driver: <lang ada>with Ada.Text_IO, Ada.Integer_Text_IO; with Integrate;
procedure Numerical_Integration is
type Scalar is digits 18; package Scalar_Text_IO is new Ada.Text_IO.Float_IO (Scalar);
generic with function F (X : Scalar) return Scalar; Name : String; From, To : Scalar; Steps : Positive; procedure Test;
procedure Test is package Integrate_Scalar_F is new Integrate (Scalar, F); use Ada.Text_IO, Ada.Integer_Text_IO, Integrate_Scalar_F, Scalar_Text_IO; begin Put (Name & " integrated from "); Put (From); Put (" to "); Put (To); Put (" in "); Put (Steps); Put_Line (" steps:");
Put ("Rectangular (left): "); Put (Left_Rectangular (From, To, Steps)); New_Line;
Put ("Rectangular (right): "); Put (Right_Rectangular (From, To, Steps)); New_Line;
Put ("Rectangular (midpoint): "); Put (Midpoint_Rectangular (From, To, Steps)); New_Line;
Put ("Trapezium: "); Put (Trapezium (From, To, Steps)); New_Line;
Put ("Simpson's: "); Put (Simpsons (From, To, Steps)); New_Line;
New_Line; end Test;
begin
Ada.Integer_Text_IO.Default_Width := 0; Scalar_Text_IO.Default_Fore := 0; Scalar_Text_IO.Default_Exp := 0;
Cubed:
declare function F (X : Scalar) return Scalar is begin return X ** 3; end F; procedure Run is new Test (F => F, Name => "x^3", From => 0.0, To => 1.0, Steps => 100); begin Run; end Cubed;
One_Over_X:
declare function F (X : Scalar) return Scalar is begin return 1.0 / X; end F; procedure Run is new Test (F => F, Name => "1/x", From => 1.0, To => 100.0, Steps => 1_000); begin Run; end One_Over_X;
X:
declare function F (X : Scalar) return Scalar is begin return X; end F; procedure Run_1 is new Test (F => F, Name => "x", From => 0.0, To => 5_000.0, Steps => 5_000_000); procedure Run_2 is new Test (F => F, Name => "x", From => 0.0, To => 6_000.0, Steps => 6_000_000); begin Run_1; Run_2; end X;
end Numerical_Integration; </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 <stdio.h>
- include <stdlib.h>
- include <math.h>
double int_leftrect(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); return h*sum;
}
double int_rightrect(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+h); 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+h/2.0); return h*sum;
}
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 += 2.0*func(from + i * h); return h * sum / 2.0;
}
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;
double x; 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>
<lang c>/* test */ double f3(double x) {
return x;
}
double f3a(double x) {
return x*x/2.0;
}
double f2(double x) {
return 1.0/x;
}
double f2a(double x) {
return log(x);
}
double f1(double x) {
return x*x*x;
}
double f1a(double x) {
return x*x*x*x/4.0;
}
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[] = { f1, f2, f3, f3 }; rfunc If[] = { f1a, f2a, f3a, f3a }; double ivals[] = { 0.0, 1.0, 1.0, 100.0, 0.0, 5000.0, 0.0, 6000.0 }; double approx[] = { 100.0, 1000.0, 5000000.0, 6000000.0 }; for(j=0; j < (sizeof(rf) / sizeof(rfunc)); j++) { for(i=0; i < 5 ; i++) { ic = (*f[i])(ivals[2*j], ivals[2*j+1], approx[j], rf[j]); printf("%10s [ 0,1] num: %+lf, an: %lf\n", names[i], ic, INTG((*If[j]), ivals[2*j], ivals[2*j+1])); } printf("\n"); }
}</lang>
C#
<lang csharp>using System; using System.Collections.Generic; using System.Linq;
public class Interval {
public Interval(double leftEndpoint, double size) { LeftEndpoint = leftEndpoint; RightEndpoint = leftEndpoint + size; }
public double LeftEndpoint { get; set; }
public double RightEndpoint { get; set; }
public double Size { get { return RightEndpoint - LeftEndpoint; } }
public double Center { get { return (LeftEndpoint + RightEndpoint) / 2; } }
public IEnumerable<Interval> Subdivide(int subintervalCount) { double subintervalSize = Size / subintervalCount; return Enumerable.Range(0, subintervalCount).Select(index => new Interval(LeftEndpoint + index * subintervalSize, subintervalSize)); }
}
public class DefiniteIntegral {
public DefiniteIntegral(Func<double, double> integrand, Interval domain) { Integrand = integrand; Domain = domain; }
public Func<double, double> Integrand { get; set; }
public Interval Domain { get; set; }
public double SampleIntegrand(ApproximationMethod approximationMethod, Interval subdomain) { switch (approximationMethod) { case ApproximationMethod.RectangleLeft: return Integrand(subdomain.LeftEndpoint); case ApproximationMethod.RectangleMidpoint: return Integrand(subdomain.Center); case ApproximationMethod.RectangleRight: return Integrand(subdomain.RightEndpoint); case ApproximationMethod.Trapezium: return (Integrand(subdomain.LeftEndpoint) + Integrand(subdomain.RightEndpoint)) / 2; case ApproximationMethod.Simpson: return (Integrand(subdomain.LeftEndpoint) + 4 * Integrand(subdomain.Center) + Integrand(subdomain.RightEndpoint)) / 6; default: throw new NotImplementedException(); } }
public double Approximate(ApproximationMethod approximationMethod, int subdomainCount) { return Domain.Size * Domain.Subdivide(subdomainCount).Sum(subdomain => SampleIntegrand(approximationMethod, subdomain)) / subdomainCount; }
public enum ApproximationMethod { RectangleLeft, RectangleMidpoint, RectangleRight, Trapezium, Simpson }
}
public class Program {
private static void TestApproximationMethods(DefiniteIntegral integral, int subdomainCount) { foreach (DefiniteIntegral.ApproximationMethod approximationMethod in Enum.GetValues(typeof(DefiniteIntegral.ApproximationMethod))) { Console.WriteLine(integral.Approximate(approximationMethod, subdomainCount)); } }
public static void Main() { TestApproximationMethods(new DefiniteIntegral(x => x * x * x, new Interval(0, 1)), 10000); TestApproximationMethods(new DefiniteIntegral(x => 1 / x, new Interval(1, 99)), 1000); TestApproximationMethods(new DefiniteIntegral(x => x, new Interval(0, 5000)), 500000); TestApproximationMethods(new DefiniteIntegral(x => x, new Interval(0, 6000)), 6000000); }
}</lang> Output: <lang>0.2499500025 0.24999999875 0.2500500025 0.250000002499999 0.25 4.65499105751468 4.60476254867838 4.55698105751468 4.60598605751468 4.60517038495713 12499975 12500000 12500025 12500000 12500000 17999997 18000000 18000003 18000000 18000000</lang>
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
<lang d>import std.stdio, std.typecons, std.typetuple;
template integrate(alias method) {
double integrate(F, Float)(F f, Float a, Float b, int steps) { double s = 0.0; double h = (b - a) / steps; foreach (i; 0 .. steps) s += method(f, a + h*i, h); return h * s; }
}
double rectangularLeft(F, Float)(F f, Float x, Float h) {
return f(x);
}
double rectangularMiddle(F, Float)(F f, Float x, Float h) {
return f(x + h/2);
}
double rectangularRight(F, Float)(F f, Float x, Float h) {
return f(x + h);
}
double trapezium(F, Float)(F f, Float x, Float h) {
return (f(x) + f(x + h)) / 2;
}
double simpson(F, Float)(F f, Float x, Float h) {
return (f(x) + 4*f(x + h/2) + f(x + h)) / 6;
}
void main() {
auto args = [ tuple((double x){ return x ^^ 3; }, 0.0, 1.0, 10), tuple((double x){ return 1 / x; }, 1.0, 100.0, 1000), tuple((double x){ return x; }, 0.0, 5_000.0, 5_000_000), tuple((double x){ return x; }, 0.0, 6_000.0, 6_000_000)];
alias TypeTuple!(integrate!rectangularLeft, integrate!rectangularMiddle, integrate!rectangularRight, integrate!trapezium, integrate!simpson) ints;
alias TypeTuple!("rectangular left: ", "rectangular middle: ", "rectangular right: ", "trapezium: ", "simpson: ") names;
foreach (a; args) { foreach (i, n; names) writefln("%s %f", n, ints[i](a.tupleof)); writeln(); }
} </lang> Output:
rectangular left: 0.202500 rectangular middle: 0.248750 rectangular right: 0.302500 trapezium: 0.252500 simpson: 0.250000 rectangular left: 4.654991 rectangular middle: 4.604763 rectangular right: 4.556981 trapezium: 4.605986 simpson: 4.605170 rectangular left: 12499997.500000 rectangular middle: 12500000.000000 rectangular right: 12500002.500000 trapezium: 12500000.000000 simpson: 12500000.000000 rectangular left: 17999997.000000 rectangular middle: 18000000.000000 rectangular right: 18000003.000000 trapezium: 18000000.000000 simpson: 18000000.000000
A faster version
This version avoids function pointers and delegates (same output, 0.45 seconds run time): <lang d>import std.stdio, std.typecons, std.typetuple;
template integrate(alias method) {
template integrate(alias f) { double integrate(Float)(Float a, Float b, int steps) { Float s = 0.0; Float h = (b - a) / steps; foreach (i; 0 .. steps) s += method!(f, Float)(a + h*i, h); return h * s; } }
}
double rectangularLeft(alias f, Float)(Float x, Float h) {
return f(x);
}
double rectangularMiddle(alias f, Float)(Float x, Float h) {
return f(x + h/2);
}
double rectangularRight(alias f, Float)(Float x, Float h) {
return f(x + h);
}
double trapezium(alias f, Float)(Float x, Float h) {
return (f(x) + f(x + h)) / 2;
}
double simpson(alias f, Float)(Float x, Float h) {
return (f(x) + 4*f(x + h/2) + f(x + h)) / 6;
}
double f1(double x) { return x ^^ 3; } double f2(double x) { return 1 / x; } double f3(double x) { return x; } alias TypeTuple!(f1, f2, f3, f3) funcs;
void main() {
alias TypeTuple!("rectangular left: ", "rectangular middle: ", "rectangular right: ", "trapezium: ", "simpson: ") names;
alias TypeTuple!(integrate!rectangularLeft, integrate!rectangularMiddle, integrate!rectangularRight, integrate!trapezium, integrate!simpson) ints;
auto args = [tuple(0.0, 1.0, 10), tuple(1.0, 100.0, 1_000), tuple(0.0, 5_000.0, 5_000_000), tuple(0.0, 6_000.0, 6_000_000)];
foreach (i, f; funcs) { foreach (j, n; names) { alias ints[j] integ; writefln("%s %f", n, integ!f(args[i].tupleof)); } writeln(); }
}</lang>
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>
Go
<lang go>package main
import (
"fmt" "math" "sort"
)
// specification for an integration type spec struct {
lower, upper float64 // bounds for integration n int // number of parts exact float64 // expected answer fs string // mathematical description of function f func(float64) float64 // function to integrate
}
// test cases per task description var data = []spec{
spec{0, 1, 100, .25, "x^3", func(x float64) float64 { return x * x * x }}, spec{1, 100, 1000, float64(math.Log(100)), "1/x", func(x float64) float64 { return 1 / x }}, spec{0, 5000, 5e5, 12.5e6, "x", func(x float64) float64 { return x }}, spec{0, 6000, 6e6, 18e6, "x", func(x float64) float64 { return x }},
}
// object for associating a printable function name with an integration method type method struct {
name string integrate func(spec) float64
}
// integration methods implemented per task description var methods = []method{
method{"Rectangular (left) ", rectLeft}, method{"Rectangular (right) ", rectRight}, method{"Rectangular (midpoint)", rectMid}, method{"Trapezium ", trap}, method{"Simpson's ", simpson},
}
func rectLeft(t spec) float64 {
parts := make([]float64, t.n) r := t.upper - t.lower nf := float64(t.n) x0 := t.lower for i := range parts { x1 := t.lower + float64(i+1)*r/nf // x1-x0 better than r/nf. // (with r/nf, the represenation error accumulates) parts[i] = t.f(x0) * (x1 - x0) x0 = x1 } return sum(parts)
}
func rectRight(t spec) float64 {
parts := make([]float64, t.n) r := t.upper - t.lower nf := float64(t.n) x0 := t.lower for i := range parts { x1 := t.lower + float64(i+1)*r/nf parts[i] = t.f(x1) * (x1 - x0) x0 = x1 } return sum(parts)
}
func rectMid(t spec) float64 {
parts := make([]float64, t.n) r := t.upper - t.lower nf := float64(t.n) // there's a tiny gloss in the x1-x0 trick here. the correct way // would be to compute x's at division boundaries, but we don't need // those x's for anything else. (the function is evaluated on x's // at division midpoints rather than division boundaries.) so, we // reuse the midpoint x's, knowing that they will average out just // as well. we just need one extra point, so we use lower-.5. x0 := t.lower - .5*r/nf for i := range parts { x1 := t.lower + (float64(i)+.5)*r/nf parts[i] = t.f(x1) * (x1 - x0) x0 = x1 } return sum(parts)
}
func trap(t spec) float64 {
parts := make([]float64, t.n) r := t.upper - t.lower nf := float64(t.n) x0 := t.lower f0 := t.f(x0) for i := range parts { x1 := t.lower + float64(i+1)*r/nf f1 := t.f(x1) parts[i] = (f0 + f1) * .5 * (x1 - x0) x0, f0 = x1, f1 } return sum(parts)
}
func simpson(t spec) float64 {
parts := make([]float64, 2*t.n+1) r := t.upper - t.lower nf := float64(t.n) // similar to the rectangle midpoint logic explained above, // we play a little loose with the values used for dx and dx0. dx0 := r / nf parts[0] = t.f(t.lower) * dx0 parts[1] = t.f(t.lower+dx0*.5) * dx0 * 4 x0 := t.lower + dx0 for i := 1; i < t.n; i++ { x1 := t.lower + float64(i+1)*r/nf xmid := (x0 + x1) * .5 dx := x1 - x0 parts[2*i] = t.f(x0) * dx * 2 parts[2*i+1] = t.f(xmid) * dx * 4 x0 = x1 } parts[2*t.n] = t.f(t.upper) * dx0 return sum(parts) / 6
}
// sum a list of numbers avoiding loss of precision // (there might be better ways...) func sum(parts []float64) float64 {
sort.SortFloat64s(parts) ip := sort.SearchFloat64s(parts, 0) in := ip - 1 sum := parts[ip] ip++ for { if sum < 0 { if ip < len(parts) { // add the next biggest positive part sum += parts[ip] ip++ } else { // no more positive parts. // add remaining negatives and exit loop for ; in >= 0; in-- { sum += parts[in] } break } } else { if in >= 0 { // add the next biggest negative part sum += parts[in] in-- } else { // no more negative parts. // add remaining positives and exit loop for ; ip < len(parts); ip++ { sum += parts[ip] } break } } } return sum
}
func main() {
for _, t := range data { fmt.Println("Test case: f(x) =", t.fs) fmt.Println("Integration from", t.lower, "to", t.upper, "in", t.n, "parts") fmt.Printf("Exact result %.7e Error\n", t.exact) for _, m := range methods { a := m.integrate(t) e := a - t.exact if e < 0 { e = -e } fmt.Printf("%s %.7e %.7e\n", m.name, a, e) } fmt.Println("") }
}</lang> Output:
Test case: f(x) = x^3 Integration from 0 to 1 in 100 parts Exact result 2.5000000e-01 Error Rectangular (left) 2.4502500e-01 4.9750000e-03 Rectangular (right) 2.5502500e-01 5.0250000e-03 Rectangular (midpoint) 2.4998750e-01 1.2500000e-05 Trapezium 2.5002500e-01 2.5000000e-05 Simpson's 2.5000000e-01 5.5511151e-17 Test case: f(x) = 1/x Integration from 1 to 100 in 1000 parts Exact result 4.6051702e+00 Error Rectangular (left) 4.6549911e+00 4.9820872e-02 Rectangular (right) 4.5569811e+00 4.8189128e-02 Rectangular (midpoint) 4.6047625e+00 4.0763731e-04 Trapezium 4.6059861e+00 8.1587153e-04 Simpson's 4.6051704e+00 1.9896905e-07 Test case: f(x) = x Integration from 0 to 5000 in 500000 parts Exact result 1.2500000e+07 Error Rectangular (left) 1.2499975e+07 2.5000000e+01 Rectangular (right) 1.2500025e+07 2.5000000e+01 Rectangular (midpoint) 1.2500000e+07 1.6763806e-08 Trapezium 1.2500000e+07 3.7252903e-08 Simpson's 1.2500000e+07 2.9802322e-08 Test case: f(x) = x Integration from 0 to 6000 in 6000000 parts Exact result 1.8000000e+07 Error Rectangular (left) 1.7999997e+07 3.0000000e+00 Rectangular (right) 1.8000003e+07 3.0000000e+00 Rectangular (midpoint) 1.8000000e+07 1.4901161e-08 Trapezium 1.8000000e+07 1.8626451e-08 Simpson's 1.8000000e+07 7.4505806e-09
Groovy
Solution: <lang groovy>def assertBounds = { List bounds, int nRect ->
assert (bounds.size() == 2) && (bounds[0] instanceof Double) && (bounds[1] instanceof Double) && (nRect > 0)
}
def integral = { List bounds, int nRectangles, Closure f, List pointGuide, Closure integralCalculator->
double a = bounds[0], b = bounds[1], h = (b - a)/nRectangles def xPoints = pointGuide.collect { double it -> a + it*h } def fPoints = xPoints.collect { x -> f(x) } integralCalculator(h, fPoints)
}
def leftRectIntegral = { List bounds, int nRect, Closure f ->
assertBounds(bounds, nRect) integral(bounds, nRect, f, (0..<nRect)) { h, fPoints -> h*fPoints.sum() }
}
def rightRectIntegral = { List bounds, int nRect, Closure f ->
assertBounds(bounds, nRect) integral(bounds, nRect, f, (1..nRect)) { h, fPoints -> h*fPoints.sum() }
}
def midRectIntegral = { List bounds, int nRect, Closure f ->
assertBounds(bounds, nRect) integral(bounds, nRect, f, ((0.5d)..nRect)) { h, fPoints -> h*fPoints.sum() }
}
def trapezoidIntegral = { List bounds, int nRect, Closure f ->
assertBounds(bounds, nRect) integral(bounds, nRect, f, (0..nRect)) { h, fPoints -> def fLeft = fPoints[0..<nRect] def fRight = fPoints[1..nRect] h/2*(fLeft + fRight).sum() }
}
def simpsonsIntegral = { List bounds, int nSimpRect, Closure f ->
assertBounds(bounds, nSimpRect) integral(bounds, nSimpRect*2, f, (0..(nSimpRect*2))) { h, fPoints -> def fLeft = fPoints[(0..<nSimpRect*2).step(2)] def fMid = fPoints[(1..<nSimpRect*2).step(2)] def fRight = fPoints[(2..nSimpRect*2).step(2)] h/3*((fLeft + fRight).sum() + 4*(fMid.sum())) }
}</lang>
Test:
Each "nRect" (number of rectangles) value given below is the minimum value that meets the tolerance condition for the given circumstances (function-to-integrate, integral-type and integral-bounds). <lang groovy>double tolerance = 0.0001 // allowable "wrongness", ensures accuracy to 1 in 10,000
double sinIntegralCalculated = -(Math.cos(Math.PI) - Math.cos(0d)) assert (leftRectIntegral([0d, Math.PI], 129, Math.&sin) - sinIntegralCalculated).abs() < tolerance assert (rightRectIntegral([0d, Math.PI], 129, Math.&sin) - sinIntegralCalculated).abs() < tolerance assert (midRectIntegral([0d, Math.PI], 91, Math.&sin) - sinIntegralCalculated).abs() < tolerance assert (trapezoidIntegral([0d, Math.PI], 129, Math.&sin) - sinIntegralCalculated).abs() < tolerance assert (simpsonsIntegral([0d, Math.PI], 6, Math.&sin) - sinIntegralCalculated).abs() < tolerance
double cubeIntegralCalculated = 1d/4d *(10d**4 - 0d**4) assert ((leftRectIntegral([0d, 10d], 20000) { it**3 } - cubeIntegralCalculated)/cubeIntegralCalculated).abs() < tolerance assert ((rightRectIntegral([0d, 10d], 20001) { it**3 } - cubeIntegralCalculated)/cubeIntegralCalculated).abs() < tolerance assert ((midRectIntegral([0d, 10d], 71) { it**3 } - cubeIntegralCalculated)/cubeIntegralCalculated).abs() < tolerance assert ((trapezoidIntegral([0d, 10d], 101) { it**3 } - cubeIntegralCalculated)/cubeIntegralCalculated).abs() < tolerance // I can name that tune in one note! assert (simpsonsIntegral([0d, 10d], 1) { it**3 } == cubeIntegralCalculated) assert (simpsonsIntegral([0d, Math.PI], 1) { it**3 } == (1d/4d *(Math.PI**4 - 0d**4))) assert (simpsonsIntegral([-7.23d, Math.PI], 1) { it**3 } == (1d/4d *(Math.PI**4 - (-7.23d)**4)))
double quarticIntegralCalculated = 1d/5d *(10d**5 - 0d**5) assert ((leftRectIntegral([0d, 10d], 25000) { it**4 } - quarticIntegralCalculated)/quarticIntegralCalculated).abs() < tolerance assert ((rightRectIntegral([0d, 10d], 25001) { it**4 } - quarticIntegralCalculated)/quarticIntegralCalculated).abs() < tolerance assert ((midRectIntegral([0d, 10d], 92) { it**4 } - quarticIntegralCalculated)/quarticIntegralCalculated).abs() < tolerance assert ((trapezoidIntegral([0d, 10d], 130) { it**4 } - quarticIntegralCalculated)/quarticIntegralCalculated).abs() < tolerance assert ((simpsonsIntegral([0d, 10d], 5) { it**4 } - quarticIntegralCalculated)/quarticIntegralCalculated).abs() < tolerance
def cubicPoly = { it**3 + 2*it**2 + 7*it + 12d } def cubicPolyAntiDeriv = { 1/4*it**4 + 2/3*it**3 + 7/2*it**2 + 12*it } double cubicPolyIntegralCalculated = (cubicPolyAntiDeriv(10d) - cubicPolyAntiDeriv(0d)) assert ((leftRectIntegral([0d, 10d], 20000, cubicPoly) - cubicPolyIntegralCalculated)/cubicPolyIntegralCalculated).abs() < tolerance assert ((rightRectIntegral([0d, 10d], 20001, cubicPoly) - cubicPolyIntegralCalculated)/cubicPolyIntegralCalculated).abs() < tolerance assert ((midRectIntegral([0d, 10d], 71, cubicPoly) - cubicPolyIntegralCalculated)/cubicPolyIntegralCalculated).abs() < tolerance assert ((trapezoidIntegral([0d, 10d], 101, cubicPoly) - cubicPolyIntegralCalculated)/cubicPolyIntegralCalculated).abs() < tolerance // I can name that tune in one note! assert ((simpsonsIntegral([0d, 10d], 1, cubicPoly) - cubicPolyIntegralCalculated)/cubicPolyIntegralCalculated).abs() < tolerance**2.75 // 1 in 100 billion
double cpIntegralCalc0ToPI = (cubicPolyAntiDeriv(Math.PI) - cubicPolyAntiDeriv(0d)) assert ((simpsonsIntegral([0d, Math.PI], 1, cubicPoly) - cpIntegralCalc0ToPI)/ cpIntegralCalc0ToPI).abs() < tolerance**2.75 // 1 in 100 billion double cpIntegralCalcMinusEToPI = (cubicPolyAntiDeriv(Math.PI) - cubicPolyAntiDeriv(-Math.E)) assert ((simpsonsIntegral([-Math.E, Math.PI], 1, cubicPoly) - cpIntegralCalcMinusEToPI)/ cpIntegralCalcMinusEToPI).abs() < tolerance**2.5 // 1 in 10 billion</lang>
Requested Demonstrations: <lang groovy>println "f(x) = x**3, where x is [0,1], with 100 approximations. The exact result is 1/4, or 0.25." println ([" LeftRect": leftRectIntegral([0d, 1d], 100) { it**3 }]) println (["RightRect": rightRectIntegral([0d, 1d], 100) { it**3 }]) println ([" MidRect": midRectIntegral([0d, 1d], 100) { it**3 }]) println (["Trapezoid": trapezoidIntegral([0d, 1d], 100) { it**3 }]) println ([" Simpsons": simpsonsIntegral([0d, 1d], 100) { it**3 }]) println ()
println "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." println ([" LeftRect": leftRectIntegral([1d, 100d], 1000) { 1/it }]) println (["RightRect": rightRectIntegral([1d, 100d], 1000) { 1/it }]) println ([" MidRect": midRectIntegral([1d, 100d], 1000) { 1/it }]) println (["Trapezoid": trapezoidIntegral([1d, 100d], 1000) { 1/it }]) println ([" Simpsons": simpsonsIntegral([1d, 100d], 1000) { 1/it }]) println ()
println "f(x) = x, where x is [0,5000], with 5,000,000 approximations. The exact result is 12,500,000." println ([" LeftRect": leftRectIntegral([0d, 5000d], 5000000) { it }]) println (["RightRect": rightRectIntegral([0d, 5000d], 5000000) { it }]) println ([" MidRect": midRectIntegral([0d, 5000d], 5000000) { it }]) println (["Trapezoid": trapezoidIntegral([0d, 5000d], 5000000) { it }]) println ([" Simpsons": simpsonsIntegral([0d, 5000d], 5000000) { it }]) println ()
println "f(x) = x, where x is [0,6000], with 6,000,000 approximations. The exact result is 18,000,000." println ([" LeftRect": leftRectIntegral([0d, 6000d], 6000000) { it }]) println (["RightRect": rightRectIntegral([0d, 6000d], 6000000) { it }]) println ([" MidRect": midRectIntegral([0d, 6000d], 6000000) { it }]) println (["Trapezoid": trapezoidIntegral([0d, 6000d], 6000000) { it }]) println ([" Simpsons": simpsonsIntegral([0d, 6000d], 6000000) { it }]) println ()</lang>
Output:
f(x) = x**3, where x is [0,1], with 100 approximations. The exact result is 1/4, or 0.25. [ LeftRect:0.24502500000000002] [RightRect:0.255025] [ MidRect:0.24998750000000008] [Trapezoid:0.250025] [ Simpsons:0.25000000000000006] 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. [ LeftRect:4.65499105751468] [RightRect:4.55698105751468] [ MidRect:4.604762548678376] [Trapezoid:4.605986057514673] [ Simpsons:4.605170384957142] f(x) = x, where x is [0,5000], with 5,000,000 approximations. The exact result is 12,500,000. [ LeftRect:1.24999975E7] [RightRect:1.25000025E7] [ MidRect:1.25E7] [Trapezoid:1.25E7] [ Simpsons:1.25E7] f(x) = x, where x is [0,6000], with 6,000,000 approximations. The exact result is 18,000,000. [ LeftRect:1.7999997000000004E7] [RightRect:1.8000003000000004E7] [ MidRect:1.7999999999999993E7] [Trapezoid:1.7999999999999996E7] [ Simpsons:1.7999999999999993E7]
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
The whole program:
<lang haskell>approx f xs ws = sum [w * f x | (x,w) <- zip xs ws]
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]]
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
intLeftRect = integrateClosed 1 [1,0] intMidRect = integrateOpen 1 [1] intRightRect = integrateClosed 1 [0,1] intTrapezium = integrateClosed 2 [1,1] intSimpson = integrateClosed 3 [1,4,1]
uncurry4 f ~(a, b, c, d) = f a b c d
main = do
let m1 = "rectangular left: " let m2 = "rectangular middle: " let m3 = "rectangular right: " let m4 = "trapezium: " let m5 = "simpson: "
let arg1 = ((\x -> x ^ 3), 0, 1, 100) putStrLn $ m1 ++ (show $ uncurry4 intLeftRect arg1) putStrLn $ m2 ++ (show $ uncurry4 intMidRect arg1) putStrLn $ m3 ++ (show $ uncurry4 intRightRect arg1) putStrLn $ m4 ++ (show $ uncurry4 intTrapezium arg1) putStrLn $ m5 ++ (show $ uncurry4 intSimpson arg1) putStrLn ""
let arg2 = ((\x -> 1 / x), 1, 100, 1000) putStrLn $ m1 ++ (show $ uncurry4 intLeftRect arg2) putStrLn $ m2 ++ (show $ uncurry4 intMidRect arg2) putStrLn $ m3 ++ (show $ uncurry4 intRightRect arg2) putStrLn $ m4 ++ (show $ uncurry4 intTrapezium arg2) putStrLn $ m5 ++ (show $ uncurry4 intSimpson arg2) putStrLn ""
let arg3 = ((\x -> x), 0, 5000, 5000000) putStrLn $ m1 ++ (show $ uncurry4 intLeftRect arg3) putStrLn $ m2 ++ (show $ uncurry4 intMidRect arg3) putStrLn $ m3 ++ (show $ uncurry4 intRightRect arg3) putStrLn $ m4 ++ (show $ uncurry4 intTrapezium arg3) putStrLn $ m5 ++ (show $ uncurry4 intSimpson arg3) putStrLn ""
let arg4 = ((\x -> x), 0, 6000, 6000000) putStrLn $ m1 ++ (show $ uncurry4 intLeftRect arg4) putStrLn $ m2 ++ (show $ uncurry4 intMidRect arg4) putStrLn $ m3 ++ (show $ uncurry4 intRightRect arg4) putStrLn $ m4 ++ (show $ uncurry4 intTrapezium arg4) putStrLn $ m5 ++ (show $ uncurry4 intSimpson arg4)</lang>
Output:
rectangular left: 0.24502500000000005 rectangular middle: 0.24998750000000006 rectangular right: 0.25502500000000006 trapezium: 0.25002500000000005 simpson: 0.25000000000000006 rectangular left: 4.65499105751468 rectangular middle: 4.604762548678376 rectangular right: 4.55698105751468 trapezium: 4.605986057514681 simpson: 4.605170384957134 rectangular left: 1.24999975e7 rectangular middle: 1.25e7 rectangular right: 1.25000025e7 trapezium: 1.25e7 simpson: 1.2499999999999993e7 rectangular left: 1.7999997000000004e7 rectangular middle: 1.7999999999999993e7 rectangular right: 1.8000003000000004e7 trapezium: 1.8000000000000004e7 simpson: 1.7999999999999993e7
Runtime: about 7 seconds.
J
Solution:
<lang j>integrate=: adverb define
'a b steps'=. 3{.y,128 size=. (b - a)%steps size * +/ u |:2 ]\ 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
Required Examples
<lang j> Ir=: rectangle integrate
It=: trapezium integrate Is=: simpson integrate ^&3 Ir 0 1 100
0.249987
^&3 It 0 1 100
0.250025
^&3 Is 0 1 100
0.25
% Ir 1 100 1000
4.60476
% It 1 100 1000
4.60599
% Is 1 100 1000
4.60517
] Ir 0 5000 5e6
1.25e7
] It 0 5000 5e6
1.25e7
] Is 0 5000 5e6
1.25e7
] Ir 0 6000 6e6
1.8e7
] It 0 6000 6e6
1.8e7
] Is 0 6000 6e6
1.8e7</lang>
Older Examples
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>
Aside
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
<lang java5>class NumericalIntegration {
interface FPFunction { double eval(double n); } public static double rectangularLeft(double a, double b, int n, FPFunction f) { return rectangular(a, b, n, f, 0); } public static double rectangularMidpoint(double a, double b, int n, FPFunction f) { return rectangular(a, b, n, f, 1); } public static double rectangularRight(double a, double b, int n, FPFunction f) { return rectangular(a, b, n, f, 2); } public static double trapezium(double a, double b, int n, FPFunction f) { double range = checkParamsGetRange(a, b, n); double nFloat = (double)n; double sum = 0.0; for (int i = 1; i < n; i++) { double x = a + range * (double)i / nFloat; sum += f.eval(x); } sum += (f.eval(a) + f.eval(b)) / 2.0; return sum * range / nFloat; } public static double simpsons(double a, double b, int n, FPFunction f) { double range = checkParamsGetRange(a, b, n); double nFloat = (double)n; double sum1 = f.eval(a + range / (nFloat * 2.0)); double sum2 = 0.0; for (int i = 1; i < n; i++) { double x1 = a + range * ((double)i + 0.5) / nFloat; sum1 += f.eval(x1); double x2 = a + range * (double)i / nFloat; sum2 += f.eval(x2); } return (f.eval(a) + f.eval(b) + sum1 * 4.0 + sum2 * 2.0) * range / (nFloat * 6.0); } private static double rectangular(double a, double b, int n, FPFunction f, int mode) { double range = checkParamsGetRange(a, b, n); double modeOffset = (double)mode / 2.0; double nFloat = (double)n; double sum = 0.0; for (int i = 0; i < n; i++) { double x = a + range * ((double)i + modeOffset) / nFloat; sum += f.eval(x); } return sum * range / nFloat; } private static double checkParamsGetRange(double a, double b, int n) { if (n <= 0) throw new IllegalArgumentException("Invalid value of n"); double range = b - a; if (range <= 0) throw new IllegalArgumentException("Invalid range"); return range; } private static void testFunction(String fname, double a, double b, int n, FPFunction f) { System.out.println("Testing function \"" + fname + "\", a=" + a + ", b=" + b + ", n=" + n); System.out.println("rectangularLeft: " + rectangularLeft(a, b, n, f)); System.out.println("rectangularMidpoint: " + rectangularMidpoint(a, b, n, f)); System.out.println("rectangularRight: " + rectangularRight(a, b, n, f)); System.out.println("trapezium: " + trapezium(a, b, n, f)); System.out.println("simpsons: " + simpsons(a, b, n, f)); System.out.println(); return; } public static void main(String[] args) { testFunction("x^3", 0.0, 1.0, 100, new FPFunction() { public double eval(double n) { return n * n * n; } } ); testFunction("1/x", 1.0, 100.0, 1000, new FPFunction() { public double eval(double n) { return 1.0 / n; } } ); testFunction("x", 0.0, 5000.0, 5000000, new FPFunction() { public double eval(double n) { return n; } } ); testFunction("x", 0.0, 6000.0, 6000000, new FPFunction() { public double eval(double n) { return n; } } ); return; }
} </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>
Lua
<lang lua>function leftRect( f, a, b, n )
local h = (b - a) / n local x = a local sum = 0 for i = 1, 100 do sum = sum + a + f(x) x = x + h end
return sum * h
end
function rightRect( f, a, b, n )
local h = (b - a) / n local x = b local sum = 0 for i = 1, 100 do sum = sum + a + f(x) x = x - h end return sum * h
end
function midRect( f, a, b, n )
local h = (b - a) / n local x = a + h/2 local sum = 0 for i = 1, 100 do sum = sum + a + f(x) x = x + h end
return sum * h
end
function trapezium( f, a, b, n )
local h = (b - a) / n local x = a local sum = 0 for i = 1, 100 do sum = sum + f(x)*2 x = x + h end
return (b - a) * sum / (2 * n)
end
function simpson( f, a, b, n )
local h = (b - a) / n local sum1 = f(a + h/2) local sum2 = 0
for i = 1, n-1 do sum1 = sum1 + f(a + h * i + h/2) sum2 = sum2 + f(a + h * i) end
return (h/6) * (f(a) + f(b) + 4*sum1 + 2*sum2)
end
int_methods = { leftRect, rightRect, midRect, trapezium, simpson }
for i = 1, 5 do
print( int_methods[i]( function(x) return x^3 end, 0, 1, 100 ) ) print( int_methods[i]( function(x) return 1/x end, 1, 100, 1000 ) ) print( int_methods[i]( function(x) return x end, 0, 5000, 5000000 ) ) print( int_methods[i]( function(x) return x end, 0, 6000, 6000000 ) )
end</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 $h2 = $h/2; my $sum1 = f($a + $h2); my $sum2 = 0; for $a+$h, *+$h ... $b-$h { $sum1 += f($_ + $h2); $sum2 += f($_); } ($h / 6) * (f($a) + f($b) + 4*$sum1 + 2*$sum2);
}
sub tryem($f, $a, $b, $n, $exact) {
say "\n$f\n in [$a..$b] / $n"; eval "my &f = $f; say ' exact result: ', $exact; 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;"
}
tryem '{ $_ ** 3 }', 0, 1, 100, 0.25;
tryem '1 / *', 1, 100, 1000, log(100);
tryem '{$_}', 0, 5_000, 10_000, 12_500_000;
tryem '{$_}', 0, 6_000, 12_000, 18_000_000;</lang> (We run the last two tests with fewer iterations to avoid burning 60 hours of CPU, since rakudo hasn't been optimized yet.)
Output: <lang>{ $_ ** 3 }
in [0..1] / 100 exact result: 0.25 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 / *
in [1..100] / 1000 exact result: 4.60517018598809 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
{$_}
in [0..5000] / 10000 exact result: 12500000 rectangle method left: 12498750 rectangle method right: 12501250 rectangle method mid: 12500000
composite trapezoidal rule: 12500000
quadratic simpsons rule: 12500000
{$_}
in [0..6000] / 12000 exact result: 18000000 rectangle method left: 17998500 rectangle method right: 18001500 rectangle method mid: 18000000
composite trapezoidal rule: 18000000
quadratic simpsons rule: 18000000</lang>
Note that these integrations are done with rationals rather than floats, so should be fairly precise (though of course with so few iterations they are not terribly accurate (except when they are)). 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)))</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
- Tests
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 ((0., 5000., 5000000, identity),
(0., 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>
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 0.0 to 5000.0 (5000000 steps) = 12499997.5 identity integrated using mid_rect from 0.0 to 5000.0 (5000000 steps) = 12500000.0 identity integrated using right_rect from 0.0 to 5000.0 (5000000 steps) = 12500002.5 identity integrated using trapezium from 0.0 to 5000.0 (5000000 steps) = 12500000.0 identity integrated using simpson from 0.0 to 5000.0 (5000000 steps) = 12500000.0 identity integrated using left_rect from Fraction(0, 1) to Fraction(5000, 1) (5000000 steps and fractions) = 12499997.5 identity integrated using mid_rect from Fraction(0, 1) to Fraction(5000, 1) (5000000 steps and fractions) = 12500000.0 identity integrated using right_rect from Fraction(0, 1) to Fraction(5000, 1) (5000000 steps and fractions) = 12500002.5 identity integrated using trapezium from Fraction(0, 1) to Fraction(5000, 1) (5000000 steps and fractions) = 12500000.0 identity integrated using simpson from Fraction(0, 1) to Fraction(5000, 1) (5000000 steps and fractions) = 12500000.0 identity integrated using left_rect from 0.0 to 6000.0 (6000000 steps) = 17999997.000000004 identity integrated using mid_rect from 0.0 to 6000.0 (6000000 steps) = 17999999.999999993 identity integrated using right_rect from 0.0 to 6000.0 (6000000 steps) = 18000003.000000004 identity integrated using trapezium from 0.0 to 6000.0 (6000000 steps) = 17999999.999999993 identity integrated using simpson from 0.0 to 6000.0 (6000000 steps) = 17999999.999999993 identity integrated using left_rect from Fraction(0, 1) to Fraction(6000, 1) (6000000 steps and fractions) = 17999997.0 identity integrated using mid_rect from Fraction(0, 1) to Fraction(6000, 1) (6000000 steps and fractions) = 18000000.0 identity integrated using right_rect from Fraction(0, 1) to Fraction(6000, 1) (6000000 steps and fractions) = 18000003.0 identity integrated using trapezium from Fraction(0, 1) to Fraction(6000, 1) (6000000 steps and fractions) = 17999999.999999993 identity integrated using simpson from Fraction(0, 1) to Fraction(6000, 1) (6000000 steps and fractions) = 17999999.999999993
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>
R
<lang R>integrate_simpsons <- function(f,a,b,n){
h = (b-a)/n sum1 = f(a + h/2) sum2 = 0 for (i in 1:(n-1)) { sum1 = sum1 + f(a + h * i + h/2) sum2 = sum2 + f(a + h * i) } (h / 6) * (f(a) + f(b) + (4*sum1) + (2*sum2))
}
integrate_rect <- function(f,a,b,n,k=0) {
h = (b-a)/n x = 0 for (i in 1:n) { x = x + f(a + ((i-k)*h)) * h } x
}
integrate_trapez <- function(f,a,b,n) {
h = (b-a)/n x=0 for (k in 1:(n-1)) { x = x + f(a + k*h) } (h/2) * (f(a) + f(b)) + (h*x)
}
f1 <- (function(x) {x^3}) f2 <- (function(x) {1/x}) f3 <- (function(x) {x}) f4 <- (function(x) {x})
integrate_simpsons(f1,0,1,100) #0.25 integrate_simpsons(f2,1,100,1000) # 4.60517 integrate_simpsons(f3,0,5000,5000000) # 12500000 integrate_simpsons(f4,0,6000,6000000) # 1.8e+07
integrate_rect(f1,0,1,100,1) #TopLeft 0.245025 integrate_rect(f1,0,1,100,0.5) #Mid 0.2499875 integrate_rect(f1,0,1,100,0) #TopRight 0.255025
integrate_trapez(f1,0,1,100) # 0.250025</lang>
REXX
Note: there was virtually no difference between
numeric digits 9 (the default0
and
numeric digits 20.
<lang rexx>
/*REXX program numerically integrates using five different methods. */
numeric digits 20
do test=1 for 4
select when test==1 then do; L=0; H= 1; i= 100; end when test==2 then do; L=1; H= 100; i= 1000; end when test==3 then do; L=0; H=5000; i=5000000; end when test==4 then do; L=0; H=6000; i=5000000; end end
say center('test' test,79,'─') say ' left_rectangular('L","H','i")=" left_rectangular(L,H,i) say ' midpoint_rectangular('L","H','i")=" midpoint_rectangular(L,H,i) say ' right_rectangular('L","H','i")=" right_rectangular(L,H,i) say ' simpson('L","H','i")=" simpson(L,H,i) say ' trapezoid('L","H','i")=" trapezoid(L,H,i) end /*test*/
exit
/*─────────────────────────────────────LEFT_RECTANGULAR procedure───────*/
left_rectangular: procedure expose test; parse arg a,b,n
h=(b-a)/n
sum=0
do x=a by h for n sum=sum+f(x) end
return sum*h/1
/*─────────────────────────────────────MIDPOINT_RECTANGULAR procedure───*/
midpoint_rectangular: procedure expose test; parse arg a,b,n
h=(b-a)/n
sum=0
do x=a+h/2 by h for n sum=sum+f(x) end
return sum*h/1
/*─────────────────────────────────────RIGHT_RECTANGULAR procedure──────*/
right_rectangular: procedure expose test; parse arg a,b,n
h=(b-a)/n
sum=0
do x=a+h by h for n sum=sum+f(x) end
return sum*h/1
/*─────────────────────────────────────SIMPSON procedure────────────────*/
simpson: procedure expose test; parse arg a,b,n
h=(b-a)/n
sum1=f(a+h/2)
sum2=0
do x=1 for n-1 sum1=sum1+f(a+h*x+h*.5) sum2=sum2+f(a+x*h) end
return h*(f(a)+f(b)+4*sum1+2*sum2)/6
/*─────────────────────────────────────TRAPEZOID procedure──────────────*/
trapezoid: procedure expose test; parse arg a,b,n
h=(b-a)/n
sum=0
do x=a to b by h sum=sum+h*(f(x)+f(x+h))*.5 end
return sum
/*─────────────────────────────────────F function (depends on TEST)─────*/
f: procedure expose test; parse arg z
select when test==1 then return z**3 when test==2 then return 1/z otherwise return z end
</lang> Output:
────────────────────────────────────test 1───────────────────────────────────── left_rectangular(0,1,100)= 0.245025 midpoint_rectangular(0,1,100)= 0.2499875 right_rectangular(0,1,100)= 0.255025 simpson(0,1,100)= 0.25 trapezoid(0,1,100)= 0.260176505 ────────────────────────────────────test 2───────────────────────────────────── left_rectangular(1,100,1000)= 4.6549910575146761473 midpoint_rectangular(1,100,1000)= 4.604762548678375185 right_rectangular(1,100,1000)= 4.5569810575146761472 simpson(1,100,1000)= 4.6051703849571421725 trapezoid(1,100,1000)= 4.6069755679493458225 ────────────────────────────────────test 3───────────────────────────────────── left_rectangular(0,5000,5000000)= 12499997.5 midpoint_rectangular(0,5000,5000000)= 12500000 right_rectangular(0,5000,5000000)= 12500002.5 simpson(0,5000,5000000)= 12500000 trapezoid(0,5000,5000000)= 12500005.0000005 ────────────────────────────────────test 4───────────────────────────────────── left_rectangular(0,6000,5000000)= 17999996.4 midpoint_rectangular(0,6000,5000000)= 18000000 right_rectangular(0,6000,5000000)= 18000003.6 simpson(0,6000,5000000)= 18000000 trapezoid(0,6000,5000000)= 18000007.20000072
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
- ALGOL 68
- AutoHotkey
- BASIC
- BASIC examples needing attention
- Examples needing attention
- C
- C sharp
- C++
- Common Lisp
- D
- E
- Forth
- Fortran
- Go
- Groovy
- Haskell
- J
- Java
- Logo
- Lua
- MATLAB
- Sample inputs
- OCaml
- Pascal
- Perl 6
- PL/I
- PicoLisp
- PureBasic
- Python
- R
- R examples needing attention
- REXX
- Ruby
- Scheme
- Standard ML
- Tcl
- TI-89 BASIC examples needing attention
- Ursala
- M4/Omit
- Arithmetic
- Mathematics