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

15,130 bytes added ,  1 month ago
Added FreeBASIC
(Added PostScript in front of →‎{{header|Prolog}})
(Added FreeBASIC)
(14 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 530 ⟶ 691:
{{works with|GnuCOBOL|3.1.2.0}}
 
{{improve|COBOL|All the numeric fields are unsigned in this program, so dalta < 0 will never be true - the pictures should have a leading s.}}
Rather than do recursions, I compute the x-intervals iteratively. No attempt is made to avoid recomputing values of sin(x).
 
Line 545 ⟶ 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 591 ⟶ 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)
Line 667 ⟶ 823:
{{out}}
<pre>$ cobc -std=cobol2014 -x adaptive_simpson_task.cob && ./adaptive_simpson_task
0.459697694</pre>
 
=={{header|Common Lisp}}==
Line 945 ⟶ 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,200 ⟶ 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,305 ⟶ 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,312 ⟶ 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,318 ⟶ 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,330 ⟶ 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,341 ⟶ 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 2,196 ⟶ 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,282 ⟶ 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/342791...363437"