Linear Time Discrete System

We will start by estimating the underlying dynamical system of a time discrete process based on some measurements via Dynamic Mode Decomposition on a simple linear system of the form $u(k+1) = A u(k)$.

At first, we simulate the corresponding system using OrdinaryDiffEq.jl and generate a DiscreteDataDrivenProblem from the simulated data.

using DataDrivenDiffEq
using LinearAlgebra
using OrdinaryDiffEq
using DataDrivenDMD
using Plots

A = [0.9 -0.2; 0.0 0.2]
u0 = [10.0; -10.0]
tspan = (0.0, 11.0)

f(u, p, t) = A * u

sys = DiscreteProblem(f, u0, tspan)
sol = solve(sys, FunctionMap());

Next, we transform our simulated solution into a DataDrivenProblem. Given that the solution knows it is a discrete solution, we can simply write

prob = DataDrivenProblem(sol)

Discrete DataDrivenProblem{Float64} ##DDProblem#264 in 2 dimensions and 12 samples

And plot the solution and the problem

plot(sol, label = string.([:x₁ :x₂]))
scatter!(prob)

To estimate the underlying operator in the states $x_1, x_2$, we solve the estimation problem using the DMDSVD algorithm for approximating the operator. First, we will have a look at the DataDrivenSolution

res = solve(prob, DMDSVD(), digits = 1)

"DataDrivenSolution{Float64}"

We see that the system has been recovered correctly, indicated by the small error and high AIC score of the result. We can confirm this by looking at the resulting Basis

get_basis(res)

\[ \begin{align} Difference(t; dt=1.0, update=false)\left( \mathtt{x{_1}} \right) &= \mathtt{p{_1}} \mathtt{x{_1}} + \mathtt{p{_2}} \mathtt{x{_2}} \\ Difference(t; dt=1.0, update=false)\left( \mathtt{x{_2}} \right) &= \mathtt{p{_3}} \mathtt{x{_2}} \end{align} \]

And also plot the prediction of the recovered dynamics

plot(res)

Copy-Pasteable Code

using DataDrivenDiffEq
using LinearAlgebra
using OrdinaryDiffEq
using DataDrivenDMD

A = [0.9 -0.2; 0.0 0.2]
u0 = [10.0; -10.0]
tspan = (0.0, 11.0)

f(u, p, t) = A * u

sys = DiscreteProblem(f, u0, tspan)
sol = solve(sys, FunctionMap());

prob = DataDrivenProblem(sol)

res = solve(prob, DMDSVD(), digits = 1)

get_basis(res)

# This file was generated using Literate.jl, https://github.com/fredrikekre/Literate.jl

This page was generated using Literate.jl.