Finding Lyapunov function to control the Duffing Oscillator

Question

I'm trying to do this exercise but I'm struggling to find the Lyapunov function.
We have the state version of the Duffing oscillator below:

$\dot{x_1} = x_2$
$\dot{x_2} = x_1 - x_1^3 - \delta x_2$

which results in the unforced Duffing oscillator if δ = 0, and the damped Duffing oscillator for δ > 0.

We want:

1. Identify the fixed points and their local stability properties.
2. Find a Lyapunov or Lyapunov-like function for the stable fixed point(s).
3. Use the function to characterize the size of the region of attraction for the damped version at each fixed point.

What i did is:

1. I find the fixed points : (0,0);(1,0);(-1,0) and I compute the jacobian for each fixed points to look at their local stability.
   I find that the points (1,0);(-1,0) have the same jacobian and are stable when $\delta \gt 0$ and the point (0,0) is not stable.
2. I tried with the lyapunov function:
   $V(x_1,x_2)=\frac 12 (x_1^2+x_2^2)$
   And I got:
   $\dot V(x_1,x_2) =2x_1x_2-x_2x_1^3-\delta x_2^2$
   Which is not $\lt 0$
   So I was thinking trying something like:
   $V(x_1,x_2)=\alpha x_1^m+\beta x_2^n$
   But after computing the derivative it seems to give nothing too.
   Do you have any suggestion for the lyapunov function?
3. For this part honestly I don't really know what to do maybe I was thinking of finding a lyapunov function at the question 2 which is negative for a region on $x_1$ or $x_2$.

Answer 1

From the Jacobian that you obtained for the stable equilibrium points $(-1, 0)$, and $(1, 0)$, the linearized system becomes

$$ \left[\matrix{\dot{x}_{1} \cr \dot{x}_{2}}\right] = \left[\matrix{0 & 1 \cr -2 & -\delta}\right] \left[\matrix{x_{1} \cr x_{2}}\right] .$$

Consider the Lyapunov function in the following quadratic form:

$$ V(\mathbf{x}) = \frac{1}{2} \mathbf{x}^{T} \mathbf{P} \mathbf{x} = \frac{1}{2} \left[\matrix{x_{1} & x_{2}}\right] \left[\matrix{\frac{\delta^{2} + 2 \delta + 2}{\delta} & 1 \cr 1 & \frac{\delta + 1}{\delta}}\right] \left[\matrix{x_{1} \cr x_{2}}\right] .$$

The positive definite $\mathbf{P}$ is found by solving the Lyapunov equation $\mathbf{P} \mathbf{A} + \mathbf{A}^{T} \mathbf{P} = - \mathbf{Q}$. In this case, $\mathbf{Q}$ is chosen as

$$ \mathbf{Q} = \left[\matrix{2 & 0 \cr 0 & \delta}\right] .$$

The derivative of $V(\mathbf{x})$ along the trajectories of the linear system is given by

$$ \dot{V}(\mathbf{x}) = - 2 x_{1}^{2} - \delta x_{2}^{2} < 0 \quad \text{for} \quad \mathbf{x} \neq \mathbf{0} .$$

Remember that this result, however, is local.

Region of Attraction

To estimate the region of attraction, the previous Lyapunov function can be used to determine a domain $D$ about the stable equilibrium by solving the derivative of $V(\mathbf{x})$ along the trajectories of the nonlinear system (Duffing Oscillator) where $\dot{V}(\mathbf{x})$ is negative definite.

$$ \dot{V} = \left(\frac{\delta^{2} + 2 \delta + 2}{\delta}\right) x_{1} \dot{x}_{1} + \dot{x}_{1} x_{2} + x_{1} \dot{x}_{2} + \left(\frac{\delta + 1}{\delta}\right) x_{2} \dot{x}_{2} $$

$$ \dot{V} = \left(\frac{\delta^{2} + 2 \delta + 2}{\delta}\right) x_{1} x_{2} + x_{2}^{2} + x_{1} \left(x_{1} - x_{1}^{3} - \delta x_{2}\right) + \left(\frac{\delta + 1}{\delta}\right) x_{2} \left(x_{1} - x_{1}^{3} - \delta x_{2}\right) $$

$$ \dot{V} = - x_{1}^{4} - \left[\left(\frac{\delta + 1}{\delta}\right) x_{1} x_{2} - 1\right] x_{1}^{2} + \left(\frac{3 \delta + 3}{\delta}\right) x_{1} x_{2} - \delta x_{2}^{2} $$

For example, if $\delta = 1$, the region of attraction is found to be satisfied in $D_{(1, 0)} = \left\{\sqrt{2} < x_{1} < \sqrt{5}, x_{2} \in \mathbb{R}\right\}$, and $D_{(-1, 0)} = \left\{-\sqrt{5} < x_{1} < -\sqrt{2}, x_{2} \in \mathbb{R}\right\}$.

The streamplot for the Duffing Oscillator with $\delta = 1$ is shown below. Note that $\pm \sqrt{5} \approx \pm 2.236$.