Numerical integration/Adaptive Simpson's method: Difference between revisions
Numerical integration/Adaptive Simpson's method (view source)
Revision as of 21:33, 19 May 2023
, 1 year agoAdded PostScript in front of →{{header|Prolog}}
(→{{header|COBOL}}: Bug fix.) |
(Added PostScript in front of →{{header|Prolog}}) |
||
Line 1,661:
Simpson's integration of sine from 0 to 1 = 0.459698
</pre>
=={{header|PostScript}}==
{{trans|COBOL}}
The following program makes little effort to be efficient. It iteratively computes the intervals, rather than do recursions.
<syntaxhighlight lang="postscript">
%!PS-Adobe-3.0 EPSF-3.0
%%BoundingBox: 0 0 280 30
/origstate save def
% We have to convert radians to degrees.
/f { 57.29577951308232 mul sin } def
/a 0.0 def
/b 1.0 def
/tol 0.000000001 def
/bisect-current-interval {
numer 2 mul /numer exch def
denom 2 mul /denom exch def
} def
/go-to-next-interval {
numer 1 add /numer exch def
{
numer 2 mod 1 eq { exit } if
numer 2 idiv /numer exch def
denom 2 idiv /denom exch def
} loop
} def
/there-are-no-more-intervals {
numer denom eq
denom 1 eq
and
} def
/x0 {
numer denom div
dup 1 exch sub
a mul
exch b mul
add
} def
/x1 {
x0 0.75 mul
x4 0.25 mul
add
} def
/x2 {
x0 0.5 mul
x4 0.5 mul
add
} def
/x3 {
x0 0.25 mul
x4 0.75 mul
add
} def
/x4 {
numer 1 add denom div
dup 1 exch sub
a mul
exch b mul
add
} def
/simpson-rule-thrice {
x0 f /y0 exch def
x1 f /y1 exch def
x2 f /y2 exch def
x3 f /y3 exch def
x4 f /y4 exch def
x4 x0 sub 6.0 div /h04 exch def
x2 x0 sub 6.0 div /h02 exch def
x4 x2 sub 6.0 div /h24 exch def
y2 4.0 mul y0 add y4 add h04 mul /whole exch def
y1 4.0 mul y0 add y2 add h02 mul /left exch def
y3 4.0 mul y2 add y4 add h24 mul /right exch def
} def
/adaptive-simpson {
/numer 0 def
/denom 1 def
/sum 0.0 def
{
there-are-no-more-intervals { exit } if
simpson-rule-thrice
left right add whole sub /delta exch def
x4 x0 sub tol mul /tol0 exch def
x2 x0 sub tol mul /tol1 exch def
tol1 tol0 eq delta abs 15.0 tol0 mul le or
{
left right add delta 15.0 div add
sum add /sum exch def
go-to-next-interval
}
{
bisect-current-interval
} ifelse
} loop
} def
/Times-Roman findfont 12 scalefont setfont
10 10 moveto (estimate of integral of sin x dx from 0 to 1: ) show
adaptive-simpson
sum 20 string cvs show
showpage
origstate restore
%%EOF
</syntaxhighlight>
{{out}}
I ran the program thus:
<pre>magick convert -flatten adaptive_simpson_task.eps adaptive_simpson_task.png</pre>
And obtained:
[[File:Adaptive simpson task EPS.png|none|alt=Output from the PostScript for Adaptive Simpson. It shows 0.459698]]
=={{header|Prolog}}==
| |||