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

10,263 bytes added ,  1 month ago
Added FreeBASIC
(Awk before C)
(Added FreeBASIC)
(11 intermediate revisions by 4 users not shown)
Line 561:
estimate of ∫ sin x dx from 0 to 1, by call to the template: 0.459698</pre>
 
=={{header|AwkAWK}}==
<syntaxhighlight lang="awk">
BEGIN {
Line 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 706 ⟶ 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 752 ⟶ 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 828 ⟶ 823:
{{out}}
<pre>$ cobc -std=cobol2014 -x adaptive_simpson_task.cob && ./adaptive_simpson_task
0.459697694</pre>
 
=={{header|Common Lisp}}==
Line 1,106 ⟶ 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,483 ⟶ 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,490 ⟶ 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,496 ⟶ 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,508 ⟶ 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,519 ⟶ 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,374 ⟶ 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,460 ⟶ 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/342820...363437"