DAE / Stiff ODE Examples

Example Code

Representative script:

"""IMSL DASPG example: Brusselator reaction-diffusion system.

Demonstrates daspg on the Brusselator ODE system (no spatial dimension):

    x' = 1 - (b+1)*x + a*x^2*y
    y' = b*x - a*x^2*y

with a=1, b=3, x(0)=1.5, y(0)=3.0.

The Brusselator exhibits limit-cycle behaviour.

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


def run_demo_imsl_daspg() -> Dict[str, object]:
    """Run IMSL DASPG example: Brusselator limit-cycle ODE.

    Args:
        None

    Returns:
        Dict[str, object]: Result dict with keys ``t`` (ndarray),
            ``y`` (ndarray), and ``plot_path`` (str).
    """
    a, b = 1.0, 3.0
    x0, y0 = 1.5, 3.0
    tf = 20.0

    def brusselator(t, y):
        """Brusselator right-hand side.

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

        Returns:
            list: Derivatives [x', y'].
        """
        return [
            1.0 - (b + 1.0) * y[0] + a * y[0] ** 2 * y[1],
            b * y[0] - a * y[0] ** 2 * y[1],
        ]

    yprime0 = brusselator(0.0, [x0, y0])
    result = daspg(brusselator, (0.0, tf), [x0, y0], yprime0,
                   rtol=1e-7, atol=1e-9)

    print("\nIMSL DASPG Example: Brusselator System (a=1, b=3)")
    print("-" * 55)
    print(f"{'Parameter':<28} {'Value':>20}")
    print("-" * 55)
    print(f"{'Integration interval':<28} {f'[0, {tf}]':>20}")
    print(f"{'a':<28} {a:>20}")
    print(f"{'b':<28} {b:>20}")
    print(f"{'x(0)':<28} {x0:>20}")
    print(f"{'y(0)':<28} {y0:>20}")
    print(f"{'Converged':<28} {str(result.success):>20}")
    print(f"{'Steps taken':<28} {result.n_steps:>20}")
    print(f"{'x final':<28} {result.y[0, -1]:>20.6f}")
    print(f"{'y final':<28} {result.y[1, -1]:>20.6f}")
    print("-" * 55)

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

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

    axes[0].plot(result.t, result.y[0], color="#0369a1", linewidth=2.0, label="x(t)")
    axes[0].plot(result.t, result.y[1], color="#0891b2", linewidth=1.5,
                 linestyle="--", label="y(t)")
    axes[0].set_xlabel("t")
    axes[0].set_ylabel("Concentration")
    axes[0].set_title("Brusselator: Time Evolution")
    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].scatter([x0], [y0], color="#d62728", s=60, zorder=5, label="Initial condition")
    axes[1].set_xlabel("x")
    axes[1].set_ylabel("y")
    axes[1].set_title("Brusselator: Phase Portrait (Limit Cycle)")
    axes[1].legend()
    axes[1].grid(True, alpha=0.3)

    fig.suptitle("IMSL DASPG: Brusselator Reaction System", 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_daspg()

Input (Console)

Run the DAE/stiff solver script from the package root:

python examples/example_imsl_daspg.py

Console Output

IMSL DASPG Example: Brusselator System (a=1, b=3)
-------------------------------------------------------
Parameter                                   Value
-------------------------------------------------------
Integration interval                    [0, 20.0]
a                                             1.0
b                                             3.0
x(0)                                          1.5
y(0)                                          3.0
Converged                                    True
Steps taken                                   578
x final                                  0.498637
y final                                  4.596780
-------------------------------------------------------

Plot Output

Generated SVG plot:

Description

This example solves the Brusselator ODE system (a=1, b=3):

\[ \begin{align}\begin{aligned}x' = 1 - (b+1)x + ax^2 y\\y' = bx - ax^2 y\end{aligned}\end{align} \]

using daspg (backed by the Radau implicit solver). The system exhibits a stable limit cycle. The left panel shows the time evolution and the right panel shows the phase portrait converging to the limit cycle attractor.