ER3.2. Método directo de Lyapunov aplicado a un modelo no lineal

Introducción (planteamiento)

Dado el siguiente modelo, determinar la estabilidad de uno de sus puntos de equilibrio aplicando el método directo de Lyapunov:

$\left\{ \begin{aligned} &\dot{x}_1=-x_1+x_{1}^{2}x_2\\ &\dot{x}_2=x_1-x_2\\ \end{aligned} \right.$

Método (plan de solución)

Pasos: (i) obtener con la función pplane o plotpp de MATLAB el retrato de fase con los puntos de equilibrio (para tener un mejor panorama del problema, dado que el modelo es de segundo orden); (ii) calcular los puntos de equilibrio; (iii) proponer una función definida positiva y aplicar el método de Lyapunov para uno de los puntos de equilibrio; (iv) analizar los resultados.

Resultados (solución)

Retrato de fase (ver código de MATLAB en este enlace):

Se observan al menos tres puntos de equilibrio, cuyo cálculo matemático se realiza de la siguiente manera:

$\left\{ \begin{array}{l}-x_1+x_{1}^{2}x_2=0\\ x_1-x_2=0\\ \end{array} \right.$

De la segunda ecuación:

$x_1=x_2$

Reemplazando en la segunda ecuación:

$-x_1+x_{1}^{3}=0$

Factorizando:

$x_1(1-x_{1}^{2})=0$

Solución de la última ecuación: $x_1=\{0,-1,1\}$. De esta manera, los puntos de equilibrio son: (0,0), (-1,-1), (1,1). Se analiza ahora la estabilidad del punto en el origen, dado que para las otros dos es necesario hacer un cambio de variables para ubicarlos en el origen, pero al ser puntos de silla (inestables) no existe una función de Lyapunov. Se toma como función de energía la siguiente función (la más obvia en todos los casos):

$V(\mathbf{x})=x_{1}^{2}+x_{2}^{2}$

El gráfico de la función anterior se puede obtener con MATLAB.

x1 = [-5:0.01:5]; 

x2 = x1; 

[x1,x2] = meshgrid(x1,x2); 

V = x1.^2 + x2.^2; 

meshc(x1,x2,V)

Al derivar esta función y reemplazar con el modelo específico de este caso se obtiene:

$\dot{V}(\mathbf{x})=-2(x_{1}^{2}+x_{2}^{2})+2x_1x_2(x_{1}^{2}+1)$

De esta expresión no se puede deducir directamente si es definida negativa, por lo que es necesario realizar ciertas operaciones para determinar si existe un dominio D en el que la función es definida negativa:

$-2(x_{1}^{2}+x_{2}^{2})+2x_1x_2(x_{1}^{2}+1)<0$

$x_1x_2(x_{1}^{2}+1)<x_{1}^{2}+x_{2}^{2}$

$x_1x_2<\frac{x_{1}^{2}+x_{2}^{2}}{x_{1}^{2}+1}<1$   si   $\left| x_2 \right|<1$

Por lo tanto, existe una región D que contiene el punto de equilibrio (0,0) en la cual la derivada de la función de energía es definida negativa, dado que $x_1x_2$ es menor que uno si ambas variables son menores que uno con cualquier signo. Por lo tanto, el punto de equilibrio es asintóticamente estable. 

Consejo: antes de probar si  $\dot{V}(\mathbf{x})$  es definida negativa se puede graficar o programar con MATLAB una función para probar si lo es en una región cercana al origen:

x1 = -5:0.05:5; 

x2 = x1; 

[x1, x2] = meshgrid(x1, x2); 

V = -2*(x1^2 + x2^2) + 2*x1*x2*(x1^2 + 1); 

meshc(x1, x2, V) ; 

max(V(.)) % Para ver si max = 0 solo en el origen

 

Discusión y verificación

El tipo de estabilidad del punto de equilibrio coincide con el tipo dado por el retrato de fase, como se ve en el gráfico, por lo que el resultado es correcto. En este caso se tuvo suerte con la selección de la función definida positiva, pero en caso contrario existen diversos métodos para buscar una forma de dicha función.