Numerical integration: Difference between revisions
(C++ version) |
m (→{{header|ALGOL 68}}: Should probably drop the n in favour of an error requirement, also from left and right rect methods) |
||
Line 114: | Line 114: | ||
end Integrate; |
end Integrate; |
||
=={{header|ALGOL 68}}== |
|||
MODE F = PROC(LONG REAL)LONG REAL; |
|||
############### |
|||
## left rect ## |
|||
############### |
|||
PROC left rect = (F f, LONG REAL a, b, INT n) LONG REAL: |
|||
BEGIN |
|||
LONG REAL h= (b - a) / n; |
|||
LONG REAL sum:= 0; |
|||
LONG REAL x:= a; |
|||
WHILE x <= b - h DO |
|||
sum := sum + (h * f(x)); |
|||
x +:= h |
|||
OD; |
|||
sum |
|||
END # left rect #; |
|||
################# |
|||
## right rect ## |
|||
################# |
|||
PROC right rect = (F f, LONG REAL a, b, INT n) LONG REAL: |
|||
BEGIN |
|||
LONG REAL h= (b - a) / n; |
|||
LONG REAL sum:= 0; |
|||
LONG REAL x:= a + h; |
|||
WHILE x <= b - h DO |
|||
sum := sum + (h * f(x)); |
|||
x +:= h |
|||
OD; |
|||
sum |
|||
END # right rect #; |
|||
############### |
|||
## mid rect ## |
|||
############### |
|||
PROC mid rect = (F f, LONG REAL a, b, INT n) LONG REAL: |
|||
BEGIN |
|||
LONG REAL h= (b - a) / n; |
|||
LONG REAL sum:= 0; |
|||
LONG REAL x:= a; |
|||
WHILE x <= b - h DO |
|||
sum := sum + (h / 2) * (f(x) + f(x + h)); |
|||
x +:= h |
|||
OD; |
|||
sum |
|||
END # mid rect #; |
|||
############### |
|||
## trapezium ## |
|||
############### |
|||
PROC trapezium = (F f, LONG REAL a, b, INT n) LONG REAL: |
|||
BEGIN |
|||
LONG REAL h= (b - a) / n; |
|||
LONG REAL sum:= f(a) + f(b); |
|||
LONG REAL x:= 1; |
|||
WHILE x <= n - 1 DO |
|||
sum := sum + 2 * f(a + x * h ); |
|||
x +:= 1 |
|||
OD; |
|||
(b - a) / (2 * n) * sum |
|||
END # trapezium #; |
|||
############# |
|||
## simpson ## |
|||
############# |
|||
PROC simpson = (F f, LONG REAL a, b, INT n) LONG REAL: |
|||
BEGIN |
|||
LONG REAL h= (b - a) / n; |
|||
LONG REAL sum1:= 0; |
|||
LONG REAL sum2:= 0; |
|||
INT limit:= n - 1; |
|||
FOR i FROM 0 TO limit DO |
|||
sum1 := sum1 + f(a + h * LONG REAL(i) + h / 2) |
|||
OD; |
|||
FOR i FROM 1 TO limit DO |
|||
sum2 +:= f(a + h * LONG REAL(i)) |
|||
OD; |
|||
h / 6 * (f(a) + f(b) + 4 * sum1 + 2 * sum2) |
|||
END # simpson #; |
|||
SKIP |
|||
=={{header|BASIC}}== |
=={{header|BASIC}}== |
Revision as of 22:02, 14 February 2008
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 the wiki) 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.
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;
ALGOL 68
MODE F = PROC(LONG REAL)LONG REAL; ############### ## left rect ## ############### PROC left rect = (F f, LONG REAL a, b, INT n) LONG REAL: BEGIN LONG REAL h= (b - a) / n; LONG REAL sum:= 0; LONG REAL x:= a; WHILE x <= b - h DO sum := sum + (h * f(x)); x +:= h OD; sum END # left rect #; ################# ## right rect ## ################# PROC right rect = (F f, LONG REAL a, b, INT n) LONG REAL: BEGIN LONG REAL h= (b - a) / n; LONG REAL sum:= 0; LONG REAL x:= a + h; WHILE x <= b - h DO sum := sum + (h * f(x)); x +:= h OD; sum END # right rect #; ############### ## mid rect ## ############### PROC mid rect = (F f, LONG REAL a, b, INT n) LONG REAL: BEGIN LONG REAL h= (b - a) / n; LONG REAL sum:= 0; LONG REAL x:= a; WHILE x <= b - h DO sum := sum + (h / 2) * (f(x) + f(x + h)); x +:= h OD; sum END # mid rect #; ############### ## trapezium ## ############### PROC trapezium = (F f, LONG REAL a, b, INT n) LONG REAL: BEGIN LONG REAL h= (b - a) / n; LONG REAL sum:= f(a) + f(b); LONG REAL x:= 1; WHILE x <= n - 1 DO sum := sum + 2 * f(a + x * h ); x +:= 1 OD; (b - a) / (2 * n) * sum END # trapezium #; ############# ## simpson ## ############# PROC simpson = (F f, LONG REAL a, b, INT n) LONG REAL: BEGIN LONG REAL h= (b - a) / n; LONG REAL sum1:= 0; LONG REAL sum2:= 0; INT limit:= n - 1; FOR i FROM 0 TO limit DO sum1 := sum1 + f(a + h * LONG REAL(i) + h / 2) OD; FOR i FROM 1 TO limit DO sum2 +:= f(a + h * LONG REAL(i)) OD; h / 6 * (f(a) + f(b) + 4 * sum1 + 2 * sum2) END # simpson #; SKIP
BASIC
Compiler: QuickBasic 4.5
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 * (b - a) / n)) NEXT i trap = (b - a) / (2 * n) * 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
C++
Due to their similarity, it makes sense to make the integration method a policy.
// the integration routine template<typename Method, typename F, typename Float) double integrate(F f, Float a, Float b, int steps, Method m) { double s = 0; double h = (b-a)/steps; for (int i = 0; i < steps; ++i) s += m(f, a + h*steps, h); return h*s; } // methods class rectangular { public: enum position_type { left, middle, right }; rect(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; ) 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());
Java
The function in this example is assumed to be f(double x).
public class Integrate{ public static double leftRect(double a, double b, double n){ double h = (b - a) / n; double sum = 0; for(double x = a;x <= b - h;x += h) sum += h * (f(x)); return 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 += h * (f(x)); return 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 += (h / 2) * (f(x) + f(x + h)); return 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-1;i++) sum += 2 * f((a + i * (b - a) / n)); return (b - a) / (2 * n) * 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 - 1;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 }