naive implementation of adaptive Simpson integration

Task #128 - implement simple numerical integration

CMakeLists.txt

@@ -193,6 +193,8 @@ set (YACAS_SCRIPTS
     scripts/multivar.rep/sparsetree.ys.def
     scripts/newly.rep/code.ys
     scripts/newly.rep/code.ys.def
+    scripts/nintegrate.rep/code.ys
+    scripts/nintegrate.rep/code.ys.def
     scripts/numbers.rep/GaussianIntegers.ys
     scripts/numbers.rep/GaussianIntegers.ys.def
     scripts/numbers.rep/NumberTheory.ys

scripts/Makefile.am

@@ -19,7 +19,8 @@ SCRIPTFILES1 = 	        										    \
 	functional.rep/code.ys functional.rep/code.ys.def functional.rep/om.ys				    \
 	graph.rep/code.ys graph.rep/code.ys.def								    \
 	html.rep/code.ys html.rep/code.ys.def 								    \
-	integrate.rep/code.ys integrate.rep/code.ys.def integrate.rep/om.ys   \
+	integrate.rep/code.ys integrate.rep/code.ys.def integrate.rep/om.ys				    \
+	nintegrate.rep/code.ys nintegrate.rep/code.ys.def						    \
 	io.rep/code.ys io.rep/code.ys.def io.rep/print.ys io.rep/formula.ys io.rep/errors.ys 		    \
 	io.rep/defaultprint.ys io.rep/defaultprint.ys.def 						    \
 	limit.rep/code.ys limit.rep/code.ys.def limit.rep/om.ys						    \

scripts/nintegrate.rep/code.ys

@@ -0,0 +1,42 @@
+100 # AdaptiveSimpson(_f, var_IsVariable, a_IsRationalOrNumber, b_IsRationalOrNumber, tol_IsRationalOrNumber) <-- [
+    Local(c, fa, fb, fc, h, q, q1, q2, stp);
+
+    c := (a + b) / 2;
+
+    fa := N(`WithValue(@var, a, f));
+    fc := N(`WithValue(@var, c, f));
+    fb := N(`WithValue(@var, b, f));
+
+    h := (b - a);
+    q := (b - a) / 6 * (fa + 4 * fc + fb);
+
+    stp(f, var, a, b, c, fa, fb, fc, q0, tol) := [
+        Local(d, e, fd, fe, q, q1, q2);
+
+        d := (a + c) / 2;
+        e := (c + b) / 2;
+
+        fd := N(`WithValue(@var, d, f));
+        fe := N(`WithValue(@var, e, f));
+
+        q1 := (c - a) / 6 * (fa + 4 * fd + fc);
+        q2 := (b - c) / 6 * (fc + 4 * fe + fb);
+
+        q := N(q1 + q2);
+
+        If (Abs(q - q0) > tol, [
+            q1 := stp (f, var, a, c, d, fa, fc, fd, q1, tol);
+            q2 := stp (f, var, c, b, e, fc, fb, fe, q2, tol);
+            q := N(q1 + q2);
+        ]);
+
+        q;
+    ];
+
+    stp(f, var, a, b, c, fa, fb, fc, q, tol);
+];
+
+
+100 # (NIntegrate(var_IsVariable, a_IsRationalOrNumber, b_IsRationalOrNumber)(_f)) <-- [
+    AdaptiveSimpson(f, var, a, b, 1e-8);
+];

scripts/nintegrate.rep/code.ys.def

@@ -0,0 +1,2 @@
+NIntegrate
+}

scripts/packages.ys

@@ -28,6 +28,7 @@ Defun(DefFileList,{}) {
     "multivar.rep/code.ys",
     "multivar.rep/sparsetree.ys",
     "newly.rep/code.ys",
+    "nintegrate.rep/code.ys",
     "numbers.rep/GaussianIntegers.ys",
     "numbers.rep/NumberTheory.ys",
     "numbers.rep/code.ys",

scripts/stdopers.ys

@@ -45,6 +45,7 @@ Postfix("!", 30);
 Postfix("!!", 30);
 Infix("***", 50);
 Bodied("Integrate",60000);
+Bodied("NIntegrate",60000);
 
 Bodied("Limit",60000);