# Nonlinear Time Continuous System

Similarly, we can use the Extended Dynamic Mode Decomposition via a nonlinear Basis of observables. Here, we will look at a rather famous example with a finite dimensional solution.

using DataDrivenDiffEq
using LinearAlgebra
using ModelingToolkit
using OrdinaryDiffEq
using Plots

function slow_manifold(du, u, p, t)
du = p * u
du = p * (u - u^2)
end

u0 = [3.0; -2.0]
tspan = (0.0, 5.0)
p = [-0.8; -0.7]

problem = ODEProblem(slow_manifold, u0, tspan, p)
solution = solve(problem, Tsit5(), saveat = 0.01);
plot(solution)

Since we are dealing with a continuous system in time, we define the associated DataDrivenProblem accordingly using the measured states X, their derivatives DX and the time t.

prob = ContinuousDataDrivenProblem(solution)
Continuous DataDrivenProblem{Float64} ##DDProblem#466 in 2 dimensions and 501 samples

Additionally, we need to define the Basis for our lifting, before we solve the problem in the lifted space.

@parameters t
@variables u[1:2](t)
Ψ = Basis([u; u^2], u, independent_variable = t)
res = solve(prob, Ψ, DMDPINV(), digits = 1)
system = result(res)
Linear Solution with 2 equations and 3 parameters.
Returncode: solved
L₂ Norm error : [0.0, 1.0114581578656041e-29]
AIC : [-Inf, -33442.5509179254]
R² : [1.0, 1.0]

Model ##Koopman#475 with 2 equations
States : u(t) u(t)
Parameters : p₁ p₂ p₃
Independent variable: t
Equations
Differential(t)(u(t)) = p₁*u(t)
Differential(t)(u(t)) = p₂*u(t) + p₃*(u(t)^2)
[-0.8, -0.7, 0.7]

The underlying dynamics have been recovered correctly by the algorithm! Similarly we could use sparse identification to solve the problem

sparse_res = solve(prob, Ψ, STLSQ(), digits = 1)
println(sparse_res)
Linear Solution with 2 equations and 3 parameters.
Returncode: solved
L₂ Norm error : [0.0, 1.0114581578656041e-29]
AIC : [-Inf, -33442.5509179254]
R² : [1.0, 1.0]

And the resulting system

sparse_system = result(sparse_res)
println(sparse_system)
Model ##Basis#478 with 2 equations
States : u(t) u(t)
Parameters : p₁ p₂ p₃
Independent variable: t
Equations
Differential(t)(u(t)) = p₁*u(t)
Differential(t)(u(t)) = p₂*u(t) + p₃*(u(t)^2)

We can also directly look at the parameters of each result

parameter_map(res)
3-element Vector{Pair{SymbolicUtils.Sym{Real, Base.ImmutableDict{DataType, Any}}, Float64}}:
p₁ => -0.8
p₂ => -0.7
p₃ => 0.7

Note that we are using parameter_map instead of just parameters, which returns a vector suitable to use with ModelingToolkit.

parameter_map(sparse_res)
3-element Vector{Pair{SymbolicUtils.Sym{Real, Base.ImmutableDict{DataType, Any}}, Float64}}:
p₁ => -0.8
p₂ => -0.7
p₃ => 0.7

To simulate the system, we create an ODESystem from the result

Both results can be converted into an ODESystem

@named sys = ODESystem(
equations(sparse_system),
get_iv(sparse_system),
states(sparse_system),
parameters(sparse_system)
);

x0 = [u => u0, u => u0]
ps = parameter_map(sparse_res)
3-element Vector{Pair{SymbolicUtils.Sym{Real, Base.ImmutableDict{DataType, Any}}, Float64}}:
p₁ => -0.8
p₂ => -0.7
p₃ => 0.7

And simulated using OrdinaryDiffEq.jl using the (known) initial conditions and the parameter mapping of the estimation.

ode_prob = ODEProblem(sys, x0, tspan, ps)
estimate = solve(ode_prob, Tsit5(), saveat = prob.t);

And look at the result

plot(solution, color = :black)
plot!(estimate, color = :red, linestyle = :dash)

## Copy-Pasteable Code

using DataDrivenDiffEq
using LinearAlgebra
using ModelingToolkit
using OrdinaryDiffEq

function slow_manifold(du, u, p, t)
du = p * u
du = p * (u - u^2)
end

u0 = [3.0; -2.0]
tspan = (0.0, 5.0)
p = [-0.8; -0.7]

problem = ODEProblem(slow_manifold, u0, tspan, p)
solution = solve(problem, Tsit5(), saveat = 0.01);

@parameters t
@variables u[1:2](t)
Ψ = Basis([u; u^2], u, independent_variable = t)
res = solve(prob, Ψ, DMDPINV(), digits = 1)
system = result(res)

sparse_res = solve(prob, Ψ, STLSQ(), digits = 1)

sparse_system = result(sparse_res)

parameter_map(res)

parameter_map(sparse_res)

@named sys = ODESystem(
equations(sparse_system),
get_iv(sparse_system),
states(sparse_system),
parameters(sparse_system)
);

x0 = [u => u0, u => u0]
ps = parameter_map(sparse_res)

ode_prob = ODEProblem(sys, x0, tspan, ps)
estimate = solve(ode_prob, Tsit5(), saveat = prob.t);

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