Rosetta Code

Numerical integration/Adaptive Simpson's method: Difference between revisions

Page Discussion

Numerical integration/Adaptive Simpson's method (view source)

Revision as of 16:59, 10 April 2024

19,216 bytes added ,  1 month ago
Added FreeBASIC
(Added ObjectIcon ahead of OCaml →‎{{header|OCaml}})
(Added FreeBASIC)
(18 intermediate revisions by 5 users not shown)
Line 83:
Simpson's integration of sine from 0 to 1 = 0.459697694
</pre>
 
=={{header|Ada}}==
 
<syntaxhighlight lang="ada">
with ada.text_io; use ada.text_io;
with ada.numerics; use ada.numerics;
 
with ada.numerics.elementary_functions;
use ada.numerics.elementary_functions;
 
procedure adaptive_simpson_task is
 
subtype flt is float;
type function_flt_to_flt is
access function (x : in flt) return flt;
 
procedure simpson_rule (f : in function_flt_to_flt;
a, fa, b, fb : in flt;
m, fm, quadval : out flt) is
begin
m := 0.5 * (a + b);
fm := f(m);
quadval := ((b - a) / 6.0) * (fa + (4.0 * fm) + fb);
end;
 
function recursive_simpson (f : in function_flt_to_flt;
a, fa, b, fb : in flt;
tol, whole, m, fm : in flt;
depth : in integer) return flt is
lm, flm, left : flt;
rm, frm, right : flt;
diff, tol2, quadval : flt;
begin
simpson_rule (f, a, fa, m, fm, lm, flm, left);
simpson_rule (f, m, fm, b, fb, rm, frm, right);
diff := left + right - whole;
tol2 := 0.5 * tol;
if depth <= 0 or tol2 = tol or abs (diff) <= 15.0 * tol then
quadval := left + right + (diff / 15.0);
else
quadval := recursive_simpson (f, a, fa, m, fm, tol2,
left, lm, flm, depth - 1)
+ recursive_simpson (f, m, fm, b, fb, tol2,
right, rm, frm, depth - 1);
end if;
return quadval;
end;
 
function quad_asr (f : in function_flt_to_flt;
a, b, tol : in flt;
depth : in integer) return flt is
fa, fb, m, fm, whole : flt;
begin
fa := f(a);
fb := f(b);
simpson_rule (f, a, fa, b, fb, m, fm, whole);
return recursive_simpson (f, a, fa, b, fb, tol,
whole, m, fm, depth);
end;
 
function sine (x : in flt) return flt is
begin
return sin (x);
end;
 
quadval : flt;
 
begin
quadval := quad_asr (sine'access, 0.0, 1.0, 1.0e-5, 100);
put ("estimate of ∫ sin x dx from 0 to 1:");
put_line (quadval'image);
end;
 
-- local variables:
-- mode: indented-text
-- tab-width: 2
-- end:
</syntaxhighlight>
 
{{out}}
<pre>$ gnatmake -q adaptive_simpson_task.adb && ./adaptive_simpson_task
estimate of ∫ sin x dx from 0 to 1: 4.59698E-01</pre>
 
=={{header|ALGOL 60}}==
Line 478 ⟶ 560:
estimate of ∫ sin x dx from 0 to 1, by call to a compiled function: 0.459698
estimate of ∫ sin x dx from 0 to 1, by call to the template: 0.459698</pre>
 
=={{header|AWK}}==
<syntaxhighlight lang="awk">
BEGIN {
printf ("estimate of ∫ sin x dx from 0 to 1: %f\n",
quad_asr(0.0, 1.0, 1.0e-9, 100))
}
 
function f(x)
{
return sin(x)
}
 
function midpoint(a, b)
{
return 0.5 * (a + b)
}
 
function simpson_rule(a, fa, b, fb, fm)
{
return ((b - a) / 6.0) * (fa + (4.0 * fm) + fb)
}
 
function recursive_simpson(a, fa, b, fb, tol, whole, m, fm, depth,
# Local variables:
lm, flm, left,
rm, frm, right,
delta, tol_,
leftval, rightval, quadval)
{
lm = midpoint(a, m)
flm = f(lm)
left = simpson_rule(a, fa, m, fm, flm)
 
rm = midpoint(m, b)
frm = f(rm)
right = simpson_rule(m, fm, b, fb, frm)
 
delta = left + right - whole
tol_ = 0.5 * tol
 
if (depth <= 0 || tol_ == tol || abs(delta) <= 15.0 * tol)
quadval = left + right + (delta / 15.0)
else
{
leftval = recursive_simpson(a, fa, m, fm, tol_,
left, lm, flm, depth - 1)
rightval = recursive_simpson(m, fm, b, fb, tol_,
right, rm, frm, depth - 1)
quadval = leftval + rightval
}
return quadval
}
 
function quad_asr(a, b, tol, depth,
# Local variables:
fa, fb, whole, m, fm)
{
fa = f(a)
fb = f(b)
 
m = midpoint(a, b)
fm = f(m)
whole = simpson_rule(a, fa, b, fb, fm)
 
return recursive_simpson(a, fa, b, fb, tol,
whole, m, fm, depth)
}
 
function abs(x)
{
return (x < 0 ? -x : x)
}
</syntaxhighlight>
 
{{out}}
Just about any modern Awk should work: nawk, mawk, gawk, busybox awk, goawk, etc. But the ancient Old Awk does not work.
<pre>$ nawk -f adaptive_simpson_task.awk
estimate of ∫ sin x dx from 0 to 1: 0.459698</pre>
 
=={{header|C}}==
Line 544 ⟶ 705:
data division.
working-storage section.
01 numer picture 9s9(25). *> 25 decimal digits.
01 denom picture 9s9(25).
01 numer2 picture 9s9(25).
01 a picture 9s9(5)V9(20). *> 5+20 decimal digits.
01 b picture 9s9(5)V9(20).
01 tol picture 9s9(5)V9(20).
01 x picture 9s9(5)V9(20).
01 y picture 9s9(5)V9(20).
01 x0 picture 9s9(5)V9(20).
01 y0 picture 9s9(5)V9(20).
01 x1 picture 9s9(5)V9(20).
01 y1 picture 9s9(5)V9(20).
01 x2 picture 9s9(5)V9(20).
01 y2 picture 9s9(5)V9(20).
01 x3 picture 9s9(5)V9(20).
01 y3 picture 9s9(5)V9(20).
01 x4 picture 9s9(5)V9(20).
01 y4 picture 9s9(5)V9(20).
01 ruleval picture 9s9(5)V9(20).
01 ruleval0 picture 9s9(5)V9(20).
01 ruleval1 picture 9s9(5)V9(20).
01 delta picture 9s9(5)V9(20).
01 abs-delta picture 9(5)V9(20). *> unsigned.
01 tol0 picture 9s9(5)V9(20).
01 tol1 picture 9s9(5)V9(20).
01 delta15 picture 9s9(5)V9(20).
01 stepval picture 9s9(5)V9(20).
01 quadval picture 9s9(5)V9(20).
01 result picture 9V9-9.9(9).
procedure division.
Line 590 ⟶ 751:
perform simpson-rule-thrice
compute delta = ruleval0 + ruleval1 - ruleval
ifcompute abs-delta <= 0 thendelta
compute abs-delta = -delta
else
compute abs-delta = delta
end-if
compute tol0 = tol * (x4 - x0)
compute tol1 = tol * (x2 - x0)
if (tol1 = tol0) or (abs-delta <= 15 * tol0) then
compute delta15 = delta / 15
compute stepval = ruleval0 + ruleval1 + delta15
Line 666 ⟶ 823:
{{out}}
<pre>$ cobc -std=cobol2014 -x adaptive_simpson_task.cob && ./adaptive_simpson_task
0.459697694</pre>
 
=={{header|Common Lisp}}==
Line 793 ⟶ 950:
Simpson's rule integration of sin from 0 to 1 is: 0.4596976941573994
</pre>
 
=={{header|Forth}}==
{{works with|Gforth|0.7.3}}
 
<syntaxhighlight lang="forth">
\ -*- mode: forth -*-
\
\ The following was written for Gforth, using its implementation of
\ local variables and of f=
\
 
Defer f
 
: simpson-rule { F: a F: fa F: b F: fb -- m fm quadval }
a b f+ 0.5e0 f* \ -- m
fdup f \ -- m fm
fdup 4.0e0 f* fa f+ fb f+
b a f- 6.0e0 f/ f* \ -- m fm quadval
;
 
: recursive-simpson
{ F: a F: fa F: b F: fb F: tol F: whole F: m F: fm depth -- quadval }
a fa m fm simpson-rule { F: lm F: flm F: left }
m fm b fb simpson-rule { F: rm F: frm F: right }
left right f+ whole f- { F: delta }
tol 0.5e0 f* { F: tol_ }
depth 0 <=
tol_ tol f= or
delta fabs 15.0e0 tol f* f<= or
if
left right f+ delta 15.0e0 f/ f+
else
a fa m fm tol_ left lm flm depth 1 - recurse
m fm b fb tol_ right rm frm depth 1 - recurse
f+
then
;
 
: quad-asr { F: a F: b F: tol depth -- quadval }
a f { F: fa }
b f { F: fb }
a fa b fb simpson-rule { F: m F: fm F: whole }
a fa b fb tol whole m fm depth recursive-simpson
;
 
:noname fsin ; IS f
 
." estimate of ∫ sin x dx from 0 to 1: "
0.0e0 1.0e0 1.0e-12 100 quad-asr f.
cr
</syntaxhighlight>
 
{{out}}
<pre>$ gforth adaptive_simpson_task.fs -e bye
estimate of ∫ sin x dx from 0 to 1: 0.45969769413186</pre>
 
=={{header|Fortran}}==
Line 889 ⟶ 1,101:
 
Such a bit of executable code on the stack often is not allowed on systems hardened against hackers. Otherwise you probably needn't be concerned, even though warnings are being generated both by the compiler ''and'' the linker. :)
 
=={{header|FreeBASIC}}==
{{trans|Wren}}
<syntaxhighlight lang="vbnet">Type QuadSimpsonMem
As Double m, fm, simp
End Type
 
Function quadSimpsonMem(f As Function(As Double) As Double, a As Double, fa As Double, b As Double, fb As Double) As QuadSimpsonMem
Dim As QuadSimpsonMem r
r.m = (a + b) / 2
r.fm = f(r.m)
r.simp = Abs(b - a) / 6 * (fa + 4 * r.fm + fb)
Return r
End Function
 
Function quadAsrRec(f As Function(As Double) As Double, a As Double, fa As Double, b As Double, fb As Double, eps As Double, whole As Double, m As Double, fm As Double) As Double
Dim As QuadSimpsonMem r1 = quadSimpsonMem(f, a, fa, m, fm)
Dim As QuadSimpsonMem r2 = quadSimpsonMem(f, m, fm, b, fb)
Dim As Double delta = r1.simp + r2.simp - whole
If Abs(delta) < eps * 15 Then Return r1.simp + r2.simp + delta / 15
Return quadAsrRec(f, a, fa, m, fm, eps / 2, r1.simp, r1.m, r1.fm) + quadAsrRec(f, m, fm, b, fb, eps / 2, r2.simp, r2.m, r2.fm)
End Function
 
Function quadAsr(f As Function(As Double) As Double, a As Double, b As Double, eps As Double) As Double
Dim As Double fa = f(a), fb = f(b)
Dim As QuadSimpsonMem r = quadSimpsonMem(f, a, fa, b, fb)
Return quadAsrRec(f, a, fa, b, fb, eps, r.simp, r.m, r.fm)
End Function
 
Function sinx(x As Double) As Double
Return Sin(x)
End Function
 
Dim As Double a = 0, b = 1
Dim As Double result = quadAsr(@sinx, a, b, 1e-09)
Print "Simpson's integration of sine from"; a; " to"; b; " ="; result
 
Sleep</syntaxhighlight>
{{out}}
<pre>Simpson's integration of sine from 0 to 1 = 0.4596976941317858</pre>
 
=={{header|Go}}==
Line 1,144 ⟶ 1,396:
Simpson's rule integration of sin from 0 to 1 is: 0.45969769415739936
</pre>
 
=={{header|K}}==
{{works with|ngn/k}}<syntaxhighlight lang=K>simpsonsRule:{[f;a;b] m: (a+b)%2; h: b-a; m, (h%6)*f[a]+f[b]+4*f[m]}
 
recursiveSimpson:{[f;a;b;tol;whole;m;depth]
(lm;left): simpsonsRule[f;a;m]
(rm;right): simpsonsRule[f;m;b]
delta: left+right-whole
tol2: tol%2
$[(1>depth)|(tol2=tol)|~(15*tol)<delta|-delta
left+right+delta%15
recursiveSimpson[f;a;m;tol2;left;lm;depth-1]+recursiveSimpson[f;m;b;tol2;right;rm;depth-1]]}
 
quadAsr:{[f;a;b;tol;depth]; (m;whole): simpsonsRule[f;a;b]; recursiveSimpson[f;a;b;tol;whole;m;depth]}</syntaxhighlight>
 
<syntaxhighlight lang=K>quadAsr[`sin;0;1;1e-12;100]
0.45969769413186023</syntaxhighlight>
 
=={{header|Kotlin}}==
Line 1,249 ⟶ 1,518:
(0.45969769413186023) // reference value
</syntaxhighlight>
 
=={{header|Mercury}}==
<syntaxhighlight lang="mercury">
%%% -*- mode: mercury; prolog-indent-width: 2; -*-
 
:- module adaptive_simpson_task_mercury.
 
:- interface.
:- import_module io.
:- pred main(io::di, io::uo) is det.
 
:- implementation.
:- import_module float.
:- import_module int.
:- import_module math.
:- import_module string.
 
:- pred simpson_rule((func(float) = float)::in,
float::in, float::in, float::in, float::in,
float::out, float::out, float::out) is det.
simpson_rule(F, A, FA, B, FB, M, FM, QuadVal) :-
FM = F(M), (M = 0.5 * (A + B)),
(QuadVal = ((B - A) / 6.0) * (FA + (4.0 * FM) + FB)).
 
:- func recursive_simpson(func(float) = float,
float, float, float, float,
float, float, float, float, int) = float.
recursive_simpson(F, A, FA, B, FB, Tol,
Whole, M, FM, Depth) = QuadVal :-
simpson_rule(F, A, FA, M, FM, LM, FLM, Left),
simpson_rule(F, M, FM, B, FB, RM, FRM, Right),
(Left + Right - Whole = Delta),
(0.5 * Tol = Tol_),
(if (Depth =< 0; Tol_ = Tol; abs(Delta) =< 15.0 * Tol)
then (QuadVal = Left + Right + (Delta / 15.0))
else (QuadVal = recursive_simpson(F, A, FA, M, FM, Tol_,
Left, LM, FLM, Depth - 1)
+ recursive_simpson(F, M, FM, B, FB, Tol_,
Right, RM, FRM, Depth - 1))).
 
:- func quad_asr(func(float) = float,
float, float, float, int) = float.
quad_asr(F, A, B, Tol, Depth) = QuadVal :-
F(A) = FA, F(B) = FB,
simpson_rule(F, A, FA, B, FB, M, FM, Whole),
(QuadVal = recursive_simpson(F, A, FA, B, FB, Tol,
Whole, M, FM, Depth)).
 
main(!IO) :-
print("estimate of ∫ sin x dx from 0 to 1: ", !IO),
(QuadVal = quad_asr(sin, 0.0, 1.0, 1.0e-9, 100)),
print(from_float(QuadVal), !IO),
nl(!IO).
 
:- end_module adaptive_simpson_task_mercury.
</syntaxhighlight>
 
{{out}}
<pre>$ mmc --no-verbose-make --make --use-subdirs adaptive_simpson_task_mercury && ./adaptive_simpson_task_mercury
estimate of ∫ sin x dx from 0 to 1: 0.4596976941317858</pre>
 
=={{header|Modula-2}}==
{{works with|GCC|13.1.0}}
 
<syntaxhighlight lang="modula2">
MODULE adaptive_simpson_task_Modula2;
 
(* ISO Modula-2 libraries. *)
IMPORT RealMath;
IMPORT SRealIO;
IMPORT STextIO;
 
TYPE functionReal2Real = PROCEDURE (REAL) : REAL;
 
PROCEDURE simpsonRule (f : functionReal2Real;
a, fa, b, fb : REAL;
VAR m, fm, quadVal : REAL);
BEGIN
m := 0.5 * (a + b);
fm := f(m);
quadVal := ((b - a) / 6.0) * (fa + (4.0 * fm) + fb)
END simpsonRule;
 
PROCEDURE recursiveSimpson (f : functionReal2Real;
a, fa, b, fb, tol, whole, m, fm : REAL;
depth : INTEGER) : REAL;
VAR lm, flm, left : REAL;
rm, frm, right : REAL;
delta, tol2, quadVal : REAL;
BEGIN
simpsonRule (f, a, fa, m, fm, lm, flm, left);
simpsonRule (f, m, fm, b, fb, rm, frm, right);
delta := left + right - whole;
tol2 := 0.5 * tol;
IF (depth <= 0) OR (tol2 = tol) OR (ABS (delta) <= 15.0 * tol) THEN
quadVal := left + right + (delta / 15.0)
ELSE
quadVal := (recursiveSimpson (f, a, fa, m, fm, tol2,
left, lm, flm, depth - 1)
+ recursiveSimpson (f, m, fm, b, fb, tol2,
right, rm, frm, depth - 1))
END;
RETURN quadVal;
END recursiveSimpson;
 
PROCEDURE quadASR (f : functionReal2Real;
a, b, tol : REAL;
depth : INTEGER) : REAL;
VAR fa, fb, m, fm, whole : REAL;
BEGIN
fa := f(a); fb := f(b);
simpsonRule (f, a, fa, b, fb, m, fm, whole);
RETURN recursiveSimpson (f, a, fa, b, fb, tol,
whole, m, fm, depth)
END quadASR;
 
BEGIN
STextIO.WriteString ('estimate of ∫ sin x dx from 0 to 1: ');
SRealIO.WriteReal (quadASR (RealMath.sin, 0.0, 1.0,
0.000000001, 100),
10);
STextIO.WriteLn;
END adaptive_simpson_task_Modula2.
</syntaxhighlight>
 
{{out}}
<pre>$ gm2 -fiso adaptive_simpson_task_Modula2.mod && ./a.out
estimate of ∫ sin x dx from 0 to 1: 0.45969769</pre>
 
=={{header|Nim}}==
Line 1,256 ⟶ 1,653:
type Func = float -> float
 
funcproc quadSimpsonsMem(f: Func; a, fa, b, fb: float): tuple[m, fm, val: float] =
## Evaluates the Simpson's Rule, also returning m and f(m) to reuse
result.m = (a + b) / 2
Line 1,262 ⟶ 1,659:
result.val = abs(b - a) * (fa + 4 * result.fm + fb) / 6
 
funcproc quadAsr(f: Func; a, fa, b, fb, eps, whole, m, fm: float): float =
## Efficient recursive implementation of adaptive Simpson's rule.
## Function values at the start, middle, end of the intervals are retained.
Line 1,274 ⟶ 1,671:
f.quadAsr(m, fm, b, fb, eps / 2, right, rm, frm)
 
funcproc quadAsr(f: Func; a, b, eps: float): float =
## Integrate f from a to b using Adaptive Simpson's Rule with max error of eps.
let fa = f(a)
Line 1,285 ⟶ 1,682:
{{out}}
<pre>Simpson's integration of sine from 0 to 1 = 0.4596976941317859</pre>
 
=={{header|Oberon-2}}==
{{trans|Modula-2}}
{{works with|Oxford Oberon-2 Compiler|3.3.0}}
 
<syntaxhighlight lang="modula2">
(* -*- mode: Modula2 -*- *)
 
MODULE adaptiveSimpsonTask;
 
IMPORT Math, Out;
 
TYPE functionReal2Real = PROCEDURE (x : REAL) : REAL;
 
PROCEDURE simpsonRule (f : functionReal2Real;
a, fa, b, fb : REAL;
VAR m, fm, quadVal : REAL);
BEGIN
m := 0.5 * (a + b);
fm := f(m);
quadVal := ((b - a) / 6.0) * (fa + (4.0 * fm) + fb)
END simpsonRule;
 
PROCEDURE recursiveSimpson (f : functionReal2Real;
a, fa, b, fb, tol, whole, m, fm : REAL;
depth : INTEGER) : REAL;
VAR lm, flm, left : REAL;
rm, frm, right : REAL;
delta, tol2, quadVal : REAL;
BEGIN
simpsonRule (f, a, fa, m, fm, lm, flm, left);
simpsonRule (f, m, fm, b, fb, rm, frm, right);
delta := left + right - whole;
tol2 := 0.5 * tol;
IF (depth <= 0) OR (tol2 = tol) OR (ABS (delta) <= 15.0 * tol) THEN
quadVal := left + right + (delta / 15.0)
ELSE
quadVal := (recursiveSimpson (f, a, fa, m, fm, tol2,
left, lm, flm, depth - 1)
+ recursiveSimpson (f, m, fm, b, fb, tol2,
right, rm, frm, depth - 1))
END;
RETURN quadVal
END recursiveSimpson;
 
PROCEDURE quadASR (f : functionReal2Real;
a, b, tol : REAL;
depth : INTEGER) : REAL;
VAR fa, fb, m, fm, whole : REAL;
BEGIN
fa := f(a); fb := f(b);
simpsonRule (f, a, fa, b, fb, m, fm, whole);
RETURN recursiveSimpson (f, a, fa, b, fb, tol,
whole, m, fm, depth)
END quadASR;
 
BEGIN
Out.String ('estimate of ∫ sin x dx from 0 to 1: ');
Out.Real (quadASR (Math.Sin, 0.0, 1.0, 0.000001, 100), 7);
Out.Ln
END adaptiveSimpsonTask.
</syntaxhighlight>
 
{{out}}
<pre>$ obc adaptiveSimpsonTask.Mod && ./a.out
estimate of ∫ sin x dx from 0 to 1: 0.459698
</pre>
 
=={{header|ObjectIcon}}==
Line 1,605 ⟶ 2,069:
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}}==
Line 2,014 ⟶ 2,604:
<pre>$ gsi adaptive_simpson_task.scm
estimated definite integral of sin(x) for x from 0 to 1: .4596976941317858</pre>
 
=={{header|Shen}}==
<syntaxhighlight lang="shen">
(define simpson-rule
F A FA B FB ->
(let M (* 0.5 (+ A B))
FM (F M)
QuadVal (* (/ (- B A) 6.0) (+ FA (* 4.0 FM) FB))
[M FM QuadVal]))
 
(define recursive-simpson
F A FA B FB Tol Whole M FM Depth ->
(let TheLeftStuff (simpson-rule F A FA M FM)
TheRightStuff (simpson-rule F M FM B FB)
(recursive-simpson-united F A FA B FB Tol Whole M FM Depth
TheLeftStuff TheRightStuff)))
 
(define recursive-simpson-united
F A FA B FB Tol Whole M FM Depth
[LM FLM Left] [RM FRM Right] ->
(let Delta (- (+ Left Right) Whole)
Tol_ (* 0.5 Tol)
(if (or (<= Depth 0)
(= Tol_ Tol)
(<= (abs Delta) (* 15.0 Tol)))
(+ Left Right (/ Delta 15.0))
(+ (recursive-simpson F A FA M FM Tol_
Left LM FLM (- Depth 1))
(recursive-simpson F M FM B FB Tol_
Right RM FRM (- Depth 1))))))
 
(define quad-asr
F A B Tol Depth ->
(let FA (F A)
FB (F B)
TheWholeStuff (simpson-rule F A FA B FB)
(quad-asr-part-two F A FA B FB Tol Depth TheWholeStuff)))
 
(define quad-asr-part-two
F A FA B FB Tol Depth [M FM Whole] ->
(recursive-simpson F A FA B FB Tol Whole M FM Depth))
 
\\ The Shen standard library contains only an unserious
\\ implementation of sin(x), so I might as well toss together
\\ my own! Thus: an uneconomized Maclaurin expansion. Gotten
\\ via "taylor(sin(x),x,0,20);" in Maxima. Using up to x**13
\\ should give an error bound on the order of 1e-12, for
\\ 0 <= x <= 1. That is much better than we need.
(define sin
X -> (compute-sin X (* X X)))
(define compute-sin
X XX ->
(* X (+ 1.0
(* XX (+ (/ -1.0 6.0)
(* XX (+ (/ 1.0 120.0)
(* XX (+ (/ -1.0 5040.0)
(* XX (+ (/ 1.0 362880.0)
(* XX (+ (/ -1.0 39916800.0)
(* XX (/ 1.0 6227020800.0)))))))))))))))
 
(define abs
\\ Like sin, abs is in the standard library, but it is very
\\ easy to write.
X -> (if (< X 0) (- 0 X) X))
 
(output "estimate of the definite integral of sin x dx from 0 to 1: ")
(print (quad-asr (lambda x (sin x)) 0.0 1.0 1.0e-9 100))
(nl)
 
\\ local variables:
\\ mode: indented-text
\\ tab-width: 2
\\ end:
</syntaxhighlight>
 
{{out}}
<pre>$ shen-scheme script adaptive_simpson_task.shen
estimate of the definite integral of sin x dx from 0 to 1: 0.4596976941318334</pre>
 
=={{header|Sidef}}==
Line 2,100 ⟶ 2,768:
{{trans|Go}}
{{libheader|Wren-fmt}}
<syntaxhighlight lang="ecmascriptwren">import "./fmt" for Fmt
 
/* "structured" adaptive version, translated from Racket */
Jjuanhdez
2,122

edits

Retrieved from "https://rosettacode.org/wiki/Special:MobileDiff/342774...363437"