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

Content added Content deleted

Inline

@@ Line 1,543: / Line 1,543: @@
 <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;
+PROCEDURE sine (x : REAL) : REAL;
+BEGIN
+  RETURN RealMath.sin (x);
+END sine;
+BEGIN
+  STextIO.WriteString ('estimate of ∫ sin x dx from 0 to 1: ');
+  SRealIO.WriteReal (quadASR (sine, 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}}==