IVP Solver Examples

Example Code

Representative script:

"""IMSL IVPAG example: Van der Pol oscillator.

Demonstrates the use of ivpag with method='Radau' for a mildly stiff
Van der Pol oscillator (mu=1):

    y1' = y2
    y2' = mu*(1 - y1^2)*y2 - y1

with initial conditions y1(0) = 2, y2(0) = 0.

Outputs:
- Table printed to stdout
- SVG plot saved to test_output/example_imsl_ivpag.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 ivpag


def run_demo_imsl_ivpag() -> Dict[str, object]:
    """Run IMSL IVPAG example: Van der Pol oscillator with mu=1.

    Args:
        None

    Returns:
        Dict[str, object]: Result dict with keys ``t`` (ndarray),
            ``y`` (ndarray), and ``plot_path`` (str).
    """
    mu = 1.0

    def van_der_pol(t, y):
        """Van der Pol right-hand side.

        Args:
            t (float): Time.
            y (list): State [y1, y2].

        Returns:
            list: Derivatives [y1', y2'].
        """
        return [y[1], mu * (1.0 - y[0] ** 2) * y[1] - y[0]]

    result = ivpag(van_der_pol, (0.0, 20.0), [2.0, 0.0],
                   method="Radau", rtol=1e-6, atol=1e-9)

    print("\nIMSL IVPAG Example: Van der Pol Oscillator (mu=1)")
    print("-" * 55)
    print(f"{'Parameter':<25} {'Value':>20}")
    print("-" * 55)
    print(f"{'Integration interval':<25} {'[0, 20]':>20}")
    print(f"{'Method':<25} {'Radau':>20}")
    print(f"{'y1(0)':<25} {'2.0':>20}")
    print(f"{'y2(0)':<25} {'0.0':>20}")
    print(f"{'Converged':<25} {str(result.success):>20}")
    print(f"{'Steps taken':<25} {result.n_steps:>20}")
    print(f"{'Function evaluations':<25} {result.n_eval:>20}")
    print(f"{'y1 final':<25} {result.y[0, -1]:>20.6f}")
    print("-" * 55)

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

    fig, axes = plt.subplots(1, 2, figsize=(12, 5))

    axes[0].plot(result.t, result.y[0], color="#0369a1", linewidth=2.0, label=r"$y_1(t)$")
    axes[0].plot(result.t, result.y[1], color="#0891b2", linewidth=1.5,
                 linestyle="--", label=r"$y_2(t)$")
    axes[0].set_xlabel("t")
    axes[0].set_ylabel("y")
    axes[0].set_title("Van der Pol Oscillator: Time Domain")
    axes[0].legend()
    axes[0].grid(True, alpha=0.3)

    axes[1].plot(result.y[0], result.y[1], color="#0369a1", linewidth=1.5)
    axes[1].set_xlabel(r"$y_1$")
    axes[1].set_ylabel(r"$y_2$")
    axes[1].set_title("Van der Pol: Phase Portrait")
    axes[1].grid(True, alpha=0.3)

    fig.suptitle("IMSL IVPAG: Van der Pol Oscillator (μ=1)", fontsize=13, fontweight="bold")
    fig.tight_layout()
    fig.savefig(plot_path, format="svg")
    plt.close(fig)

    return {"t": result.t, "y": result.y, "plot_path": str(plot_path)}


if __name__ == "__main__":
    run_demo_imsl_ivpag()

Input (Console)

Run the IVP solver script from the package root:

python examples/example_imsl_ivpag.py

Console Output

IMSL IVPAG Example: Van der Pol Oscillator (mu=1)
-------------------------------------------------------
Parameter                                Value
-------------------------------------------------------
Integration interval                   [0, 20]
Method                                   Radau
y1(0)                                      2.0
y2(0)                                      0.0
Converged                                 True
Steps taken                                496
Function evaluations                      3767
y1 final                              2.008150
-------------------------------------------------------

Plot Output

Generated SVG plot:

Description

This example solves the Van der Pol oscillator (mu=1) using ivpag with the Radau implicit Runge-Kutta method. The left panel shows the time evolution of both state variables, and the right panel shows the phase portrait converging to the limit cycle.