Overview

I provide a procedure to construct a specified invariant manifold for a specified system of ordinary differential equations or delay 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. 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, such as to derive spectral sub-manifold models of forced nonlinear normal modes.

Execute on your computer? The procedure uses computer algebra, the package Reduce, to construct approximations to the invariant manifolds.

  1. So download and install Reduce, and then download InvariantManifold.zip.
  2. Startup Reduce in the unzipped  InvariantManifold  folder.
  3. Execute  in_tex "invariantManifold.tex"$  to load the procedure.
  4. Test by executing  exampleslowman();  and confirm the output is as in  invariantManifold.pdf 
  5. 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

Non-autonomous ODEs?

Systems with slowly-varying time-dependence, or with sinusoidal time-dependence may be analysed here. Systems with more general time-dependence are significantly more difficult, but are analysed via the web page Normal form of stochastic or deterministic multiscale differential equations.

Submit your system of ODE/DDEs for analysis

Fill in the fields below for your ODE/DDE 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/DDEs 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/DDEs must be multinomial in form;
  • Use the syntax of Reduce for the algebraic expressions

Enter the magic word "a z a l e a" into

Description Specify your ODEs
RHS function the nonlinear vector function: parentheses around a comma separated list of multinomial expressions. (However, due to a glitch in the current web-site's system do not use "0" for the RHS component of a 'variable' that is to be constant in time, instead use "small^9" as this is effectively zero)
Invariant eigenvalues a comma separated list of m eigenvalues, real or complex, zero for a slow mode.
Invariant subspace m basis vectors: must be eigenvectors of\(~L\) corresponding to the m eigenvalues: comma separated list of eigenvectors; each eigenvector is parentheses around comma separated numbers.
Adjoint basis m basis vectors: comma separated list of vectors; each vector is parentheses around comma separated numbers; for each of the m invariant eigenvalues, in order, may be complex conjugate of a corresponding eigenvector of the adjoint matrix.
Order of error of the analysis in the `nonlinear' terms on the RHS.
Print expressions with the following variables factored---this does not affect the analysis, but must not be empty.

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.
  • For explanations and relevant theory, see my book Modelling emergent dynamics in complex systems, or Low- dimensional modelling of dynamics via computer algebra, or the classic Simple examples of the derivation of amplitude equations for systems of equations possessing bifurcations.
  • 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). 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.
    For explanations and relevant theory, see Computer algebra derives correct initial conditions for low-dimensional dynamical models, or the classic Appropriate initial conditions for asymptotic descriptions of the long term evolution of dynamical systems.

Other example systems of ODEs/DDEs

Double Hopf bifurcation in a delay DE

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.
Description Delay ODE example
RHS function (-4*(1+small*delta)^2*(5/8*u2 +3/8*u2(pi)), +u1*(1+u2))
Centre eigenvalues i,2*i,-i,-2*i
Centre subspace (1,-i), (1,-i/2), (1,+i), (1,+i/2)
Adjoint subspace (1,-i), (1,-2*i), (1,+i), (1,+2*i)
Order of error 3
Print small,delta
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.
Description ODE example
RHS function (u2,-u2,-u3,u3) +small*(-u1, +u1-u2+u3, +u2-u3+u4, -u4)
Slow eigenvalues 0,0
Slow subspace (0,0,0,1), (1,0,0,0)
Adjoint subspace (0,0,1,1), (1,1,0,0)
Order of error 5
Print small
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.
Description ODE example
RHS function ( u3, u4, -2*u1 +u2 -small*u1^3/2 +small*3/10*(-u3+u4), u1 -2*u2 +small*3/10*(u3 -2*u4) )
All eigenvalues i, -i, i*sqrt(3), -i*sqrt(3)
All space basis (1,1,+i,+i), (1,1,-i,-i), (1,-1,i*sqrt(3),-i*sqrt(3)), (1,-1,-i*sqrt(3),i*sqrt(3))
Adjoint subspace (1,1,+i,+i), (1,1,-i,-i), (-i*sqrt(3),+i*sqrt(3),1,-1), (+i*sqrt(3),-i*sqrt(3),1,-1)
Order of error 3
Print small
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.
Description ODE example
RHS function ( u3, u4, -2*u1 +u2 -small*u1^3/2 +small*3/10*(-u3+u4)+small*f_1*sin(t), u1 -2*u2 +small*3/10*(u3 -2*u4)*f_2*cos(t/2) )
Modal eigenvalues i, -i
Modal space basis (1,1,+i,+i), (1,1,-i,-i)
Adjoint subspace (1,1,+i,+i), (1,1,-i,-i)
Order of error 5
Print f_1,f_2,small
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\).
Description ODE example
RHS function (-u2*u3+b*u2*u5, +u1*u3 -b*u1*u5, -u1*u2, -u5, u4 +b*u1*u2)
Slow eigenvalues 0,0,0
Slow subspace (1,0,0,0,0), (0,1,0,0,0), (0,0,1,0,0)
Adjoint subspace (1,0,0,0,0), (0,1,0,0,0), (0,0,1,0,0)
Order of error 5
Print small,b
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\).
Description ODE example
RHS function (-u2*u3 +b*u2*u5, +u1*u3 -b*u1*u5, -u1*u2, -u5, u4 +b*u1*u2)
All eigenvalues 0, 0, 0, i, -i
All space basis (1,0,0,0,0), (0,1,0,0,0), (0,0,1,0,0), (0,0,0,1,-i), (0,0,0,1,+i)
Adjoint subspace (1,0,0,0,0), (0,1,0,0,0), (0,0,1,0,0), (0,0,0,1,-i), (0,0,0,1,+i)
Order of error 4
Print small,b

If you like this web page, please link to it so others can find it more easily.