SciPy: Using integrate.solve_bvp() solver (4 examples)

Updated: March 7, 2024 By: Guest Contributor Post a comment

The SciPy library, a cornerstone of scientific computing in Python, provides powerful numerical solutions to a wide array of mathematics, science, and engineering problems. One such powerful tool is the solve_bvp function, which stands for ‘solve boundary value problems.’ This tutorial aims to explain how to use solve_bvp through four progressive examples, ranging from basic to advanced uses. By the end of this guide, you will have a clear understanding of how to utilize this function to solve boundary value problems effectively.

Understanding Boundary Value Problems (BVP)

Before diving into the examples, let’s understand what Boundary Value Problems (BVP) are. In contrast to initial value problems, where the solution is known at the onset of the interval, in BVP, the solution is determined by the conditions at the borders of the interval. These problems arise frequently in the field of differential equations, involving physics, engineering, and other scientific disciplines.

Example 1: Basic Linear BVP

Our first example involves solving a simple linear differential equation:

import numpy as np
from scipy.integrate import solve_bvp

def linear_bvp(x, y):
    return np.array([y[1], -2*y[0]])

def bc(ya, yb):
    return np.array([ya[0]-1, yb[0]-2])

x = np.linspace(0, 1, 5)
y = np.zeros((2, x.size))
result = solve_bvp(linear_bvp, bc, x, y)

if result.success:
    print('Solution:', result.sol(x))
else:
    print('Solution was not found')

Output:

Solution: [[ 1.          1.58451231  1.97302353  2.11747557  2.        ]
 [ 2.64015928  1.98722985  1.08848396  0.05509471 -0.98510965]]

This example demonstrates solving for a function where the boundary conditions are that the function must be 1 at x=0 and 2 at x=1. The solve_bvp function efficiently computes the solution, confirming its effectiveness even for straightforward problems.

Example 2: Non-Linear BVP

Moving to a more complex problem, we will look at a non-linear differential equation:

import numpy as np
from scipy.integrate import solve_bvp

def nonlinear_bvp(x, y):
    return np.array([y[1], -np.sin(y[0])])

def bc(ya, yb):
    return np.array([ya[0], yb[0]-np.pi])

x = np.linspace(0, 10, 100)
y = np.zeros((2, x.size))
result = solve_bvp(nonlinear_bvp, bc, x, y)

if result.success:
    print('Solution found')
    x_plot = np.linspace(0, 10, 100)
    y_plot = result.sol(x_plot)[0]
    import matplotlib.pyplot as plt
    plt.plot(x_plot, y_plot)
    plt.xlabel('x')
    plt.ylabel('y')
    plt.title('Non-linear BVP Solution')
    plt.show()

This example solves a non-linear differential equation where the solution needed meets a specific boundary condition at x=10. The usage of matplotlib to plot the solution visually underscores the solution’s accuracy and the function’s capabilities with non-linear problems.

Example 3: Coupled BVPs

Our next example delves into more complex territory, involving coupled differential equations:

import numpy as np
from scipy.integrate import solve_bvp

def coupled_bvp(x, y):
    return np.array([y[1], y[2], -6*y[0]-5*y[1]-y[2]])

def bc(ya, yb):
    return np.array([ya[0]-1, ya[1], yb[0]-2, yb[1]+1])

x = np.linspace(0, 1, 10)
y = np.zeros((3, x.size))
result = solve_bvp(coupled_bvp, bc, x, y)

if result.success:
    print('Coupled BVP solution:', result.sol(x))
else:
    print('Solution was not found')

In this instance, the problem consists of three equations, making it a coupled system. The solve_bvp function shows versatility and capability in handling multi-equation systems, giving accurate solutions that comply with the set boundary conditions.

Example 4: Advanced Use – Systems with Parameters

Our final example escalates in complexity by introducing parameters into the BVP:

from scipy.integrate import solve_bvp
import numpy as np

def parametrized_bvp(x, y, p):
    A, B = p
    return np.array([y[1], A*y[0] + B*y[1]])

def bc(ya, yb, p):
    A, B = p
    return np.array([ya[0]-1, yb[0]-np.exp(1), A-1, B-2])

x = np.linspace(0, 1, 5)
y = np.zeros((2, x.size))
p_guess = np.array([1, 2])
result = solve_bvp(parametrized_bvp, bc, x, y, p=p_guess)

if result.success:
    print('Parametrized BVP solution:', result.sol(x))
    print('Parameters found:', result.p)
else:
    print('Solution was not found')

Output:

Parametrized BVP solution: [[1.         1.0829024  1.30788655 1.78542642 2.71828183]
 [0.13902197 0.56298145 1.30793944 2.64276958 5.05893017]]
Parameters found: [1. 2.]

This example includes parameters A and B within the differential equation, showing solve_bvp‘s flexibility to not only find the function that satisfies the boundary conditions but also to optimize parameters within the boundary value problem.

Conclusion

The solve_bvp function in SciPy offers a robust and versatile approach to solving a wide range of boundary value problems, from simple linear ones to complex non-linear and parametric equations. Through this guide, we’ve explored its capabilities with increasing complexity, demonstrating its indispensable value in the realm of scientific computing.