Lecture 6 - Center Manifold Theory

Based on notes created by Sam Coogan and Murat Arcak. Licensed under a “Creative Commons Attribution-NonCommercial-ShareAlike 4.0 International License”

Additional Reading

• Khalil, Chapter 8.1
• Sastry, Chapter 7.6.1

Overview

• Center Manifold Theory

Motivation for Center Manifold Theory

Center manifold theory is used to study the stability of equilibrium points where linearization fails (for example, when the linearization has purely imaginary or zero eigenvalues).

Theorem (4.7 from Khalil).
Let x=0 be an equilibrium point for the nonlinear system \dot{x} = f(x) where f: D \to \mathbb{R}^n is continuously differentiable and D \subset \mathbb{R}^n is a neighborhood of the origin. Let A = \left.\frac{\partial f}{\partial x}(x) \right|_{x=0} Then,

1. x^* = 0 is asymptotically stable if all eigenvalues of A have negative real parts.
2. x^* = 0 is unstable if at least one eigenvalue of A has a positive real part.

Note: If A has some eigenvalues with zero real part and no eigenvalues with positive real part, then the theorem is inconclusive.

Let’s assume that A has k eigenvalues with zero real parts and m=n-k eigenvalues with negative real parts. Then, we have two options:

1. Analyze a n-dimensional nonlinear system directly (difficult).
2. Use center manifold theory to reduce the system to a k-dimensional system

Mathematical Preliminaries

A k-dimensional manifold in \mathbb{R}^n (1 \leq k < n) is informally the solution to \eta(x) = 0 with \eta: \mathbb{R}^n \to \mathbb{R}^{n-k} sufficiently smooth.

Example 1: The unit circle \{x \in \mathbb{R}^2 \mid x_1^2 + x_2^2 = 1\} is a one-dimensional manifold in \mathbb{R}^2.

Example 2: The unit sphere \{x \in \mathbb{R}^n \text{ s.t. } \sum_{i=1}^n x_i^2 = 1\} is an n-1 dimensional manifold in \mathbb{R}^n.

A manifold is an invariant manifold if: \eta(x(0)) = 0 \implies \eta(x(t)) \equiv 0 \quad \forall t \in [0, t_1) \subset \mathbb{R} where [0, t_1) is any time interval over which x(t) is defined.

Center Manifold Theory

Consider the general nonlinear system with an equilibrium point at the origin:

\dot x = f(x) \qquad f(0) = 0

Suppose A \triangleq \left.\frac{\partial f}{\partial x}\right|_{x = 0} has k eigenvalues with zero real parts, and m = n-k eigenvalues with negative real parts.

Define a change of coordinates \begin{bmatrix} y \\ z \end{bmatrix} = T x such that TAT^{-1} = \begin{bmatrix} A_1 & 0 \\ 0 & A_2 \end{bmatrix} where the eigenvalues of A_1 have zero real parts and the eigenvalues of A_2 have negative real parts.

Rewrite \dot x = f(x) in the new coordinates:

\begin{aligned} \dot y &= A_1 y + g_1(y, z) \\ \dot z &= A_2 z + g_2(y, z) \end{aligned} where g_i(0, 0) = 0, \frac{\partial g_i}{\partial y}(0, 0) = 0, \frac{\partial g_i}{\partial z}(0, 0) = 0, i = 1, 2.

g_1 and g_2 inherit the properties of \tilde{f} in the equation: \dot{x} = f(x) = A x + \tilde{f}(x)

with \tilde{f}(x) = f(x) - \frac{\partial f}{\partial x}(x)\mid_{x=0}, which has the properties \tilde{f}(0) = 0 and \frac{\partial f}{\partial x}(0) = 0

Theorem 1:
There exists an invariant manifold z = h(y) defined in a neighborhood of the origin such that h(0) = 0 \qquad \frac{\partial h}{\partial y}(0) = 0.

z = h(y) is called a center manifold in this case.

Reduced System: \dot y = A_1 y + g_1(y, h(y)) \qquad y \in \mathbb{R}^k

Theorem 2:
If y = 0 is asymptotically stable (resp., unstable) for the reduced system, then x=0 is asymptotically stable (resp., unstable) for the full system \dot x = f(x).

Characterizing the Center Manifold

Define w \triangleq z - h(y) and note that it satisfies \begin{aligned} \dot w &= \dot{z} - \frac{\partial h}{\partial y} \dot{y} \\ &= A_2 z + g_2(y, z) - \frac{\partial h}{\partial y}\Big( A_1 y + g_1(y, z)\Big). \end{aligned}

The invariance of z = h(y) means that w = 0 implies \dot w = 0. Thus, the expression above must vanish when we substitute z = h(y):

A_2 h(y) + g_2(y, h(y)) - \frac{\partial h}{\partial y}\Big( A_1 y + g_1(y, h(y))\Big) = 0

To find h(y) solve this partial differential equation for h as a function on y.

If the exact solution is unavailable, an approximation might be sufficient.

For scalar y, expand h(y) as h(y) = h_2 y^2 + \dots + h_p y^p + O(y^{p+1}) where h_1 = h_0 = 0 because h(0) = \frac{\partial h}{\partial y}(0) = 0. The notation O(y^{p+1}) refers to the higher order terms of power p+1 and above.

Example (8.2 from Khalil): \begin{aligned} \dot y &= y z \\ \dot z &= -z + a y^2 \qquad a \ne 0 \end{aligned} This is of the form above with g_1(y, z) = y z, g_2(y, z) = a y^2, A_2 = -1.

Thus h(y) must satisfy -h(y) + a y^2 - \frac{\partial h}{\partial y} y h(y) = 0. Try h(y) = h_2 y^2 + O(y^3): \begin{aligned} 0 &= -h_2 y^2 + O(y^3) + a y^2 - (2 h_2 y + O(y^2)) y (h_2 y^2 + O(y^3)) \\ &= (a - h_2) y^2 + O(y^3) \\ &\implies h_2 = a \end{aligned} Reduced System: \dot y = y(a y^2 + O(y^3)) = a y^3 + O(y^4).

If a < 0, the full system is asymptotically stable in a neighborhood around the origin. If a > 0, unstable.

Example 2

Consider the nonlinear system:

\dot{x}_1 = x_1 x_2^3 \\ \dot{x}_2 = -x_2 - x_1^2 + 2 x_2^8

import numpy as np
import matplotlib.pyplot as plt

# Define the system
def f(x1, x2):
  dx1 = x1 * x2**3
  dx2 = -x2 - x1**2 + 2 * x2**8
  return dx1, dx2

# Create a grid
x1 = np.linspace(-2, 2, 20)
x2 = np.linspace(-2, 2, 20)
X1, X2 = np.meshgrid(x1, x2)
DX1, DX2 = f(X1, X2)

# Plot streamplot
plt.figure(figsize=(6, 6))
plt.streamplot(X1, X2, DX1, DX2, color='b', density=1.2, linewidth=1)
plt.xlabel('$x_1$')
plt.ylabel('$x_2$')
plt.title('Streamplot of the Nonlinear System')
plt.grid(True)
plt.show()

The equilibrium point is x^* = (0,0) which has the linearization:

A = \left.\frac{\partial f}{\partial x}\right|_{x=0} = \begin{bmatrix} 0 & 0 \\ 0 & -1 \end{bmatrix}

Since the eigenvalues of A are \lambda_1 = 0 and \lambda_2 = -1, the linearization is inconclusive about the stability of the equilibrium point so instead we will use center manifold theory.

This system is already in the form required for center manifold theory. But for clarity, we will define our change of coordinates as:

\begin{bmatrix} y \\ z \end{bmatrix} = \underbrace{\begin{bmatrix} 1 & 0 \\ 0 & 1 \end{bmatrix}}_{P^{-1}} \begin{bmatrix} x_1 \\ x_2 \end{bmatrix} for some similarity transformation \begin{bmatrix} A_1 & 0 \\ 0 & A_2 \end{bmatrix} = P^{-1}AP.

Using the change of coordinates, our transformed system is:

\dot{y} = y z^3 \\ \dot{z} = -z - y^2 + 2 z^8

To solve for the center manifold, we assume that z = h(y) where h(y) is a function we need to determine. We use the expression derived earlier for the characterization of the center manifold w \triangleq z - h(y) = 0:

\begin{align*} \dot{w} &= \dot{z} - \frac{\partial h (y)}{\partial y} \dot{y} = 0 \\ &= (-z - y^2 + 2 z^8) - \frac{\partial h (y)}{\partial y} (y z^3) \\ &= (-h(y) - y^2 + 2 h(y)^8) - \frac{\partial h (y)}{\partial y} (y h(y)^3) \\ \end{align*}

Let’s approximate h(y) as the power expansion h(y) = h_2 y^2 + O(3). Plugging this in gives us the expression:

\begin{align*} \dot{w} &= (-(h_2 y^2 + O(3)) - y^2 + 2 (h_2 y^2 + O(3))^8) - \frac{\partial (h_2 y^2 + O(3))}{\partial y} (y (h_2 y^2 + O(3))^3) \\ &= (-(h_2 + 1) y^2 + O(3)) \\ &= (-h_2 - 1) y^2 + O(3) \end{align*}

This gives us h_2 = -1 which results in the center manifold approximation: h(y) = -y^2 + O(3).

Plugging z = h(y) = -y^2 + O(3) back into the y-dynamics gives us the reduced system on the center manifold: \dot{y} = y (-y^2 + O(3))^3 = -y^7 + O(9)

Therefore, for solutions near the origin (y \approx 0), we can conclude that the equilibrium point is stable since \dot{y} < 0 for small nonzero y.

# Overlay the center manifold approximation x2 = -x1^2
plt.figure(figsize=(6, 6))
plt.streamplot(X1, X2, DX1, DX2, color='b', density=1.2, linewidth=1)
plt.xlabel('$x_1$')
plt.ylabel('$x_2$')
plt.title('Streamplot with Center Manifold Approximation')
plt.grid(True)

# Plot the center manifold: x2 = -x1^2
x1_curve = np.linspace(-2, 2, 200)
x2_curve = -x1_curve**2
plt.plot(x1_curve, x2_curve, 'r', linewidth=2, label='$x_2 = -x_1^2$')
plt.ylim([-2, 2])
plt.legend()
plt.show()

Higher order manifold

We can show that a higher order power expansion of the center manifold would give us a better approximation. We can use code to systematically solve for an 8^{th} order approximation:

from sympy import symbols, Eq, solve, expand, collect, series

# Define the symbols
y = symbols('y')
h2, h3, h4, h5, h6, h7, h8 = symbols('h2 h3 h4 h5 h6 h7 h8')

# Define h(y)
h_y = h2 * y**2 + h3 * y**3 + h4 * y**4 + h5 * y**5 + h6 * y**6 + h7 * y**7 + h8 * y**8

# Compute the derivative
h_y_prime = h_y.diff(y)

# Define the invariance condition
invariance_condition = (-h_y - y**2 + 2 * (h_y)**8) - (h_y_prime * y * (h_y)**3)

# Expand and truncate the series at y^3
invariance_condition = expand(invariance_condition).series(y, 0, 9).removeO()

# Collect terms in y
coeff_y0 = invariance_condition.coeff(y, 0)
coeff_y1 = invariance_condition.coeff(y, 1)
coeff_y2 = invariance_condition.coeff(y, 2)
coeff_y3 = invariance_condition.coeff(y, 3)
coeff_y4 = invariance_condition.coeff(y, 4)
coeff_y5 = invariance_condition.coeff(y, 5)
coeff_y6 = invariance_condition.coeff(y, 6)
coeff_y7 = invariance_condition.coeff(y, 7)
coeff_y8 = invariance_condition.coeff(y, 8)

# Solve for h2, h3, ..., h8 such that all coefficients vanish
solutions = solve([Eq(coeff_y0, 0), Eq(coeff_y1, 0), Eq(coeff_y2, 0), Eq(coeff_y3, 0), Eq(coeff_y4, 0), Eq(coeff_y5, 0), Eq(coeff_y6, 0), Eq(coeff_y7, 0), Eq(coeff_y8, 0)], (h2, h3, h4, h5, h6, h7, h8))

# Print the solutions
print("coefficients for 8th order expansion:")
print(solutions)

coefficients for 8th order expansion:
[(-1, 0, 0, 0, 0, 0, -2)]

This gives us the manifold z = h(y) = -y^2 - 2 y^8. Plugging this into our y-dynamics gives us the reduced system on the center manifold: \dot{y} = y z^3 = y( -y^2 - 2y^8 ) = -y^3 - 2 y^9 which again concludes local stability at the origin.

# Overlay the center manifold approximation x2 = -x1^2 - 2 x1^8
plt.figure(figsize=(6, 6))
plt.streamplot(X1, X2, DX1, DX2, color='b', density=1.2, linewidth=1)
plt.xlabel('$x_1$')
plt.ylabel('$x_2$')
plt.title('Streamplot with Center Manifold Approximation')
plt.grid(True)

# Plot the center manifold: x2 = -x1^2
x1_curve = np.linspace(-2, 2, 200)
x2_curve = -x1_curve**2 - 2 * x1_curve**8
plt.plot(x1_curve, x2_curve, 'r', linewidth=2, label='$x_2 = -x_1^2 - 2 x_1^8$')
plt.ylim([-2, 2])
plt.legend()
plt.show()