numerical-methods-guide

Numerical Methods Guide

A skill for applying numerical methods in scientific computing and research. Covers root finding, numerical integration, ODE solvers, optimization, interpolation, and error analysis with practical implementations in Python.

Root Finding

Newton's Method and Alternatives

import numpy as np


def newton_method(f, df, x0: float, tol: float = 1e-10,
                  max_iter: int = 100) -> dict:
    """
    Newton's method for finding roots of f(x) = 0.

    Args:
        f: Function whose root we seek
        df: Derivative of f
        x0: Initial guess
        tol: Convergence tolerance
        max_iter: Maximum iterations
    """
    x = x0
    history = [x]

    for i in range(max_iter):
        fx = f(x)
        dfx = df(x)

        if abs(dfx) < 1e-15:
            return {"root": x, "converged": False,
                    "reason": "Zero derivative encountered"}

        x_new = x - fx / dfx
        history.append(x_new)

        if abs(x_new - x) < tol:
            return {
                "root": x_new,
                "converged": True,
                "iterations": i + 1,
                "f_at_root": f(x_new),
                "convergence": "quadratic"
            }

        x = x_new

    return {"root": x, "converged": False, "reason": "Max iterations reached"}

Method Selection Guide

Numerical Integration

Quadrature Methods

from scipy import integrate


def numerical_integration_comparison(f, a: float, b: float) -> dict:
    """
    Compare numerical integration methods.

    Args:
        f: Function to integrate
        a: Lower bound
        b: Upper bound
    """
    # Adaptive Gaussian quadrature (recommended default)
    quad_result, quad_error = integrate.quad(f, a, b)

    # Simpson's rule (fixed-point)
    n_points = 101
    x = np.linspace(a, b, n_points)
    simps_result = integrate.simpson(f(x), x=x)

    # Romberg integration
    romb_result = integrate.romberg(f, a, b)

    return {
        "quad": {"value": quad_result, "error_estimate": quad_error},
        "simpson": {"value": simps_result, "n_points": n_points},
        "romberg": {"value": romb_result},
        "recommendation": (
            "Use scipy.integrate.quad for most cases. "
            "It adaptively chooses points for accuracy."
        )
    }