Introduction
Numerical methods are computational techniques used to solve mathematical problems that cannot be solved analytically. As a developer, you’ll encounter situations where you need to find roots of equations, optimize functions, integrate data, or solve systems of linear equations. This guide covers essential numerical methods with practical implementations in Python.
Root Finding Methods
Bisection Method
The bisection method finds roots by repeatedly dividing an interval in half:
def bisection(f, a, b, tol=1e-10, max_iter=100):
"""
Find root of function f in interval [a, b] using bisection method.
Args:
f: Continuous function with f(a) and f(b) having opposite signs
a, b: Interval endpoints
tol: Tolerance for convergence
max_iter: Maximum iterations
Returns:
Root approximation and number of iterations
"""
if f(a) * f(b) >= 0:
raise ValueError("f(a) and f(b) must have opposite signs")
iteration = 0
while (b - a) > tol and iteration < max_iter:
c = (a + b) / 2
fc = f(c)
if abs(fc) < tol:
return c, iteration
if f(a) * fc < 0:
b = c
else:
a = c
iteration += 1
return (a + b) / 2, iteration
def test_bisection():
# Find root of f(x) = x^2 - 2 (sqrt(2) โ 1.414)
f = lambda x: x**2 - 2
root, iterations = bisection(f, 0, 2)
print(f"Root: {root:.10f}, Iterations: {iterations}")
print(f"Verification: f(root) = {f(root):.15f}")
# Find root of f(x) = cos(x) - x
f = lambda x: math.cos(x) - x
root, iterations = bisection(f, 0, 1)
print(f"Root: {root:.10f}, Iterations: {iterations}")
Newton’s Method
Newton’s method uses the derivative for faster convergence:
def newton(f, df, x0, tol=1e-10, max_iter=100):
"""
Find root using Newton's method.
Args:
f: Function to find root for
df: Derivative of f
x0: Initial guess
tol: Tolerance
max_iter: Maximum iterations
Returns:
Root approximation and number of iterations
"""
x = x0
for iteration in range(max_iter):
fx = f(x)
if abs(fx) < tol:
return dfx = x, iteration
df(x)
if abs(dfx) < 1e-15:
raise ValueError("Derivative too small")
x = x - fx / dfx
return x, max_iter
def newton_raphson(f, x0, tol=1e-10, max_iter=100):
"""
Newton-Raphson method with numerical derivative.
"""
x = x0
h = 1e-10
for iteration in range(max_iter):
fx = f(x)
if abs(fx) < tol:
return x, iteration
# Numerical derivative
dfx = (f(x + h) - f(x - h)) / (2 * h)
if abs(dfx) < 1e-15:
raise ValueError("Derivative too small")
x = x - fx / dfx
return x, max_iter
def test_newton():
# f(x) = x^2 - 2, df(x) = 2x
f = lambda x: x**2 - 2
df = lambda x: 2*x
root, iterations = newton(f, df, 1)
print(f"Newton - Root: {root:.10f}, Iterations: {iterations}")
root, iterations = newton_raphson(f, 1)
print(f"Newton-Raphson - Root: {root:.10f}, Iterations: {iterations}")
Secant Method
When derivatives are unavailable, the secant method uses secant lines:
def secant(f, x0, x1, tol=1e-10, max_iter=100):
"""
Find root using secant method.
"""
x_prev = x0
x_curr = x1
for iteration in range(max_iter):
f_prev = f(x_prev)
f_curr = f(x_curr)
if abs(f_curr) < tol:
return x_curr, iteration
if abs(f_curr - f_prev) < 1e-15:
raise ValueError("Division by zero in secant method")
# Secant formula
x_next = x_curr - f_curr * (x_curr - x_prev) / (f_curr - f_prev)
x_prev = x_curr
x_curr = x_next
return x_curr, max_iter
Numerical Integration
Trapezoidal Rule
def trapezoidal(f, a, b, n=1000):
"""
Approximate integral using trapezoidal rule.
Args:
f: Function to integrate
a, b: Integration bounds
n: Number of subintervals
Returns:
Approximate integral value
"""
h = (b - a) / n
result = 0.5 * (f(a) + f(b))
for i in range(1, n):
result += f(a + i * h)
return result * h
def simpson(f, a, b, n=1000):
"""
Approximate integral using Simpson's rule.
Requires n to be even.
"""
if n % 2 == 1:
n += 1
h = (b - a) / n
result = f(a) + f(b)
for i in range(1, n):
x = a + i * h
coefficient = 4 if i % 2 == 1 else 2
result += coefficient * f(x)
return result * h / 3
def adaptive_integration(f, a, b, tol=1e-10, max_depth=50):
"""
Adaptive Simpson's rule for automatic error control.
"""
def simpson_simple(f, a, b):
c = (a + b) / 3
h = b - a
return h / 6 * (f(a) + 4*f(c) + f(b))
def recursive_asr(f, a, b, tol, whole, depth):
c = (a + b) / 2
left = simpson_simple(f, a, c)
right = simpson_simple(f, c, b)
if depth >= max_depth or abs(left + right - whole) <= 15 * tol:
return left + right + (left + right - whole) / 15
return (recursive_asr(f, a, c, tol/2, left, depth+1) +
recursive_asr(f, c, b, tol/2, right, depth+1))
whole = simpson_simple(f, a, b)
return recursive_asr(f, a, b, tol, whole, 0)
def test_integration():
import math
# Integrate sin(x) from 0 to pi
f = lambda x: math.sin(x)
exact = 2.0 # โซโ^ฯ sin(x) dx = 2
result_trap = trapezoidal(f, 0, math.pi, 1000)
result_simp = simpson(f, 0, math.pi, 1000)
result_adap = adaptive_integration(f, 0, math.pi)
print(f"Exact: {exact:.10f}")
print(f"Trapezoidal: {result_trap:.10f} (error: {abs(exact - result_trap):.2e})")
print(f"Simpson: {result_simp:.10f} (error: {abs(exact - result_simp):.2e})")
print(f"Adaptive: {result_adap:.10f} (error: {abs(exact - result_adap):.2e})")
Gaussian Quadrature
import numpy as np
def gaussian_quadrature(f, n):
"""
Gauss-Legendre quadrature for accurate integration.
Args:
f: Function to integrate
n: Number of quadrature points
Returns:
Approximate integral over [-1, 1]
"""
# Get nodes and weights
nodes, weights = np.polynomial.legendre.leggauss(n)
# Evaluate
result = sum(w * f(x) for w, x in zip(weights, nodes))
return result
def gaussian_quadrature_general(f, a, b, n):
"""
Gauss-Legendre quadrature over arbitrary interval [a, b].
"""
nodes, weights = np.polynomial.legendre.leggauss(n)
# Transform to [a, b]
transformed_nodes = 0.5 * (b - a) * nodes + 0.5 * (b + a)
result = 0.5 * (b - a) * sum(w * f(x) for w, x in zip(weights, transformed_nodes))
return result
Numerical Differentiation
def central_difference(f, x, h=1e-10):
"""
Central difference approximation for derivative.
"""
return (f(x + h) - f(x - h)) / (2 * h)
def second_derivative(f, x, h=1e-10):
"""
Second derivative using central differences.
"""
return (f(x + h) - 2*f(x) + f(x - h)) / (h**2)
def richardson_extrapolation(f, x, h=0.1):
"""
Richardson extrapolation for more accurate derivative.
"""
D1 = central_difference(f, x, h)
D2 = central_difference(f, x, h/2)
# Extrapolated value
return (4*D2 - D1) / 3
def gradient(f, x, h=1e-10):
"""
Compute gradient of multivariate function.
"""
n = len(x)
grad = np.zeros(n)
for i in range(n):
x_plus = x.copy()
x_minus = x.copy()
x_plus[i] += h
x_minus[i] -= h
grad[i] = (f(x_plus) - f(x_minus)) / (2*h)
return grad
Solving Linear Systems
Gaussian Elimination
import numpy as np
def gaussian_elimination(A, b):
"""
Solve Ax = b using Gaussian elimination with partial pivoting.
Args:
A: Coefficient matrix (n x n)
b: Right-hand side vector (n)
Returns:
Solution vector x
"""
n = len(b)
# Augment matrix
M = np.column_stack([A, b]).astype(float)
# Forward elimination
for i in range(n):
# Partial pivoting
max_row = i + np.argmax(abs(M[i:, i]))
M[[i, max_row]] = M[[max_row, i]]
# Eliminate column
for j in range(i + 1, n):
factor = M[j, i] / M[i, i]
M[j, i:] -= factor * M[i, i:]
# Back substitution
x = np.zeros(n)
for i in range(n - 1, -1, -1):
x[i] = (M[i, -1] - np.dot(M[i, :-1], x)) / M[i, i]
return x
def lu_decomposition(A):
"""
LU decomposition with pivoting: PA = LU
"""
n = len(A)
U = A.astype(float).copy()
L = np.eye(n)
P = np.eye(n)
for k in range(n - 1):
# Find pivot
max_idx = k + np.argmax(abs(U[k:, k]))
# Swap rows
U[[k, max_idx]] = U[[max_idx, k]]
P[[k, max_idx]] = P[[max_idx, k]]
if k > 0:
L[[k, max_idx], :k] = L[[max_idx, k], :k]
# Eliminate
for i in range(k + 1, n):
factor = U[i, k] / U[k, k]
L[i, k] = factor
U[i, k:] -= factor * U[k, k:]
return P, L, U
def solve_lu(P, L, U, b):
"""Solve Ax = b using LU decomposition."""
# Pb = LUx
Pb = np.dot(P, b)
# Ly = Pb
y = np.linalg.solve(L, Pb)
# Ux = y
x = np.linalg.solve(U, y)
return x
Iterative Methods
def jacobi(A, b, x0=None, tol=1e-10, max_iter=1000):
"""
Jacobi iteration for solving linear systems.
"""
n = len(b)
if x0 is None:
x = np.zeros(n)
else:
x = x0.copy()
D = np.diag(np.diag(A))
R = A - D
for iteration in range(max_iter):
x_new = np.dot(np.linalg.inv(D), b - np.dot(R, x))
if np.linalg.norm(x_new - x) < tol:
return x_new, iteration
x = x_new
return x, max_iter
def gauss_seidel(A, b, x0=None, tol=1e-10, max_iter=1000):
"""
Gauss-Seidel iteration.
"""
n = len(b)
if x0 is None:
x = np.zeros(n)
else:
x = x0.copy()
for iteration in range(max_iter):
x_old = x.copy()
for i in range(n):
x[i] = (b[i] - np.dot(A[i, :i], x[:i]) - np.dot(A[i, i+1:], x_old[i+1:])) / A[i, i]
if np.linalg.norm(x - x_old) < tol:
return x, iteration
return x, max_iter
def conjugate_gradient(A, b, x0=None, tol=1e-10, max_iter=1000):
"""
Conjugate gradient method for symmetric positive-definite matrices.
"""
n = len(b)
if x0 is None:
x = np.zeros(n)
else:
x = x0.copy()
r = b - np.dot(A, x)
p = r.copy()
rsold = np.dot(r, r)
for iteration in range(max_iter):
Ap = np.dot(A, p)
alpha = rsold / np.dot(p, Ap)
x = x + alpha * p
r = r - alpha * Ap
rsnew = np.dot(r, r)
if np.sqrt(rsnew) < tol:
return x, iteration
p = r + (rsnew / rsold) * p
rsold = rsnew
return x, max_iter
Optimization Methods
Gradient Descent
def gradient_descent(f, gradient, x0, learning_rate=0.01, tol=1e-10, max_iter=10000):
"""
Basic gradient descent optimization.
"""
x = x0.copy()
for iteration in range(max_iter):
grad = gradient(x)
x = x - learning_rate * grad
if np.linalg.norm(grad) < tol:
return x, iteration
return x, max_iter
def gradient_descent_momentum(f, gradient, x0, learning_rate=0.01,
momentum=0.9, tol=1e-10, max_iter=10000):
"""
Gradient descent with momentum.
"""
x = x0.copy()
v = np.zeros_like(x)
for iteration in range(max_iter):
grad = gradient(x)
v = momentum * v - learning_rate * grad
x = x + v
if np.linalg.norm(grad) < tol:
return x, iteration
return x, max_iter
def adam(f, gradient, x0, learning_rate=0.001, beta1=0.9,
beta2=0.999, eps=1e-8, tol=1e-10, max_iter=10000):
"""
Adam optimizer.
"""
x = x0.copy()
m = np.zeros_like(x)
v = np.zeros_like(x)
for iteration in range(max_iter):
grad = gradient(x)
# Update biased first moment estimate
m = beta1 * m + (1 - beta1) * grad
# Update biased second raw moment estimate
v = beta2 * v + (1 - beta2) * (grad ** 2)
# Compute bias-corrected estimates
m_hat = m / (1 - beta1 ** (iteration + 1))
v_hat = v / (1 - beta2 ** (iteration + 1))
# Update
x = x - learning_rate * m_hat / (np.sqrt(v_hat) + eps)
if np.linalg.norm(grad) < tol:
return x, iteration
return x, max_iter
def test_optimization():
# Rosenbrock function: f(x, y) = (a-x)^2 + b(y-x^2)^2
# Minimum at (a, a^2)
a, b = 1, 100
rosenbrock = lambda x: (a - x[0])**2 + b * (x[1] - x[0]**2)**2
rosenbrock_grad = lambda x: np.array([
-2*(a - x[0]) - 4*b*x[0]*(x[1] - x[0]**2),
2*b*(x[1] - x[0]**2)
])
x0 = np.array([0.0, 0.0])
result, iterations = gradient_descent(rosenbrock, rosenbrock_grad, x0)
print(f"Gradient Descent: {result}, iterations: {iterations}")
result, iterations = gradient_descent_momentum(rosenbrock, rosenbrock_grad, x0)
print(f"Momentum: {result}, iterations: {iterations}")
result, iterations = adam(rosenbrock, rosenbrock_grad, x0)
print(f"Adam: {result}, iterations: {iterations}")
Eigenvalue Computation
Power Iteration
def power_iteration(A, num_iterations=1000, tol=1e-10):
"""
Find dominant eigenvalue using power iteration.
"""
n = A.shape[0]
v = np.random.rand(n)
v = v / np.linalg.norm(v)
for iteration in range(num_iterations):
Av = np.dot(A, v)
eigenvalue = np.dot(v, Av)
v_new = Av / np.linalg.norm(Av)
if np.linalg.norm(v_new - v) < tol:
return eigenvalue, v_new, iteration
v = v_new
return eigenvalue, v, num_iterations
def inverse_power_iteration(A, shift=0, num_iterations=1000, tol=1e-10):
"""
Find eigenvalue closest to shift using inverse power iteration.
"""
n = A.shape[0]
A_shifted = A - shift * np.eye(n)
v = np.random.rand(n)
v = v / np.linalg.norm(v)
for iteration in range(num_iterations):
v_new = np.linalg.solve(A_shifted, v)
eigenvalue = np.dot(v, v_new)
v_new = v_new / np.linalg.norm(v_new)
if np.linalg.norm(v_new - v) < tol:
return eigenvalue + shift, v_new, iteration
v = v_new
return eigenvalue + shift, v_new, num_iterations
def qr_algorithm(A, num_iterations=100, tol=1e-10):
"""
QR algorithm for computing all eigenvalues.
"""
A_k = A.astype(float).copy()
n = A.shape[0]
for iteration in range(num_iterations):
# QR decomposition
Q, R = np.linalg.qr(A_k)
# Update
A_k = np.dot(R, Q)
# Check convergence
off_diagonal = np.sum(np.abs(np.tril(A_k, -1)))
if off_diagonal < tol:
break
eigenvalues = np.diag(A_k)
return eigenvalues
Conclusion
Numerical methods are essential tools for developers working on scientific computing, data analysis, optimization, and many other domains. The implementations provided in this guide give you starting points for building robust numerical solutions. For production use, consider using established libraries like NumPy, SciPy, or specialized numerical computing packages that provide optimized, well-tested implementations.
Key takeaways:
- Choose appropriate methods based on problem characteristics
- Understand convergence properties and limitations
- Implement proper error handling and tolerance checking
- Consider using established libraries for production code
Comments