| Title: | Koopman Operator and Dynamic Mode Decomposition for Dynamical Systems |
|---|---|
| Description: | Dynamic Mode Decomposition (DMD) with Koopman operator theory extensions, powered by a Rust backend via 'extendr'. Provides standard DMD as described in Schmid (2010) <doi:10.1017/S0022112010001217>, DMD with control for forced linear systems following Proctor, Brunton, and Kutz (2016) <doi:10.1137/15M1013857>, Extended DMD with lifting functions, Hankel-DMD via time-delay embedding, Generalized Laplace Analysis for direct eigenfunction computation, and harmonic time averages and mesochronic harmonic plots for phase space analysis as developed in Mezic (2020) <doi:10.48550/arXiv.2009.05883>. Includes built-in area-preserving and chaotic maps for experimentation. |
| Authors: | James Harris [aut, cre, cph], The authors of the dependency Rust crates [ctb] (see inst/AUTHORS file for details) |
| Maintainer: | James Harris <[email protected]> |
| License: | MIT + file LICENSE |
| Version: | 0.2.2 |
| Built: | 2026-10-15 09:49:30 UTC |
| Source: | https://github.com/cran/koopman.dmd |
Dynamic Mode Decomposition (DMD) with Koopman operator theory extensions, powered by a Rust backend via 'extendr'. Provides standard DMD as described in Schmid (2010) doi:10.1017/S0022112010001217, DMD with control for forced linear systems following Proctor, Brunton, and Kutz (2016) doi:10.1137/15M1013857, Extended DMD with lifting functions, Hankel-DMD via time-delay embedding, Generalized Laplace Analysis for direct eigenfunction computation, and harmonic time averages and mesochronic harmonic plots for phase space analysis as developed in Mezic (2020) doi:10.48550/arXiv.2009.05883. Includes built-in area-preserving and chaotic maps for experimentation.
Maintainer: James Harris [email protected] [copyright holder]
Other contributors:
The authors of the dependency Rust crates (see inst/AUTHORS file for details) [contributor]
Useful links:
Report bugs at https://github.com/jimeharrisjr/rust-dmd/issues
Classify phase space points by HTA magnitude
classify_phase_space( hta_magnitudes, resonating_threshold = 0.01, chaotic_threshold = 1e-04 )classify_phase_space( hta_magnitudes, resonating_threshold = 0.01, chaotic_threshold = 1e-04 )
hta_magnitudes |
Numeric vector of |HTA| values. |
resonating_threshold |
Threshold for resonating. Default 0.01. |
chaotic_threshold |
Threshold for chaotic. Default 0.0001. |
Integer vector (1=resonating, 2=chaotic, 3=non-resonating).
# Classify orbits from their HTA magnitudes mags <- c(0.5, 0.001, 1e-06, 0.1) classify_phase_space(mags) # Thresholds are adjustable classify_phase_space(mags, resonating_threshold = 0.1, chaotic_threshold = 0.001)# Classify orbits from their HTA magnitudes mags <- c(0.5, 0.001, 1e-06, 0.1) classify_phase_space(mags) # Thresholds are adjustable classify_phase_space(mags, resonating_threshold = 0.1, chaotic_threshold = 0.001)
Perform Dynamic Mode Decomposition on time-series data.
dmd( X, rank = NULL, center = FALSE, dt = 1, lifting = NULL, lifting_param = NULL )dmd( X, rank = NULL, center = FALSE, dt = 1, lifting = NULL, lifting_param = NULL )
X |
Numeric matrix (n_vars x n_time). |
rank |
Integer truncation rank, or NULL for automatic. |
center |
Logical, center data by subtracting row means. |
dt |
Numeric time step. |
lifting |
Character lifting type or NULL. |
lifting_param |
Integer lifting parameter or NULL. |
An S3 object of class "dmd".
# One row per variable, one column per time step t <- seq(0, 10, length.out = 100) X <- rbind(sin(t), cos(t)) d <- dmd(X, rank = 2, dt = t[2] - t[1]) print(d) summary(d) # Extended DMD: lift into a higher-dimensional space where nonlinear # dynamics become approximately linear d2 <- dmd(X, lifting = "polynomial", lifting_param = 2)# One row per variable, one column per time step t <- seq(0, 10, length.out = 100) X <- rbind(sin(t), cos(t)) d <- dmd(X, rank = 2, dt = t[2] - t[1]) print(d) summary(d) # Extended DMD: lift into a higher-dimensional space where nonlinear # dynamics become approximately linear d2 <- dmd(X, lifting = "polynomial", lifting_param = 2)
DMD dominant modes
dmd_dominant_modes( object, n = 3, criterion = c("amplitude", "energy", "stability") )dmd_dominant_modes( object, n = 3, criterion = c("amplitude", "energy", "stability") )
object |
A dmd object. |
n |
Number of modes. |
criterion |
"amplitude", "energy", or "stability". |
Integer vector of 1-based mode indices.
t <- seq(0, 10, length.out = 100) X <- rbind(sin(t), cos(t)) d <- dmd(X, rank = 2, dt = t[2] - t[1]) # Indices of the most significant modes dmd_dominant_modes(d, n = 1)t <- seq(0, 10, length.out = 100) X <- rbind(sin(t), cos(t)) d <- dmd(X, rank = 2, dt = t[2] - t[1]) # Indices of the most significant modes dmd_dominant_modes(d, n = 1)
DMD reconstruction error
dmd_error(object)dmd_error(object)
object |
A dmd object. |
List with error metrics (rmse, mae, mape, relative_error).
t <- seq(0, 10, length.out = 100) X <- rbind(sin(t), cos(t)) d <- dmd(X, rank = 2, dt = t[2] - t[1]) # Reconstruction error metrics dmd_error(d)t <- seq(0, 10, length.out = 100) X <- rbind(sin(t), cos(t)) d <- dmd(X, rank = 2, dt = t[2] - t[1]) # Reconstruction error metrics dmd_error(d)
Reconstruct data from DMD modes
dmd_reconstruct(object, n_steps = NULL, modes = NULL)dmd_reconstruct(object, n_steps = NULL, modes = NULL)
object |
A dmd object. |
n_steps |
Number of time steps. Defaults to original length. |
modes |
Integer vector of mode indices (1-based), or NULL for all. |
Numeric matrix.
t <- seq(0, 10, length.out = 100) X <- rbind(sin(t), cos(t)) d <- dmd(X, rank = 2, dt = t[2] - t[1]) recon <- dmd_reconstruct(d, n_steps = ncol(X)) dim(recon)t <- seq(0, 10, length.out = 100) X <- rbind(sin(t), cos(t)) d <- dmd(X, rank = 2, dt = t[2] - t[1]) recon <- dmd_reconstruct(d, n_steps = ncol(X)) dim(recon)
DMD residual analysis
dmd_residual(object)dmd_residual(object)
object |
A dmd object. |
List with residual_norm and residual_relative.
t <- seq(0, 10, length.out = 100) X <- rbind(sin(t), cos(t)) d <- dmd(X, rank = 2, dt = t[2] - t[1]) dmd_residual(d)t <- seq(0, 10, length.out = 100) X <- rbind(sin(t), cos(t)) d <- dmd(X, rank = 2, dt = t[2] - t[1]) dmd_residual(d)
DMD spectrum analysis
dmd_spectrum(object, dt = NULL)dmd_spectrum(object, dt = NULL)
object |
A dmd object. |
dt |
Time step. Uses stored dt by default. |
Data frame with mode information.
t <- seq(0, 10, length.out = 100) X <- rbind(sin(t), cos(t)) d <- dmd(X, rank = 2, dt = t[2] - t[1]) # Frequency, growth rate, amplitude and stability of each mode dmd_spectrum(d)t <- seq(0, 10, length.out = 100) X <- rbind(sin(t), cos(t)) d <- dmd(X, rank = 2, dt = t[2] - t[1]) # Frequency, growth rate, amplitude and stability of each mode dmd_spectrum(d)
DMD stability analysis
dmd_stability(object, tol = 1e-06)dmd_stability(object, tol = 1e-06)
object |
A dmd object. |
tol |
Tolerance for marginal classification. |
List with stability information.
t <- seq(0, 10, length.out = 100) X <- rbind(sin(t), cos(t)) d <- dmd(X, rank = 2, dt = t[2] - t[1]) dmd_stability(d)t <- seq(0, 10, length.out = 100) X <- rbind(sin(t), cos(t)) d <- dmd(X, rank = 2, dt = t[2] - t[1]) dmd_stability(d)
Identify the forced linear system from
snapshot pairs and control inputs, following Proctor, Brunton and Kutz
(2016). Unlike dmd, which takes one contiguous trajectory,
dmdc takes explicit pair matrices: X1 holds states at time
, X2 the states one step later, and U the control
input applied during each transition. Columns may therefore come from many
concatenated trajectories.
dmdc( X1, X2, U = NULL, rank_input = NULL, rank_output = NULL, dt = 1, known_B = NULL )dmdc( X1, X2, U = NULL, rank_input = NULL, rank_output = NULL, dt = 1, known_B = NULL )
X1 |
Numeric matrix of states at time t (n_states x n_pairs). |
X2 |
Numeric matrix of states at time t+1 (n_states x n_pairs). |
U |
Numeric matrix of control inputs during each transition (n_inputs x n_pairs), or NULL for an autonomous multi-trajectory fit. A vector is taken as a single input row. |
rank_input |
Integer truncation rank for the regression-input SVD, or NULL for automatic (99 percent cumulative variance). |
rank_output |
Integer rank of the output basis (SVD of |
dt |
Numeric time step between snapshot pairs. |
known_B |
Known input matrix B (n_states x n_inputs), or NULL to estimate B jointly with A. |
Two identification modes are available. With known_B = NULL (the
default), A and B are solved jointly from the stacked
regression ; this requires the input to be
persistently exciting and exogenous (not state feedback). When the input
coupling is known by construction, pass it as known_B and only
A is estimated.
An S3 object of class "dmdc" with components a, b,
a_tilde, b_tilde, basis, eigenvalues_re,
eigenvalues_im, singular_values, rank_input,
rank_output, dt, n_states, and n_inputs.
Proctor, J. L., Brunton, S. L., and Kutz, J. N. (2016). Dynamic Mode Decomposition with Control. SIAM Journal on Applied Dynamical Systems, 15(1), 142-161. doi:10.1137/15M1013857
dmd for autonomous systems,
predict.dmdc to simulate the identified system.
# Simulate x_{t+1} = A0 x_t + B0 u_t with a persistently exciting input A0 <- matrix(c(0.9, 0, 0.1, 0.8), 2, 2) B0 <- matrix(c(0.5, 1), 2, 1) m <- 120 X1 <- matrix(0, 2, m) X2 <- matrix(0, 2, m) U <- matrix(0, 1, m) x <- c(1, -0.5) for (t in seq_len(m)) { u_t <- sin(0.7 * (t - 1)) + 0.5 * cos(2.3 * (t - 1) + 1) X1[, t] <- x U[, t] <- u_t x <- as.numeric(A0 %*% x + B0 * u_t) X2[, t] <- x } # Identify A and B jointly; both are recovered to machine precision fit <- dmdc(X1, X2, U, rank_input = 3) round(fit$a, 6) round(fit$b, 6) # Known B: pin the input matrix and estimate only A fit2 <- dmdc(X1, X2, U, rank_input = 2, known_B = B0) round(fit2$a, 6)# Simulate x_{t+1} = A0 x_t + B0 u_t with a persistently exciting input A0 <- matrix(c(0.9, 0, 0.1, 0.8), 2, 2) B0 <- matrix(c(0.5, 1), 2, 1) m <- 120 X1 <- matrix(0, 2, m) X2 <- matrix(0, 2, m) U <- matrix(0, 1, m) x <- c(1, -0.5) for (t in seq_len(m)) { u_t <- sin(0.7 * (t - 1)) + 0.5 * cos(2.3 * (t - 1) + 1) X1[, t] <- x U[, t] <- u_t x <- as.numeric(A0 %*% x + B0 * u_t) X2[, t] <- x } # Identify A and B jointly; both are recovered to machine precision fit <- dmdc(X1, X2, U, rank_input = 3) round(fit$a, 6) round(fit$b, 6) # Known B: pin the input matrix and estimate only A fit2 <- dmdc(X1, X2, U, rank_input = 2, known_B = B0) round(fit2$a, 6)
Per-mode frequency, growth rate and stability for the operator identified
by dmdc. DMDc has no mode amplitudes, so the
amplitude column is reported as 0.
dmdc_spectrum(object, dt = NULL)dmdc_spectrum(object, dt = NULL)
object |
A dmdc object. |
dt |
Time step. Uses the stored dt by default. |
Data frame with mode information.
A0 <- matrix(c(0.9, 0, 0.1, 0.8), 2, 2) B0 <- matrix(c(0.5, 1), 2, 1) m <- 60 X1 <- matrix(0, 2, m); X2 <- matrix(0, 2, m); U <- matrix(0, 1, m) x <- c(1, -0.5) for (t in seq_len(m)) { u_t <- sin(0.7 * (t - 1)) + 0.5 * cos(2.3 * (t - 1) + 1) X1[, t] <- x U[, t] <- u_t x <- as.numeric(A0 %*% x + B0 * u_t) X2[, t] <- x } fit <- dmdc(X1, X2, U, rank_input = 3) dmdc_spectrum(fit)A0 <- matrix(c(0.9, 0, 0.1, 0.8), 2, 2) B0 <- matrix(c(0.5, 1), 2, 1) m <- 60 X1 <- matrix(0, 2, m); X2 <- matrix(0, 2, m); U <- matrix(0, 1, m) x <- c(1, -0.5) for (t in seq_len(m)) { u_t <- sin(0.7 * (t - 1)) + 0.5 * cos(2.3 * (t - 1) + 1) X1[, t] <- x U[, t] <- u_t x <- as.numeric(A0 %*% x + B0 * u_t) X2[, t] <- x } fit <- dmdc(X1, X2, U, rank_input = 3) dmdc_spectrum(fit)
Classify the stability of the operator identified by dmdc
from the eigenvalues of .
dmdc_stability(object, tol = 1e-06)dmdc_stability(object, tol = 1e-06)
object |
A dmdc object. |
tol |
Tolerance for marginal classification. |
List with stability information (is_stable,
is_unstable, is_marginal, spectral_radius).
A0 <- matrix(c(0.9, 0, 0.1, 0.8), 2, 2) B0 <- matrix(c(0.5, 1), 2, 1) m <- 60 X1 <- matrix(0, 2, m); X2 <- matrix(0, 2, m); U <- matrix(0, 1, m) x <- c(1, -0.5) for (t in seq_len(m)) { u_t <- sin(0.7 * (t - 1)) + 0.5 * cos(2.3 * (t - 1) + 1) X1[, t] <- x U[, t] <- u_t x <- as.numeric(A0 %*% x + B0 * u_t) X2[, t] <- x } fit <- dmdc(X1, X2, U, rank_input = 3) dmdc_stability(fit)A0 <- matrix(c(0.9, 0, 0.1, 0.8), 2, 2) B0 <- matrix(c(0.5, 1), 2, 1) m <- 60 X1 <- matrix(0, 2, m); X2 <- matrix(0, 2, m); U <- matrix(0, 1, m) x <- c(1, -0.5) for (t in seq_len(m)) { u_t <- sin(0.7 * (t - 1)) + 0.5 * cos(2.3 * (t - 1) + 1) X1[, t] <- x U[, t] <- u_t x <- as.numeric(A0 %*% x + B0 * u_t) X2[, t] <- x } fit <- dmdc(X1, X2, U, rank_input = 3) dmdc_stability(fit)
Extended standard map (3D)
extended_standard_map(state, epsilon = 0.01, delta = 0.001)extended_standard_map(state, epsilon = 0.01, delta = 0.001)
state |
Numeric vector of length 3. |
epsilon |
Perturbation parameter. |
delta |
Coupling parameter. |
Updated state vector.
extended_standard_map(c(0.1, 0.2, 0.3))extended_standard_map(c(0.1, 0.2, 0.3))
Froeschle map (4D coupled standard maps)
froeschle_map(state, epsilon1 = 0.02, epsilon2 = 0.02, eta = 0.01)froeschle_map(state, epsilon1 = 0.02, epsilon2 = 0.02, eta = 0.01)
state |
Numeric vector of length 4. |
epsilon1 |
First perturbation. |
epsilon2 |
Second perturbation. |
eta |
Coupling parameter. |
Updated state vector.
froeschle_map(c(0.1, 0.2, 0.3, 0.4))froeschle_map(c(0.1, 0.2, 0.3, 0.4))
Generate a trajectory from a built-in map
generate_trajectory(map_name, initial_condition, n_iter, ...)generate_trajectory(map_name, initial_condition, n_iter, ...)
map_name |
Character: "standard", "froeschle", "extended_standard", "henon", or "logistic". |
initial_condition |
Numeric vector. |
n_iter |
Integer number of iterations. |
... |
Map parameters passed as named arguments. |
Numeric matrix (n_dim x n_iter+1).
# Chirikov standard map with default parameters traj <- generate_trajectory("standard", c(0.1, 0.2), 100) dim(traj) # Map parameters are passed through `...` henon <- generate_trajectory("henon", c(0, 0), 100, a = 1.4, b = 0.3) logistic <- generate_trajectory("logistic", 0.5, 100, r = 3.9)# Chirikov standard map with default parameters traj <- generate_trajectory("standard", c(0.1, 0.2), 100) dim(traj) # Map parameters are passed through `...` henon <- generate_trajectory("henon", c(0, 0), 100, a = 1.4, b = 0.3) logistic <- generate_trajectory("logistic", 0.5, 100, r = 3.9)
Generalized Laplace Analysis
gla(y, eigenvalues = NULL, n_eigenvalues = 5, tol = 1e-06, max_iter = NULL)gla(y, eigenvalues = NULL, n_eigenvalues = 5, tol = 1e-06, max_iter = NULL)
y |
Numeric matrix (n_obs x n_time). |
eigenvalues |
Complex vector of known eigenvalues, or NULL. |
n_eigenvalues |
Number of eigenvalues to estimate. |
tol |
Convergence tolerance. |
max_iter |
Maximum iterations, or NULL. |
An S3 object of class "gla".
t <- seq(0, 10, length.out = 200) y <- rbind(sin(t), cos(t)) g <- gla(y, n_eigenvalues = 2) print(g)t <- seq(0, 10, length.out = 200) y <- rbind(sin(t), cos(t)) g <- gla(y, n_eigenvalues = 2) print(g)
Reconstruct from GLA
gla_reconstruct(object, modes_to_use = NULL)gla_reconstruct(object, modes_to_use = NULL)
object |
A gla object. |
modes_to_use |
Integer vector of mode indices (1-based), or NULL. |
Numeric matrix.
t <- seq(0, 10, length.out = 200) y <- rbind(sin(t), cos(t)) g <- gla(y, n_eigenvalues = 2) recon <- gla_reconstruct(g) dim(recon)t <- seq(0, 10, length.out = 200) y <- rbind(sin(t), cos(t)) g <- gla(y, n_eigenvalues = 2) recon <- gla_reconstruct(g) dim(recon)
Hankel-DMD (Time-Delay Embedding DMD)
hankel_dmd(y, delays = NULL, rank = NULL, dt = 1)hankel_dmd(y, delays = NULL, rank = NULL, dt = 1)
y |
Numeric matrix (n_obs x n_time). |
delays |
Integer number of delays, or NULL for automatic. |
rank |
Integer truncation rank, or NULL for automatic. |
dt |
Numeric time step. |
An S3 object of class "hankel_dmd".
# A scalar signal, as a 1-row matrix t <- seq(0, 4 * pi, length.out = 200) y <- matrix(sin(t), nrow = 1) h <- hankel_dmd(y, delays = 10) print(h)# A scalar signal, as a 1-row matrix t <- seq(0, 4 * pi, length.out = 200) y <- matrix(sin(t), nrow = 1) h <- hankel_dmd(y, delays = 10) print(h)
Reconstruct from Hankel-DMD
hankel_reconstruct(object, n_steps)hankel_reconstruct(object, n_steps)
object |
A hankel_dmd object. |
n_steps |
Number of time steps. |
Numeric matrix.
t <- seq(0, 4 * pi, length.out = 200) y <- matrix(sin(t), nrow = 1) h <- hankel_dmd(y, delays = 10) recon <- hankel_reconstruct(h, 50)t <- seq(0, 4 * pi, length.out = 200) y <- matrix(sin(t), nrow = 1) h <- hankel_dmd(y, delays = 10) recon <- hankel_reconstruct(h, 50)
Harmonic Time Average
harmonic_time_average( map_name, initial_condition, observable = "sin_pi", omega = 0.1, n_iter = 10000, ... )harmonic_time_average( map_name, initial_condition, observable = "sin_pi", omega = 0.1, n_iter = 10000, ... )
map_name |
Character map name. |
initial_condition |
Numeric vector. |
observable |
Character observable name. |
omega |
Numeric frequency. |
n_iter |
Integer iterations. |
... |
Map parameters. |
List with magnitude, phase, hta_re, hta_im.
# Harmonic time average of an orbit of the Chirikov standard map harmonic_time_average("standard", c(0.1, 0.2), "sin_pi", 0.1, 500)# Harmonic time average of an orbit of the Chirikov standard map harmonic_time_average("standard", c(0.1, 0.2), "sin_pi", 0.1, 500)
Henon map (2D dissipative)
henon_map(state, a = 1.4, b = 0.3)henon_map(state, a = 1.4, b = 0.3)
state |
Numeric vector c(x, y). |
a |
Parameter a. |
b |
Parameter b. |
Updated state vector.
henon_map(c(0, 0)) # Classic chaotic parameters henon_map(c(0, 0), a = 1.4, b = 0.3)henon_map(c(0, 0)) # Classic chaotic parameters henon_map(c(0, 0), a = 1.4, b = 0.3)
HTA convergence analysis
hta_convergence( map_name, initial_condition, observable = "sin_pi", omega = 0.1, n_iter = 10000, ... )hta_convergence( map_name, initial_condition, observable = "sin_pi", omega = 0.1, n_iter = 10000, ... )
map_name |
Character map name. |
initial_condition |
Numeric vector. |
observable |
Character observable name. |
omega |
Numeric frequency. |
n_iter |
Integer iterations. |
... |
Map parameters. |
List with times, hta_magnitudes, convergence_rate, dynamics_type.
# How the time average converges along the orbit conv <- hta_convergence("standard", c(0.1, 0.2), "sin_pi", 0.1, 500) str(conv)# How the time average converges along the orbit conv <- hta_convergence("standard", c(0.1, 0.2), "sin_pi", 0.1, 500) str(conv)
Logistic map (1D)
logistic_map(state, r = 3.9)logistic_map(state, r = 3.9)
state |
Numeric scalar. |
r |
Growth rate parameter. |
Updated state value.
logistic_map(0.5) # In the chaotic regime logistic_map(0.5, r = 3.9)logistic_map(0.5) # In the chaotic regime logistic_map(0.5, r = 3.9)
Mesochronic harmonic plot computation
mesochronic_compute( map_name, x_range = c(0, 1), y_range = c(0, 1), resolution = 100, observable = "sin_pi", omega = 0.1, n_iter = 30000, ... )mesochronic_compute( map_name, x_range = c(0, 1), y_range = c(0, 1), resolution = 100, observable = "sin_pi", omega = 0.1, n_iter = 30000, ... )
map_name |
Character map name. |
x_range |
Numeric vector c(min, max). |
y_range |
Numeric vector c(min, max). |
resolution |
Integer grid resolution. |
observable |
Character observable name. |
omega |
Numeric frequency. |
n_iter |
Integer iterations. |
... |
Map parameters. |
List with hta_matrix, phase_matrix, x_coords, y_coords.
# Mesochronic plot over a coarse grid. Raise `resolution` and `n_iter` # for publication-quality figures; both cost time roughly linearly. mhp <- mesochronic_compute("standard", c(0, 1), c(0, 1), 10, "sin_pi", 0.1, 100) str(mhp)# Mesochronic plot over a coarse grid. Raise `resolution` and `n_iter` # for publication-quality figures; both cost time roughly linearly. mhp <- mesochronic_compute("standard", c(0, 1), c(0, 1), 10, "sin_pi", 0.1, 100) str(mhp)
Predict from DMD model
## S3 method for class 'dmd' predict(object, n_ahead = 10, x0 = NULL, method = c("modes", "matrix"), ...)## S3 method for class 'dmd' predict(object, n_ahead = 10, x0 = NULL, method = c("modes", "matrix"), ...)
object |
A dmd object. |
n_ahead |
Number of steps to predict. |
x0 |
Optional initial condition vector. |
method |
"modes" (default) or "matrix". |
... |
Additional arguments (ignored). |
Numeric matrix of predictions.
t <- seq(0, 10, length.out = 100) X <- rbind(sin(t), cos(t)) d <- dmd(X, rank = 2, dt = t[2] - t[1]) # Forecast 10 steps beyond the input pred <- predict(d, n_ahead = 10) dim(pred)t <- seq(0, 10, length.out = 100) X <- rbind(sin(t), cos(t)) d <- dmd(X, rank = 2, dt = t[2] - t[1]) # Forecast 10 steps beyond the input pred <- predict(d, n_ahead = 10) dim(pred)
Simulate the identified system forward from
an initial state under a given control input sequence.
## S3 method for class 'dmdc' predict(object, U = NULL, x0 = NULL, n_ahead = NULL, ...)## S3 method for class 'dmdc' predict(object, U = NULL, x0 = NULL, n_ahead = NULL, ...)
object |
A dmdc object. |
U |
Numeric matrix of control inputs (n_inputs x n_steps); the number
of columns sets the prediction horizon. A vector is taken as a single
input row. NULL applies zero input for |
x0 |
Initial state vector. Defaults to the first stored snapshot. |
n_ahead |
Number of steps when |
... |
Additional arguments (ignored). |
Numeric matrix of predicted states
(n_states x k).
A0 <- matrix(c(0.9, 0, 0.1, 0.8), 2, 2) B0 <- matrix(c(0.5, 1), 2, 1) m <- 120 X1 <- matrix(0, 2, m) X2 <- matrix(0, 2, m) U <- matrix(0, 1, m) x <- c(1, -0.5) for (t in seq_len(m)) { u_t <- sin(0.7 * (t - 1)) + 0.5 * cos(2.3 * (t - 1) + 1) X1[, t] <- x U[, t] <- u_t x <- as.numeric(A0 %*% x + B0 * u_t) X2[, t] <- x } fit <- dmdc(X1, X2, U, rank_input = 3) # Replaying the training inputs reproduces the observed successors pred <- predict(fit, U = U) max(abs(pred - X2)) # Zero-input (free) response from a chosen state free <- predict(fit, x0 = c(1, 1), n_ahead = 10) dim(free)A0 <- matrix(c(0.9, 0, 0.1, 0.8), 2, 2) B0 <- matrix(c(0.5, 1), 2, 1) m <- 120 X1 <- matrix(0, 2, m) X2 <- matrix(0, 2, m) U <- matrix(0, 1, m) x <- c(1, -0.5) for (t in seq_len(m)) { u_t <- sin(0.7 * (t - 1)) + 0.5 * cos(2.3 * (t - 1) + 1) X1[, t] <- x U[, t] <- u_t x <- as.numeric(A0 %*% x + B0 * u_t) X2[, t] <- x } fit <- dmdc(X1, X2, U, rank_input = 3) # Replaying the training inputs reproduces the observed successors pred <- predict(fit, U = U) max(abs(pred - X2)) # Zero-input (free) response from a chosen state free <- predict(fit, x0 = c(1, 1), n_ahead = 10) dim(free)
Predict from GLA
## S3 method for class 'gla' predict(object, n_ahead = 10, ...)## S3 method for class 'gla' predict(object, n_ahead = 10, ...)
object |
A gla object. |
n_ahead |
Number of steps to predict. |
... |
Additional arguments (ignored). |
Numeric matrix.
t <- seq(0, 10, length.out = 200) y <- rbind(sin(t), cos(t)) g <- gla(y, n_eigenvalues = 2) pred <- predict(g, n_ahead = 5)t <- seq(0, 10, length.out = 200) y <- rbind(sin(t), cos(t)) g <- gla(y, n_eigenvalues = 2) pred <- predict(g, n_ahead = 5)
Predict from Hankel-DMD
## S3 method for class 'hankel_dmd' predict(object, n_ahead = 10, ...)## S3 method for class 'hankel_dmd' predict(object, n_ahead = 10, ...)
object |
A hankel_dmd object. |
n_ahead |
Number of steps to predict. |
... |
Additional arguments (ignored). |
Numeric matrix.
t <- seq(0, 4 * pi, length.out = 200) y <- matrix(sin(t), nrow = 1) h <- hankel_dmd(y, delays = 10) pred <- predict(h, n_ahead = 10)t <- seq(0, 4 * pi, length.out = 200) y <- matrix(sin(t), nrow = 1) h <- hankel_dmd(y, delays = 10) pred <- predict(h, n_ahead = 10)
Standard map (Chirikov)
standard_map(state, epsilon = 0.12)standard_map(state, epsilon = 0.12)
state |
Numeric vector c(x, y). |
epsilon |
Perturbation parameter. |
Updated state vector.
# One iteration from a given state standard_map(c(0.1, 0.2)) # Stronger nonlinearity standard_map(c(0.1, 0.2), epsilon = 0.5)# One iteration from a given state standard_map(c(0.1, 0.2)) # Stronger nonlinearity standard_map(c(0.1, 0.2), epsilon = 0.5)