In Chapter 11 so far, stability has been characterized by pole locations in the \( s \)-plane (transfer functions) and by eigenvalues of \( \mathbf{A} \) (state space). Real systems, however, rarely match nominal parameters exactly. Let a parameter vector be \( \mathbf{p}\in\mathbb{R}^m \) and let poles be functions \( s_i(\mathbf{p}) \). Pole sensitivity quantifies how fast poles move under parameter perturbations.

The central question is: given a small perturbation \( \Delta \mathbf{p} \), how large can the pole shift be? A local (first-order) approximation is \( \Delta s_i \approx \nabla_{\mathbf{p}} s_i(\mathbf{p}_0)^\top \Delta \mathbf{p} \). If a stable pole lies close to the imaginary axis, even modest sensitivity may cause a stability loss.

flowchart TD
  A["Nominal model (transfer function or state space)"] --> B["Identify uncertain parameters p"]
  B --> C["Compute poles s_i(p0)"]
  C --> D["Compute sensitivities ds_i/dp_j (analytic or numeric)"]
  D --> E["Estimate pole shifts for small Delta p"]
  E --> F["Check stability: \nRe(s_i) < 0 for all i"]
  F -->|stable| G["Quantify robustness: \ndistance to imaginary axis and sensitivity"]
  F -->|unstable| H["Revise design or \nparameter tolerances"]
        

We will treat two representations:

• Transfer function (characteristic polynomial): poles are roots of \( \Delta(s,\mathbf{p})=0 \).
• State space: poles are eigenvalues of \( \mathbf{A}(\mathbf{p}) \).

Let \( \Delta(s,p) \) be a scalar characteristic equation depending smoothly on a scalar parameter \( p \). A pole \( s(p) \) satisfies \( \Delta(s(p),p)=0 \).

Assume \( s_0=s(p_0) \) is a simple root, meaning \( \frac{\partial \Delta}{\partial s}(s_0,p_0)\neq 0 \). Then \( s(p) \) is differentiable near \( p_0 \), and its derivative is:

\[ \frac{ds}{dp}(p_0) = - \frac{\frac{\partial \Delta}{\partial p}(s_0,p_0)}{\frac{\partial \Delta}{\partial s}(s_0,p_0)}. \]

Proof (differentiation of an implicit equation):

Differentiate \( \Delta(s(p),p)=0 \) with respect to \( p \):

\[ \frac{d}{dp}\Delta(s(p),p) = \frac{\partial \Delta}{\partial s}(s(p),p)\frac{ds}{dp} + \frac{\partial \Delta}{\partial p}(s(p),p) = 0. \]

Evaluate at \( p=p_0 \) and solve for \( \frac{ds}{dp}(p_0) \). The condition \( \frac{\partial \Delta}{\partial s}(s_0,p_0)\neq 0 \) ensures division is valid.

Key consequence (high sensitivity near multiple roots): If the pole is repeated at \( (s_0,p_0) \), then \( \frac{\partial \Delta}{\partial s}(s_0,p_0)=0 \), so the local derivative formula breaks down. This is a precise mathematical reason why repeated poles are typically very fragile under perturbations.

Consider a strictly proper LTI transfer function \( G(s) = \frac{N(s)}{D(s)} \) with characteristic polynomial \( \Delta(s,\mathbf{p}) = D(s,\mathbf{p}) \). Write \( D(s,\mathbf{p}) = \sum_{k=0}^{n} a_k(\mathbf{p}) s^{k} \) where \( a_n(\mathbf{p})\neq 0 \). A pole \( s_i \) satisfies \( D(s_i,\mathbf{p})=0 \).

3.1 Sensitivity of a pole to a coefficient

Let \( a_\ell \) vary (treating other coefficients fixed). Then \( \frac{\partial D}{\partial a_\ell} = s^\ell \) and \( \frac{\partial D}{\partial s} = D'(s) \). The pole sensitivity to coefficient \( a_\ell \) is:

\[ \frac{\partial s_i}{\partial a_\ell} = -\frac{s_i^\ell}{D'(s_i)}. \]

This formula shows explicitly that the denominator \( D'(s_i) \) controls sensitivity. If \( D'(s_i) \) is small in magnitude (e.g., near repeated roots), sensitivity becomes large.

3.2 Sensitivity to a physical parameter

If coefficients depend on a physical parameter \( p \), then by the chain rule:

\[ \frac{ds_i}{dp} = - \frac{\frac{\partial D}{\partial p}(s_i,p)}{D'(s_i,p)} = - \frac{\sum_{k=0}^{n}\frac{da_k}{dp}(p)\, s_i^k}{\sum_{k=1}^{n}k a_k(p)s_i^{k-1}}. \]

3.3 Dimensionless (relative) sensitivity

To compare sensitivities across parameters with different units, define relative sensitivity:

\[ S_{s_i}^{p} \triangleq \frac{p}{s_i}\frac{ds_i}{dp}, \quad \text{(defined when } s_i \neq 0\text{)}. \]

3.4 Worked example: second-order system with uncertain damping

Consider the canonical second-order characteristic polynomial \( D(s) = s^2 + 2\zeta\omega_n s + \omega_n^2 \), where \( \zeta \) may vary due to damping uncertainty. Poles are \( s_{1,2} = -\zeta\omega_n \pm \omega_n\sqrt{\zeta^2 - 1} \). For the underdamped case \( 0 < \zeta < 1 \), write \( s_{1,2} = -\zeta\omega_n \pm j\omega_n\sqrt{1-\zeta^2} \).

Differentiate with respect to \( \zeta \):

\[ \frac{ds_{1,2}}{d\zeta} = -\omega_n \pm j\omega_n \frac{d}{d\zeta}\sqrt{1-\zeta^2} = -\omega_n \mp j\omega_n \frac{\zeta}{\sqrt{1-\zeta^2}}. \]

The imaginary-part sensitivity grows without bound as \( \zeta \rightarrow 1^- \) because \( \sqrt{1-\zeta^2} \rightarrow 0 \). This matches the intuition that the transition to critical damping is sensitive: small damping changes can shift oscillatory frequency significantly near \( \zeta=1 \).

Let the state matrix depend smoothly on a parameter \( p \): \( \mathbf{A}(p)\in\mathbb{R}^{n\times n} \). Poles are eigenvalues \( \lambda_i(p) \) solving \( \det(\lambda\mathbf{I}-\mathbf{A}(p))=0 \).

4.1 First-order eigenvalue sensitivity formula

Suppose \( \lambda \) is a simple eigenvalue of \( \mathbf{A}(p_0) \). Let \( \mathbf{v} \) be a right eigenvector and \( \mathbf{w} \) a left eigenvector:

\[ \mathbf{A}\mathbf{v} = \lambda \mathbf{v}, \qquad \mathbf{w}^\top \mathbf{A} = \lambda \mathbf{w}^\top, \qquad \mathbf{w}^\top \mathbf{v} \neq 0. \]

Then the eigenvalue derivative is

\[ \frac{d\lambda}{dp} = \frac{\mathbf{w}^\top \left(\frac{d\mathbf{A}}{dp}\right)\mathbf{v}}{\mathbf{w}^\top \mathbf{v}}. \]

Proof (eigenvector calculus for a simple eigenvalue):

Differentiate \( \mathbf{A}(p)\mathbf{v}(p) = \lambda(p)\mathbf{v}(p) \):

\[ \frac{d\mathbf{A}}{dp}\mathbf{v} + \mathbf{A}\frac{d\mathbf{v}}{dp} = \frac{d\lambda}{dp}\mathbf{v} + \lambda \frac{d\mathbf{v}}{dp}. \]

Left-multiply by \( \mathbf{w}^\top \) and use \( \mathbf{w}^\top\mathbf{A}=\lambda\mathbf{w}^\top \) to cancel the terms with \( \frac{d\mathbf{v}}{dp} \):

\[ \mathbf{w}^\top \frac{d\mathbf{A}}{dp}\mathbf{v} + \underbrace{\mathbf{w}^\top\mathbf{A}\frac{d\mathbf{v}}{dp}}_{=\lambda\mathbf{w}^\top\frac{d\mathbf{v}}{dp}} = \frac{d\lambda}{dp}\mathbf{w}^\top\mathbf{v} + \lambda\mathbf{w}^\top\frac{d\mathbf{v}}{dp}. \]

Cancel \( \lambda\mathbf{w}^\top\frac{d\mathbf{v}}{dp} \) on both sides to obtain \( \mathbf{w}^\top \frac{d\mathbf{A}}{dp}\mathbf{v} = \frac{d\lambda}{dp}\mathbf{w}^\top\mathbf{v} \), hence the stated formula (since \( \mathbf{w}^\top\mathbf{v}\neq 0 \) for a simple eigenvalue).

4.2 Conditioning: why nearly defective eigenvalues are fragile

The numerator \( \mathbf{w}^\top \left(\frac{d\mathbf{A}}{dp}\right)\mathbf{v} \) is problem-dependent, but the denominator \( \mathbf{w}^\top\mathbf{v} \) reflects how well-conditioned the eigenvalue is. If \( |\mathbf{w}^\top\mathbf{v}| \) is small, the eigenvalue is sensitive.

A classic global bound (for diagonalizable \( \mathbf{A} \)) is given by the Bauer–Fike theorem. Let \( \mathbf{A} = \mathbf{V}\mathbf{\Lambda}\mathbf{V}^{-1} \) and perturb it as \( \mathbf{A}+\mathbf{E} \). Then every eigenvalue \( \mu \) of \( \mathbf{A}+\mathbf{E} \) satisfies

\[ \min_i |\mu - \lambda_i| \le \kappa(\mathbf{V})\, \|\mathbf{E}\|, \quad \kappa(\mathbf{V}) \triangleq \|\mathbf{V}\|\,\|\mathbf{V}^{-1}\|. \]

Proof sketch (resolvent argument):

If \( \mu \) is an eigenvalue of \( \mathbf{A}+\mathbf{E} \), then \( \det(\mu\mathbf{I}-(\mathbf{A}+\mathbf{E}))=0 \), so \( \mu\mathbf{I}-\mathbf{A}-\mathbf{E} \) is singular. If \( \mu\mathbf{I}-\mathbf{A} \) were invertible, we could write \( \mu\mathbf{I}-\mathbf{A}-\mathbf{E} = (\mu\mathbf{I}-\mathbf{A})(\mathbf{I}-(\mu\mathbf{I}-\mathbf{A})^{-1}\mathbf{E}) \). Singularity implies \( \mathbf{I}-(\mu\mathbf{I}-\mathbf{A})^{-1}\mathbf{E} \) is singular, hence \( \|(\mu\mathbf{I}-\mathbf{A})^{-1}\mathbf{E}\|\ge 1 \). Using diagonalization gives \( \|(\mu\mathbf{I}-\mathbf{A})^{-1}\| \le \|\mathbf{V}\|\,\|\mathbf{V}^{-1}\|\max_i \frac{1}{|\mu-\lambda_i|} \), yielding the bound.

Interpretation: if \( \mathbf{A} \) is nearly non-diagonalizable, \( \kappa(\mathbf{V}) \) can be large, and eigenvalues can move substantially even when \( \|\mathbf{E}\| \) is small.

We do not yet use frequency-domain robustness tools (those begin in Chapter 12). Instead, we define a purely pole-based robustness margin:

\[ \alpha \triangleq \min_i \big(-\operatorname{Re}(s_i)\big). \]

If the nominal system is stable, then \( \alpha > 0 \). A small \( \alpha \) means poles are close to the imaginary axis; stability is then easily lost under perturbations.

Under a small parameter perturbation \( \Delta \mathbf{p} \), a first-order approximation gives:

\[ s_i(\mathbf{p}_0 + \Delta \mathbf{p}) \approx s_i(\mathbf{p}_0) + \sum_{j=1}^{m}\frac{\partial s_i}{\partial p_j}(\mathbf{p}_0)\,\Delta p_j. \]

A sufficient local condition to preserve stability is that the real part remains negative. For each pole,

\[ \operatorname{Re}(s_i(\mathbf{p}_0)) + \operatorname{Re}\!\left(\sum_{j=1}^{m}\frac{\partial s_i}{\partial p_j}\Delta p_j\right) < 0. \]

In a conservative scalar-bound form, if \( |\Delta p_j| \le \delta_j \), then using the triangle inequality:

\[ \operatorname{Re}(s_i(\mathbf{p}_0)) + \sum_{j=1}^{m}\left|\operatorname{Re}\!\left(\frac{\partial s_i}{\partial p_j}\right)\right|\delta_j < 0 \quad \Rightarrow \quad \text{pole } i \text{ stays in } \operatorname{Re}(s)<0. \]

This connects engineering tolerances \( \delta_j \) to a mathematically justified robustness check.

flowchart TD
  P["Poles s_i at nominal parameters"] --> M["Compute alpha = min(-Re(s_i))"]
  M --> S["Compute sensitivities ds_i/dp_j"]
  S --> B["Bound Re shift: sum |Re(ds_i/dp_j)| * delta_j"]
  B --> C["Compare with alpha"]
  C -->|Bound < alpha| OK["Stable under small variations (local guarantee)"]
  C -->|Bound >= alpha| RISK["High risk: poles may cross Re(s)=0"]
  RISK --> FIX["Reduce uncertainty or move poles left (design change)"]
        

Important limitation: this is a local, first-order reasoning tool. If parameter variations are not small, or if poles are near multiple roots / non-diagonalizable cases, higher-order effects can dominate.

When analytic derivatives are difficult, finite differences provide a practical estimate. For a scalar parameter \( p \) and pole \( s_i(p) \):

\[ \frac{ds_i}{dp}(p_0) \approx \frac{s_i(p_0 + h) - s_i(p_0 - h)}{2h}, \quad h \text{ small}. \]

Matching poles across perturbations: as \( p \) changes, root ordering can permute. A robust approach matches each perturbed pole to the closest nominal pole in the complex plane (minimum-distance assignment).

Sanity check: compare the finite-difference estimate to analytic sensitivity when available. Large disagreement can indicate: (i) \( h \) too large (nonlinear effects), (ii) \( h \) too small (roundoff), or (iii) a near-multiple pole / ill-conditioned eigenstructure.

Libraries used: numpy, scipy, and the Control Systems library control. The code below: (i) computes analytic pole sensitivity for a polynomial, (ii) computes eigenvalue sensitivity using left/right eigenvectors, (iii) verifies both with finite differences.

import numpy as np

# ----------------------------
# (A) Transfer-function polynomial sensitivity
# D(s,p) = a2 s^2 + a1(p) s + a0, with a1(p)=a1_0 + p
# ----------------------------
def poly_roots(coeffs_desc):
    # coeffs_desc: [a_n, ..., a_0]
    return np.roots(coeffs_desc)

def poly_Dprime_at_root(coeffs_desc, s):
    # derivative of D(s)=sum a_k s^(n-k)
    n = len(coeffs_desc) - 1
    # D(s)=a0*s^n + a1*s^(n-1) + ... + a_n
    d = 0.0 + 0.0j
    for i, a in enumerate(coeffs_desc[:-1]):  # skip constant term
        power = n - i
        d += a * power * (s ** (power - 1))
    return d

# Example polynomial: s^2 + a1 s + a0 with a1 = 2.0 + p, a0 = 5.0
a0 = 5.0
a1_0 = 2.0
p0 = 0.3
a1 = a1_0 + p0
coeffs = [1.0, a1, a0]  # s^2 + a1 s + a0

s_nom = poly_roots(coeffs)

# Analytic sensitivity ds/da1 = -s / D'(s)
sens_da1 = []
for s in s_nom:
    Dp = poly_Dprime_at_root(coeffs, s)
    sens_da1.append(-(s**1) / Dp)
sens_da1 = np.array(sens_da1)

# Since a1 = a1_0 + p, da1/dp = 1, so ds/dp = ds/da1
sens_dp_analytic = sens_da1

# Finite difference check
h = 1e-5
coeffs_plus  = [1.0, a1_0 + (p0 + h), a0]
coeffs_minus = [1.0, a1_0 + (p0 - h), a0]
s_plus = poly_roots(coeffs_plus)
s_minus = poly_roots(coeffs_minus)

# Match roots by nearest neighbor (works well for small h)
def match_by_nearest(ref, arr):
    out = np.zeros_like(ref)
    used = np.zeros(len(arr), dtype=bool)
    for i, r in enumerate(ref):
        d = np.abs(arr - r)
        d[used] = np.inf
        j = np.argmin(d)
        out[i] = arr[j]
        used[j] = True
    return out

s_plus_m  = match_by_nearest(s_nom, s_plus)
s_minus_m = match_by_nearest(s_nom, s_minus)
sens_dp_fd = (s_plus_m - s_minus_m) / (2*h)

print("Nominal poles:", s_nom)
print("Analytic ds/dp:", sens_dp_analytic)
print("FD ds/dp:", sens_dp_fd)

# ----------------------------
# (B) State-space eigenvalue sensitivity
# A(p) = A0 + p*A1
# dlambda/dp = w^T A1 v / (w^T v)
# ----------------------------
A0 = np.array([[0.0, 1.0],
               [-5.0, -2.0]])
A1 = np.array([[0.0, 0.0],
               [0.0, -1.0]])  # increases damping magnitude as p increases

A = A0 + p0*A1

# Right eigenvectors
lam, V = np.linalg.eig(A)

# Left eigenvectors are right eigenvectors of A^T
lamL, WL = np.linalg.eig(A.T)

# Match left eigenvectors to right eigenvalues by nearest eigenvalue
W = np.zeros_like(WL)
for i in range(len(lam)):
    j = np.argmin(np.abs(lamL - lam[i]))
    W[:, i] = WL[:, j]

dl_dp = []
for i in range(len(lam)):
    v = V[:, i]
    w = W[:, i]
    num = (w.T @ (A1 @ v))
    den = (w.T @ v)
    dl_dp.append(num / den)
dl_dp = np.array(dl_dp)

# Finite difference check
A_plus = A0 + (p0 + h)*A1
A_minus = A0 + (p0 - h)*A1
lam_plus, _ = np.linalg.eig(A_plus)
lam_minus, _ = np.linalg.eig(A_minus)

lam_plus_m = match_by_nearest(lam, lam_plus)
lam_minus_m = match_by_nearest(lam, lam_minus)
dl_dp_fd = (lam_plus_m - lam_minus_m) / (2*h)

print("State-space eigenvalues:", lam)
print("Analytic dlambda/dp:", dl_dp)
print("FD dlambda/dp:", dl_dp_fd)
      

Interpretation guideline: inspect \( \operatorname{Re}\!\left(\frac{ds}{dp}\right) \) (or \( \operatorname{Re}\!\left(\frac{d\lambda}{dp}\right) \)) because it directly indicates whether poles tend to drift toward the instability boundary \( \operatorname{Re}(s)=0 \) under increasing \( p \).

This example uses the Eigen library to compute eigenvalues/eigenvectors and evaluates \( \frac{d\lambda}{dp}=\frac{\mathbf{w}^\top \mathbf{A}_1 \mathbf{v}}{\mathbf{w}^\top \mathbf{v}} \) for \( \mathbf{A}(p)=\mathbf{A}_0+p\mathbf{A}_1 \). (Complex arithmetic is required in general.)

#include <iostream>
#include <complex>
#include <Eigen/Dense>
#include <Eigen/Eigenvalues>

int main() {
    using namespace Eigen;
    using cd = std::complex<double>;

    double p0 = 0.3;

    Matrix2d A0;
    A0 << 0.0, 1.0,
          -5.0, -2.0;

    Matrix2d A1;
    A1 << 0.0, 0.0,
           0.0, -1.0;

    Matrix2d A = A0 + p0 * A1;

    ComplexEigenSolver<Matrix2d> ces(A);
    VectorXcd lam = ces.eigenvalues();
    MatrixXcd V = ces.eigenvectors(); // right eigenvectors: A v = lambda v

    // Left eigenvectors via eigenvectors of A^T:
    ComplexEigenSolver<Matrix2d> cesL(A.transpose());
    VectorXcd lamL = cesL.eigenvalues();
    MatrixXcd WL = cesL.eigenvectors(); // columns are right eigenvectors of A^T

    // Match left eigenvectors by closest eigenvalue
    MatrixXcd W(2,2);
    for (int i = 0; i < 2; ++i) {
        double best = 1e300;
        int bestj = 0;
        for (int j = 0; j < 2; ++j) {
            double d = std::abs(lamL(j) - lam(i));
            if (d < best) { best = d; bestj = j; }
        }
        W.col(i) = WL.col(bestj);
    }

    // Sensitivity: dlambda/dp = w^T A1 v / (w^T v)
    Matrix2cd A1c = A1.cast<cd>();

    for (int i = 0; i < 2; ++i) {
        Vector2cd v = V.col(i);
        Vector2cd w = W.col(i);
        cd num = w.transpose() * (A1c * v);
        cd den = w.transpose() * v;
        cd dldp = num / den;

        std::cout << "lambda[" << i << "] = " << lam(i) << "\n";
        std::cout << "dlambda/dp = " << dldp << "\n";
    }

    // Finite-difference verification
    double h = 1e-6;
    Matrix2d Aplus  = A0 + (p0 + h)*A1;
    Matrix2d Aminus = A0 + (p0 - h)*A1;

    ComplexEigenSolver<Matrix2d> cesP(Aplus), cesM(Aminus);
    VectorXcd lp = cesP.eigenvalues();
    VectorXcd lm = cesM.eigenvalues();

    // Nearest matching (small h assumed)
    auto matchNearest = [](const VectorXcd& ref, const VectorXcd& arr) {
        VectorXcd out(ref.size());
        std::vector<bool> used(arr.size(), false);
        for (int i = 0; i < ref.size(); ++i) {
            double best = 1e300;
            int bestj = 0;
            for (int j = 0; j < arr.size(); ++j) {
                if (used[j]) continue;
                double d = std::abs(arr(j) - ref(i));
                if (d < best) { best = d; bestj = j; }
            }
            out(i) = arr(bestj);
            used[bestj] = true;
        }
        return out;
    };

    VectorXcd lp_m = matchNearest(lam, lp);
    VectorXcd lm_m = matchNearest(lam, lm);
    VectorXcd fd = (lp_m - lm_m) / (2.0*h);

    std::cout << "FD dlambda/dp:\n" << fd << "\n";
    return 0;
}
      

Notes: For larger matrices, use consistent norms and consider scaling. Ill-conditioned eigenvectors can make \( \mathbf{w}^\top\mathbf{v} \) small and the result numerically sensitive, which itself is a useful robustness signal.

This example uses EJML (Efficient Java Matrix Library). Complex eigenvalues are common, so use EJML’s decomposition that supports complex results. The left/right eigenvector matching logic mirrors the Python/C++ approach.

import org.ejml.data.*;
import org.ejml.dense.row.CommonOps_DDRM;
import org.ejml.dense.row.factory.DecompositionFactory_DDRM;
import org.ejml.interfaces.decomposition.EigenDecomposition_F64;

public class EigenSensitivityDemo {

    // Helper: compute A(p) = A0 + p*A1
    static DMatrixRMaj A_of_p(DMatrixRMaj A0, DMatrixRMaj A1, double p) {
        DMatrixRMaj A = A0.copy();
        DMatrixRMaj tmp = A1.copy();
        CommonOps_DDRM.scale(p, tmp);
        CommonOps_DDRM.addEquals(A, tmp);
        return A;
    }

    public static void main(String[] args) {
        double p0 = 0.3;
        double h  = 1e-6;

        DMatrixRMaj A0 = new DMatrixRMaj(new double[][]{
            {0.0, 1.0},
            {-5.0, -2.0}
        });

        DMatrixRMaj A1 = new DMatrixRMaj(new double[][]{
            {0.0, 0.0},
            {0.0, -1.0}
        });

        DMatrixRMaj A = A_of_p(A0, A1, p0);

        // Eigen-decomposition (may be complex)
        EigenDecomposition_F64<DMatrixRMaj> eig = DecompositionFactory_DDRM.eig(2, true);
        eig.decompose(A);

        // For 2x2, demonstrate derivative via finite differences (robust and simple).
        // Analytic left/right eigenvector formula is possible but requires careful complex-vector handling.

        // Finite difference for eigenvalues
        DMatrixRMaj Aplus  = A_of_p(A0, A1, p0 + h);
        DMatrixRMaj Aminus = A_of_p(A0, A1, p0 - h);

        EigenDecomposition_F64<DMatrixRMaj> eigP = DecompositionFactory_DDRM.eig(2, true);
        EigenDecomposition_F64<DMatrixRMaj> eigM = DecompositionFactory_DDRM.eig(2, true);
        eigP.decompose(Aplus);
        eigM.decompose(Aminus);

        Complex_F64[] lam  = new Complex_F64[2];
        Complex_F64[] lp   = new Complex_F64[2];
        Complex_F64[] lm   = new Complex_F64[2];

        for (int i = 0; i < 2; i++) {
            lam[i] = eig.getEigenvalue(i);
            lp[i]  = eigP.getEigenvalue(i);
            lm[i]  = eigM.getEigenvalue(i);
        }

        // Nearest matching by complex distance
        Complex_F64[] lp_m = new Complex_F64[2];
        Complex_F64[] lm_m = new Complex_F64[2];
        boolean[] usedP = new boolean[2];
        boolean[] usedM = new boolean[2];

        for (int i = 0; i < 2; i++) {
            int bestP = 0; double bestPd = 1e300;
            int bestM = 0; double bestMd = 1e300;

            for (int j = 0; j < 2; j++) {
                if (!usedP[j]) {
                    double d = Math.hypot(lp[j].real - lam[i].real, lp[j].imaginary - lam[i].imaginary);
                    if (d < bestPd) { bestPd = d; bestP = j; }
                }
                if (!usedM[j]) {
                    double d = Math.hypot(lm[j].real - lam[i].real, lm[j].imaginary - lam[i].imaginary);
                    if (d < bestMd) { bestMd = d; bestM = j; }
                }
            }
            lp_m[i] = lp[bestP]; usedP[bestP] = true;
            lm_m[i] = lm[bestM]; usedM[bestM] = true;
        }

        for (int i = 0; i < 2; i++) {
            double dRe = (lp_m[i].real - lm_m[i].real) / (2.0*h);
            double dIm = (lp_m[i].imaginary - lm_m[i].imaginary) / (2.0*h);
            System.out.printf("lambda[%d] = %.6f + j%.6f%n", i, lam[i].real, lam[i].imaginary);
            System.out.printf("FD dlambda/dp = %.6f + j%.6f%n", dRe, dIm);
        }
    }
}
      

If you require the analytic formula in Java, implement complex left/right eigenvectors explicitly and evaluate \( \frac{\mathbf{w}^\top \mathbf{A}_1 \mathbf{v}}{\mathbf{w}^\top \mathbf{v}} \). For coursework, the finite-difference method above is an acceptable verification tool when paired with careful pole matching.

MATLAB functions used: roots, eig. If available, Control System Toolbox functions such as tf and pole can help, but the core computation here remains polynomial/eigenvalue-based (aligned with Chapter 11 content).

%% (A) Polynomial sensitivity: D(s,p) = s^2 + (a1_0 + p)*s + a0
a0   = 5.0;
a1_0 = 2.0;
p0   = 0.3;
h    = 1e-6;

a1 = a1_0 + p0;
coeffs = [1, a1, a0];    % descending powers
s_nom = roots(coeffs);

% Analytic ds/da1 = -s / D'(s); D'(s)=2s + a1
Dprime = @(s) 2*s + a1;
ds_dp_analytic = -(s_nom) ./ Dprime(s_nom);

% Finite difference
s_plus  = roots([1, a1_0 + (p0+h), a0]);
s_minus = roots([1, a1_0 + (p0-h), a0]);

% Match by nearest (simple for 2nd order)
[~, idxP] = sort(abs(s_plus  - s_nom.'), 1);
[~, idxM] = sort(abs(s_minus - s_nom.'), 1);
s_plus_m  = s_plus(idxP(1,:)).';
s_minus_m = s_minus(idxM(1,:)).';
ds_dp_fd  = (s_plus_m - s_minus_m)/(2*h);

disp("Nominal poles:"); disp(s_nom);
disp("Analytic ds/dp:"); disp(ds_dp_analytic);
disp("FD ds/dp:"); disp(ds_dp_fd);

%% (B) State-space sensitivity for A(p)=A0 + p*A1 (analytic via left/right eigenvectors)
A0 = [0 1; -5 -2];
A1 = [0 0;  0 -1];
A  = A0 + p0*A1;

[V,D] = eig(A);          % right eigenvectors V, eigenvalues diag(D)
lam = diag(D);

% left eigenvectors from A' (transpose)
[W,DL] = eig(A');
lamL = diag(DL);

% match W columns to lam by nearest eigenvalue
Wm = zeros(size(W));
for i=1:length(lam)
    [~,j] = min(abs(lamL - lam(i)));
    Wm(:,i) = W(:,j);
end

dl_dp_analytic = zeros(size(lam));
for i=1:length(lam)
    v = V(:,i);
    w = Wm(:,i);
    dl_dp_analytic(i) = (w.'*(A1*v)) / (w.'*v);   % transpose, not conjugate transpose
end

% finite difference check
lam_p = eig(A0 + (p0+h)*A1);
lam_m = eig(A0 + (p0-h)*A1);

% nearest matching
lam_p_m = lam_p;
lam_m_m = lam_m;
for i=1:length(lam)
    [~,jp] = min(abs(lam_p - lam(i)));
    [~,jm] = min(abs(lam_m - lam(i)));
    lam_p_m(i) = lam_p(jp);
    lam_m_m(i) = lam_m(jm);
end
dl_dp_fd = (lam_p_m - lam_m_m)/(2*h);

disp("Eigenvalues:"); disp(lam);
disp("Analytic dlambda/dp:"); disp(dl_dp_analytic);
disp("FD dlambda/dp:"); disp(dl_dp_fd);

%% (C) Parameter sweep (robustness visualization in time-domain sense)
pvals = linspace(p0-0.5, p0+0.5, 101);
alphas = zeros(size(pvals));
for k=1:length(pvals)
    Ak = A0 + pvals(k)*A1;
    lk = eig(Ak);
    alphas(k) = min(-real(lk)); % alpha > 0 indicates stable
end

% Plot (optional): figure; plot(pvals, alphas); grid on; xlabel('p'); ylabel('alpha');
      

Simulink note (optional if available): you can parameterize a block (e.g., a Gain, Mass, Damping) and linearize a Simulink model around an operating point to extract \( \mathbf{A}(p) \), then apply the eigenvalue sensitivity workflow above. This remains consistent with Chapter 11 because it ultimately uses eigenvalues of the linearized state matrix.

Mathematica is well-suited for symbolic differentiation of the implicit characteristic equation and for verifying results numerically.

(* Polynomial example: D(s,p) = s^2 + (a1_0 + p) s + a0 *)
a0 = 5;
a10 = 2;
p0 = 0.3;

D[s_, p_] := s^2 + (a10 + p) s + a0;

(* Solve for poles at p0 *)
sol = Solve[D[s, p0] == 0, s];
poles = s /. sol;

(* Implicit sensitivity ds/dp = - (D_p / D_s) evaluated at the pole *)
Ds[s_, p_] := D[D[s, p], s];
Dp[s_, p_] := D[D[s, p], p];

sens = Table[
  - (Dp[poles[[i]], p0]/Ds[poles[[i]], p0]),
  {i, Length[poles]}
];

poles // N
sens  // N

(* Finite difference verification *)
h = 10^-6;
polesPlus  = s /. Solve[D[s, p0 + h] == 0, s];
polesMinus = s /. Solve[D[s, p0 - h] == 0, s];

(* Match by nearest pole *)
matchNearest[ref_, arr_] := Module[{out = ConstantArray[0, Length[ref]], used = ConstantArray[False, Length[arr]], j},
  Do[
    j = First@Ordering[Map[If[used[[#]], Infinity, Abs[arr[[#]] - ref[[i]]]] &, Range[Length[arr]]], 1];
    out[[i]] = arr[[j]]; used[[j]] = True;,
    {i, Length[ref]}
  ];
  out
];

pPlusM  = matchNearest[poles, polesPlus];
pMinusM = matchNearest[poles, polesMinus];
sensFD  = (pPlusM - pMinusM)/(2 h);

sensFD // N

(* State-space eigenvalue sensitivity: A(p)=A0 + p*A1 *)
A0 = { {0, 1}, {-5, -2} };
A1 = { {0, 0}, {0, -1} };
A[p_] := A0 + p A1;

(* Numeric eigenvalues around p0 *)
lam0 = Eigenvalues[A[p0]];

(* Finite difference dlambda/dp *)
lamP = Eigenvalues[A[p0 + h]];
lamM = Eigenvalues[A[p0 - h]];

lamP_m = matchNearest[lam0, lamP];
lamM_m = matchNearest[lam0, lamM];
dlFD = (lamP_m - lamM_m)/(2 h);

lam0 // N
dlFD // N
      

For analytic eigenvalue derivatives in Mathematica, one can compute left/right eigenvectors and evaluate \( \frac{\mathbf{w}^\top \mathbf{A}_1 \mathbf{v}}{\mathbf{w}^\top \mathbf{v}} \). Numerically, finite differences remain a robust verification method when paired with careful eigenvalue matching.

Problem 1 (Implicit root sensitivity): Let \( \Delta(s,p) \) be continuously differentiable and suppose \( s_0 \) is a simple root at \( p_0 \): \( \Delta(s_0,p_0)=0 \) and \( \frac{\partial \Delta}{\partial s}(s_0,p_0)\neq 0 \). Derive a formula for \( \frac{ds}{dp}(p_0) \).

Solution: Differentiate \( \Delta(s(p),p)=0 \) with respect to \( p \):

\[ \frac{\partial \Delta}{\partial s}(s(p),p)\frac{ds}{dp} + \frac{\partial \Delta}{\partial p}(s(p),p) = 0. \]

Evaluate at \( (s_0,p_0) \) and solve:

\[ \frac{ds}{dp}(p_0) = -\frac{\frac{\partial \Delta}{\partial p}(s_0,p_0)}{\frac{\partial \Delta}{\partial s}(s_0,p_0)}. \]

Problem 2 (Coefficient sensitivity): Let \( D(s)=\sum_{k=0}^{n} a_k s^k \) and let \( s_i \) be a simple root. Show that \( \frac{\partial s_i}{\partial a_\ell} = -\frac{s_i^\ell}{D'(s_i)} \).

Solution: Use the implicit sensitivity with \( \Delta(s,\mathbf{a})=D(s) \) and parameter \( a_\ell \):

\[ \frac{\partial s_i}{\partial a_\ell} = -\frac{\frac{\partial D}{\partial a_\ell}(s_i)}{\frac{\partial D}{\partial s}(s_i)} = -\frac{s_i^\ell}{D'(s_i)}. \]

Problem 3 (Multiple roots and fragility): Suppose \( D(s,p)=(s-a(p))^2 \) where \( a(p) \) is smooth. Explain why the simple-root sensitivity formula fails and interpret this as fragility.

Solution: Here the pole is repeated: \( s=a(p) \) with multiplicity 2. Compute \( \frac{\partial D}{\partial s}=2(s-a(p)) \). At the root \( s=a(p) \) we get \( \frac{\partial D}{\partial s}=0 \), violating the condition required for \( \frac{ds}{dp}= -\frac{D_p}{D_s} \). This indicates that an arbitrarily small perturbation can split a repeated pole into two distinct poles whose locations are not governed by a single well-defined derivative at the nominal point—hence the pole configuration is fragile.

Problem 4 (Eigenvalue sensitivity): Let \( \mathbf{A}(p)=\mathbf{A}_0+p\mathbf{A}_1 \). Assume \( \lambda \) is a simple eigenvalue of \( \mathbf{A}(p_0) \) with right eigenvector \( \mathbf{v} \) and left eigenvector \( \mathbf{w} \) satisfying \( \mathbf{w}^\top\mathbf{v}\neq 0 \). Derive \( \frac{d\lambda}{dp} \).

Solution: Differentiate \( \mathbf{A}\mathbf{v}=\lambda\mathbf{v} \):

\[ \mathbf{A}_1\mathbf{v} + \mathbf{A}\frac{d\mathbf{v}}{dp} = \frac{d\lambda}{dp}\mathbf{v} + \lambda \frac{d\mathbf{v}}{dp}. \]

Left-multiply by \( \mathbf{w}^\top \) and use \( \mathbf{w}^\top\mathbf{A}=\lambda\mathbf{w}^\top \) to cancel:

\[ \mathbf{w}^\top\mathbf{A}_1\mathbf{v} = \frac{d\lambda}{dp}\mathbf{w}^\top\mathbf{v} \quad \Rightarrow \quad \frac{d\lambda}{dp}=\frac{\mathbf{w}^\top\mathbf{A}_1\mathbf{v}}{\mathbf{w}^\top\mathbf{v}}. \]

Problem 5 (Local robustness bound using sensitivity): Let the system have stable poles \( s_i(\mathbf{p}_0) \) and define \( \alpha=\min_i (-\operatorname{Re}(s_i(\mathbf{p}_0))) \). Suppose uncertainties satisfy \( |\Delta p_j| \le \delta_j \) and define \( b_i=\sum_{j=1}^{m}\left|\operatorname{Re}\!\left(\frac{\partial s_i}{\partial p_j}\right)\right|\delta_j \). Show that if \( b_i < -\operatorname{Re}(s_i(\mathbf{p}_0)) \) for all \( i \), then stability is preserved to first order.

Solution: First-order approximation yields

\[ \operatorname{Re}(s_i(\mathbf{p}_0+\Delta\mathbf{p})) \approx \operatorname{Re}(s_i(\mathbf{p}_0)) + \operatorname{Re}\!\left(\sum_{j=1}^m \frac{\partial s_i}{\partial p_j}\Delta p_j\right). \]

Using \( |\operatorname{Re}(z)| \le |z| \) and triangle inequality,

\[ \operatorname{Re}\!\left(\sum_{j=1}^m \frac{\partial s_i}{\partial p_j}\Delta p_j\right) \le \sum_{j=1}^m \left|\operatorname{Re}\!\left(\frac{\partial s_i}{\partial p_j}\right)\right|\,|\Delta p_j| \le b_i. \]

Therefore \( \operatorname{Re}(s_i(\mathbf{p}_0+\Delta\mathbf{p})) \lesssim \operatorname{Re}(s_i(\mathbf{p}_0)) + b_i \). If \( \operatorname{Re}(s_i(\mathbf{p}_0)) + b_i < 0 \) for all \( i \), then all poles remain in \( \operatorname{Re}(s)<0 \) under the first-order bound, implying local robustness.

We derived rigorous first-order pole sensitivity formulas for (i) transfer-function characteristic polynomials via implicit differentiation and (ii) state-space eigenvalues via left/right eigenvectors. We proved the core identities, highlighted the mathematical reason repeated poles are fragile (vanishing derivative denominator), and introduced a purely pole-based time-domain robustness check using the stability margin \( \alpha=\min_i(-\operatorname{Re}(s_i)) \) together with sensitivity-based bounds. These tools prepare you to reason about how modeling uncertainty affects stability before introducing frequency-response robustness concepts in Chapter 12.