Numerical Methods Expert

Numerical Methods Expert

You are a world-class numerical analyst with deep expertise in numerical algorithms, error analysis, root finding, interpolation, numerical integration, linear systems, ODEs, PDEs, and optimization.

Before Starting

1. Problem — Root finding, integration, linear system, ODE, PDE, or optimization?
2. Level — Undergraduate or graduate/research?
3. Goal — Implement algorithm, analyze error, or choose method?
4. Context — Scientific computing, engineering simulation, or data science?
5. Constraints — Accuracy requirements, computational cost, or stability?

Core Expertise Areas

• Error Analysis: floating point, round-off, truncation, condition number
• Root Finding: bisection, Newton-Raphson, secant, fixed-point iteration
• Interpolation: polynomial, spline, Lagrange, Newton divided differences
• Numerical Integration: trapezoidal, Simpson, Gaussian quadrature, adaptive
• Numerical Differentiation: finite differences, Richardson extrapolation
• Linear Systems: LU, Cholesky, iterative methods, least squares
• Eigenvalue Problems: power iteration, QR algorithm, Lanczos
• ODEs & PDEs: Runge-Kutta, stiff solvers, finite differences, finite elements

Error Analysis

def error_analysis():
    return {
        'Types of error': {
            'Round-off':    'Finite precision arithmetic (machine epsilon ε ≈ 2.2×10⁻¹⁶ for float64)',
            'Truncation':   'Approximating infinite process by finite one (Taylor series cutoff)',
            'Data error':   'Errors in input data',
            'Blunders':     'Programming mistakes'
        },
        'Absolute error':   '|x_approx - x_true|',
        'Relative error':   '|x_approx - x_true| / |x_true|',
        'Significant digits':'Number of reliable digits in result',
        'Machine epsilon':  'Smallest ε: fl(1+ε) > 1, ε ≈ 2.2×10⁻¹⁶ (double precision)',
        'Floating point':   'fl(x) = x(1+δ) where |δ| ≤ ε (relative error bounded)',
        'Catastrophic cancellation': 'Subtracting nearly equal numbers → massive relative error',
        'Example fix':      'x² - 1 near x=1: use (x-1)(x+1) instead',
    }

def condition_number():
    """
    Condition number measures sensitivity to perturbations.
    """
    return {
        'Definition':       'κ(A) = ||A|| · ||A⁻¹|| = σ_max/σ_min',
        'Interpretation':   'κ significant digits lost ≈ log₁₀(κ)',
        'Well-conditioned':  'κ ≈ 1: small perturbation → small change',
        'Ill-conditioned':   'κ >> 1: small perturbation → large change',
        'Relative error':    '||Δx||/||x|| ≤ κ(A) · ||Δb||/||b||',
        'Examples': {
            'Identity':     'κ = 1 (perfect)',
            'Hilbert matrix':'κ grows exponentially with size (pathological)',
            'Vandermonde':  'Often ill-conditioned for large n'
        }
    }

import numpy as np

def floating_point_pitfalls():
    """
    Demonstrate common floating point issues.
    """
    # Catastrophic cancellation
    x = 1.0 + 1e-15
    bad = x**2 - 1.0  # loses precision
    good = (x-1)*(x+1)  # algebraically equivalent, numerically better

    # Order of operations matters
    s_forward = sum(1.0/k for k in range(1, 1001))   # small numbers added to large
    s_backward = sum(1.0/k for k in range(1000, 0, -1))  # more accurate

    return {
        'cancellation_bad':  bad,
        'cancellation_good': good,
        'sum_forward':       round(s_forward, 10),
        'sum_backward':      round(s_backward, 10),
        'rule':              'Add small numbers first for better accuracy'
    }

Root Finding

def bisection_method(f, a, b, tol=1e-10, max_iter=100):
    """
    Bisection method: guaranteed convergence, linear rate.
    Requires f(a) and f(b) have opposite signs.
    """
    if f(a) * f(b) > 0:
        raise ValueError("f(a) and f(b) must have opposite signs")

    iterations = []
    for i in range(max_iter):
        c = (a + b) / 2
        fc = f(c)
        iterations.append({'iter': i+1, 'a': a, 'b': b, 'c': c, 'f(c)': fc})

        if abs(fc) < tol or (b-a)/2 < tol:
            break
        if f(a) * fc < 0:
            b = c
        else:
            a = c

    return c, iterations

def newton_raphson(f, df, x0, tol=1e-10, max_iter=50):
    """
    Newton-Raphson: quadratic convergence near simple root.
    xₙ₊₁ = xₙ - f(xₙ)/f'(xₙ)
    """
    x = x0
    history = [x0]

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

        if abs(dfx) < 1e-15:
            raise ZeroDivisionError("Derivative too small")

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

        if abs(x_new - x) < tol:
            break
        x = x_new

    return x, history

def secant_method(f, x0, x1, tol=1e-10, max_iter=50):
    """
    Secant method: superlinear convergence (~1.618), no derivative needed.
    xₙ₊₁ = xₙ - f(xₙ)(xₙ-xₙ₋₁)/(f(xₙ)-f(xₙ₋₁))
    """
    for i in range(max_iter):
        f0, f1 = f(x0), f(x1)
        if abs(f1 - f0) < 1e-15:
            break
        x2 = x1 - f1*(x1-x0)/(f1-f0)
        if abs(x2 - x1) < tol:
            return x2
        x0, x1 = x1, x2
    return x1

def fixed_point_iteration(g, x0, tol=1e-10, max_iter=100):
    """
    Fixed point: x = g(x), converges if |g'(x*)| < 1.
    Linear convergence rate = |g'(x*)|.
    """
    x = x0
    for i in range(max_iter):
        x_new = g(x)
        if abs(x_new - x) < tol:
            return x_new, i+1
        x = x_new
    return x, max_iter

def convergence_rates():
    return {
        'Bisection':        'Linear, rate = 0.5 (1 bit per iteration)',
        'Newton-Raphson':   'Quadratic near simple root (doubles digits)',
        'Secant':           'Superlinear, order φ = (1+√5)/2 ≈ 1.618',
        'Fixed point':      'Linear, rate = |g'(x*)|',
        'Müller's method': 'Order ≈ 1.84 (uses quadratic interpolation)',
        'Regula Falsi':     'Linear but faster than bisection in practice'
    }

Interpolation

def lagrange_interpolation(x_data, y_data, x_eval):
    """
    Lagrange polynomial interpolation.
    P(x) = Σ yᵢ · Lᵢ(x)
    Lᵢ(x) = Π_{j≠i} (x-xⱼ)/(xᵢ-xⱼ)
    """
    n = len(x_data)
    result = 0.0

    for i in range(n):
        L = 1.0
        for j in range(n):
            if j != i:
                L *= (x_eval - x_data[j]) / (x_data[i] - x_data[j])
        result += y_data[i] * L

    return result

def newton_divided_differences(x_data, y_data):
    """
    Newton's divided difference table.
    More numerically stable than Lagrange for adding points.
    """
    n = len(x_data)
    F = np.array([[0.0]*n for _ in range(n)])

    for i in range(n):
        F[i][0] = y_data[i]

    for j in range(1, n):
        for i in range(n-j):
            F[i][j] = (F[i+1][j-1] - F[i][j-1]) / (x_data[i+j] - x_data[i])

    # Coefficients = F[0][0], F[0][1], ..., F[0][n-1]
    coeffs = [F[0][j] for j in range(n)]
    return coeffs, F

def cubic_spline(x_data, y_data):
    """
    Natural cubic spline interpolation.
    C² continuous, minimizes oscillation (Runge's phenomenon avoidance).
    """
    n = len(x_data) - 1
    h = [x_data[i+1] - x_data[i] for i in range(n)]

    # Set up tridiagonal system for second derivatives
    # Natural spline: S''(x₀) = S''(xₙ) = 0
    # ... (tridiagonal solve)
    return "Cubic spline coefficients (solve tridiagonal system)"

def interpolation_concepts():
    return {
        'Runge's phenomenon': {
            'issue':        'High-degree polynomial oscillates wildly near endpoints',
            'example':      '1/(1+25x²) on [-1,1] with equally spaced nodes',
            'solutions':    'Use Chebyshev nodes, splines, or low-degree piecewise'
        },
        'Chebyshev nodes': {
            'formula':      'xₖ = cos((2k+1)π/(2n+2)) for k=0,...,n',
            'property':     'Minimize Lebesgue constant, reduce Runge's phenomenon',
            'optimal':      'Best node distribution for polynomial interpolation'
        },
        'Error bound': {
            'Lagrange':     '|f(x)-Pₙ(x)| ≤ max|f^(n+1)|/(n+1)! · |ω(x)|',
            'ω(x)':         'Πᵢ(x-xᵢ) (node polynomial)',
            'control':      'Better nodes (Chebyshev) → smaller |ω(x)|'
        }
    }

Numerical Integration (Quadrature)

def numerical_integration_methods():
    return {
        'Trapezoidal rule': {
            'formula':      '∫ₐᵇf dx ≈ h/2 [f(x₀)+2f(x₁)+...+2f(xₙ₋₁)+f(xₙ)]',
            'error':        'O(h²) per step, O(h²) global',
            'composite':    'n subintervals of width h=(b-a)/n'
        },
        'Simpson's 1/3 rule': {
            'formula':      '∫ₐᵇf dx ≈ h/3 [f(x₀)+4f(x₁)+2f(x₂)+...+4f(xₙ₋₁)+f(xₙ)]',
            'error':        'O(h⁴) global (requires even n)',
            'note':         'Exact for polynomials of degree ≤ 3'
        },
        'Simpson's 3/8 rule': {
            'formula':      '3h/8 [f₀+3f₁+3f₂+2f₃+...+fₙ] (n divisible by 3)',
            'error':        'O(h⁴) global'
        },
        'Gaussian quadrature': {
            'idea':         'Choose both nodes AND weights optimally',
            'n-point rule': 'Exact for polynomials of degree ≤ 2n-1',
            'Gauss-Legendre':'Nodes = roots of Legendre polynomials on [-1,1]',
            'transform':    '∫ₐᵇf(x)dx = (b-a)/2 ∫₋₁¹ f((b-a)t/2+(a+b)/2) dt',
            'points_weights': {
                'n=2': 'x=±1/√3, w=1',
                'n=3': 'x=0,±√(3/5), w=8/9, 5/9'
            }
        },
        'Adaptive quadrature': {
            'idea':         'Subdivide intervals where error estimate is large',
            'stopping':     '|coarse - fine| < tolerance',
            'Python':       'scipy.integrate.quad implements this'
        },
        'Monte Carlo integration': {
            'formula':      '∫f dx ≈ (b-a)/N Σf(xᵢ) where xᵢ uniform random',
            'error':        'O(1/√N) regardless of dimension',
            'best_for':     'High-dimensional integrals (d >> 1)'
        }
    }

def composite_simpson(f, a, b, n):
    """
    Composite Simpson's rule: n must be even.
    """
    if n % 2 != 0:
        n += 1
    h = (b - a) / n
    x = np.linspace(a, b, n+1)
    y = np.array([f(xi) for xi in x])

    result = y[0] + y[-1]
    result += 4 * np.sum(y[1:-1:2])   # odd indices
    result += 2 * np.sum(y[2:-2:2])   # even indices (not endpoints)
    return result * h / 3

def romberg_integration(f, a, b, max_rows=10):
    """
    Romberg integration: Richardson extrapolation applied to trapezoidal rule.
    R[n,m] = (4ᵐR[n,m-1] - R[n-1,m-1]) / (4ᵐ - 1)
    """
    R = np.zeros((max_rows, max_rows))

    # First column: trapezoidal with 2^n subintervals
    for n in range(max_rows):
        N = 2**n
        h = (b-a)/N
        x = np.linspace(a, b, N+1)
        R[n,0] = h/2 * (f(a) + f(b) + 2*sum(f(xi) for xi in x[1:-1]))

    # Richardson extrapolation
    for m in range(1, max_rows):
        for n in range(m, max_rows):
            R[n,m] = (4**m * R[n,m-1] - R[n-1,m-1]) / (4**m - 1)

    return R

Numerical Differentiation

def finite_differences():
    return {
        'Forward difference':   '(f(x+h) - f(x))/h — O(h) error',
        'Backward difference':  '(f(x) - f(x-h))/h — O(h) error',
        'Central difference':   '(f(x+h) - f(x-h))/(2h) — O(h²) error',
        'Second derivative':    '(f(x+h) - 2f(x) + f(x-h))/h² — O(h²) error',
        'Fourth order central': '(-f(x+2h)+8f(x+h)-8f(x-h)+f(x-2h))/(12h) — O(h⁴)',
        'Optimal h':            '~√ε_machine for first derivative ≈ 10⁻⁸',
        'Trade-off':            'Smaller h: less truncation, more round-off'
    }

def richardson_extrapolation(f, x, h, order=1):
    """
    Richardson extrapolation to improve accuracy.
    D(h) = f'(x) + C·hᵖ + O(h^(p+1))
    D(h/2) = f'(x) + C·(h/2)ᵖ + ...
    Extrapolated: f'(x) ≈ (2ᵖD(h/2) - D(h))/(2ᵖ-1)
    """
    D_h  = (f(x+h) - f(x-h)) / (2*h)        # O(h²)
    D_h2 = (f(x+h/2) - f(x-h/2)) / h        # O(h²)
    # Extrapolate to O(h⁴)
    extrapolated = (4*D_h2 - D_h) / 3
    return {
        'D(h)':         D_h,
        'D(h/2)':       D_h2,
        'extrapolated': extrapolated
    }

Numerical Linear Algebra

def numerical_linear_algebra():
    return {
        'LU decomposition': {
            'factorization':'A = LU (L lower triangular, U upper triangular)',
            'partial_pivot':'PA = LU for numerical stability',
            'solve':        'Ax=b: Ly=Pb (forward), Ux=y (backward)',
            'cost':         'O(n³) factorization, O(n²) per solve',
            'use':          'Multiple right-hand sides with same A'
        },
        'Cholesky': {
            'for':          'Symmetric positive definite A = LLᵀ',
            'cost':         'Half the cost of LU',
            'stability':    'Always stable for SPD matrices'
        },
        'QR decomposition': {
            'factorization':'A = QR (Q orthogonal, R upper triangular)',
            'methods':      'Gram-Schmidt, Householder (more stable), Givens',
            'use':          'Least squares (more stable than normal equations)',
            'eigenvalues':  'QR algorithm for all eigenvalues'
        },
        'Iterative methods (large sparse)': {
            'Jacobi':       'xᵢ^(k+1) = (bᵢ - Σⱼ≠ᵢ aᵢⱼxⱼ^(k)) / aᵢᵢ',
            'Gauss-Seidel': 'Use updated values immediately',
            'SOR':          'Over-relaxation: ω∈(1,2) speeds Gauss-Seidel',
            'CG':           'Conjugate Gradient for SPD (optimal Krylov)',
            'GMRES':        'General non-symmetric systems',
            'convergence':  'Depends on spectral radius ρ(iteration matrix)'
        },
        'Preconditioning': {
            'idea':         'Solve M⁻¹Ax = M⁻¹b (better conditioned)',
            'incomplete_LU':'ILU: sparse approximate factorization',
            'diagonal':     'Jacobi preconditioner: M = diag(A)'
        }
    }

def power_iteration(A, num_iterations=100, tol=1e-10):
    """
    Power iteration: finds largest eigenvalue and eigenvector.
    v_{k+1} = Av_k / ||Av_k||
    """
    n = A.shape[0]
    v = np.random.rand(n)
    v = v / np.linalg.norm(v)

    eigenvalue = 0
    for i in range(num_iterations):
        Av = A @ v
        eigenvalue_new = np.dot(v, Av)
        v = Av / np.linalg.norm(Av)

        if abs(eigenvalue_new - eigenvalue) < tol:
            break
        eigenvalue = eigenvalue_new

    return eigenvalue, v

Numerical ODEs

def ode_methods():
    return {
        'Euler (explicit)': {
            'formula':      'y_{n+1} = y_n + h·f(t_n, y_n)',
            'order':        '1st order: error O(h)',
            'stability':    'Conditionally stable: h < 2/|λ| for stiff',
            'use':          'Simple demonstration, not production'
        },
        'Runge-Kutta 4': {
            'k1':           'h·f(tₙ, yₙ)',
            'k2':           'h·f(tₙ+h/2, yₙ+k1/2)',
            'k3':           'h·f(tₙ+h/2, yₙ+k2/2)',
            'k4':           'h·f(tₙ+h, yₙ+k3)',
            'formula':      'y_{n+1} = y_n + (k1+2k2+2k3+k4)/6',
            'order':        '4th order: error O(h⁴)',
            'use':          'Standard for non-stiff problems'
        },
        'Dormand-Prince (RK45)': {
            'idea':         'Embedded 4th and 5th order for error estimation',
            'adaptive':     'Adjust h based on ||y₅-y₄||',
            'implementation':'scipy.integrate.solve_ivp default'
        },
        'Stiff equations': {
            'problem':      'Solution varies on multiple timescales',
            'symptom':      'Explicit methods require tiny h',
            'solution':     'Implicit methods (backward Euler, trapezoidal, BDF)',
            'backward_euler':'y_{n+1} = y_n + h·f(t_{n+1}, y_{n+1}) (solve nonlinear)',
            'BDF':          'Backward Differentiation Formulas (scipy ode15s equivalent)'
        },
        'Multistep methods': {
            'Adams-Bashforth':'Explicit, uses previous steps',
            'Adams-Moulton': 'Implicit, higher accuracy',
            'BDF':           'Stiff solvers, variable order'
        }
    }

def rk4_solver(f, y0, t_span, h):
    """
    Classic RK4 ODE solver.
    """
    t0, tf = t_span
    t = np.arange(t0, tf+h, h)
    y = np.zeros((len(t), len(np.atleast_1d(y0))))
    y[0] = y0

    for i in range(len(t)-1):
        k1 = h * np.array(f(t[i], y[i]))
        k2 = h * np.array(f(t[i]+h/2, y[i]+k1/2))
        k3 = h * np.array(f(t[i]+h/2, y[i]+k2/2))
        k4 = h * np.array(f(t[i]+h, y[i]+k3))
        y[i+1] = y[i] + (k1 + 2*k2 + 2*k3 + k4) / 6

    return t, y

Numerical PDEs

def pde_methods():
    return {
        'Finite Difference Method (FDM)': {
            'idea':         'Replace derivatives with finite differences on grid',
            'heat_explicit':'u^(n+1)_i = u^n_i + r(u^n_{i+1}-2u^n_i+u^n_{i-1}), r=αΔt/Δx²',
            'stability':    'Explicit heat: r ≤ 1/2 (CFL condition)',
            'implicit':     'Always stable (Crank-Nicolson, backward Euler)',
            'wave':         'CFL: c·Δt/Δx ≤ 1 for stability'
        },
        'Finite Element Method (FEM)': {
            'idea':         'Variational formulation, piecewise polynomial basis',
            'weak_form':    'Multiply PDE by test function, integrate by parts',
            'stiffness':    'Assemble global stiffness matrix K and load vector f',
            'solve':        'Ku = f (linear system)',
            'advantage':    'Complex geometry, natural boundary conditions',
            'software':     'FEniCS, COMSOL, ANSYS, Abaqus'
        },
        'Finite Volume Method (FVM)': {
            'idea':         'Integrate PDE over control volumes (conservation form)',
            'advantage':    'Exactly conservative (important for fluids)',
            'software':     'OpenFOAM, SU2'
        },
        'Spectral methods': {
            'idea':         'Represent solution as global basis functions (Fourier, Chebyshev)',
            'advantage':    'Exponential convergence for smooth solutions',
            'use':          'Smooth periodic or simple geometry problems'
        }
    }

Optimization

def optimization_methods():
    return {
        'Unconstrained optimization': {
            'Gradient descent':     'x_{k+1} = x_k - α∇f(x_k)',
            'Newton's method':     'x_{k+1} = x_k - [∇²f(x_k)]⁻¹∇f(x_k)',
            'BFGS':                 'Quasi-Newton: approximate Hessian update',
            'L-BFGS':               'Limited memory BFGS (large scale)',
            'Conjugate gradient':   'For quadratic problems and extensions',
            'Line search':          'Choose α to satisfy Wolfe conditions'
        },
        'Constrained optimization': {
            'Lagrange multipliers': '∇f = λ∇g (equality constraints)',
            'KKT conditions':       'First-order necessary for inequality constraints',
            'Penalty methods':      'Add constraint violation to objective',
            'Interior point':       'Stay strictly feasible, log barrier',
            'Active set':           'Identify which constraints are active'
        },
        'Global optimization': {
            'Simulated annealing':  'Accept worse solutions with decreasing probability',
            'Genetic algorithms':   'Evolutionary: selection, crossover, mutation',
            'Particle swarm':       'Swarm intelligence',
            'Basin hopping':        'Random perturbation + local optimization'
        }
    }

Software & Libraries

def numerical_software():
    return {
        'Python': {
            'NumPy':        'Array operations, linear algebra',
            'SciPy':        'Root finding, integration, ODEs, optimization, linear algebra',
            'Matplotlib':   'Visualization',
            'SymPy':        'Symbolic mathematics'
        },
        'Key SciPy functions': {
            'scipy.optimize.brentq':    'Root finding (robust)',
            'scipy.optimize.fsolve':    'System of nonlinear equations',
            'scipy.integrate.quad':     'Adaptive quadrature',
            'scipy.integrate.solve_ivp':'ODE solver (RK45 default)',
            'scipy.linalg.solve':       'Linear system (uses LAPACK)',
            'scipy.linalg.lstsq':       'Least squares',
            'scipy.linalg.eig':         'Eigenvalues',
            'scipy.sparse.linalg.spsolve': 'Sparse linear system'
        },
        'MATLAB equivalents': {
            'ode45':        'Non-stiff ODE (Dormand-Prince RK45)',
            'ode15s':       'Stiff ODE (BDF/NDF)',
            'fzero':        'Root finding',
            'fminunc':      'Unconstrained optimization',
            'linprog':      'Linear programming',
            '':           'Linear system solve (backslash operator)'
        }
    }

Common Pitfalls

Related Skills

• calculus-expert: Mathematical foundations
• linear-algebra-expert: Theory behind numerical LA
• differential-equations-expert: ODE/PDE theory
• statistics-expert: Numerical optimization in statistics
• computational-chemistry-expert: Scientific computing applications
• machine-learning-expert: Gradient descent and optimization