Finding the stability boundary "or attraction region" of a nonlinear differential equation

Hi everyone.

I have been stuck into this problem for weeks. How do I find the stability (or attraction ) region of a nonlinear differential equation using Matlab.

Let's say I have this equation:

if true
x' = y;
y' = -10*sin(x)  - y + 9;
end

The equilibrium point for this equation is [x , y] = [1.1198 , 0]. I wanted to draw the stability boundary of this nonlinear differential equation. By that I mean, I want to find the region where any initial point will converge to the equilibrium point, any any point outside that region will diverge. Please see attached image

Right now, I run the following Matlab code:

if true
f = @(t , x)[x(2) ; -10 * sin(x(1)) - x(2) + 9];
[T , X] = ode45(f , Tint , X0);
end

For some Tint, I plot the result in a phase portrait figure (i.e., x vs y), and I vary the initial condition (X0) till it works (i.e., some educated trial and error).

I need to find the stability region for many different variations of this differential equation. My question is: How do I find this region automatically?

Thanks for the help

Answer 1

[Yu Jiang]

Hi Abdallah

It is more practical to find an estimate of the region of attraction (ROA), instead of the exact region of attraction as you defined. Any initial point within an estimate of the ROA will converge to the equilibrium point. But it does not necessarily mean that the initial point outside of it will result in a divergent trajectory (the definition of estimate region of attraction can be found in the book “Nonlinear Systems” by Hassan Khalil).

To find an estimate of ROA, one popular way is to find a Lyapunov function for the nonlinear system and compute its contour map, from which the estimate of ROA can be observed. However, Lyapunov function is different from one system to another, and it is not trivial to find the correct Lyapunov function for a given system. Therefore, I don’t think there is an automatic way to do it.

However, I think it will be useful if you plot the vector field diagram of the system. For example, try execute the following code

x = linspace(0,4,20);
y = linspace(-5,3,20);
[x,y] = meshgrid(x,y);
u = y;
v = -10*sin(x) - y +9;
quiver(x,y,u,v,5);
axis([0 4 -5 3])

For more details about the function quiver, see the documentation .

-Yu