feat(math-engine): implement NumericalMethods with derivative, integral, and solve functions

src/core/math/NumericalMethods.js

@@ -0,0 +1,115 @@
+// src/core/math/NumericalMethods.js
+/**
+ * Numerical methods for the scientific calculator
+ * Provides derivative, integral, and solver functions
+ */
+
+const DEFAULT_TOLERANCE = 1e-6;
+const MAX_ITERATIONS = 100;
+const DEFAULT_H = 1e-5;
+
+/**
+ * Compute the derivative of an expression with respect to a variable using central difference.
+ * @param {Function} exprFunc - A function that takes a variable value and returns the expression value
+ * @param {number} x - The point at which to evaluate the derivative
+ * @param {number} h - The step size (optional, default 1e-5)
+ * @returns {number} The derivative value
+ */
+function derivative(exprFunc, x, h = DEFAULT_H) {
+  return (exprFunc(x + h) - exprFunc(x - h)) / (2 * h);
+}
+
+/**
+ * Compute the definite integral of an expression from a to b using adaptive Simpson's rule.
+ * @param {Function} exprFunc - A function that takes a variable value and returns the expression value
+ * @param {number} a - Lower limit
+ * @param {number} b - Upper limit
+ * @param {number} tolerance - Error tolerance (optional, default 1e-6)
+ * @param {number} maxRecursionDepth - Maximum recursion depth (optional, default 10)
+ * @returns {number} The integral value
+ */
+function integral(exprFunc, a, b, tolerance = DEFAULT_TOLERANCE, maxRecursionDepth = 10) {
+  // Helper function for Simpson's rule
+  function simpson(fa, fm, fb, width) {
+    return (width * (fa + 4 * fm + fb)) / 6;
+  }
+
+  // Recursive adaptive function
+  function adaptive(l, r, fl, fr, fmid, whole, depth) {
+    const leftMid = (l + (l + r) / 2) / 2;
+    const rightMid = ((l + r) / 2 + r) / 2;
+    const fleftmid = exprFunc(leftMid);
+    const frightmid = exprFunc(rightMid);
+    const left = simpson(fl, fleftmid, fmid, (r - l) / 4);
+    const right = simpson(fmid, frightmid, fr, (r - l) / 4);
+    if (depth >= maxRecursionDepth || Math.abs(left + right - whole) <= 15 * tolerance) {
+      return left + right + (left + right - whole) / 15;
+    }
+    return adaptive(l, (l + r) / 2, fl, fmid, fleftmid, left, depth + 1) +
+           adaptive((l + r) / 2, r, fmid, fr, frightmid, right, depth + 1);
+  }
+
+  const fa = exprFunc(a);
+  const fb = exprFunc(b);
+  const fmid = exprFunc((a + b) / 2);
+  const whole = simpson(fa, fmid, fb, b - a);
+  return adaptive(a, b, fa, fb, fmid, whole, 0);
+}
+
+/**
+ * Solve for a root of the expression f(x) = 0 using Newton-Raphson with secant fallback.
+ * @param {Function} exprFunc - A function that takes a variable value and returns the expression value
+ * @param {number} initialGuess - Initial guess for the root (optional, default 0)
+ * @param {number} tolerance - Convergence tolerance (optional, default 1e-10)
+ * @param {number} maxIterations - Maximum iterations (optional, default 100)
+ * @returns {number} The root value
+ */
+function solve(exprFunc, initialGuess = 0, tolerance = 1e-10, maxIterations = 100) {
+  let x0 = initialGuess;
+  let x1 = x0 + 0.1; // Second point for secant method
+  let f0 = exprFunc(x0);
+  let f1 = exprFunc(x1);
+
+  for (let i = 0; i < maxIterations; i++) {
+    // Check for convergence
+    if (Math.abs(f1) < tolerance) {
+      return x1;
+    }
+
+    // Avoid division by zero
+    if (Math.abs(f1 - f0) < Number.EPSILON) {
+      // Fallback to bisection if derivative is zero
+      // For simplicity, we'll just perturb and try again
+      x0 = x0 + 0.1;
+      x1 = x0 + 0.1;
+      f0 = exprFunc(x0);
+      f1 = exprFunc(x1);
+      continue;
+    }
+
+    // Newton-Raphson step (requires derivative, which we don't have directly)
+    // Instead, we'll use the secant method
+    const x2 = x1 - f1 * (x1 - x0) / (f1 - f0);
+    const f2 = exprFunc(x2);
+
+    // Check for convergence
+    if (Math.abs(f2) < tolerance) {
+      return x2;
+    }
+
+    // Shift for next iteration
+    x0 = x1;
+    x1 = x2;
+    f0 = f1;
+    f1 = f2;
+  }
+
+  // If we reach here, we didn't converge
+  throw new Error('Solver did not converge within the maximum number of iterations');
+}
+
+module.exports = {
+  derivative,
+  integral,
+  solve
+};