BVP Solver Examples

Example Code

Representative script:

"""IMSL BVPFD example: BVP y'' + y = 0, y(0) = 0, y(pi/2) = 1.

Demonstrates use of bvpfd (finite-difference BVP solver):

    y'' + y = 0
    y(0) = 0
    y(pi/2) = 1

Exact solution: y(x) = sin(x).

Outputs:
- Table printed to stdout
- SVG plot saved to test_output/example_imsl_bvpfd.svg
"""

from __future__ import annotations

from pathlib import Path
from typing import Dict

import matplotlib.pyplot as plt
import numpy as np

from differentialequations import bvpfd


def run_demo_imsl_bvpfd() -> Dict[str, object]:
    """Run IMSL BVPFD example: y'' + y = 0 with Dirichlet BCs.

    Args:
        None

    Returns:
        Dict[str, object]: Result dict with keys ``x`` (ndarray),
            ``y`` (ndarray), ``error_max`` (float), and ``plot_path`` (str).
    """
    def ode(x, y):
        """ODE system: [y1' = y2, y2' = -y1].

        Args:
            x (float | np.ndarray): Independent variable.
            y (np.ndarray): State vector of shape (2, N) or (2,).

        Returns:
            np.ndarray: Derivatives stacked vertically.
        """
        return np.vstack([y[1], -y[0]])

    def bc(ya, yb):
        """Boundary conditions: y(0)=0, y(pi/2)=1.

        Args:
            ya (np.ndarray): Solution at left boundary x=0.
            yb (np.ndarray): Solution at right boundary x=pi/2.

        Returns:
            np.ndarray: BC residuals of shape (2,).
        """
        return np.array([ya[0], yb[0] - 1.0])

    x_mesh = np.linspace(0, np.pi / 2, 30)
    y_guess = np.zeros((2, x_mesh.size))
    y_guess[0] = np.sin(x_mesh)

    result = bvpfd(ode, bc, x_mesh, y_guess, tol=1e-5)

    exact = np.sin(result.t)
    error_max = float(np.max(np.abs(result.y[0] - exact)))

    print("\nIMSL BVPFD Example: y'' + y = 0, y(0)=0, y(π/2)=1")
    print("-" * 55)
    print(f"{'Parameter':<30} {'Value':>20}")
    print("-" * 55)
    print(f"{'Domain':<30} {'[0, π/2]':>20}")
    print(f"{'Exact solution':<30} {'sin(x)':>20}")
    print(f"{'Converged':<30} {str(result.success):>20}")
    print(f"{'Mesh points':<30} {result.t.size:>20}")
    print(f"{'Max absolute error':<30} {error_max:>20.2e}")
    print("-" * 55)

    output_dir = Path("test_output")
    output_dir.mkdir(parents=True, exist_ok=True)
    plot_path = output_dir / "example_imsl_bvpfd.svg"

    x_fine = np.linspace(0, np.pi / 2, 200)
    fig, axes = plt.subplots(1, 2, figsize=(12, 5))

    axes[0].plot(x_fine, np.sin(x_fine), color="#999999", linewidth=2.5,
                 label="Exact: sin(x)", linestyle="--")
    axes[0].plot(result.t, result.y[0], color="#0369a1", linewidth=1.5,
                 marker="o", markersize=4, label="BVPFD solution")
    axes[0].set_xlabel("x")
    axes[0].set_ylabel("y(x)")
    axes[0].set_title("BVP Solution: y'' + y = 0")
    axes[0].legend()
    axes[0].grid(True, alpha=0.3)

    axes[1].semilogy(result.t, np.abs(result.y[0] - np.sin(result.t)) + 1e-16,
                     color="#0891b2", linewidth=1.5, marker="s", markersize=4)
    axes[1].set_xlabel("x")
    axes[1].set_ylabel("|error|")
    axes[1].set_title("Absolute Error vs. Exact Solution")
    axes[1].grid(True, alpha=0.3)

    fig.suptitle("IMSL BVPFD: Two-Point BVP y'' + y = 0", fontsize=13, fontweight="bold")
    fig.tight_layout()
    fig.savefig(plot_path, format="svg")
    plt.close(fig)

    return {"x": result.t, "y": result.y[0], "error_max": error_max, "plot_path": str(plot_path)}


if __name__ == "__main__":
    run_demo_imsl_bvpfd()

Input (Console)

Run the BVP solver script from the package root:

python examples/example_imsl_bvpfd.py

Console Output

IMSL BVPFD Example: y'' + y = 0, y(0)=0, y(π/2)=1
-------------------------------------------------------
Parameter                                     Value
-------------------------------------------------------
Domain                                     [0, π/2]
Exact solution                               sin(x)
Converged                                      True
Mesh points                                      30
Max absolute error                         6.71e-09
-------------------------------------------------------

Plot Output

Generated SVG plot:

Description

This example solves the two-point boundary value problem:

\[y'' + y = 0, \quad y(0) = 0, \quad y(\pi/2) = 1\]

whose exact solution is \(y(x) = \sin(x)\). bvpfd (finite-difference BVP solver backed by scipy.integrate.solve_bvp) is used. The left panel compares the numerical solution against the exact result, and the right panel shows the absolute error on a log scale.