45.4 Numerical Methods

Implement fundamental algorithms for root finding, integration, and interpolation.

45.4.1 Newton-Raphson Root Finding

Find x where f(x) = 0 using iteration: x_{n+1} = x_n - f(x_n)/f'(x_n)

interface Func { double f(double x); double df(double x); }

double newtonRaphson(Func func, double x0, double tol, int maxIter) {
  double x = x0;
  for (int i = 0; i < maxIter; i++) {
    double fx = func.f(x);
    if (Math.abs(fx) < tol) return x;
    double dfx = func.df(x);
    if (Math.abs(dfx) < 1e-14) throw new ArithmeticException("Derivative too small");
    x = x - fx / dfx;
  }
  throw new RuntimeException("Did not converge");
}

// Example: find sqrt(2) via f(x) = x² - 2
Func sqrt2 = new Func() {
  public double f(double x) { return x*x - 2; }
  public double df(double x) { return 2*x; }
};
double root = newtonRaphson(sqrt2, 1.0, 1e-10, 100); // ~1.414

45.4.2 Bisection Method

Guaranteed convergence if f continuous and f(a)·f(b) < 0.

double bisection(java.util.function.DoubleUnaryOperator f, double a, double b, double tol) {
  double fa = f.applyAsDouble(a);
  double fb = f.applyAsDouble(b);
  if (fa * fb > 0) throw new IllegalArgumentException("f(a) and f(b) must have opposite signs");

  while (b - a > tol) {
    double c = (a + b) / 2;
    double fc = f.applyAsDouble(c);
    if (Math.abs(fc) < tol) return c;
    if (fa * fc < 0) { b = c; fb = fc; }
    else { a = c; fa = fc; }
  }
  return (a + b) / 2;
}

45.4.3 Numerical Integration (Trapezoidal Rule)

∫ₐᵇ f(x) dx ≈ h/2 · [f(a) + 2·Σf(xᵢ) + f(b)]

double trapezoid(java.util.function.DoubleUnaryOperator f, double a, double b, int n) {
  double h = (b - a) / n;
  double sum = 0.5 * (f.applyAsDouble(a) + f.applyAsDouble(b));
  for (int i = 1; i < n; i++) {
    sum += f.applyAsDouble(a + i * h);
  }
  return h * sum;
}

// Example: ∫₀¹ x² dx = 1/3
double integral = trapezoid(x -> x*x, 0, 1, 1000); // ~0.333

45.4.4 Simpson's Rule

Higher accuracy: ∫ₐᵇ f(x) dx ≈ h/3 · [f(a) + 4·Σf(odd) + 2·Σf(even) + f(b)]

double simpson(java.util.function.DoubleUnaryOperator f, double a, double b, int n) {
  if (n % 2 != 0) throw new IllegalArgumentException("n must be even");
  double h = (b - a) / n;
  double sum = f.applyAsDouble(a) + f.applyAsDouble(b);
  for (int i = 1; i < n; i++) {
    double x = a + i * h;
    sum += (i % 2 == 0 ? 2 : 4) * f.applyAsDouble(x);
  }
  return h * sum / 3;
}

45.4.5 Linear Interpolation

double lerp(double x0, double y0, double x1, double y1, double x) {
  return y0 + (y1 - y0) * (x - x0) / (x1 - x0);
}

// For tabulated data
double interpolate(double[] xs, double[] ys, double x) {
  for (int i = 0; i < xs.length - 1; i++) {
    if (x >= xs[i] && x <= xs[i+1]) {
      return lerp(xs[i], ys[i], xs[i+1], ys[i+1], x);
    }
  }
  throw new IllegalArgumentException("x out of range");
}