Numerical integration: Difference between revisions

<lang PL/I>integrals: procedure options (main); /* 1 September 2019 */

f: procedure (x, function) returns (float(18));

declare x float(18), function fixed binary;

select (function);

when (1) return (x**3);

when (2) return (1/x);

when (3) return (x);

when (4) return (x);

end;

declare (a, b) fixed decimal (10);

declare (rect_area, trap_area, Simpson) float(18);

declare (d, dx) float(18);

declare (S1, S2) float(18);

declare N fixed decimal (15), function fixed binary;

declare k fixed decimal (7,2);

put (' Rectangle-left Rectangle-mid Rectangle-right' ||

' Trapezoid Simpson');

do function = 1 to 4;

select(function);

when (1) do; N = 100; a = 0; b = 1; end;

when (2) do; N = 1000; a = 1; b = 100; end;

when (3) do; N = 5000000; a = 0; b = 5000; end;

when (4) do; N = 6000000; a = 0; b = 6000; end;

end;

dx = (b-a)/float(N);

/* Rectangle method, left-side */

rect_area = 0;

do d = 0 to N-1;

rect_area = rect_area + dx*f(a + d*dx, function);

end;

put skip edit (rect_area) (E(25, 15));

/* Rectangle method, mid-point */

rect_area = 0;

do d = 0 to N-1;

rect_area = rect_area + dx*f(a + d*dx + dx/2, function);

end;

put edit (rect_area) (E(25, 15));

/* Rectangle method, right-side */

rect_area = 0;

do d = 1 to N;

rect_area = rect_area + dx*f(a + d*dx, function);

end;

put edit (rect_area) (E(25, 15));

/* Trapezoid method */

trap_area = 0;

do d = 0 to N-1;

trap_area = trap_area + dx*(f(a+d*dx, function) + f(a+(d+1)*dx, function))/2;

end;

put edit (trap_area) (X(1), E(25, 15));

/* Simpson's Rule */

S1 = f(a+dx/2, function);

S2 = 0;

do d = 1 to N-1;

S1 = S1 + f(a+d*dx+dx/2, function);

S2 = S2 + f(a+d*dx, function);

end;

Simpson = dx * (f(a, function) + f(b, function) + 4*S1 + 2*S2) / 6;

put edit (Simpson) (X(1), E(25, 15));

<pre>

Rectangle-left Rectangle-mid Rectangle-right Trapezoid Simpson

2.450250000000000E-0001 2.499875000000000E-0001 2.550250000000000E-0001 2.500250000000000E-0001 2.500000000000000E-0001

4.654991057514676E+0000 4.604762548678375E+0000 4.556981057514676E+0000 4.605986057514676E+0000 4.605170384957142E+0000

1.249999750000000E+0007 1.250000000000000E+0007 1.250000250000000E+0007 1.250000000000000E+0007 1.250000000000000E+0007

1.799999700000000E+0007 1.800000000000000E+0007 1.800000300000000E+0007 1.800000000000000E+0007 1.800000000000000E+0007

</pre>