Numerical integration: Difference between revisions

← Older editNewer edit →
Content deleted Content added
→‎{{header|Ada}}: Fixed incorrect integration. Split up package specification and body. Added a test driver.
Line 90: Line 90:
</lang>
</lang>
=={{header|Ada}}==
=={{header|Ada}}==
Specification of a generic package implementing the five specified kinds of numerical integration:
{{incorrect|Ada|Right_Rect stops one H too soon; Mid_Rect doesn't sample at midpoints but merely recreates Trapezium differently.}}
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
<lang ada>generic
type Scalar is digits <>;
with function F(X : Long_Float) return Long_Float;
with function F (X : Scalar) return Scalar;
package Integrate is
package Integrate is
function Left_Rect(A, B, N : Long_Float) return Long_Float;
function Left_Rectangular (A, B : Scalar; N : Positive) return Scalar;
function Right_Rect(A, B, N : Long_Float) return Long_Float;
function Right_Rectangular (A, B : Scalar; N : Positive) return Scalar;
function Mid_Rect(A, B, N : Long_Float) return Long_Float;
function Midpoint_Rectangular (A, B : Scalar; N : Positive) return Scalar;
function Trapezium(A, B, N : Long_Float) return Long_Float;
function Trapezium (A, B : Scalar; N : Positive) return Scalar;
function Simpson(A, B, N : Long_Float) return Long_Float;
function Simpsons (A, B : Scalar; N : Positive) return Scalar;
end Integrate;
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.
package body Integrate is


Body of the package implementing numerical integration:
---------------
<lang ada>package body Integrate is
-- Left_Rect --
function Left_Rectangular (A, B : Scalar; N : Positive) return Scalar is
---------------
H : constant Scalar := (B - A) / Scalar (N);

Sum : Scalar := 0.0;
function Left_Rect (A, B, N : Long_Float) return Long_Float is
H : constant Long_Float := (B - A) / N;
X : Scalar;
Sum : Long_Float := 0.0;
X : Long_Float := A;
begin
begin
while X <= B - H loop
for I in 0 .. N - 1 loop
Sum := Sum + (H * F(X));
X := A + Scalar (I) * H;
X := X + H;
Sum := Sum + H * F (X);
end loop;
end loop;
return Sum;
return Sum;
end Left_Rect;
end Left_Rectangular;


function Right_Rectangular (A, B : Scalar; N : Positive) return Scalar is
----------------
H : constant Scalar := (B - A) / Scalar (N);
-- Right_Rect --
Sum : Scalar := 0.0;
----------------
X : Scalar;

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
begin
while X <= B - H loop
for I in 1 .. N loop
Sum := Sum + (H * F(X));
X := A + Scalar (I) * H;
X := X + H;
Sum := Sum + H * F (X);
end loop;
end loop;
return Sum;
return Sum;
end Right_Rect;
end Right_Rectangular;


function Midpoint_Rectangular (A, B : Scalar; N : Positive) return Scalar is
--------------
H : constant Scalar := (B - A) / Scalar (N);
-- Mid_Rect --
Sum : Scalar := 0.0;
--------------
X : Scalar;

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
begin
while X <= B - H loop
for I in 1 .. N loop
Sum := Sum + (H / 2.0) * (F(X) + F(X + H));
X := A + Scalar (I) * H - 0.5 * H;
X := X + H;
Sum := Sum + H * F (X);
end loop;
end loop;
return Sum;
return Sum;
end Mid_Rect;
end Midpoint_Rectangular;


function Trapezium (A, B : Scalar; N : Positive) return Scalar is
---------------
H : constant Scalar := (B - A) / Scalar (N);
-- Trapezium --
Sum : Scalar := F(A) + F(B);
---------------
X : Scalar := 1.0;

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
begin
while X <= N - 1.0 loop
while X <= Scalar (N) - 1.0 loop
Sum := Sum + 2.0 * F(A + X * (B - A) / N);
Sum := Sum + 2.0 * F (A + X * (B - A) / Scalar (N));
X := X + 1.0;
X := X + 1.0;
end loop;
end loop;
return (B - A) / (2.0 * N) * Sum;
return (B - A) / (2.0 * Scalar (N)) * Sum;
end Trapezium;
end Trapezium;


function Simpsons (A, B : Scalar; N : Positive) return Scalar is
-------------
H : constant Scalar := (B - A) / Scalar (N);
-- Simpson --
Sum_1 : Scalar := 0.0;
-------------
Sum_2 : Scalar := 0.0;

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
begin
for I in 0..Limit loop
for I in 0 .. N - 1 loop
Sum1 := Sum1 + F(A + H * Long_Float(I) + H / 2.0);
Sum_1 := Sum_1 + F (A + H * Scalar (I) + 0.5 * H);
Sum_2 := Sum_2 + F (A + H * Scalar (I));
end loop;
end loop;
return H / 6.0 * (F (A) + F (B) + 4.0 * Sum_1 + 2.0 * Sum_2);
for I in 1..Limit loop
end Simpsons;
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>
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>


=={{header|ALGOL 68}}==
=={{header|ALGOL 68}}==
Retrieved from "https://rosettacode.org/wiki/Numerical_integration"