Construct invariant manifolds such as centre
manifolds, spectral submanifolds, nonlinear normal modes, of
ordinary, partial, or delay differential equations
(autonomous or periodically forced)
I provide a procedure to construct a specified invariant
manifold for a specified system of ordinary differential
equations or delay differential equations or partial
differential equations. The invariant manifold may be any of
a centre manifold, a slow manifold, an un/stable manifold, a
sub-centre manifold, a nonlinear normal form, or any
spectral submanifold. Thus the procedure may be used to
analyse pitchfork bifurcations, or oscillatory Hopf
bifurcations, or any more complicated superposition, or
multiscale patterns in PDEs. In the cases
when the neglected spectral modes all decay, the constructed
invariant manifold supplies a faithful, large time,
emergent, model of the dynamics of the differential
equations. Further, in the case of a slow manifold, this
procedure now derives vectors defining the projection onto
the invariant manifold along the isochrons: this projection
is needed for initial conditions, forcing, system
modifications, and uncertainty quantification.
The procedure now also empowers one to account for
sinusoidal time dependence in the ODEs or
PDEs, such as to derive spectral sub-manifold
models of forced nonlinear normal modes.
In the case of PDEs
The analysis assumes
that the 'spatial' gradients are small, that is, the spatial
variations are gradual in space. The manifold is then
multiscale. Currently caters for up to three `spatial'
variables.
Execute on your computer?
The procedure uses computer
algebra, the package
Reduce, to construct approximations to the invariant
manifolds.
Startup Reduce in the unzipped
InvariantManifold folder.
Execute in_tex
"invariantManifold.tex"$ to load the procedure.
Test by executing
exampleslowman(); and confirm the
output is as in invariantManifold.pdf
See examples in
diverseExamples.pdf and then try for
systems of your interest.
Maybe use the web service
Alternatively, click this link to expand this
page. Then via the web form below you may obtain an
invariant manifold of your specified system of ordinary
differential equations (ODE) or delay
differential equations (DDE), when the
ODE/DDE has fast and centre modes. The
modes may be slow, as in a pitchfork bifurcation, or
oscillatory, as in a Hopf bifurcation, or some more
complicated superposition. In the case when the fast modes
all decay, the centre manifold supplies you with a faithful
large time model of the dynamics. For example, this web
page could help you analyse the long time dynamics of the
system \[ \frac{d\vec u}{dt} =\left[\begin{array}{ccc}
2&1&2\\ 1&-1&1\\ -3&-1&-3 \end{array}\right] \vec u
+\varepsilon \left[\begin{array}{c}u_2u_3\\ -u_1u_3\\
-u_1u_2\end{array}\right]. \] As this system is already
entered for you, just enter the magic word, then click on
the Submit button to see.
The bottom of the web form lists further examples.
For example, you can obtain the equivalent modulation
equations corresponding to a given set of
ODE/DDEs that have one or more oscillatory
modes (a Hopf bifurcation for example). The analysis
provides you with equations for the evolution of the complex
amplitudes of the oscillators. This approach is better than
averaging/homogenisation.
FYI: the source code is now available for collaborative
development via the folder CentreManifold of a
Github repository
Fill in the fields below for your ODE/DDE/PDE system:
the n variables of the system must be denoted u1,u2,...,un;
any time tau delayed variables must be denoted
u1(tau),...,un(tau) for any delay tau;
specify the ODE/DDE/PDEs of the dynamical
system \(\frac{d\vec u}{dt}=\vec f(\vec u)\) by specifying
the `nonlinear' function\(~\vec f(\vec u)\) which should
have an equilibrium at the origin;
the invariant manifold is constructed to be tangential to
the invariant subspace of the origin which is defined by you
specifying the m eigenvalues and m eigenvectors,\(~\vec
e_j\), of the m zero and/or pure imaginary eigenvalues of
the matrix\(~L\) of the linearisation about the origin;
the invariant manifold is parametrised by you also
specifying m vectors,\(~\vec z_j\), which are to be
orthogonal to updates to the invariant manifold---these vectors
MAY be the eigenvectors of the m zero and/or pure
imaginary eigenvalues of the adjoint of the matrix\(~L\);
for the moment, the ODE/DDE/PDEs must be
multinomial in form;
Except that for PDEs use diff()
to input the `spatial' derivatives. For example,
diff(u2,x)\(=\partial u2/\partial x\), etc.
Wait a minute or two
The analysis may take minutes after submitting. Be patient.
Read the following. Please inform me of any problems.
In the results
The centre manifold is parametrised by m slow variables
s(j).
Each such slow variable s(j) is either
a 'real' slow variable (specified eigenvalue zero),
or
a complex amplitude of an oscillatory
mode\(~e^{i\omega t}\) for some frequency\(~\omega\)
corresponding to specified eigenvalue \(\lambda=i\omega\),
or
most generally, a (complex) amplitude of the
exponential mode\(~e^{\lambda t}\) for each specified
eigenvalue\(~\lambda\).
In a real system of ODE/DDEs, the complex
amplitudes occur in complex conjugate pairs.
I use the variable small (also appearing
as\(~\varepsilon\)) to control and order the asymptotic
expansion: introducing small into your definition
of the `nonlinear' function empowers finer control of the
asymptotics. For example, in the delay ODE
example the parameter\(~\delta\) is made 'small'.
Analogously, cubic terms may best be made `small' when added
to quadratic terms.
The code does cater for degenerate cases involving
generalised eigenvectors. But the code does this by
modifying\(~L\), the matrix of the linearisation at the
origin. The code attempts to make the smallest modification
it can to remove the degeneracy, and flags the change with
variable 'small' so you recover the original with
\(\verb|small|=1\). Be careful that the results are
relevant to what you want.
Similarly, the code tries to make all, except the
nominated eigenvalues, of\(~\vec z^*L\vec e\) to be 'small',
and also makes\(~\vec f(\vec 0)\) `small'. Be very
careful that the results are relevant to what you want.
In the case of a slow manifold (all specified
eigenvalues are zero), the code also computes a basis of
normal vectors to the isochrons (actually normal to their
tangent space at the slow manifold). For
PDEs, the bare operators such as
diff(x) represent the partial derivative
\(\partial/\partial x\). Use these basis vectors:
to project initial conditions onto the slow
manifold;
to project non-autonomous forcing onto the slow
evolution;
to predict the consequences of modifying the
original system; and
in uncertainty quantification to quantify effects
on the model of uncertainties in the original system.
This example models the double Hopf bifurcation that occurs
in the coupled delay differential equations
\(\frac{dx}{dt}=-4(1+\delta)^2 \left[\frac38y(t) +\frac58
y(t-\pi) \right]\), and \(\frac{dy}{dt} =\left[1+y(t)\right]
x(t)\) as parameter\(~\delta\) crosses zero. Define
\(u_1(t)=x(t)\) and \(u_2(t)=y(t)\); denoted as u1
and u2. The time delayed variable\(~y(t-\pi)\) is
denoted u2(pi). Copy and paste the following
entries.
The web service finds that in terms of complex
amplitudes\(~s_1(t)\), \(s_2(t)\) and their complex
conjugates\(~\bar s_1(t)\) and\(~\bar s_2(t)\), the centre
manifold \[ \vec u=\left[\begin{array}{c} e^{it}s_1
+e^{i2t}s_2 +e^{-it}\bar s_1 +e^{-i2t}\bar s_2\\ \textstyle
-ie^{it}s_1 -\frac12ie^{i2t}s_2 +ie^{-it}\bar s_1
+\frac12ie^{-i2t}\bar s_2 \end{array}\right] \] to error
\({\cal O}(\varepsilon)\). The corresponding evolution on
the centre manifold is \[\begin{array}{l} ds_1/dt=
\varepsilon(-0.04-0.09i)s_2\bar s_1 +{\cal
O}(\varepsilon^2)\\ ds_2/dt= \varepsilon(0.42-0.49i)s_1^2
+{\cal O}(\varepsilon^2). \end{array}\]
Metastability in a four state Markov chain
Variable\(~\varepsilon\) characterises the rate of exchange
between metastable states.
\[\begin{array}{l}
&\dot u_{1}=-\varepsilon u_{1}+u_{2}
\\&\dot u_{2}=\varepsilon \big(u_{3}-u_{2}+u_{1}\big)-u_{2}
\\&\dot u_{3}=\varepsilon \big(u_{4}-u_{3}+u_{2}\big)-u_{3}
\\&\dot u_{4}=-\varepsilon u_{4}+u_{3}
\end{array}
\]
The linear perturbation terms gets multiplied by small again.
Copy and paste the following.
Then\(~s_1\) and\(~s_2\) represent the probabilities of
being in states one and two, and in three and four,
respectively.
Nonlinear normal modes
Renson (2012) explored finite element construction of the
nonlinear normal modes of a pair of coupled oscillators.
Defining two new variables one of their example systems is
\[\begin{array}{rcl}
&&\dot x_1=x_3\,,
\\&&\dot x_2=x_4\,,
\\&&\dot x_3=-2x_1+x_2-\frac12x_1^3+\frac3{10}(-x_3+x_4),
\\&&\dot x_4=x_1-2x_2+\frac3{10}(x_3-2x_4).
\end{array}\]
Copy and paste the following code which makes the linear
damping to be effectively small (which then makes it small
squared); consequently scale the smallness of the cubic
nonlinearity.
The square root frequencies do not cause any trouble
(although one may need to reformat the LaTeX of the exp
operator). In the model, observe that \(s_1=s_2=0\) is
invariant, as is \(s_3=s_4=0\). These are the nonlinear
normal modes.
Harmonically forced nonlinear normal mode
Let's apply periodic forcing to the previous example, both
direct and parametric. For example, here derive the effect
on the mode with frequency one. Defining two new variables
one of their example systems is
\[\begin{array}{rcl}
&&\dot x_1=x_3\,,
\\&&\dot x_2=x_4\,,
\\&&\dot x_3=-2x_1+x_2-\frac12x_1^3+\frac3{10}(-x_3+x_4)+f_1\cos(t),
\\&&\dot x_4=x_1-2x_2+\frac3{10}(x_3-2x_4)f_2\sin(t/2),
\end{array}\]
where \(f_1\) is the strength of the direct forcing, and
\(f_2~\)is the strength of the parametric oscillation in the
last ODE. Copy and paste the following code
which makes both the linear damping and the direct forcing
to be effectively small (which then makes it small squared);
consequently scale the smallness of the cubic nonlinearity.
In the derived ODEs for the modulation of the
frequency one mode, see that the direct forcing drives
effects first seen in terms linear in\(~f_1\). However, the
parametric forcing drives effects quadratic in\(~f_2\) and
so our higher-order, systematic, analysis is required.
Slow manifold among fast oscillations
Lorenz (1986) explored a five equation toy model to
illustrate the quasi-geostrophic approximation and its slow
manifold.
\[\begin{array}{rcl}
&&\dot u=-vw+bvz\,,
\\&&\dot v=uw-buz\,,
\\&&\dot w=-uv\,,
\\&&\dot x=-z
\\&&\dot z=x+buv\,.
\end{array}\]
Copy and paste the following code to find the 3D slow
manifold among the rapid oscillations of\(~x,z\).
Alternatively, find that there is not a coordinate transform
that separates the system into a normal form of separated
fast and slow modes (Cox and Roberts. Initialisation and the
quasi-geostrophic slow manifold.
http://arXiv.org/abs/nlin.CD/0303011, 1994.) Copy and
paste the following to find that, as well as the slow modes
modulating the fast oscillations, the fast oscillations
drive a slow response in\(~w\).
Gradual space variation in a three component system
This is an example from "Slowly varying, macroscale models
emerge from microscale dynamics over multiscale domains"
(2005). The PDEs model a random walker who
walks in any one of three directions in the \(xy\)-plane,
and randomly in time swaps between the three directions.
What is a long-term advection-diffusion model for the
walkers probability distribution?
Long waves in three component Grad's Moment system
This is Example B from "Dynamically Optimal Projection onto
Slow Spectral Manifolds for Linear Systems" by Kogelbauer
and Karlin (2025). The PDEs is apparently a
classical system in kinetic theory. Variables are the
pressure \(u_1\), velocity \(u_2\), and
stress \(u_3\). What are PDEs for the
long-wave pressure-velocity waves?