Numerical integration/Adaptive Simpson's method: Difference between revisions
Content added Content deleted
(→{{header|COBOL}}: Bug fix.) |
(Added PostScript in front of →{{header|Prolog}}) |
||
Line 1,661: | Line 1,661: | ||
Simpson's integration of sine from 0 to 1 = 0.459698 |
Simpson's integration of sine from 0 to 1 = 0.459698 |
||
</pre> |
</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}}== |
=={{header|Prolog}}== |