# Stability analysis of transcendental equation (stability crossing curves)

I am working with a non-linear delay system with three scalar delays. After taking the Laplace transform of the linearized system, the characteristic function is a transcendental equation with three different complex exponentials:

$$F(s)=(s+\lambda)-\left(\frac{\lambda}{e^{-\tau_1\lambda}-e^{-\tau_2\lambda}}\right)(e^{-\tau_1(s+\lambda)}-e^{-\tau_2(s+\lambda)})+\left(\frac{e^{-\tau_1\lambda}-e^{-\tau_2\lambda}-\lambda}{1-e^{-\tau_3\lambda}}\right)(1-e^{-\tau_3(s+\lambda)})=0$$

where $0<\tau_1\leq\tau_2\leq\tau_3$ and $0<\lambda\ll1$.

One extra condition that should be satisfied for stability is that

$$\tau_2>\frac{1}{\lambda}\ln\left(\frac{1}{e^{-\lambda\tau_1}-\lambda}\right)$$

This transcendental equation has infinitely many roots (eigenvalues), but they usually lie on the left half plane, and keep the system stable. However, sometimes one or more conjugate pairs of roots cross to the right half plane (positive real part), and the function is well-behaved enough such that this crossing of the roots occur in a continuous fashion. Thus, the boundary of stability in the parameter space is given when the root with the largest real part is purely imaginary, a Hopf bifurcation.

The problem would be, then, to plot such stability crossing curves, that is, a plot in the 2D parameter space (I should fix $\lambda$ and $\tau_3$ and let $\tau_1$ and $\tau_2$ vary, for example) of the curves where the root with largest real part are purely imaginary. The difficulty of analysing this system is that the constant coefficients multiplying the complex exponentials actually depend on the delays, making it impossible to use the formalism of the above-mentioned paper.

I have already managed to plot all the roots with positive real part given some fixed parameters using Reduce, that is relatively simple.

f = (s + l) - (l/(Exp[-l*T1] - Exp[-l*T2]))*(Exp[-T1*(s + l)] -
Exp[-T2*(s + l)]) + ((Exp[-l*T1] - Exp[-l*T2] - l)/(1 -
Exp[-l*T3]))*(1 - Exp[-T3 (s + l)]) /. {l -> 0.003, T1 -> 0,
T2 -> 3, T3 -> 200};
rts = Reduce[f == 0 && Abs[s] < 2 && Re[s] >= 0, s]
ContourPlot[{Re[f /. s -> x + I y] == 0,
Im[f /. s -> x + I y] == 0}, {x, 0, 0.005}, {y, 0.07, 0.14},
RegionFunction -> (Abs[#1 + I #2] < 1 &),
Epilog -> {Red, PointSize[0.025],
Point[{Re[s], Im[s]}] /. {ToRules[rts]}}]

This is zoomed around two of the roots, the complex conjugate are always solutions as well.

The advantage is that the absolute value of the positive roots cannot be greater than $4$, so it suffices to check only inside the half circle for these roots, $(r,\theta), r \in [0,4]$ and $\theta \in [-\frac{\pi}{2},\frac{\pi}{2}]$.

The difficulty now would be to find the stability crossing curves. A first step would be to calculate the largest eigenvalue (the root with largest real part), given a choice of parameters. A step further is to plot the contours of the crossing of the second largest eigenvalue (or conjugate pairs), the third largest and so on.

Another useful plot would be the the imaginary part of the largest root (if they lie in the right half plane).

Any ideas how that could be done in a reasonable amount of time?

• Your RegionFunction option seems incomplete – Dr. belisarius Jun 4 '15 at 19:08
• Corrected, sorry! – Bruno Pace Jun 4 '15 at 21:25
• Can you make use of RootLocusPlot? – bill s Jun 5 '15 at 12:43
• @bills, I could but RootLocusPlot still doesn't deal with transcendental equations. The problem is that, with the inclusion of delays there are infinitely many branches in the root locus. However, only the branches that might potentially cross to the RHP are interesting for stability purposes. I could use Padé approximants to the exponentials using polynomial quotients, but that would still be an approximation. – Bruno Pace Jun 7 '15 at 11:00