In previous lessons of this chapter, we introduced sampling, difference equations, the \( z \)-transform, and continuous-to-discrete conversion. We now study what happens after a discrete-time model is obtained: whether trajectories decay, remain bounded, oscillate, or diverge, and how fast the transient disappears.

For a discrete-time state-space model \( x[k+1] = A x[k] + B u[k] \), \( y[k] = C x[k] + D u[k] \), the central questions are:

• Internal (state) stability: do states converge to zero for zero input?
• BIBO stability: does every bounded input produce a bounded output?
• Time response shaping: how do poles determine speed, oscillation, and overshoot?

flowchart TD
  A["Discrete model from prior lessons"] --> B["State form / Transfer function"]
  B --> C["Internal stability test"]
  B --> D["BIBO stability test"]
  C --> E["Eigenvalues / Lyapunov"]
  D --> F["Impulse response summability"]
  E --> G["Natural response decay rate"]
  F --> H["Forced response boundedness"]
  G --> I["Time-response metrics"]
  H --> I
        

Repeated substitution gives the exact solution of the discrete-time state equation for \( k \ge 0 \):

\[ x[k] = A^k x[0] + \sum_{j=0}^{k-1} A^{k-1-j} B u[j]. \]

The first term \( A^k x[0] \) is the natural response, and the summation term is the forced response. Therefore, internal stability is entirely governed by the asymptotic behavior of \( A^k \).

For a SISO LTI transfer function \( G(z) = Y(z)/U(z) \) with impulse response \( g[k] \), the convolution form is

\[ y[k] = \sum_{j=0}^{k} g[j]u[k-j]. \]

This representation is the discrete-time counterpart of continuous-time convolution and is the basis of BIBO stability proofs.

Scalar prototype. Consider \( x[k+1] = a x[k] \). The solution is \( x[k] = a^k x[0] \), so asymptotic stability holds if and only if \( |a| < 1 \).

\[ |x[k]| = |a|^k |x[0]| \xrightarrow[k\to\infty]{} 0 \quad \Longleftrightarrow \quad |a| < 1. \]

If \( |a| = 1 \), the response is marginal (constant or oscillatory); if \( |a| > 1 \), the response diverges.

Matrix case. For \( x[k+1]=Ax[k] \), the system is asymptotically stable if and only if every eigenvalue \( \lambda_i(A) \) lies strictly inside the unit disk:

\[ \rho(A) = \max_i |\lambda_i(A)| < 1. \]

This is the discrete-time analog of the continuous-time criterion \( \Re\{\lambda_i\} < 0 \). In the z-plane, stability means all poles are inside the unit circle.

Proof sketch (Jordan-form argument):

Write \( A = VJV^{-1} \) (or real Jordan/Schur form). Then \( A^k = VJ^kV^{-1} \). Each Jordan block associated with eigenvalue \( \lambda \) contributes terms of the form \( k^m \lambda^k \). If \( |\lambda| < 1 \) for all blocks, then \( k^m |\lambda|^k \to 0 \), hence \( A^k \to 0 \). Conversely, if any \( |\lambda| > 1 \) (or a defective block on \( |\lambda|=1 \)), \( A^k \) cannot converge to zero.

A powerful matrix test avoids explicit eigenvalue computation. Let \( Q = Q^T \succ 0 \). If there exists \( P = P^T \succ 0 \) such that

\[ A^T P A - P = -Q, \]

then the origin of \( x[k+1]=Ax[k] \) is asymptotically stable.

Define the Lyapunov function \( V[k] = x[k]^T P x[k] \). Along trajectories:

\[ \Delta V = V[k+1]-V[k] = x[k]^T (A^T P A - P)x[k] = -x[k]^T Q x[k] < 0 \]

for all nonzero \( x[k] \). Therefore \( V[k] \) decreases strictly and \( x[k] \to 0 \).

When \( \rho(A) < 1 \), the unique positive definite solution is given by the convergent series (also called the discrete Lyapunov or Stein series):

\[ P = \sum_{j=0}^{\infty} (A^j)^T Q A^j. \]

Substituting into the left-hand side gives a telescoping cancellation, which proves the equation exactly.

For a causal SISO LTI system with impulse response \( g[k] \), a sufficient and necessary condition for BIBO stability is absolute summability:

\[ \sum_{k=0}^{\infty} |g[k]| < \infty. \]

Proof (bounded input \(\Rightarrow\) bounded output): If \( |u[k]| \le M \) for all \( k \), then

\[ |y[k]| = \left|\sum_{j=0}^{k} g[j]u[k-j]\right| \le \sum_{j=0}^{k} |g[j]|\,|u[k-j]| \le M\sum_{j=0}^{\infty} |g[j]|. \]

Hence \( y[k] \) is bounded uniformly in \( k \). For rational LTI systems with no pole-zero cancellation in the external transfer function, this is equivalent to all poles satisfying \( |p_i| < 1 \).

Final Value Theorem (discrete-time). Under standard conditions (all poles of \( (1-z^{-1})Y(z) \) strictly inside the unit circle), the steady-state value is

\[ \lim_{k\to\infty} y[k] = \lim_{z\to 1} (1-z^{-1})Y(z). \]

This theorem is very useful for step-response steady-state calculations and error checks.

Consider a complex-conjugate pole pair \( z_{1,2} = r e^{\pm j\theta} \) with \( 0 < r < 1 \). The modal contribution to the natural response has the form

\[ y_n[k] = r^k \left( A_c \cos(k\theta) + B_c \sin(k\theta) \right). \]

Therefore:

• The oscillation period in samples is approximately \( N_p = 2\pi/\theta \).
• The envelope decays geometrically as \( r^k \).
• For sampling period \( T \), an equivalent decay time constant is \( \tau_{eq} = -T / \ln(r) \).
• The damped oscillation frequency in rad/s is \( \omega_d = \theta/T \).

For a first-order difference equation \( y[k+1] = a y[k] + b u[k] \) with unit-step input and \( y[0]=0 \), the step response is

\[ y[k] = \frac{b}{1-a}\left(1-a^k\right), \quad |a| < 1. \]

If the steady-state value is \( y_\infty = b/(1-a) \), the absolute error is \( |y_\infty-y[k]| = |y_\infty|\,|a|^k \). Thus a settling-index estimate for a relative tolerance \( \varepsilon \) is

\[ k_s \ge \frac{\ln(\varepsilon)}{\ln(|a|)}. \]

Because \( \ln(|a|) < 0 \) in the stable case, the right-hand side is positive.

From exact discretization (Lesson 4), eigenvalues map according to \( z = e^{sT} \). Let \( s = \sigma + j\omega \); then

\[ |z| = |e^{sT}| = e^{\sigma T}. \]

Hence:

\[ \sigma < 0 \; \Longleftrightarrow \; |z| < 1, \qquad \sigma = 0 \; \Longleftrightarrow \; |z| = 1. \]

This explains why the left half of the s-plane maps to the interior of the unit circle in the z-plane. It also clarifies how sampling preserves stability under exact ZOH discretization of a stable continuous-time LTI system.

In engineering practice, stability analysis and response characterization should be done systematically. The following workflow applies to both transfer-function and state-space descriptions.

flowchart TD
  A["Start with G(z) or (A,B,C,D)"] --> B["Check pole/eigenvalue locations"]
  B --> C["If all magnitudes < 1: \nstable"]
  B --> D["If any magnitude > 1: \nunstable"]
  B --> E["If magnitude = 1: \ninspect multiplicity / \nmarginal case"]
  C --> F["Compute natural + forced response"]
  F --> G["Extract rise, settling, overshoot, oscillation"]
  G --> H["Validate final value and steady-state error"]
  H --> I["Refine sampling time or redesign model"]
        

The following script tests Schur stability, constructs a Lyapunov certificate using the discrete Lyapunov series, simulates the state-space model, and computes basic step-response metrics.

Chapter16_Lesson5.py

# Chapter16_Lesson5.py
"""
Discrete-time stability and time response demo for System Dynamics (Chapter 16, Lesson 5)

Topics:
- Schur stability (eigenvalues inside unit circle)
- Natural/forced response
- Discrete-time Lyapunov equation (iterative solver)
- Step response metrics (rise time, settling time, overshoot)
"""
import numpy as np

def spectral_radius(A: np.ndarray) -> float:
    vals = np.linalg.eigvals(A)
    return float(np.max(np.abs(vals)))

def is_schur_stable(A: np.ndarray, tol: float = 1e-12) -> bool:
    vals = np.linalg.eigvals(A)
    return bool(np.all(np.abs(vals) < 1.0 - tol))

def simulate_state_space(A, B, C, D, u, x0):
    """
    x[k+1] = A x[k] + B u[k]
    y[k]   = C x[k] + D u[k]
    """
    A = np.array(A, dtype=float)
    B = np.array(B, dtype=float).reshape((-1, 1))
    C = np.array(C, dtype=float).reshape((1, -1))
    D = np.array(D, dtype=float).reshape((1, 1))
    u = np.asarray(u, dtype=float).reshape((-1,))
    x = np.array(x0, dtype=float).reshape((-1, 1))

    n = len(u)
    xs = [x.copy()]
    ys = []
    for k in range(n):
        yk = C @ x + D * u[k]
        ys.append(float(yk[0, 0]))
        x = A @ x + B * u[k]
        xs.append(x.copy())
    X = np.hstack(xs).T
    y = np.array(ys)
    return X, y

def dlyap_iterative(A, Q, max_iter=20000, tol=1e-12):
    """
    Solve P = A^T P A + Q by summation:
    P = sum_{k=0}^{inf} (A^k)^T Q A^k, valid for Schur-stable A.
    """
    A = np.array(A, dtype=float)
    Q = np.array(Q, dtype=float)
    n = A.shape[0]
    P = np.zeros((n, n))
    Ak = np.eye(n)
    for _ in range(max_iter):
        term = Ak.T @ Q @ Ak
        P_next = P + term
        if np.linalg.norm(term, ord='fro') < tol:
            return 0.5 * (P_next + P_next.T)
        P = P_next
        Ak = Ak @ A
    raise RuntimeError("dlyap_iterative did not converge. Check Schur stability or increase max_iter.")

def step_metrics(y, final_value=None, tol=0.02):
    y = np.asarray(y, dtype=float)
    N = y.size
    if final_value is None:
        final_value = float(np.mean(y[max(0, N//2):]))
    if abs(final_value) < 1e-12:
        return {"final": final_value, "rise_index": None, "settling_index": None, "overshoot_pct": 0.0}

    # rise time index (10% to 90%)
    y10 = 0.1 * final_value
    y90 = 0.9 * final_value
    idx10 = next((i for i, yi in enumerate(y) if yi >= y10), None)
    idx90 = next((i for i, yi in enumerate(y) if yi >= y90), None)
    rise_index = None if idx10 is None or idx90 is None else idx90 - idx10

    band = tol * abs(final_value)
    settling_index = None
    for i in range(N):
        if np.all(np.abs(y[i:] - final_value) <= band):
            settling_index = i
            break

    overshoot = (np.max(y) - final_value) / abs(final_value) * 100.0
    return {
        "final": final_value,
        "rise_index": rise_index,
        "settling_index": settling_index,
        "overshoot_pct": overshoot,
    }

def main():
    # A stable second-order discrete-time model
    # Poles are approximately 0.82*exp(±j0.28)
    A = np.array([[1.5770, -0.6724],
                  [1.0000,  0.0000]])
    B = np.array([[1.0],
                  [0.0]])
    C = np.array([[0.0676, 0.0604]])
    D = np.array([[0.0]])

    print("A =\n", A)
    eigA = np.linalg.eigvals(A)
    print("Eigenvalues(A) =", eigA)
    print("Spectral radius =", spectral_radius(A))
    print("Schur stable?   =", is_schur_stable(A))

    # Lyapunov certificate
    Q = np.eye(2)
    P = dlyap_iterative(A, Q)
    M = A.T @ P @ A - P
    print("\nDiscrete Lyapunov matrix P (Q=I):\n", P)
    print("Check A^T P A - P (should be -Q):\n", M)

    # Time response to step input
    N = 80
    u = np.ones(N)
    x0 = np.array([0.0, 0.0])
    X, y = simulate_state_space(A, B, C, D, u, x0)

    # Natural response from nonzero initial condition
    Xn, yn = simulate_state_space(A, B, C, D, np.zeros(N), np.array([1.0, -0.4]))

    metrics = step_metrics(y, final_value=1.0)  # DC gain chosen ~1
    print("\nStep-response metrics (sample counts):")
    for k, v in metrics.items():
        print(f"  {k}: {v}")

    print("\nFirst 12 samples of unit-step response:")
    print(np.array2string(y[:12], precision=5, separator=", "))

    print("\nFirst 12 samples of natural response (u=0, x0=[1,-0.4]):")
    print(np.array2string(yn[:12], precision=5, separator=", "))

if __name__ == "__main__":
    main()

      

The C++ and Java codes implement a compact second-order example without external libraries. They include: (i) second-order Jury stability test, (ii) eigenvalue magnitude check, and (iii) step-response simulation.

Chapter16_Lesson5.cpp

// Chapter16_Lesson5.cpp
// Discrete-time stability and time response demo (C++17, no external libraries)
#include <iostream>
#include <array>
#include <vector>
#include <cmath>
#include <iomanip>

struct ComplexPair {
    double real;
    double imag;
};

double det2(const std::array<std::array<double,2>,2>& A) {
    return A[0][0]*A[1][1] - A[0][1]*A[1][0];
}

double trace2(const std::array<std::array<double,2>,2>& A) {
    return A[0][0] + A[1][1];
}

ComplexPair eig_pair_2x2(const std::array<std::array<double,2>,2>& A, double& lambda1_abs, double& lambda2_abs) {
    double tr = trace2(A);
    double det = det2(A);
    double disc = tr*tr - 4.0*det;

    if (disc >= 0.0) {
        double s = std::sqrt(disc);
        double l1 = 0.5*(tr + s);
        double l2 = 0.5*(tr - s);
        lambda1_abs = std::fabs(l1);
        lambda2_abs = std::fabs(l2);
        return {0.0, 0.0};
    } else {
        double re = 0.5*tr;
        double im = 0.5*std::sqrt(-disc);
        double mag = std::sqrt(re*re + im*im);
        lambda1_abs = mag;
        lambda2_abs = mag;
        return {re, im};
    }
}

// Jury test for a2 z^2 + a1 z + a0 (with a2 > 0)
bool jury_second_order(double a2, double a1, double a0) {
    if (a2 <= 0.0) return false;
    // Necessary and sufficient for second-order Schur stability:
    // p(1) > 0, p(-1) > 0 for normalized orientation, and |a0| < a2
    // for polynomial a2 z^2 + a1 z + a0
    double p1 = a2 + a1 + a0;
    double pm1 = a2 - a1 + a0;
    return (std::fabs(a0) < a2) && (p1 > 0.0) && (pm1 > 0.0);
}

int main() {
    std::array<std::array<double,2>,2> A{ { {1.5770, -0.6724},
                                           {1.0000,  0.0000} } };
    std::array<double,2> B{ {1.0, 0.0} };
    std::array<double,2> C{ {0.0676, 0.0604} };
    double D = 0.0;

    // Characteristic polynomial of A: z^2 - tr(A) z + det(A)
    double a2 = 1.0;
    double a1 = -trace2(A);
    double a0 = det2(A);

    double mag1=0.0, mag2=0.0;
    ComplexPair cp = eig_pair_2x2(A, mag1, mag2);

    std::cout << std::fixed << std::setprecision(6);
    std::cout << "Characteristic polynomial: z^2 + (" << a1 << ") z + (" << a0 << ")\n";
    std::cout << "Jury stable (2nd order test)? " << (jury_second_order(a2, a1, a0) ? "true" : "false") << "\n";
    std::cout << "Eigenvalue magnitudes: " << mag1 << ", " << mag2 << "\n";
    if (cp.imag != 0.0) {
        std::cout << "Complex conjugate poles: " << cp.real << " +/- j" << cp.imag << "\n";
    }

    // Step response simulation
    std::array<double,2> x{ {0.0, 0.0} };
    const int N = 60;
    std::vector<double> y;
    y.reserve(N);

    for (int k=0; k<N; ++k) {
        double u = 1.0;
        double yk = C[0]*x[0] + C[1]*x[1] + D*u;
        y.push_back(yk);

        std::array<double,2> xn{};
        xn[0] = A[0][0]*x[0] + A[0][1]*x[1] + B[0]*u;
        xn[1] = A[1][0]*x[0] + A[1][1]*x[1] + B[1]*u;
        x = xn;
    }

    std::cout << "\nFirst 15 step-response samples:\n";
    for (int k=0; k<15; ++k) {
        std::cout << "k=" << std::setw(2) << k << "  y=" << y[k] << "\n";
    }

    // Approximate settling sample for 2% band around final value 1.0
    int k_settle = -1;
    for (int k=0; k<N; ++k) {
        bool inside = true;
        for (int j=k; j<N; ++j) {
            if (std::fabs(y[j] - 1.0) > 0.02) { inside = false; break; }
        }
        if (inside) { k_settle = k; break; }
    }
    std::cout << "\n2% settling sample index = " << k_settle << "\n";
    return 0;
}

      

Chapter16_Lesson5.java

// Chapter16_Lesson5.java
// Discrete-time stability and step response demo in Java (no external libraries)
import java.util.Arrays;

public class Chapter16_Lesson5 {

    static class Complex {
        double re, im;
        Complex(double re, double im) { this.re = re; this.im = im; }
        double abs() { return Math.hypot(re, im); }
        public String toString() {
            if (Math.abs(im) < 1e-12) return String.format("%.6f", re);
            return String.format("%.6f %s j%.6f", re, (im >= 0 ? "+" : "-"), Math.abs(im));
        }
    }

    static double trace2(double[][] A) { return A[0][0] + A[1][1]; }
    static double det2(double[][] A) { return A[0][0]*A[1][1] - A[0][1]*A[1][0]; }

    static Complex[] eig2x2(double[][] A) {
        double tr = trace2(A);
        double det = det2(A);
        double disc = tr*tr - 4.0*det;
        if (disc >= 0) {
            double s = Math.sqrt(disc);
            return new Complex[] {
                new Complex(0.5*(tr + s), 0.0),
                new Complex(0.5*(tr - s), 0.0)
            };
        } else {
            double re = 0.5*tr;
            double im = 0.5*Math.sqrt(-disc);
            return new Complex[] {
                new Complex(re, im),
                new Complex(re, -im)
            };
        }
    }

    // Jury criterion specialized to second-order polynomial z^2 + a1 z + a0
    static boolean jurySecondOrder(double a1, double a0) {
        return (Math.abs(a0) < 1.0) && (1.0 + a1 + a0 > 0.0) && (1.0 - a1 + a0 > 0.0);
    }

    public static void main(String[] args) {
        double[][] A = {
            {1.5770, -0.6724},
            {1.0000,  0.0000}
        };
        double[] B = {1.0, 0.0};
        double[] C = {0.0676, 0.0604};
        double D = 0.0;

        double a1 = -trace2(A);
        double a0 = det2(A);

        Complex[] eigs = eig2x2(A);
        System.out.printf("Characteristic polynomial: z^2 + (%.6f) z + (%.6f)%n", a1, a0);
        System.out.println("Jury stable? " + jurySecondOrder(a1, a0));
        System.out.println("Eigenvalues: " + eigs[0] + ", " + eigs[1]);
        System.out.printf("Magnitudes: %.6f, %.6f%n", eigs[0].abs(), eigs[1].abs());

        int N = 70;
        double[] x = {0.0, 0.0};
        double[] y = new double[N];

        for (int k = 0; k < N; k++) {
            double u = 1.0;
            y[k] = C[0]*x[0] + C[1]*x[1] + D*u;

            double[] xn = new double[2];
            xn[0] = A[0][0]*x[0] + A[0][1]*x[1] + B[0]*u;
            xn[1] = A[1][0]*x[0] + A[1][1]*x[1] + B[1]*u;
            x = xn;
        }

        System.out.println("\nFirst 12 step samples:");
        for (int k = 0; k < 12; k++) {
            System.out.printf("k=%2d  y=%.6f%n", k, y[k]);
        }

        // 2% settling index about final value 1
        int settle = -1;
        for (int k = 0; k < N; k++) {
            boolean ok = true;
            for (int j = k; j < N; j++) {
                if (Math.abs(y[j] - 1.0) > 0.02) { ok = false; break; }
            }
            if (ok) { settle = k; break; }
        }
        System.out.println("\n2% settling sample index = " + settle);

        // Optional: natural response from nonzero initial condition
        x = new double[]{1.0, -0.4};
        double[] yn = new double[12];
        for (int k = 0; k < 12; k++) {
            yn[k] = C[0]*x[0] + C[1]*x[1];
            double[] xn = new double[2];
            xn[0] = A[0][0]*x[0] + A[0][1]*x[1];
            xn[1] = A[1][0]*x[0] + A[1][1]*x[1];
            x = xn;
        }
        System.out.println("\nNatural response first 12 samples:");
        System.out.println(Arrays.toString(yn));
    }
}

      

The MATLAB script uses discrete transfer-function and state-space tools to check poles, solve the discrete Lyapunov equation, and compute step-response metrics. The Simulink exercise script programmatically builds a discrete transfer-function model. The Mathematica code reproduces the same ideas symbolically and numerically.

Chapter16_Lesson5.m

% Chapter16_Lesson5.m
% Discrete-time stability and time response (MATLAB / Simulink-oriented)
% Requires Control System Toolbox for tf/step; dlyap is in Control System Toolbox.

clear; clc;

Ts = 0.1;  % sampling period [s]

% A stable 2nd-order discrete transfer function with DC gain = 1
num = [0.0676 0.0604];
den = [1 -1.5770 0.6724];
Gz = tf(num, den, Ts);

disp('Discrete transfer function G(z):');
Gz

p = pole(Gz);
fprintf('Poles of G(z):\n');
disp(p.');
fprintf('Magnitudes of poles:\n');
disp(abs(p).');

if all(abs(p) < 1)
    fprintf('System is Schur stable (all poles inside unit circle).\n');
else
    fprintf('System is NOT Schur stable.\n');
end

% State-space form and discrete Lyapunov equation
sys_ss = ss(Gz);
A = sys_ss.A; B = sys_ss.B; C = sys_ss.C; D = sys_ss.D;
Q = eye(size(A));
P = dlyap(A', Q);  % solves A' P A - P + Q = 0
disp('Lyapunov matrix P from A'' P A - P = -Q:');
disp(P);

% Step response and basic metrics
N = 100;
k = 0:N-1;
t = k * Ts;
[y, t] = step(Gz, t);

figure;
stairs(t, y, 'LineWidth', 1.3); grid on;
xlabel('Time (s)'); ylabel('y[k]');
title('Unit-Step Response of Discrete-Time System');

info = stepinfo(y, t, 1.0, 'SettlingTimeThreshold', 0.02);
disp('Step response metrics:');
disp(info);

% Natural response (zero input, nonzero initial condition)
x0 = [1; -0.4];
u = zeros(N,1);
[y_nat, t_nat, x_nat] = lsim(sys_ss, u, t, x0);

figure;
stairs(t_nat, y_nat, 'LineWidth', 1.3); grid on;
xlabel('Time (s)'); ylabel('y_n[k]');
title('Natural Response (u=0, x_0=[1;-0.4])');

% Optional Simulink note:
% You can build an equivalent model with:
%   Step -> Discrete Transfer Fcn -> Scope
% and set numerator/denominator to the vectors above with sample time Ts.

      

Chapter16_Lesson5_Ex1.m

% Chapter16_Lesson5_Ex1.m
% Programmatically builds a simple Simulink model:
% Step -> Discrete Transfer Fcn -> Scope
% This is useful for sampled-data visualization in the lesson.

modelName = 'Chapter16_Lesson5_Simulink';
if bdIsLoaded(modelName)
    close_system(modelName, 0);
end
new_system(modelName);
open_system(modelName);

% Add blocks
add_block('simulink/Sources/Step', [modelName '/Step'], ...
    'Position', [40 80 70 110], ...
    'Time', '0', 'Before', '0', 'After', '1');

add_block('simulink/Discrete/Discrete Transfer Fcn', [modelName '/Gz'], ...
    'Position', [150 70 300 120]);

add_block('simulink/Sinks/Scope', [modelName '/Scope'], ...
    'Position', [380 75 430 115]);

% Parameters
set_param([modelName '/Gz'], ...
    'Numerator', '[0.0676 0.0604]', ...
    'Denominator', '[1 -1.5770 0.6724]', ...
    'SampleTime', '0.1');

% Connect lines
add_line(modelName, 'Step/1', 'Gz/1', 'autorouting', 'on');
add_line(modelName, 'Gz/1', 'Scope/1', 'autorouting', 'on');

% Simulation settings
set_param(modelName, 'StopTime', '10');
save_system(modelName);
disp(['Created model: ' modelName]);

% To run:
% sim(modelName);

      

Chapter16_Lesson5.nb

(* Chapter16_Lesson5.nb *)
(* Wolfram Language code for discrete-time stability and time response *)

ClearAll["Global`*"];

(* Define discrete-time state-space model *)
A = { {1.5770, -0.6724}, {1.0, 0.0} };
B = { {1.0}, {0.0} };
C = { {0.0676, 0.0604} };
D = { {0.0} };

eig = Eigenvalues[A];
Print["Eigenvalues(A) = ", eig];
Print["Magnitudes = ", Abs[eig]];
Print["Schur stable? ", And @@ Thread[Abs[eig] < 1]];

(* Discrete Lyapunov equation: A^T P A - P == -Q *)
Q = IdentityMatrix[2];
p11 =.; p12=.; p22=.;
P = { {p11, p12}, {p12, p22} };
eqs = Flatten[Transpose[A].P.A - P + Q == ConstantArray[0, {2, 2}]];
solP = Solve[eqs, {p11, p12, p22}][[1]];
Psol = P /. solP // N;
Print["P = ", MatrixForm[Psol]];
Print["Check A^T P A - P = ", MatrixForm[(Transpose[A].Psol.A - Psol) // N]];

(* Simulate step response *)
Nsteps = 80;
u[k_] := 1.0;
x[0] = {0.0, 0.0};
y[k_] := (C.{x[k]}[[1]] + D[[1,1]]*u[k]) // First;

Do[
  x[k + 1] = (A.{x[k]}[[1]] + B[[All,1]]*u[k]) // N
, {k, 0, Nsteps - 1}];

yvals = Table[y[k], {k, 0, Nsteps - 1}] // N;
Print["First 12 step samples = ", Take[yvals, 12]];

ListLinePlot[yvals,
  PlotRange -> All,
  AxesLabel -> {"k", "y[k]"},
  PlotLabel -> "Unit-Step Response (Discrete-Time)",
  InterpolationOrder -> 0
]

(* Natural response *)
Clear[xn, yn];
xn[0] = {1.0, -0.4};
yn[k_] := (C.{xn[k]}[[1]]) // First;
Do[
  xn[k + 1] = (A.{xn[k]}[[1]]) // N
, {k, 0, Nsteps - 1}];
ynvals = Table[yn[k], {k, 0, 20}] // N;
Print["Natural response first 10 samples = ", Take[ynvals, 10]];

      

For a second-order polynomial \( p(z) = z^2 + a_1 z + a_0 \), the Jury conditions are equivalent to all roots lying inside the unit circle:

\[ |a_0| < 1, \qquad 1 + a_1 + a_0 > 0, \qquad 1 - a_1 + a_0 > 0. \]

These inequalities are the discrete-time analog of Routh-Hurwitz tests for continuous-time second-order polynomials and are very convenient for parameter-range studies.

A common second-order discrete denominator can be parameterized directly by pole radius and angle:

\[ (z-re^{j\theta})(z-re^{-j\theta}) = z^2 - 2r\cos(\theta)z + r^2. \]

Thus the denominator coefficients are \( a_1 = -2r\cos(\theta) \) and \( a_0 = r^2 \). Stability is immediate when \( 0 < r < 1 \).

Problem 1 (Scalar Stability): Consider the difference equation \( x[k+1] = 0.92\,x[k] \) with \( x[0]=5 \). Determine asymptotic stability and estimate the smallest \( k \) such that \( |x[k]| \le 0.05 \).

Solution: Since \( |0.92| < 1 \), the origin is asymptotically stable.

\[ |x[k]| = 5(0.92)^k \le 0.05 \quad \Longleftrightarrow \quad (0.92)^k \le 0.01. \]

\[ k \ge \frac{\ln(0.01)}{\ln(0.92)} \approx 55.26. \]

Therefore the smallest integer is \( k=56 \).

Problem 2 (Matrix Stability and Lyapunov Certificate): Let \( A = \begin{bmatrix} 0.6 & 0.2 \\ -0.1 & 0.7 \end{bmatrix} \). Verify Schur stability and explain how a discrete Lyapunov equation can certify it.

Solution: The characteristic polynomial is

\[ p(z)=z^2 - 1.3z + 0.44. \]

The roots are complex conjugates with magnitude \( \sqrt{0.44} \approx 0.6633 < 1 \), so the matrix is Schur stable. For any chosen \( Q=Q^T\succ 0 \), the equation \( A^T P A - P = -Q \) has a unique positive definite solution \( P \), and the Lyapunov function \( V[k]=x[k]^TPx[k] \) proves asymptotic stability.

Problem 3 (Jury Test): Determine whether the polynomial \( p(z)=z^2-1.4z+0.52 \) is stable (all roots inside the unit circle).

Solution: Here \( a_1=-1.4 \), \( a_0=0.52 \). Check:

\[ |a_0| = 0.52 < 1, \qquad 1+a_1+a_0 = 1-1.4+0.52 = 0.12 > 0, \\ 1-a_1+a_0 = 1+1.4+0.52 = 2.92 > 0. \]

All Jury conditions hold, so the polynomial is Schur stable.

Problem 4 (Discrete Step Response Formula): For \( y[k+1]=0.8y[k]+0.1u[k] \), \( y[0]=0 \), and unit-step input \( u[k]=1 \), derive \( y[k] \) and its steady-state value.

Solution: The homogeneous part is \( y_h[k]=c(0.8)^k \). A constant particular solution \( y_p \) satisfies

\[ y_p = 0.8y_p + 0.1 \quad \Longrightarrow \quad y_p = \frac{0.1}{0.2} = 0.5. \]

Using \( y[0]=0 \), the complete response is

\[ y[k] = 0.5\left(1-0.8^k\right). \]

Therefore \( \lim_{k\to\infty} y[k] = 0.5 \).

Problem 5 (Pole Geometry and Oscillation): A discrete-time mode has poles \( z_{1,2} = 0.9 e^{\pm j0.3} \) and sampling time \( T=0.05 \) s. Estimate (i) the oscillation period in samples, (ii) the equivalent damped frequency in rad/s, and (iii) the decay time constant.

Solution:

\[ N_p = \frac{2\pi}{0.3} \approx 20.94 \text{ samples}, \qquad \omega_d = \frac{0.3}{0.05} = 6 \text{ rad/s}. \]

\[ \tau_{eq} = -\frac{T}{\ln(0.9)} = -\frac{0.05}{\ln(0.9)} \approx 0.475 \text{ s}. \]

So the response oscillates with period about 21 samples and decays with an envelope time constant of approximately 0.475 s.

We established the core stability and response results for discrete-time LTI systems: internal stability is characterized by eigenvalues (or poles) inside the unit circle, BIBO stability is characterized by absolute summability of the impulse response, and time-response features such as settling and oscillation are directly encoded in pole radius and angle. We also linked these ideas to the discrete Lyapunov equation and provided multi-language implementations for analysis and simulation workflows.