47 lines
4.1 KiB
Rexx
47 lines
4.1 KiB
Rexx
/*REXX pgm performs numerical integration using 5 different algorithms and show results.*/
|
|
numeric digits 20 /*use twenty decimal digits precision. */
|
|
|
|
do test=1 for 4 /*perform the 4 different test suites. */
|
|
if test==1 then do; L= 0; H= 1; i= 100; end
|
|
if test==2 then do; L= 1; H= 100; i= 1000; end
|
|
if test==3 then do; L= 0; H= 5000; i= 5000000; end
|
|
if test==4 then do; L= 0; H= 6000; i= 5000000; end
|
|
say
|
|
say center('test' test, 65, "─") /*display a header for the test suite. */
|
|
say ' left rectangular('L", "H', 'i") ──► " left_rect(L, H, i)
|
|
say ' midpoint rectangular('L", "H', 'i") ──► " midpoint_rect(L, H, i)
|
|
say ' right rectangular('L", "H', 'i") ──► " right_rect(L, H, i)
|
|
say ' Simpson('L", "H', 'i") ──► " Simpson(L, H, i)
|
|
say ' trapezium('L", "H', 'i") ──► " trapezium(L, H, i)
|
|
end /*test*/
|
|
exit /*stick a fork in it, we're all done. */
|
|
/*──────────────────────────────────────────────────────────────────────────────────────*/
|
|
f: if test==1 then return arg(1) **3 /*choose the cube function. */
|
|
if test==2 then return 1 / arg(1) /* " " reciprocal " */
|
|
return arg(1) /* " " "as-is" " */
|
|
/*──────────────────────────────────────────────────────────────────────────────────────*/
|
|
left_rect: procedure expose test; parse arg a,b,n; h= (b-a) / n
|
|
$= 0
|
|
do x=a by h for n; $=$+f(x); end /*x*/
|
|
return $*h/1
|
|
/*──────────────────────────────────────────────────────────────────────────────────────*/
|
|
midpoint_rect: procedure expose test; parse arg a,b,n; h= (b-a) / n
|
|
$= 0
|
|
do x=a+h/2 by h for n; $=$+f(x); end /*x*/
|
|
return $*h/1
|
|
/*──────────────────────────────────────────────────────────────────────────────────────*/
|
|
right_rect: procedure expose test; parse arg a,b,n; h= (b-a) / n
|
|
$= 0
|
|
do x=a+h by h for n; $=$+f(x); end /*x*/
|
|
return $*h/1
|
|
/*──────────────────────────────────────────────────────────────────────────────────────*/
|
|
Simpson: procedure expose test; parse arg a,b,n; h= (b-a) / n
|
|
$= f(a + h/2)
|
|
@= 0; do x=1 for n-1; $=$+f(a+h*x+h*.5); @=@+f(a+x*h); end /*x*/
|
|
|
|
return h*(f(a) + f(b) + 4*$ + 2*@) / 6
|
|
/*──────────────────────────────────────────────────────────────────────────────────────*/
|
|
trapezium: procedure expose test; parse arg a,b,n; h=(b-a)/n
|
|
$= 0
|
|
do x=a by h for n; $=$+(f(x)+f(x+h)); end /*x*/
|
|
return $*h/2
|