Numerical integration: Difference between revisions
Content added Content deleted
m (→{{header|REXX}}: optimized two subroutines.) |
(→{{header|REXX}}: optimized the ƒ function (it's about 4% faster).) |
||
Line 4,286: | Line 4,286: | ||
=={{header|REXX}}== |
=={{header|REXX}}== |
||
Note: there was virtually no difference in accuracy between '''numeric digits 9''' (the default) and '''numeric digits 20'''. |
Note: there was virtually no difference in accuracy between '''numeric digits 9''' (the default) and '''numeric digits 20'''. |
||
<lang rexx>/*REXX pgm performs numerical integration using |
<lang rexx>/*REXX pgm performs numerical integration using 5 different algorithms and show results.*/ |
||
numeric digits 20 /*use twenty decimal digits precision. */ |
numeric digits 20 /*use twenty decimal digits precision. */ |
||
Line 4,295: | Line 4,295: | ||
if test==4 then do; L= 0; H= 6000; i= 6000000; end |
if test==4 then do; L= 0; H= 6000; i= 6000000; end |
||
say center('test' test, 79, "═") /*display a header for the test suite. */ |
say center('test' test, 79, "═") /*display a header for the test suite. */ |
||
say ' left rectangular('L", "H', 'i") ──► " |
say ' left rectangular('L", "H', 'i") ──► " left_rect(L, H, i) |
||
say ' midpoint rectangular('L", "H', 'i") ──► " |
say ' midpoint rectangular('L", "H', 'i") ──► " midpoint_rect(L, H, i) |
||
say ' right rectangular('L", "H', 'i") ──► " |
say ' right rectangular('L", "H', 'i") ──► " right_rect(L, H, i) |
||
say ' Simpson('L", "H', 'i") ──► " |
say ' Simpson('L", "H', 'i") ──► " Simpson(L, H, i) |
||
say ' trapezium('L", "H', 'i") ──► " |
say ' trapezium('L", "H', 'i") ──► " trapezium(L, H, i) |
||
end /*test*/ |
end /*test*/ |
||
exit /*stick a fork in it, we're all done. */ |
exit /*stick a fork in it, we're all done. */ |
||
/*──────────────────────────────────────────────────────────────────────────────────────*/ |
/*──────────────────────────────────────────────────────────────────────────────────────*/ |
||
f: |
f: parse arg y; if test>2 then return y /*choose the "as─is" function. */ |
||
if test==1 then return y**3 /* " " cube function. */ |
|||
return 1/y /* " " reciprocal " */ |
|||
/*──────────────────────────────────────────────────────────────────────────────────────*/ |
/*──────────────────────────────────────────────────────────────────────────────────────*/ |
||
left_rect: procedure expose test; parse arg a,b,n; $= 0; h= (b-a)/n |
left_rect: procedure expose test; parse arg a,b,n; $= 0; h= (b-a)/n |
||
do x=a by h for n; |
do x=a by h for n; $= $ + f(x); end /*x*/ |
||
return $*h/1 |
return $*h/1 |
||
/*──────────────────────────────────────────────────────────────────────────────────────*/ |
/*──────────────────────────────────────────────────────────────────────────────────────*/ |
||
midpoint_rect: procedure expose test; parse arg a,b,n; $= 0; h= (b-a)/n |
|||
do x=a+h/2 by h for n; |
do x=a+h/2 by h for n; $= $ + f(x); end /*x*/ |
||
return $*h/1 |
return $*h/1 |
||
/*──────────────────────────────────────────────────────────────────────────────────────*/ |
/*──────────────────────────────────────────────────────────────────────────────────────*/ |
||
right_rect: procedure expose test; parse arg a,b,n; $= 0; h= (b-a)/n |
right_rect: procedure expose test; parse arg a,b,n; $= 0; h= (b-a)/n |
||
do x=a+h by h for n; |
do x=a+h by h for n; $= $ + f(x); end /*x*/ |
||
return $*h/1 |
return $*h/1 |
||
/*──────────────────────────────────────────────────────────────────────────────────────*/ |
/*──────────────────────────────────────────────────────────────────────────────────────*/ |
||
Simpson: procedure expose test; parse arg a,b,n; h= (b-a)/n |
Simpson: procedure expose test; parse arg a,b,n; h= (b-a)/n; hh= h/2 |
||
$= f(a + h/2) |
|||
@= 0; |
@= 0; do x=1 for n-1; hx=h*x; $=$+f(a+hx+hh); @=@+f(a+hx); end /*x*/ |
||
return h |
return h * (f(a) + f(b) + 4*$ + 2*@) / 6 |
||
/*──────────────────────────────────────────────────────────────────────────────────────*/ |
|||
/*─────────--───────────────────────────────────────────────────────────────────────────*/ |
|||
trapezium: procedure expose test; parse arg a,b,n; $= 0; h= (b-a)/n |
trapezium: procedure expose test; parse arg a,b,n; $= 0; h= (b-a) / n |
||
do x=a by h for n; |
do x=a by h for n; $= $ + (f(x) + f(x+h)); end /*x*/ |
||
return $*h/2</lang> |
return $*h/2</lang> |
||
{{out|output|text= when using the default inputs:}} |
{{out|output|text= when using the default inputs:}} |
||
<pre> |
<pre> |