How to Implement Custom Numerical Solvers with NumPy

Introduction

In computational mathematics, the design of custom numerical solvers is often essential for tackling a wide array of problems with specific features or requirements not covered by general-purpose algorithms. Python, equipped with the powerful NumPy library, is a fantastic tool for scientists and engineers looking to develop their numerical methods. This tutorial will guide you through implementing custom numerical solvers using NumPy, starting with the basics and progressing to more advanced examples.

What are Numerical Solvers?

Numerical solvers are algorithms that provide approximate solutions to mathematical problems. These problems can range from solving simple equations to complex differential equations encountered in various fields such as physics, engineering, and finance. The beauty of using NumPy lies in its array manipulation capabilities, which are conducive to numerical computations.

Getting Started

import numpy as np

# Make sure NumPy is installed
print(np.__version__)

this should output the version of NumPy which confirms that the installation is successful:

1.22.3

Basic Example: Bisection Method

The bisection method is a simple root-finding algorithm that repeatedly narrows down an interval containing a zero of a function. Without further ado, let’s implement this algorithm to find the square root of 2.

def bisection_method(f, a, b, tol):
    assert f(a) * f(b) < 0, 'Function values at the boundaries must have opposite signs'
    while (b - a) / 2.0 > tol:
        midpoint = (a + b) / 2.0
        if f(midpoint) == 0:
            return midpoint
        elif f(midpoint) * f(a) < 0:
            b = midpoint
        else:
            a = midpoint
    return (a + b) / 2.0

root = bisection_method(lambda x: x**2 - 2, 1, 2, 1e-7)
print(f'The root is: {root}') 

The bisection_method function should return an approximation of the square root of 2 within the specified tolerance:

The root is: 1.4142135623842478

Intermediate Example: Euler’s Method for Ordinary Differential Equations

Euler’s Method is a straightforward numerical procedure for solving ordinary differential equations (ODEs) with a given initial value. It is the simplest explicit method for numerical integration of ODEs and serves as the basis for more sophisticated methods such as the Runge–Kutta methods. Let’s consider solving the ODE dy/dx = x + y with initial condition y(0) = 1.

def euler_method(dy_dx, x0, y0, x1, h):
    x = np.arange(x0, x1+h, h)
    y = np.zeros(len(x))
    y[0] = y0
    for i in range(1, len(x)):
        y[i] = y[i-1] + h * dy_dx(x[i-1], y[i-1])

    return x, y

x, y = euler_method(lambda x, y: x + y, 0, 1, 5, 0.1)
for xi, yi in zip(x, y):
    print(f'({xi}, {yi})') 

Upon running the euler_method function, you will observe the discrete points approximating the solution to the ODE:

Advanced Example: Implementing a Runge–Kutta Method

Runge–Kutta methods are among the most popular methods used to solve ODEs and are known for their accuracy and efficiency over the Euler’s method in many situations. The following is an implementation of the classical fourth-order Runge–Kutta method (often referred to as the RK4 method).

def rk4(dy_dx, x0, y0, x1, n):
    x = np.linspace(x0, x1, n + 1)
    y = np.zeros(n + 1)
    y[0] = y0
    h = (x1 - x0) / n
    for i in range(n):
        k1 = dy_dx(x[i], y[i])
        k2 = dy_dx(x[i] + h/2, y[i] + h/2 * k1)
        k3 = dy_dx(x[i] + h/2, y[i] + h/2 * k2)
        k4 = dy_dx(x[i] + h, y[i] + h * k3)
        y[i + 1] = y[i] + h * (k1 + 2*k2 + 2*k3 + k4) / 6

    return x, y

x, y = rk4(lambda x, y: x + y, 0, 1, 5, 100)
for xi, yi in zip(x, y):
    print(f'({xi}, {yi})') 

Conclusion

In this tutorial, we’ve explored how to leverage NumPy to implement custom numerical solvers ranging from a simple bisection method, through the Euler method, and up to the popular RK4 method for ODEs. Embracing these solvers enables the solving of numerous mathematical problems effectively. It is clear that with practice and understanding the underlying mathematical concepts, one can tailor algorithms using Python and NumPy that are perfectly suited to specific problem domains.