R Examples

Complete worked examples with plots showing predicted vs actual results

1. Core DMD -- Predicted vs Actual

This example fits a standard DMD model to a damped oscillator and predicts beyond the training window.

library(koopman.dmd)

# Create a damped oscillator signal
dt <- 0.05
n_train <- 200; n_total <- 300
t_all <- seq(0, (n_total - 1) * dt, by = dt)
x1 <- exp(-0.05 * t_all) * sin(2 * pi * 0.5 * t_all)
x2 <- exp(-0.05 * t_all) * cos(2 * pi * 0.5 * t_all)
X_train <- rbind(x1[1:n_train], x2[1:n_train])

# Fit DMD with rank 2
result <- dmd(X_train, rank = 2, dt = dt)
summary(result)

# Predict the full time range (including out-of-sample)
pred <- predict(result, n_ahead = n_total - 1)

Core DMD predicted vs actual for a damped oscillator signal

DMD captures the damped oscillation dynamics perfectly. The dashed line extends beyond the training window (vertical line) and matches the true signal.

2. Eigenvalue Spectrum

Examining the eigenvalues in the complex plane reveals the stability and frequency content of the decomposition.

# Eigenvalues in the complex plane
spec <- dmd_spectrum(result)
print(spec)

# Stability analysis
stab <- dmd_stability(result)
cat("Spectral radius:", stab$spectral_radius, "\n")
cat("Stable:", stab$is_stable, "\n")

DMD eigenvalues plotted in the complex plane relative to the unit circle

Both eigenvalues lie just inside the unit circle (|lambda|=0.9975), confirming the slight damping. The frequency 0.5 Hz matches the input signal.

3. Reconstruction Residuals

Reconstruction error quantifies how well DMD approximates the original data.

recon <- dmd_reconstruct(result)
err <- dmd_error(result)
cat("RMSE:", err$rmse, "\n")
cat("Relative error:", err$relative_error, "\n")

DMD reconstruction residuals over time showing near-zero error

For a purely linear system with rank matching the true rank, DMD achieves near-zero reconstruction error.

4. Extended DMD with Polynomial Lifting

When the dynamics are nonlinear, polynomial lifting maps the data into a higher-dimensional space where DMD can better approximate the Koopman operator.

# Nonlinear signal: sin^2(t)
t3 <- seq(0, 9.95, by = 0.05)
X3 <- matrix(sin(t3)^2, nrow = 1)

# Standard DMD
res_std <- dmd(X3, dt = 0.05)

# Extended DMD with degree-2 polynomial lifting
res_lift <- dmd(X3, lifting = "polynomial", lifting_param = 2, dt = 0.05)

# Compare errors
cat("Standard RMSE:", dmd_error(res_std)$rmse, "\n")
cat("Lifted RMSE:  ", dmd_error(res_lift)$rmse, "\n")

Comparison of standard DMD vs extended DMD with polynomial lifting on a nonlinear signal

Polynomial lifting maps the nonlinear sin^2 signal into a space where DMD can better approximate the dynamics.

5. Hankel-DMD for Scalar Time Series

Hankel-DMD augments a scalar measurement with time-delay embeddings, enabling DMD to recover multiple frequencies from a single observable.

# Two-frequency scalar signal
dt4 <- 0.02; t4 <- seq(0, 9.98, by = dt4)
signal <- sin(2 * pi * 1.0 * t4) + 0.4 * sin(2 * pi * 3.0 * t4)
X4 <- matrix(signal, nrow = 1)

# Train on first 400 samples, predict last 100
hresult <- hankel_dmd(X4[, 1:400, drop = FALSE], delays = 30, rank = 4, dt = dt4)

# Reconstruct and predict
hrecon <- hankel_reconstruct(hresult, 400)
hpred <- predict(hresult, n_ahead = 100)

# Recovered frequencies
eig_freq <- abs(atan2(hresult$eigenvalues_im, hresult$eigenvalues_re)) / (2 * pi * dt4)
cat("Recovered frequencies:", round(eig_freq, 2), "Hz\n")

Hankel-DMD reconstruction and prediction of a two-frequency scalar signal

Hankel-DMD recovers the correct frequencies (1.0 Hz and 3.0 Hz) from a scalar measurement via time-delay embedding.

6. GLA Reconstruction and Prediction

Generalized Laplace Analysis (GLA) directly computes Koopman eigenfunctions from trajectory data using an iterative algorithm.

# Simple oscillator
t5 <- seq(0, 39.9, by = 0.1)
gla_data <- rbind(sin(t5), cos(t5))

# GLA with 2 eigenvalues
gresult <- gla(gla_data, n_eigenvalues = 2, tol = 1e-4)
grecon <- gla_reconstruct(gresult)
gpred <- predict(gresult, n_ahead = 100)

GLA reconstruction and out-of-sample prediction for a pure oscillator

GLA directly computes Koopman eigenfunctions. For a pure oscillator, it achieves exact reconstruction and accurate extrapolation.

7. Dynamical Map Phase Portraits

The library includes several classic dynamical systems for generating trajectories and studying phase space structure.

# Standard map -- regular and chaotic
traj_reg <- generate_trajectory("standard", c(0.1, 0.2), 5000, epsilon = 0.12)
traj_chaos <- generate_trajectory("standard", c(0.5, 0.5), 50000, epsilon = 0.97)

# Henon attractor
traj_henon <- generate_trajectory("henon", c(0, 0), 20000, a = 1.4, b = 0.3)

# Logistic map
traj_log <- generate_trajectory("logistic", 0.4, 200, r = 3.9)

Phase portraits of the standard map, Henon attractor, and logistic map

The library includes classic dynamical systems. The standard map transitions from regular islands to global chaos as epsilon increases. The Henon map shows a strange attractor, and the logistic map exhibits chaotic time evolution.

8. HTA Convergence

The harmonic time average (HTA) converges differently for regular and chaotic orbits, providing a diagnostic for orbit classification.

# Regular orbit converges to a non-zero value
conv_reg <- hta_convergence("standard", c(0.1, 0.2),
                             observable = "sin_pi", omega = 0.0,
                             n_iter = 50000, epsilon = 0.12)

# Chaotic orbit decays toward zero
conv_chaos <- hta_convergence("standard", c(0.5, 0.5),
                               observable = "sin_pi", omega = 0.0,
                               n_iter = 50000, epsilon = 0.97)

HTA convergence comparison between regular and chaotic orbits

The HTA magnitude converges to a stable non-zero value for regular orbits (phase-locked to the frequency), while it decays for chaotic orbits (no persistent periodicity).

9. Mesochronic Harmonic Plot

The mesochronic plot computes the HTA over a grid of initial conditions, revealing the phase space structure of a dynamical system.

meso <- mesochronic_compute("standard",
                             x_range = c(0, 1), y_range = c(0, 1),
                             resolution = 100, observable = "sin_pi",
                             omega = 0.0, n_iter = 5000, epsilon = 0.12)

# Plot magnitude heatmap
image(meso$x_coords, meso$y_coords, meso$hta_matrix,
      col = hcl.colors(256, "viridis"),
      xlab = "x", ylab = "y", main = "|HTA| Magnitude")

Mesochronic heatmap of the standard map showing regular islands and chaotic regions

The mesochronic plot reveals the phase space structure of the standard map. Bright regions correspond to regular islands where orbits resonate with the test frequency. Dark regions indicate chaotic motion.

10. Method Comparison

This example applies all three decomposition methods to the same multi-frequency signal and compares their out-of-sample prediction accuracy.

# Multi-frequency test signal
dt10 <- 0.01; t10 <- seq(0, 2.99, by = dt10)
signal <- sin(2 * pi * 2 * t10) + 0.3 * sin(2 * pi * 7 * t10)
X10 <- matrix(signal, nrow = 1)

# Train on 250 samples, predict 50
X10_train <- X10[, 1:250, drop = FALSE]

dmd_res <- dmd(X10_train, rank = 4, dt = dt10)
hdmd_res <- hankel_dmd(X10_train, delays = 30, rank = 4, dt = dt10)
gla_res <- gla(X10_train, n_eigenvalues = 4, tol = 1e-4)

# Compare prediction RMSE
dmd_pred <- predict(dmd_res, n_ahead = 50)
hdmd_pred <- predict(hdmd_res, n_ahead = 50)
gla_pred <- predict(gla_res, n_ahead = 50)

Out-of-sample prediction comparison between standard DMD, Hankel-DMD, and GLA

All three methods are compared on the same out-of-sample prediction task. Performance varies by signal characteristics -- standard DMD works well for multi-variable systems, Hankel-DMD excels at scalar frequency extraction, and GLA provides direct eigenfunction computation.