1

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:

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

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, and any point outside that region will diverge. Please see attached image at http://www.mathworks.com/matlabcentral/answers/146562-finding-the-stability-boundary-or-attraction-region-of-a-nonlinear-differential-equation

Right now, I run the following Matlab code:

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

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

Abdallah
  • 11
  • 2

1 Answers1

2

At a very high computational cost, you could do a brute force search: discretize your search space, choose a long Tint (to account for solutions that converge in a really long time), and integrate for each starting point.
In order to accelerate the process a bit, you could implement some stopping conditions for the integration, and prune the branches that get too far away from your equilibrium point, as well as the ones that go near the equilibrium point.

Of course, this way, you would only get an approximation of the stability boundary, due to 1. the discretization step chosen for your search space, and 2. possible 'mispruning' if a solution converges after getting further from the equilibrium point than your stopping condition, or diverges after getting closer to the equilibrium than your stopping condition.

For help implementing stopping conditions, look at the event functions part of the ode45 doc, as well as this example, or this one

Etienne Pellegrini
  • 738
  • 1
  • 5
  • 11