Pendulum

System Description

Driven damped pendulum:

\[ \begin{aligned} \dot{\theta} &= \omega \\ \dot{\omega} &= -\alpha \omega - K \sin(\theta) + T \end{aligned} \]

System Parameters

Parameter	Symbol	Value
Damping coefficient	\(\alpha\)	0.1
Torque	\(T\)	0.5
Stiffness	\(K\)	1.0

Sampling

Dimension: \(D = 2\)
Sample size: \(N = 10000\)
Distribution: \(\rho\) = Uniform
Region of interest: \(\mathcal{Q}(\theta, \omega) : [\psi - \pi, \psi + \pi] \times [-10, 10]\) where \(\psi = \arcsin(T/K)\)

Solver

Setting	Value
Method	Dopri5 (Diffrax)
Time span	\([0, 1000]\)
Steps	1000 (\(f_s\) = 1 Hz)
Relative tolerance	1e-08
Absolute tolerance	1e-06

Feature Extraction

Log-delta feature on angular velocity:

States: \(\omega\) (state 1)
Formula: \(\Delta = \log_{10}(|\max(\omega) - \text{mean}(\omega)| + \epsilon)\)
Transient cutoff: \(t^* = 950.0\)

Clustering

Method: k-NN (k=1)
Template ICs:
FP: \([0.4, 0.0]\) — Fixed point (stable equilibrium)
LC: \([2.7, 0.0]\) — Limit cycle (rotational motion)

Reproduction Code

Setup

def setup_pendulum_system() -> SetupProperties:
    n = 10000

    device = "cuda" if torch.cuda.is_available() else "cpu"
    print(f"Setting up pendulum system on device: {device}")

    params: PendulumParams = {"alpha": 0.1, "T": 0.5, "K": 1.0}

    ode_system = PendulumJaxODE(params)

    sampler = UniformRandomSampler(
        min_limits=[-np.pi + np.arcsin(params["T"] / params["K"]), -10.0],
        max_limits=[np.pi + np.arcsin(params["T"] / params["K"]), 10.0],
        device=device,
    )

    solver = JaxSolver(
        t_span=(0, 1000),
        t_steps=1000,
        # t_eval=(950.0, 1000.0),
        device=device,
        rtol=1e-8,
        atol=1e-6,
        cache_dir=".pybasin_cache/pendulum",
    )

    feature_extractor = TorchFeatureExtractor(
        time_steady=950.0,
        features=None,
        features_per_state={
            1: {"log_delta": None},
        },
        normalize=False,
    )

    template_y0 = [
        [0.4, 0.0],
        [2.7, 0.0],
    ]
    classifier_labels = ["FP", "LC"]

    knn = KNeighborsClassifier(n_neighbors=1)

    template_integrator = TemplateIntegrator(
        template_y0=template_y0,
        labels=classifier_labels,
        ode_params=params,
    )

    return {
        "n": n,
        "ode_system": ode_system,
        "sampler": sampler,
        "solver": solver,
        "feature_extractor": feature_extractor,
        "estimator": knn,
        "template_integrator": template_integrator,
    }

Main Estimation

def main():
    props = setup_pendulum_system()

    bse = BasinStabilityEstimator(
        n=props["n"],
        ode_system=props["ode_system"],
        sampler=props["sampler"],
        solver=props.get("solver"),
        feature_extractor=props.get("feature_extractor"),
        predictor=props.get("estimator"),
        template_integrator=props.get("template_integrator"),
        output_dir="results_case1",
        feature_selector=None,
    )

    bse.run()

    return bse

Case 1: Baseline Results

Comparison with MATLAB bSTAB

Overall Classification Quality:

Matthews Correlation Coefficient: 1.0000

Attractor	pybasin BS ± SE	bSTAB BS ± SE
FP	0.15200 ± 0.00359	0.15200 ± 0.00359
LC	0.84800 ± 0.00359	0.84800 ± 0.00359

Visualizations

Basin Stability

State Space

Feature Space

Template Trajectories

Trajectories

Case 2: Parameter Sweep

Comparison with MATLAB bSTAB

Average MCC = 1.0000

The average excludes cases where there is only a single attractor and the basin stability values are the same since MCC is 0 for single class cases, and would therefore drop the average.

Parameter	Attractor	pybasin BS ± SE	bSTAB BS ± SE	MCC
0.01	FP	1.00000 ± 0.00000	1.00000 ± 0.00000	0.0000
0.06	FP	1.00000 ± 0.00000	1.00000 ± 0.00000	0.0000
0.11	FP	1.00000 ± 0.00000	1.00000 ± 0.00000	0.0000
0.16	FP	0.48300 ± 0.00500	0.48300 ± 0.00500	1.0000
	LC	0.51700 ± 0.00500	0.51700 ± 0.00500
0.21	FP	0.39670 ± 0.00489	0.39670 ± 0.00489	1.0000
	LC	0.60330 ± 0.00489	0.60330 ± 0.00489
0.26	FP	0.32680 ± 0.00469	0.32680 ± 0.00469	1.0000
	LC	0.67320 ± 0.00469	0.67320 ± 0.00469
0.31	FP	0.27560 ± 0.00447	0.27560 ± 0.00447	1.0000
	LC	0.72440 ± 0.00447	0.72440 ± 0.00447
0.36	FP	0.23800 ± 0.00426	0.23800 ± 0.00426	1.0000
	LC	0.76200 ± 0.00426	0.76200 ± 0.00426
0.41	FP	0.20650 ± 0.00405	0.20650 ± 0.00405	1.0000
	LC	0.79350 ± 0.00405	0.79350 ± 0.00405
0.46	FP	0.17030 ± 0.00376	0.17030 ± 0.00376	1.0000
	LC	0.82970 ± 0.00376	0.82970 ± 0.00376
0.51	FP	0.14700 ± 0.00354	0.14700 ± 0.00354	1.0000
	LC	0.85300 ± 0.00354	0.85300 ± 0.00354
0.56	FP	0.12790 ± 0.00334	0.12790 ± 0.00334	1.0000
	LC	0.87210 ± 0.00334	0.87210 ± 0.00334
0.61	FP	0.10340 ± 0.00304	0.10340 ± 0.00304	1.0000
	LC	0.89660 ± 0.00304	0.89660 ± 0.00304
0.66	FP	0.08390 ± 0.00277	0.08390 ± 0.00277	1.0000
	LC	0.91610 ± 0.00277	0.91610 ± 0.00277
0.71	FP	0.06530 ± 0.00247	0.06530 ± 0.00247	1.0000
	LC	0.93470 ± 0.00247	0.93470 ± 0.00247
0.76	FP	0.05450 ± 0.00227	0.05450 ± 0.00227	1.0000
	LC	0.94550 ± 0.00227	0.94550 ± 0.00227
0.81	FP	0.03910 ± 0.00194	0.03910 ± 0.00194	1.0000
	LC	0.96090 ± 0.00194	0.96090 ± 0.00194
0.86	FP	0.02440 ± 0.00154	0.02440 ± 0.00154	1.0000
	LC	0.97560 ± 0.00154	0.97560 ± 0.00154
0.91	FP	0.01720 ± 0.00130	0.01720 ± 0.00130	1.0000
	LC	0.98280 ± 0.00130	0.98280 ± 0.00130
0.96	FP	0.00710 ± 0.00084	0.00710 ± 0.00084	1.0000
	LC	0.99290 ± 0.00084	0.99290 ± 0.00084

Visualizations

Basin Stability Variation

Bifurcation Diagram

Case 3: Sample Size Convergence Study

This hyperparameter study varies the number of initial conditions \(N\) from ~50 to ~5000 (using \(5 \times \text{logspace}(1, 3, 20)\)) to assess how basin stability estimates converge as sample size increases. The relative standard error decreases as \(\text{SE}/\mathcal{S}_{\mathcal{B}} \sim 1/\sqrt{N}\).

Comparison with MATLAB bSTAB

Average MCC = 1.0000

Parameter	Attractor	pybasin BS ± SE	bSTAB BS ± SE	MCC
50	FP	0.10000 ± 0.04243	0.10000 ± 0.04243	1.0000
	LC	0.90000 ± 0.04243	0.90000 ± 0.04243
63.7137	FP	0.09375 ± 0.03644	0.09375 ± 0.03644	1.0000
	LC	0.90625 ± 0.03644	0.90625 ± 0.03644
81.1888	FP	0.15854 ± 0.04033	0.15854 ± 0.04033	1.0000
	LC	0.84146 ± 0.04033	0.84146 ± 0.04033
103.4569	FP	0.17308 ± 0.03710	0.17308 ± 0.03710	1.0000
	LC	0.82692 ± 0.03710	0.82692 ± 0.03710
131.8325	FP	0.15152 ± 0.03121	0.15152 ± 0.03121	1.0000
	LC	0.84848 ± 0.03121	0.84848 ± 0.03121
167.9909	FP	0.16667 ± 0.02875	0.16667 ± 0.02875	1.0000
	LC	0.83333 ± 0.02875	0.83333 ± 0.02875
214.0666	FP	0.12558 ± 0.02260	0.12558 ± 0.02260	1.0000
	LC	0.87442 ± 0.02260	0.87442 ± 0.02260
272.7797	FP	0.13187 ± 0.02048	0.13187 ± 0.02048	1.0000
	LC	0.86813 ± 0.02048	0.86813 ± 0.02048
347.5964	FP	0.14080 ± 0.01865	0.14080 ± 0.01865	1.0000
	LC	0.85920 ± 0.01865	0.85920 ± 0.01865
442.9334	FP	0.16253 ± 0.01753	0.16253 ± 0.01753	1.0000
	LC	0.83747 ± 0.01753	0.83747 ± 0.01753
564.4189	FP	0.16106 ± 0.01546	0.16106 ± 0.01546	1.0000
	LC	0.83894 ± 0.01546	0.83894 ± 0.01546
719.2249	FP	0.15278 ± 0.01341	0.15278 ± 0.01341	1.0000
	LC	0.84722 ± 0.01341	0.84722 ± 0.01341
916.4904	FP	0.12650 ± 0.01098	0.12650 ± 0.01098	1.0000
	LC	0.87350 ± 0.01098	0.87350 ± 0.01098
1167.8607	FP	0.15582 ± 0.01061	0.15582 ± 0.01061	1.0000
	LC	0.84418 ± 0.01061	0.84418 ± 0.01061
1488.1757	FP	0.16588 ± 0.00964	0.16588 ± 0.00964	1.0000
	LC	0.83412 ± 0.00964	0.83412 ± 0.00964
1896.3451	FP	0.14655 ± 0.00812	0.14655 ± 0.00812	1.0000
	LC	0.85345 ± 0.00812	0.85345 ± 0.00812
2416.4651	FP	0.16177 ± 0.00749	0.16177 ± 0.00749	1.0000
	LC	0.83823 ± 0.00749	0.83823 ± 0.00749
3079.2411	FP	0.15227 ± 0.00647	0.15227 ± 0.00647	1.0000
	LC	0.84773 ± 0.00647	0.84773 ± 0.00647
3923.7999	FP	0.15851 ± 0.00583	0.15851 ± 0.00583	1.0000
	LC	0.84149 ± 0.00583	0.84149 ± 0.00583
5000	FP	0.16180 ± 0.00521	0.16180 ± 0.00521	1.0000
	LC	0.83820 ± 0.00521	0.83820 ± 0.00521

Visualizations

Basin Stability Variation