MORFE.jl / Tutorials / Building a full-order model
CORE LIBRARY < 15 S · NO FEM

Building a full-order model

Open the notebook ↗

In MORFE.jl full-order models have 3 components: the matrices of the linear terms, the nonlinear terms, and (optionally) a external system. This tutorial builds all three from scratch, draws what each one means and builds a model that holds them together.

FIG · 01 A MultilinearMap is one term of the right-hand side of the full-order dynamics. All smooth curves here are produced by evaluate_term! on a concrete MultilinearMap. The third panel is a term mixing state and external state, $F = k\,u\,(r_1+r_2)$ with $r_1+r_2 = 2\cos\Omega t$. Drag it to orbit.
Reproduce this page

Every line on this page comes from building_a_full-order_model. It needs no FEM backend and runs in a few seconds:

julia --project -e 'include("examples/internals/full_order_model/main.jl")'

What this tutorial covers

Mathematical formulation

MORFE considers $\mathrm{ORD}$-th order ODEs of the form

$$ \mathbf{B}_0\,\mathbf{u} + \mathbf{B}_1\,\dot{\mathbf{u}} + \mathbf{B}_{2} \, \ddot{\mathbf{u}} + \cdots = \sum_j \mathbf{T}_{j} (\underbrace{\mathbf{u}, \dotsc, \mathbf{u}}_{m_{j,0}}, \underbrace{\dot{\mathbf{u}}, \dotsc, \dot{\mathbf{u}}}_{m_{j,1}}, \dotsc, \underbrace{\mathbf{r}, \dotsc, \mathbf{r}}_{\widehat{m}_{j}}) \\[10pt] \dot{\mathbf{r}} = \mathbf{E}(\mathbf{r})$$

The matrices $\mathbf{B}_j$ describe the internal linear part; and the $\mathbf{T}_j$ express the remaining nonlinear and external terms as MultilinearMaps. The external state $\mathbf{r}$ satisfies its own autonomous dynamics $\dot{\mathbf{r}} = \mathbf{E}(\mathbf{r})$ and enters $\mathbf{T}_j$ to express non-autonomous terms. The first $m_{j,0}$ arguments of $\mathbf{T}_j$ have the zeroth derivative $\mathbf{u}$, the following $m_{j,1}$ arguments have the first derivative $\dot{\mathbf{u}}$, and so on. This defines the multiindex $(m_{j,0}, m_{j,1}, \dotsc, m_{j, \mathrm{ORD}-1})$. The external multiplicity $\widehat{m}_j$ counts how many copies of the external state $\mathbf{r}$ enter the term $\mathbf{T}_j$.

01

Multilinear terms

The $k$-th entry of multiindex defines $m_{j,k-1}$ which counts how many arguments take the derivative $u^{(k-1)}$. The length of a multiindex is the order of the system. So for a second-order model, (3, 0) means f!(res, u, u, u) and (1, 1) means f!(res, u, u̇). 

using MORFE

# M ü + C u̇ + K u = F(u, u̇, r), so a *hardening* Duffing spring contributes F = −k₃ u³
hard!(res, u1, u2, u3) = (res .+= -6.0 .* u1 .* u2 .* u3)
drag!(res, v1, v2) = (res .+= -0.8 .* v1 .* v2)
mixed!(res, u, r) = (res .+= -1.5 .* u .* (r[1] + r[2]))

term_hard = MultilinearMap(duff!; multiindex = (3, 0), multiplicity_external = 0, fully_asymmetric = false)
term_drag = MultilinearMap(drag!; multiindex = (0, 2), multiplicity_external = 0)
term_mixed = MultilinearMap(mixed!; multiindex = (1, 0), multiplicity_external = 1)
multiindexthe arguments aredegree
(2, 0)$\mathbf{u}$, $\mathbf{u}$2
(3, 0)$\mathbf{u}$, $\mathbf{u}$, $\mathbf{u}$3
(0, 2)$\dot{\mathbf{u}}$, $\dot{\mathbf{u}}$2
(1, 1)$\mathbf{u}$, $\dot{\mathbf{u}}$2
(1, 0), 1$\mathbf{u}$, $\mathbf{r}$2
(0, 0), 1$\mathbf{r}$1

When keyword arguments are missing, the constructor any the assumed values.

The parametrisation method requires smooth maps

MORFE assumes MultilinearMap is linear in each argument; but cannot explicitly enforce this property. The user is responsible for ensuring multilinearity. For instance, MORFE cannot handle non-smooth functions like $F(u̇) = |u̇| u̇$ of Figure 1, panel 2.

02

Periodic excitation

An ExternalSystem is an autonomous ODE $\dot{\mathbf{r}} = \mathbf{E}(\mathbf{r})$ whose state enters the model's multilinear terms as an extra argument. The classical case of harmonic forcing considers $\dot{r}_{1} = + i\Omega \, r_{1}$ and $\dot{r}_{2} = - i\Omega \, r_{2}$ with solutions $r_{1}(t) = r_{1,0} \, e^{ i\Omega t}$ and $r_{2}(t) = r_{2,0} \, e^{- i\Omega t}$. Fixing $r_{1,0} = r_{2,0} = 1$ and adding both states yields $r_1 + r_2 = 2 \cos(\Omega t)$. 

harmonic = ExternalSystem((im*Ω, -im*Ω)) # a circle
quasi = ExternalSystem((im*Ω₁, -im*Ω₁, im*Ω₂, -im*Ω₂)) # the closure is a torus

Passing a tuple of eigenvalues builds a diagonal driver. But excitation can be more general. Suppose the two signals you actually want are

$$r_1 = -\cos\Omega t, \qquad r_2 = 0.03 \sin\Omega t - \cos 4\Omega t \, .$$

We define auxiliary states, $r_3 = \sin\Omega t$ and $r_4 = \sin 4\Omega t$, to close this as a first order linear system. Differentiating and writing sines and cosines back in terms of the state yields

$$\dot r_1 = \Omega r_3, \qquad \dot r_2 = -0.03\Omega r_1 + 4\Omega r_4, \qquad \dot r_3 = -\Omega r_1, \qquad \dot r_4 = -4\Omega r_2 + 0.12\Omega r_3 \, .$$
In code, we implement this linear ODE, $\dot{\mathbf{r}} = \mathbf{A} \mathbf{r}$, as follows: 
A = [ 0.0     0.0   Ω      0.0 ;
     -0.03Ω  0.0   0.0    4Ω  ;
     -Ω       0.0   0.0    0.0 ;
      0.0    -4Ω    0.12Ω  0.0 ]

multiharmonic = ExternalSystem(linear_polynomial(A))

When applicable, MORFE performs a change of coordinates $\mathbf{r} = \mathbf{Q}\widehat{\mathbf{z}}$, which re-expresses the external system as $\dot{\widehat{\mathbf{z}}} = \mathbf{U}\widehat{\mathbf{z}} + \mathcal{O}(\|\widehat{\mathbf{z}}\|^2)$ with $\mathbf{U}$ upper triangular, and says so: 

Info: ExternalSystem: the linear matrix was not upper triangular — non-zero entries below
the diagonal at [(2, 1), (3, 1), (4, 2), (4, 3)].  It has been re-based rather than rejected:
  · route: eigenvector basis (A is real and diagonalisable).  U is diagonal,
    and the conjugate pairing of the external variables is preserved exactly.
  · the external coordinates are now r′, related to the physical r by r = Q r′;
    `external_basis(sys)` returns Q.

Nothing else has to change: nonlinear terms are still fed the physical external argument, and to_physical_external maps a state back whenever you want to read one.

In every case the dynamics are stored as a polynomial, so solving one needs nothing more than evaluate and a standard solver; e.g. Runge–Kutta for time integration.

FIG · 02 Three drivers — switch tabs to compare. Time plot on the left, phase portrait on the right. A purely imaginary pair gives a circle. The multiharmonic driver draws an M. Two incommensurate frequencies give a quasi-periodic orbit that never closes.
03

Chaotic excitation

Now we treat an external signal with nonlinear dynamics by considering the Lorenz system. MORFE demands the origin as an equilibrium; therefore we must centre the polynomial expansion at an equilibrium. The Lorenz system possesses three equilibria: the trivial origin and the nontrivial pair $C_\pm = (\pm C,\, \pm C,\, \rho-1)$, with $C = \sqrt{\beta(\rho-1)}$. The chaotic attractor, the famous butterfly, orbits these nontrivial equilibria. We choose $C_+$ and shift the origin there: $(X,Y,Z) = (x,y,z) - C_+$. The translation of coordinates transforms the system as follows:

$$\begin{aligned} \begin{cases} \dot x = \sigma (y - x) \\ \dot y = x(\rho - z) - y \\ \dot z = xy - \beta z \end{cases} \quad \iff \quad \begin{cases} \dot X = \sigma (Y - X) \\ \dot Y = X - Y - Z(X + C) \\ \dot Z = C(X + Y) + XY - \beta Z \end{cases} \end{aligned}$$

The system in $(X,Y,Z)$ has an equilibrium at the origin, satisfying the requirement.

FIG · 03 Integrating the ExternalSystem and reading back with to_physical_external yields the famous butterfly of the Lorenz system. Left, the trajectory in phase space — drag to orbit; right, each coordinate against time. The arrow is the translation to the equilibrium point $C_+$.
Complex coordinates carry a reality condition

MORFE employs complex coordinates, but the physical flow exists in $\mathbb{R}^3$. Linearisation yields one real eigenvalue and a complex-conjugate pair. Hence, the real eigenvalue's coordinate stays real, and the conjugate pair's coordinates remain complex conjugates of each other. These simple constraints guarantee the imaginary parts cancel out in physical coordinates. In practice, time integration should enforce this at each time step.

04

Assemble the model

An NthOrderModel takes the linear matrices as a tuple $(B_0, \ldots, B_{\mathrm{ORD}})$, the nonlinear terms as a tuple, and the external system. A forced Duffing oscillator is two terms — the hardening cubic from step 1 and a forcing term that reads the external state. 

forcing!(res, r) = (res .+= 2.5 * (r[1] + r[2])) # 2.5 (r₁+r₂) = 5 cos(Ωt)
term_forcing = MultilinearMap(forcing!; multiindex = (0, 0), multiplicity_external = 1)

model = NthOrderModel((K, C, M), (term_hard, term_forcing), harmonic)

A, B = linear_first_order_matrices(model)   # companion pair for the eigensolver

A term that reads the external state needs a model that has one. Leave the ExternalSystem out and construction fails immediately, naming the term — rather than surviving to fail much later, during evaluation.

evaluate_nonlinear_terms! assembles $F$ at a given state, one degree at a time, which is all the right-hand side of a time integration needs. Driving the model with the harmonic system of step 2 gives the forced response:

FIG · 04 Damping pulls the orbit onto the periodic response — a closed loop, which the hardening cubic bends away from an ellipse. Twenty-five forcing periods are integrated; purple is the transient, red the final period. The nonlinear force at every step comes from evaluate_nonlinear_terms! on the model itself.

Recap

  • A multiindex is a calling convention: multiindex[k] slots take $x^{(k-1)}$, and its length is the system order.
  • Multilinear means linear in each slot separately — which rules out non-smooth terms.
  • An ExternalSystem may be nonlinear.
  • Internally, MORFE changes the coordinates of the ExternalSystem to make its linear part upper triangular. You can read states back with to_physical_external.
  • You must centre the polynomial expansion at an equilibrium, as it cannot have a constant term.

Where to next?