Evaluating stability of fixed points in a system of differential equations,

Hi there!

Let's say I have a set of coupled ordinary differential equations for which I have found a certain number of fixed points, say, 50. Now I want to evaluate the stability of these 50 fixed points, so I have to first compute the Jacobian matrix of these differential equations, and then evaluate this Jacobian matrix at each of the 50 fixed points. Now I have 50 Jacobian matrices that I need to find eigenvalues for. Stability theory says that for any fixed point, if all of the corresponding eigenvalues have real part < 0, then the fixed point is stable. Otherwise, the fixed point is unstable. That is at least my preliminary understanding of stability.

How can I implement this in Matlab?

Should I write the differential equations symbolically (I have the Symbolic Math Toolbox), and then compute the Jacobian 50 times, and find eigenvalues for all 50 matrices using, say, a for loop? Or should I instead do this numerically? In the recent past, I wrote codes that started symbolically, and later converted to numerical codes, using, say, Matlab's EquationsToMatrix( ) function. Then, with numerical files, I created simulation codes to study a dynamical system. I wonder if I should be doing something similar -- writing both symbolic and numerical codes in Matlab -- in order to evaluate stability.

Sorry for my lack of understanding here.

I'm sort of trying to gauge the difficulty of this task, as well as trying to understand what would be considered best practice, before I proceed, so any help would be greatly appreciated.

Thanks in advance!

[Paul]

If you don't mind me asking, what kind of system are you working on that has 50 fixed points?

[Sam Chak]

Hi @Noob

If you merely want to determine the stability of a fixed point, then it may be unnecessary to simulate the dynamical system. The isstable() function, perhaps, is sufficient. However, if you wish to investigate further to find out how stable that fixed point is, then running a simulation becomes necessary.

For example, a linearized system with a fixed point characterized by a pair of complex eigenvalues with small negative real parts and very large imaginary parts is still considered 'locally exponentially stable' from a mathematical theory standpoint. But, it is, in fact, very close to being classified as 'marginally stable' from a practical perspective.

%% Eigenvalues of the system
format short g
ev  = - eps + sqrt(realmax)*1i;
evc = conj(ev);
sol = [ev; evc];
disp("The eigenvalues are: "); disp(sol)

The eigenvalues are: 
  -2.2204e-16 +1.3408e+154i
  -2.2204e-16 -1.3408e+154i

G   = zpk([], sol, prod(sol));

%% Determine stability: True (1) or False (0)
TF  = isstable(G)   

TF = logical
   1

%% Step response of the system
G   = tf(G)

G =
 
            1.798e308
  -----------------------------
  s^2 + 4.441e-16 s + 1.798e308
 
Continuous-time transfer function.

step(G), grid on

[Sam Chak]

A system that has 50 fixed points? 🤔 I'd imagine that either the system has periodic solutions or that it functions similarly to an Elevator servicing a 50-floor tall building. If the latter is the case, it is not unusual to provide a detailed stability proof in order to obtain a building permit or approval for lift installation from the local council.

[Noob]

Hi Paul and Sam!

Oops! This is embarrassing:

To clarify, I have found fixed points for a system by sweeping through parameters, so the system is changing each time. I do not have one such system with 50 fixed points. Sorry about that. Thanks for your comments! I will read all tonight and work on this again. Thanks for your time!

[Sam Chak]

Hey @Noob,

If your system is periodically time-varying, you should consider studying Floquet Theory. This theory indicates that even if the eigenvalues have negative real parts, the system may still be unstable. Therefore, your Jacobian matrix analysis for such systems may be insufficient to conclude the stability.

However, Floquet Theory is typically not covered in introductory undergraduate courses such as Elementary Differential Equations, System Dynamics, and Control Theory. You can still find it in the following books:

[Noob]

Hi Sam!

Thanks so much for this!

I have recently been running simulations! I think it's perhaps wiser to do so, to get a feel for how good my calculations of the fixed points were. Which leads me to my new question I'll ask shortly on this site!

If you're curious about what my question is, well, I am seeing my dynamical state variables, say, for "very long times" stay constant for up to 13 to 14 digits after the decimal point. To simulate, using the pretty standard ode45, I used the most stringent tolerances that Matlab was happy with: 3e-14, for both AbsTol and RelTol (both of which I need to learn more about). Is this sort of "promising"? Is this as good as it gets? Even before formal stability calculations, I'm just trying to make sure that my calculations of fixed points agree with simulations. Can we do better than the 13 to 14 digits of constant behavior?

Thanks again!

[Sam Chak]

Hi @Noob

If you mean that the values of the state variable have remained constant and relatively consistent up to 14 decimal places for a considerable period when no external "energy" is being added to the dynamical system, this may imply that the system response has reached the steady state. This set of steady-state values is mathematically referred to as the "Equilibrium Point."

Based on your description, you may be struggling to get the solver to reach the "theoretical Equilibrium Point," but are unable to do so. Consider the ordinary differential equation as an example. The theoretical Equilibrium Point of the system is at . However, no matter how large (and positive) the value you set for the parameter k, it will never truly reach in finite time.

format long g
[t, x]  = ode45(@(t, x) - 10*x, [0 100], 1);
x(end)

ans = 
      2.92615760447614e-08

Answer 1

[John D'Errico]

First, there is nothing wrong with a loop! You know how to evaluate stability for each case. Just write the loop. Could you do something more sophisticated? Possibly, but it would be more complex to understand and code. Remember that you need to be able to read the code you wrote next year, when you need to debug it. And when the loop will be just as efficient, don't freak out about not using loops. And what happens when some of the cases are marginal? (There always will be some cases like that.) You want to be able to look more carefully at them. Again, this is something you would be able to do far bettter in a loop.

Should you do exactly as you suggest? That will work. So do it. Never worry about whether some scheme is the best, most highly optimal scheme, unless the code you write becomes a computational problem. Never pre-optimize your code! Write it. Make it work. Remember that your programming time is worth far more than a few milliseconds of CPU time, even a few minutes or more.

So I would do as you say. The only thing I would plan on is using a numerical eigenvalue decomposition, thus the double version of eig, since we expect the symbolic eig will be much slower. Always watch for issues of course, but there should be nothing to fear.

There is nothing difficult in this, unless you try to make it more difficult than it needs to be.

[Noob]

Thanks so much John!