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 five different algorithms & show results*/
<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") ──► " left_rect(L, H, i)
say ' left rectangular('L", "H', 'i") ──► " left_rect(L, H, i)
say ' midpoint rectangular('L", "H', 'i") ──► " midpt_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 ' right rectangular('L", "H', 'i") ──► " right_rect(L, H, i)
say ' Simpson('L", "H', 'i") ──► " Simpson(L, H, i)
say ' Simpson('L", "H', 'i") ──► " Simpson(L, H, i)
say ' trapezium('L", "H', 'i") ──► " 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: if test==1 then return arg(1) **3 /*choose the cube function. */
f: parse arg y; if test>2 then return y /*choose the "as─is" function. */
if test==2 then return 1 / arg(1) /* " " reciprocal " */
if test==1 then return y**3 /* " " cube function. */
return arg(1) /* " " "as─is" " */
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; $= $ + f(x); end /*x*/
do x=a by h for n; $= $ + f(x); end /*x*/
return $*h/1
return $*h/1
/*──────────────────────────────────────────────────────────────────────────────────────*/
/*──────────────────────────────────────────────────────────────────────────────────────*/
midpt_rect: procedure expose test; parse arg a,b,n; $= 0; h= (b-a)/n
midpoint_rect: procedure expose test; parse arg a,b,n; $= 0; h= (b-a)/n
do x=a+h/2 by h for n; $= $ + f(x); end /*x*/
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; $= $ + f(x); end /*x*/
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
hh= h/2; $= f(a + hh)
$= f(a + h/2)
@= 0; do x=1 for n-1; _= h*x+a; $= $+f(_+hh); @= @+f(_); end /*x*/
@= 0; do x=1 for n-1; hx=h*x; $=$+f(a+hx+hh); @=@+f(a+hx); end /*x*/


return h * (f(a) + f(b) + 4*$ + 2*@) / 6
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; $= $ + (f(x) + f(x+h)); end /*x*/
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=&nbsp; when using the default inputs:}}
{{out|output|text=&nbsp; when using the default inputs:}}
<pre>
<pre>
Retrieved from "https://rosettacode.org/wiki/Numerical_integration"