Skip to Content

SOLVE_BVP

Overview

The SOLVE_BVP function numerically solves a boundary value problem for a system of second-order ODEs with two-point boundary conditions. It is a wrapper around the scipy.integrate.solve_bvp method, simplified for use in Excel.

The general form of a two-point boundary value problem is:

y′′=f(x,y)y(a)=yay(b)=yby'' = f(x, y) y(a) = y_a y(b) = y_b

This function simplifies the problem by assuming a linear system of the form y′′=A⋅yy'' = A \cdot y, where AA is a coefficient matrix. This is equivalent to solving a first-order system by augmenting with derivatives.

For a system of nn second-order ODEs, the function internally expands it to a system of 2n2n first-order ODEs:

u1′=un+1u2′=un+2⋮un+1′=A1,1u1+A1,2u2+⋯+A1,nun⋮u2n′=An,1u1+An,2u2+⋯+An,nun\begin{align*} u_1' &= u_{n+1} \\ u_2' &= u_{n+2} \\ &\vdots \\ u_{n+1}' &= A_{1,1} u_1 + A_{1,2} u_2 + \cdots + A_{1,n} u_n \\ &\vdots \\ u_{2n}' &= A_{n,1} u_1 + A_{n,2} u_2 + \cdots + A_{n,n} u_n \end{align*}

where u1,…,unu_1, \ldots, u_n represent the original functions and un+1,…,u2nu_{n+1}, \ldots, u_{2n} represent their derivatives.

For more details on the underlying SciPy function, refer to the official SciPy documentation .

This example function is provided as-is without any representation of accuracy.

Usage

To use the function in Excel:

=SOLVE_BVP(y_double_prime_coeffs, bc_ya, bc_yb, x_mesh, y_init)
  • y_double_prime_coeffs (2D list, required): Square matrix A in y'' = A @ y with numeric entries.
  • bc_ya (2D list, required): Column vector of boundary condition values for y(a); must have the same number of rows as y_double_prime_coeffs.
  • bc_yb (2D list, required): Column vector of boundary condition values for y(b); must have the same number of rows as y_double_prime_coeffs.
  • x_mesh (2D list, required): Column vector of initial mesh points; must be strictly increasing and contain at least 2 values.
  • y_init (2D list, required): Initial guess for the solution; must have 2*n rows (where n is the dimension of y_double_prime_coeffs) and columns matching the size of x_mesh.

The function returns a 2D list (matrix) where the first column contains the final mesh x-values and the remaining columns contain the corresponding values of the first n solution components. Returns an error message (string) if inputs are invalid.

Examples

Example 1: Simple Harmonic Oscillator

This example solves the simple harmonic oscillator ODE y′′=−yy'' = -y with boundary conditions y(0)=1y(0) = 1 and y(π/2)=0y(\pi/2) = 0. The expected solution is y=cos⁡(x)y = \cos(x).

Inputs:

y_double_prime_coeffsbc_yabc_ybx_meshy_init
-1.01.00.00.01.00.0
0.78540.7071-0.7071
1.57080.0-1.0

Excel formula:

=SOLVE_BVP({{-1}}, {1}, {0}, {0;0.7854;1.5708}, {1,0.7071,0;0,-0.7071,-1})

Expected output:

xy
0.0001.000
0.3930.924
0.7850.707
1.1780.383
1.5710.000

Example 2: Scaled Harmonic Oscillator

This example solves y′′=−4yy'' = -4y with boundary conditions y(0)=0y(0) = 0 and y(π/4)=1y(\pi/4) = 1. The expected solution is y=sin⁡(2x)y = \sin(2x).

Inputs:

y_double_prime_coeffsbc_yabc_ybx_meshy_init
-4.00.01.00.00.02.0
0.39270.70711.4142
0.78541.00.0

Excel formula:

=SOLVE_BVP({{-4}}, {0}, {1}, {0;0.3927;0.7854}, {0,0.7071,1;2,1.4142,0})

Expected output:

xy
0.0000.000
0.1960.383
0.3930.707
0.5890.924
0.7851.000

Example 3: Exponential Growth

This example solves y′′=0.5yy'' = 0.5y (which exhibits exponential growth rather than oscillation) with boundary conditions y(0)=1y(0) = 1 and y(1)=2y(1) = 2.

Inputs:

y_double_prime_coeffsbc_yabc_ybx_meshy_init
0.51.02.00.01.00.5
0.33331.25710.5
0.66671.58430.5
1.02.00.5

Excel formula:

=SOLVE_BVP({{0.5}}, {1}, {2}, {0;0.3333;0.6667;1}, {1,1.2571,1.5843,2;0.5,0.5,0.5,0.5})

Expected output:

xy
0.0001.000
0.3331.257
0.6671.584
1.0002.000

Example 4: Two-Component System

This example solves a 2×2 system where y1′′=−y1y_1'' = -y_1 and y2′′=−y2y_2'' = -y_2 (two decoupled harmonic oscillators) with mixed boundary conditions: y1(0)=1y_1(0) = 1, y1(1)=0y_1(1) = 0, y2(0)=0y_2(0) = 0, y2(1)=1y_2(1) = 1.

Inputs:

y_double_prime_coeffsbc_yabc_ybx_meshy_init
-1.00.01.00.00.01.00.73490.38880.0
0.0-1.00.01.00.3333-1.0-1.0-1.0-1.0
0.66670.00.38880.73491.0
1.01.01.01.01.0

Excel formula:

=SOLVE_BVP({{-1,0;0,-1}}, {1;0}, {0;1}, {0;0.3333;0.6667;1}, {1,0.7349,0.3888,0;-1,-1,-1,-1;0,0.3888,0.7349,1;1,1,1,1})

Expected output:

xy1y2
0.0001.0000.000
0.3330.7350.389
0.6670.3890.735
1.0000.0001.000

Python Code

 
import numpy as np
from scipy.integrate import solve_bvp as scipy_solve_bvp
 
 
def solve_bvp(
    y_double_prime_coeffs,
    bc_ya,
    bc_yb,
    x_mesh,
    y_init,
):
    """
    Solve a boundary value problem for a second-order system of ODEs.
 
    Wraps scipy.integrate.solve_bvp to solve systems of the form:
    y'' = A @ y with boundary conditions y(a) = bc_ya and y(b) = bc_yb.
    For more details, see https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.solve_bvp.html
 
    Args:
        y_double_prime_coeffs (2D list): Square matrix A in y'' = A @ y.
        bc_ya (2D list): Column vector of boundary condition values at x[0].
        bc_yb (2D list): Column vector of boundary condition values at x[-1].
        x_mesh (2D list): Column vector of initial mesh points (strictly increasing).
        y_init (2D list): Initial guess for y values; rows are solution components, columns are mesh points.
 
    Returns:
        2D list where first column is x-mesh, remaining columns are y values,
        or an error message (str) when inputs are invalid.
 
    This example function is provided as-is without any representation of accuracy.
    """
 
    def _ensure_matrix(data, name):
        if isinstance(data, (int, float)):
            return [[float(data)]]
        if not isinstance(data, list) or len(data) == 0:
            return f"Invalid input: {name} must be a 2D list with at least one row."
        matrix = []
        row_length = None
        for row in data:
            if not isinstance(row, list) or len(row) == 0:
                return f"Invalid input: each row in {name} must be a non-empty list."
            try:
                float_row = [float(value) for value in row]
            except Exception:
                return f"Invalid input: {name} must contain numeric values."
            if row_length is None:
                row_length = len(float_row)
            elif len(float_row) != row_length:
                return f"Invalid input: all rows in {name} must have the same length."
            matrix.append(float_row)
        return matrix
 
    def _ensure_column(data, name):
        if isinstance(data, (int, float)):
            return [float(data)]
        if not isinstance(data, list) or len(data) == 0:
            return f"Invalid input: {name} must be a 2D list with at least one row."
        column = []
        for row in data:
            if not isinstance(row, list) or len(row) != 1:
                return f"Invalid input: each row in {name} must contain exactly one value."
            try:
                column.append(float(row[0]))
            except Exception:
                return f"Invalid input: {name} must contain numeric values."
        return column
 
    matrix_result = _ensure_matrix(y_double_prime_coeffs, "y_double_prime_coeffs")
    if isinstance(matrix_result, str):
        return matrix_result
    coeff_matrix = np.array(matrix_result, dtype=float)
    if coeff_matrix.ndim != 2:
        return "Invalid input: y_double_prime_coeffs must be a 2D list (matrix)."
    if coeff_matrix.shape[0] != coeff_matrix.shape[1]:
        return "Invalid input: y_double_prime_coeffs must form a square matrix."
 
    bc_ya_result = _ensure_column(bc_ya, "bc_ya")
    if isinstance(bc_ya_result, str):
        return bc_ya_result
    bc_a = np.array(bc_ya_result, dtype=float)
 
    bc_yb_result = _ensure_column(bc_yb, "bc_yb")
    if isinstance(bc_yb_result, str):
        return bc_yb_result
    bc_b = np.array(bc_yb_result, dtype=float)
 
    if bc_a.size != coeff_matrix.shape[0] or bc_b.size != coeff_matrix.shape[0]:
        return "Invalid input: boundary condition sizes must match y_double_prime_coeffs dimension."
 
    x_mesh_result = _ensure_column(x_mesh, "x_mesh")
    if isinstance(x_mesh_result, str):
        return x_mesh_result
    if len(x_mesh_result) < 2:
        return "Invalid input: x_mesh must contain at least 2 points."
    x = np.array(x_mesh_result, dtype=float)
    if not np.all(np.diff(x) > 0):
        return "Invalid input: x_mesh must be strictly increasing."
 
    matrix_result = _ensure_matrix(y_init, "y_init")
    if isinstance(matrix_result, str):
        return matrix_result
    y = np.array(matrix_result, dtype=float)
 
    n = coeff_matrix.shape[0]
    expected_rows = 2 * n
    if y.shape[0] != expected_rows:
        return f"Invalid input: y_init must have {expected_rows} rows (2 * dimension of y_double_prime_coeffs)."
    if y.shape[1] != x.size:
        return "Invalid input: y_init column count must match x_mesh size."
 
    def rhs(x_val, y_val):
        n = y_val.shape[0] // 2
        result = np.zeros_like(y_val)
        result[:n, :] = y_val[n:, :]
        result[n:, :] = coeff_matrix @ y_val[:n, :]
        return result
 
    def bc_residuals(ya, yb):
        n_bc = bc_a.size
        residuals = np.concatenate([ya[:n_bc] - bc_a, yb[:n_bc] - bc_b])
        return residuals
 
    n = coeff_matrix.shape[0]
    try:
        result = scipy_solve_bvp(rhs, bc_residuals, x, y, tol=0.001, max_nodes=1000)
    except Exception as exc:
        return f"scipy.integrate.solve_bvp error: {exc}"
 
    if not result.success:
        return f"scipy.integrate.solve_bvp error: {result.message}"
 
    x_out = result.x
    y_out = result.y[:n, :]
 
    if not np.all(np.isfinite(x_out)) or not np.all(np.isfinite(y_out)):
        return "scipy.integrate.solve_bvp error: result is not finite."
 
    combined = np.column_stack((x_out, y_out.T))
    return combined.tolist()

Example Workbook

Link to Workbook 

Last updated on
VAN_DER_POLSOLVE_IVP