Symbolic full-order model

General usage

As in Building a full-order model, MORFE considers $\mathrm{ORD}$-th order ODEs of the form

ExternalSystem

The right-hand side $\mathbf{E}(\mathbf{r})$ of an ExternalSystem must be a polynomial function in $\mathbf{r}$ and is not allowed to be time-dependent. Zero must also be an equilibrium point, so $\mathbf{E}(\mathbf{0})=\mathbf{0}$ has to hold. There are two ways to define one symbolically.

The first defines $\mathbf{E}(\mathbf{r})$ as an expression:

# R external states, autonomous dynamics ṙ = E(r)
@variables r[1:R]
r = collect(r)

ext_exprs = [ # E(r) — polynomial, autonomous
    #   ...
]

external_1 = externalsystem_from_symbolics(ext_exprs, r)

The second is a function shaped like in DifferentialEquations.jl, either a returning function g or in-placeg!:

function g(r, params, t)
    dr = similar(r)
    dr[1] = # ...
    # ...
    dr[R] = # ...
    return dr
end
external_2 = externalsystem_from_symbolics(g, R; p = params)


function g!(dr, r, params, t)
    dr[1] = # ...
    # ...
    dr[R] = # ...
end
external_3 = externalsystem_from_symbolics(g!, R; p = params)

Full model

A code, for a system of order $N$ in $n$ degrees of freedom, coupled to $R$ external states, would look like this:

# N+1 derivative groups (u, u̇, ü, …), n dof each — extend the list for N > 2
@variables u[1:n] du[1:n] ddu[1:n]
u, du, ddu = collect(u), collect(du), collect(ddu)
groups = (u, du, ddu)   # (u, u̇, ü, …) — ascending order

# R external states, autonomous dynamics ṙ = E(r)
@variables r[1:R]
r = collect(r)

exprs = [ # B₀u + B₁u̇ + B₂ü + ⋯ − Σⱼ Tⱼ(u,…,u,u̇,…,u̇,…,r,…,r) , written as "= 0"
    #   ...
]
ext_exprs = [ # E(r) — polynomial, autonomous
    #   ...
]

model = model_from_symbolics(exprs, groups, r, ext_exprs)

Variables third, equations fourth. Drop both for an unforced system — model_from_symbolics(exprs, groups).

Both constructors also take a DifferentialEquations.jl-shaped function, with no @variables to declare. model_from_symbolics needs an in-place f! whose arguments run from the highest derivative down to $\mathbf u$:

function f!(dᴺu, #=…=# ddu, du, u, params, t)
    dᴺu[1] = # ...
    # ...
end
model = model_from_symbolics(f!, N, n; p = params)

In the coupled form, f! takes the external state as one extra argument, after $\mathbf u$ and before p. Referencing r there is how forcing enters the model:

function f!(dᴺu, #=…=# ddu, du, u, r, params, t)
    dᴺu[1] = # ... may read r ...
    # ...
end
function g!(dr, r, params_ext, t)
    dr[1] = # ...
    # ...
end

model = model_from_symbolics(f!, N, n, g!, R; p = params, p_ext = params_ext)

Two-mass oscillator

The two-degree-of-freedom nonlinear oscillator introduced by Shaw and Pierre: two coupled masses with linear stiffness and damping, the first carrying an additional cubic restoring force.

Without forcing

Drop the $-2\cos(\Omega t)$. Declare one Symbolics.jl vector per derivative order, write the equations of motion as expressions that equal zero, and hand both to model_from_symbolics.

@variables u[1:2] du[1:2] ddu[1:2]
u, du, ddu = collect(u), collect(du), collect(ddu)

k, g, c = 1.0, 6.0, 0.1

exprs = [
    ddu[1] + c*du[1] - c*du[2] + (1 + k)*u[1] - k*u[2] + g*u[1]^3,
    ddu[2] - c*du[1] + 2*c*du[2] - k*u[1] + (1 + k)*u[2],
]

groups = (u, du, ddu)   # ascending order: (u, u̇, ü)

model = model_from_symbolics(exprs, groups)

The last entry of groups is the derivative solved for. What came out:

With forcing

As in Building a full-order model, substitute $-2\cos(\Omega t)$ by $-r_1-r_2$ with $\dot r_1 = i\Omega r_1$ and $\dot r_2 = -i\Omega r_2$, so that $r_1 + r_2 = 2\cos\Omega t$ when both start at 1, and add that ODE as an external system.

@variables u[1:2] du[1:2] ddu[1:2] r[1:2]
u, du, ddu, r = collect(u), collect(du), collect(ddu), collect(r)

k, g, c = 1.0, 6.0, 0.1
Ω = 1.3

exprs = [
    ddu[1] + c*du[1] - c*du[2] + (1 + k)*u[1] - k*u[2] + g*u[1]^3 - r[1] - r[2],
    ddu[2] - c*du[1] + 2*c*du[2] - k*u[1] + (1 + k)*u[2],
]

ext_exprs = [
    im*Ω*r[1],
    -im*Ω*r[2],
]

groups = (u, du, ddu)   # ascending order: (u, u̇, ü)
model = model_from_symbolics(exprs, groups, r, ext_exprs)

Harmonic excitation

In practice the full-order model comes from a FEM code, not from symbolic expressions. The symbolic interface earns its keep on the ExternalSystem, where forcing terms and parameters enter.

In Building a full-order model the harmonic driver $\dot r_1 = i\Omega r_1,\ \dot r_2 = -i\Omega r_2$ is directly defined from its eigenvalues: ExternalSystem((im*Ω, -im*Ω)). Written symbolically instead, the same driver is just its own right-hand side:

@variables r[1:2]
r = collect(r)
Ω = 1.3

ext_exprs = [im*Ω*r[1], -im*Ω*r[2]]
harmonic = externalsystem_from_symbolics(ext_exprs, r)   # same object as ExternalSystem((im*Ω, -im*Ω))

Unlike the eigenvalue-tuple constructor, ext_exprs need not be linear — any polynomial right-hand side works, as Step 05 shows.

Quasi-periodic and multiharmonic drivers

A pair of incommensurate frequencies closes onto a torus rather than a circle — still linear, still diagonal, just twice as many states. Two independent harmonic pairs, written the same way as Step 03:

@variables r[1:4]
r = collect(r)
Ω1, Ω2 = 1.0, sqrt(2)   # incommensurate

ext_exprs = [
    im*Ω1*r[1], -im*Ω1*r[2],
    im*Ω2*r[3], -im*Ω2*r[4],
]
quasi = externalsystem_from_symbolics(ext_exprs, r)   # same object as ExternalSystem((im*Ω1,-im*Ω1,im*Ω2,-im*Ω2))

The eigenvalue-tuple constructor only reaches diagonal systems, though. Suppose the two signals you actually want are $r_1 = -\cos\Omega t$ and $r_2 = 0.03\sin\Omega t - \cos 4\Omega t$ — not diagonal in the obvious basis. Introducing auxiliary states $r_3=\sin\Omega t$, $r_4=\sin 4\Omega t$ to close this as a first-order linear system gives $\dot{\mathbf r} = \mathbf A\mathbf r$ with (see Building a full-order model, Step 02 for the derivation):

Ω = 1.3
@variables r[1:4]
r = collect(r)

ext_exprs = [
    Ω*r[3],
    -0.03*Ω*r[1] + 4*Ω*r[4],
    -Ω*r[1],
    -4*Ω*r[2] + 0.12*Ω*r[3],
]
multiharmonic = externalsystem_from_symbolics(ext_exprs, r)

This matrix is not upper triangular, so MORFE re-bases it and reports the change of coordinates (see Building a full-order model, Step 02). external_basis returns the $\mathbf Q$ relating the physical $\mathbf r$ to the stored $\mathbf r'$, and to_physical_external maps a state back.

Chaotic nonlinear excitation

A nonlinear driver works the same way — as long as it satisfies the same conditions as everything else on this page: polynomial, with an equilibrium at the origin. The Lorenz system has three equilibria, none at the origin — so Building a full-order model, Step 03 shifts coordinates to the nontrivial equilibrium $C_+=(C,C,\rho-1)$, $C=\sqrt{\beta(\rho-1)}$, before defining it:

That shift is just the coordinates you declare — and nothing checks it for you here, as the warning above explains.

σ, ρ, β = 10.0, 28.0, 8/3
C = sqrt(β*(ρ-1))

@variables X Y Z
r = [X, Y, Z]

ext_exprs = [
    σ*(Y-X),
    X - Y - Z*(X+C),
    C*(X+Y) + X*Y - β*Z,
]
lorenz = externalsystem_from_symbolics(ext_exprs, r)

Using the function layout

The same objects can be built from plain functions, with no @variables to declare. For model_from_symbolics the function is in-place and takes the derivatives from highest to lowest; p is passed straight through, so parameters can be supplied rather than closed over.

function two_mass!(a, v, x, p, t)          # a = ü, v = u̇, x = u
    kk, gg, cc = p
    a[1] = -cc*v[1] + cc*v[2] - (1 + kk)*x[1] + kk*x[2] - gg*x[1]^3
    a[2] = cc*v[1] - 2*cc*v[2] + kk*x[1] - (1 + kk)*x[2]
end

# same B₀, B₁, B₂ and the same (3, 0) cubic as Step 01
model = model_from_symbolics(two_mass!, 2, 2; p = (k, g, c))

rotation!(dr, r, p, t) = (dr[1] = p[1]*r[2]; dr[2] = -p[1]*r[1])
rotation = externalsystem_from_symbolics(rotation!, 2; p = (Ω,))

Everything above only builds objects. The example integrates a damped oscillator $\ddot y + c\,\dot y + k\,y = \mathrm{forcing}(\mathbf r)$ together with its driver, taking the driver's right-hand side from the ExternalSystem through evaluate.