In continuous-time modeling, system evolution is described by differential equations. After sampling (studied in Lesson 1), the dynamics are naturally represented on an integer index \( k \in \mathbb{Z}_{\ge 0} \). The central objects are:

\( y[k] \) (scalar output sequence), \( u[k] \) (input sequence), and \( \mathbf{x}[k] \in \mathbb{R}^n \) (state vector sequence).

A difference equation relates current and past samples, while a discrete-time state-space model updates the internal state one sample at a time. These are discrete analogs of ODEs and continuous-time state-space equations from earlier chapters.

\[ \text{Difference-equation view:}\quad y[k+n] + a_1 y[k+n-1] + \cdots + a_n y[k] = \\ b_0 u[k] + b_1 u[k-1] + \cdots + b_m u[k-m] \]

\[ \text{State-space view:}\quad \mathbf{x}[k+1] = \mathbf{A}\mathbf{x}[k] + \mathbf{B}\mathbf{u}[k],\qquad \mathbf{y}[k] = \mathbf{C}\mathbf{x}[k] + \mathbf{D}\mathbf{u}[k] \]

One of the core goals of this lesson is to show that these two descriptions are equivalent (for linear finite-dimensional systems) under suitable state definitions.

The simplest causal first-order discrete-time dynamic model is

\[ y[k+1] = a\,y[k] + b\,u[k],\qquad k \ge 0,\qquad y[0]=y_0 \]

Here \( a \) controls memory/decay and \( b \) maps the input into the state/output recursion. For constant input \( u[k]=u_0 \), the equilibrium \( y^\star \) (if it exists) satisfies

\[ y^\star = a\,y^\star + b\,u_0 \quad \Rightarrow \quad (1-a)y^\star = b\,u_0 \quad \Rightarrow \quad y^\star = \frac{b}{1-a}u_0,\;\; a \neq 1 \]

The recursion is asymptotically stable (zero-input sense) when \( |a| < 1 \), because repeated multiplication by \( a \) drives the free response to zero.

Closed-form solution (first-order case):

\[ y[k] = a^k y_0 + \sum_{i=0}^{k-1} a^{k-1-i} b\,u[i],\qquad k \ge 1 \]

Proof by iteration:

Starting from \( y[1]=a y_0 + b u[0] \), then

\[ y[2] = a y[1] + b u[1] = a(a y_0 + b u[0]) + b u[1] = a^2 y_0 + ab u[0] + b u[1] \]

Continuing this pattern yields

\[ y[k] = a^k y_0 + b\left(a^{k-1}u[0] + a^{k-2}u[1] + \cdots + u[k-1]\right) \]

which is exactly the stated summation form. This is the discrete-time analogue of convolution with an impulse response, but the formal transform machinery is deferred to Lesson 3.

A general linear constant-coefficient difference equation of order \( n \) is

\[ y[k+n] + a_1 y[k+n-1] + \cdots + a_n y[k] = b_0 u[k] + b_1 u[k-1] + \cdots + b_m u[k-m] \]

To generate a unique solution, we need \( n \) initial conditions, e.g. \( y[0],y[1],\ldots,y[n-1] \), and the input sequence values required by the right-hand side.

The structure mirrors the nth-order ODE model studied previously, but the evolution is indexed by sample number rather than continuous time.

Homogeneous and forced parts:

For \( u[k]=0 \), we obtain the homogeneous equation:

\[ y_h[k+n] + a_1 y_h[k+n-1] + \cdots + a_n y_h[k] = 0 \]

Looking for exponential-sequence solutions \( y_h[k]=\lambda^k \) gives the characteristic polynomial

\[ \lambda^n + a_1 \lambda^{n-1} + \cdots + a_n = 0 \]

The roots \( \lambda_i \) determine modal behavior, exactly as continuous-time eigenvalues determine ODE modes. For distinct roots, the homogeneous solution is a linear combination of \( \lambda_i^k \).

Stability insight:

The homogeneous response decays to zero if all characteristic roots satisfy \( |\lambda_i| < 1 \). This is the discrete-time stability region (inside the unit circle), which will be studied more systematically in later lessons.

For multi-state systems, the canonical linear discrete-time model is

\[ \mathbf{x}[k+1] = \mathbf{A}\mathbf{x}[k] + \mathbf{B}\mathbf{u}[k] \]

\[ \mathbf{y}[k] = \mathbf{C}\mathbf{x}[k] + \mathbf{D}\mathbf{u}[k] \]

Dimensions:

\[ \mathbf{x}[k] \in \mathbb{R}^n,\quad \mathbf{u}[k] \in \mathbb{R}^p,\quad \mathbf{y}[k] \in \mathbb{R}^q \]

\[ \mathbf{A} \in \mathbb{R}^{n\times n},\quad \mathbf{B} \in \mathbb{R}^{n\times p},\quad \mathbf{C} \in \mathbb{R}^{q\times n},\quad \mathbf{D} \in \mathbb{R}^{q\times p} \]

The matrix \( \mathbf{A} \) encodes one-sample state propagation. Its powers \( \mathbf{A}^k \) play the role of the state transition operator in discrete time.

A scalar difference equation can often be embedded into this vector form by choosing state variables as delayed outputs (and, if necessary, delayed inputs).

The exact solution of the linear discrete-time state equation

\[ \mathbf{x}[k+1] = \mathbf{A}\mathbf{x}[k] + \mathbf{B}\mathbf{u}[k],\qquad \mathbf{x}[0]=\mathbf{x}_0 \]

is

\[ \mathbf{x}[k] = \mathbf{A}^k \mathbf{x}_0 + \sum_{i=0}^{k-1}\mathbf{A}^{k-1-i}\mathbf{B}\mathbf{u}[i], \qquad k \ge 1 \]

and then the output follows from \( \mathbf{y}[k]=\mathbf{C}\mathbf{x}[k]+\mathbf{D}\mathbf{u}[k] \).

Proof by induction:

For \( k=1 \),

\[ \mathbf{x}[1] = \mathbf{A}\mathbf{x}_0 + \mathbf{B}\mathbf{u}[0] = \mathbf{A}^1\mathbf{x}_0 + \sum_{i=0}^{0}\mathbf{A}^{0}\mathbf{B}\mathbf{u}[i] \]

so the formula holds. Assume it is true for some \( k \). Then

\[ \mathbf{x}[k+1] = \mathbf{A}\mathbf{x}[k] + \mathbf{B}\mathbf{u}[k] \]

Substitute the induction hypothesis:

\[ \mathbf{x}[k+1] = \mathbf{A}\left( \mathbf{A}^k \mathbf{x}_0 + \sum_{i=0}^{k-1}\mathbf{A}^{k-1-i}\mathbf{B}\mathbf{u}[i] \right) + \mathbf{B}\mathbf{u}[k] \]

\[ \mathbf{x}[k+1] = \mathbf{A}^{k+1}\mathbf{x}_0 + \sum_{i=0}^{k-1}\mathbf{A}^{k-i}\mathbf{B}\mathbf{u}[i] + \mathbf{A}^{0}\mathbf{B}\mathbf{u}[k] \]

\[ \mathbf{x}[k+1] = \mathbf{A}^{k+1}\mathbf{x}_0 + \sum_{i=0}^{k}\mathbf{A}^{k-i}\mathbf{B}\mathbf{u}[i] \]

which is exactly the same formula with \( k \) replaced by \( k+1 \). Hence proved.

Consider the second-order difference equation

\[ y[k+2] - 1.5\,y[k+1] + 0.7\,y[k] = u[k] \]

Define the state vector using delayed outputs:

\[ \mathbf{x}[k] = \begin{bmatrix} x_1[k] \\ x_2[k] \end{bmatrix} = \begin{bmatrix} y[k+1] \\ y[k] \end{bmatrix} \]

Then

\[ x_1[k+1] = y[k+2] = 1.5\,y[k+1] - 0.7\,y[k] + u[k] = 1.5\,x_1[k] - 0.7\,x_2[k] + u[k] \]

\[ x_2[k+1] = y[k+1] = x_1[k] \]

Therefore the companion-form realization is

\[ \mathbf{x}[k+1] = \begin{bmatrix} 1.5 & -0.7 \\ 1 & 0 \end{bmatrix}\mathbf{x}[k] + \begin{bmatrix} 1 \\ 0 \end{bmatrix}u[k] \]

\[ y[k] = \begin{bmatrix} 0 & 1 \end{bmatrix}\mathbf{x}[k] \]

This construction generalizes to higher-order monic difference equations and is one of the most important bridges between input-output and state-space modeling.

flowchart TD
  A["Start from nth-order difference equation"] --> B["Choose state vector as delayed outputs"]
  B --> C["Rewrite highest-shift term y[k+n]"]
  C --> D["Assemble matrix A and input matrix B"]
  D --> E["Define output equation y[k] = C x[k] + D u[k]"]
  E --> F["Check dimensions and initial conditions"]
        

For constant input \( \mathbf{u}[k]=\mathbf{u}_0 \), an equilibrium state \( \mathbf{x}^\star \) satisfies

\[ \mathbf{x}^\star = \mathbf{A}\mathbf{x}^\star + \mathbf{B}\mathbf{u}_0 \]

\[ (\mathbf{I}-\mathbf{A})\mathbf{x}^\star = \mathbf{B}\mathbf{u}_0 \]

\[ \mathbf{x}^\star = (\mathbf{I}-\mathbf{A})^{-1}\mathbf{B}\mathbf{u}_0 \quad \text{if } \det(\mathbf{I}-\mathbf{A}) \neq 0 \]

Let the equilibrium error be \( \tilde{\mathbf{x} }[k] = \mathbf{x}[k]-\mathbf{x}^\star \). Subtracting equilibrium dynamics from the state equation gives

\[ \tilde{\mathbf{x} }[k+1] = \mathbf{A}\tilde{\mathbf{x} }[k] \]

Hence stability reduces to the zero-input dynamics governed by \( \mathbf{A} \).

Modal stability theorem (diagonalizable case):

If \( \mathbf{A}=\mathbf{V}\mathbf{\Lambda}\mathbf{V}^{-1} \) and every eigenvalue satisfies \( |\lambda_i| < 1 \), then \( \mathbf{A}^k → \mathbf{0} \) as \( k → \infty \), so the origin is asymptotically stable.

Proof:

\[ \mathbf{A}^k = \mathbf{V}\mathbf{\Lambda}^k\mathbf{V}^{-1} \]

Since \( \mathbf{\Lambda}^k \) is diagonal with entries \( \lambda_i^k \), and \( |\lambda_i| < 1 \Rightarrow \lambda_i^k → 0 \), we have \( \mathbf{\Lambda}^k → \mathbf{0} \). Multiplying by the constant matrices \( \mathbf{V} \) and \( \mathbf{V}^{-1} \) preserves convergence to zero, so \( \mathbf{A}^k → \mathbf{0} \).

Discrete Lyapunov equation (preview):

A deeper criterion states that if all eigenvalues of \( \mathbf{A} \) are strictly inside the unit circle, then for any symmetric positive definite \( \mathbf{Q} \), there exists a unique symmetric positive definite \( \mathbf{P} \) such that

\[ \mathbf{A}^\top \mathbf{P}\mathbf{A} - \mathbf{P} = -\mathbf{Q} \]

This is the discrete-time counterpart of continuous-time Lyapunov stability and becomes central in advanced control design.

A practical workflow for discrete-time modeling is shown below. This is useful when students are translating a physical system (or a sampled continuous-time model) into a recurrence suitable for computation.

flowchart TD
  A["Define sampled variables y[k], u[k], x[k]"] --> B["Write causal recurrence or state update"]
  B --> C["Collect coefficients into matrices A, B, C, D"]
  C --> D["Specify initial conditions x[0] or y[0...n-1]"]
  D --> E["Simulate sample-by-sample for k = 0,1,2,..."]
  E --> F["Check equilibrium and stability (eigenvalues / roots)"]
  F --> G["Validate against physical behavior and units"]
        

The following Python script simulates (i) a first-order scalar difference equation and (ii) a second-order system realized in companion-form state-space.

Chapter16_Lesson2.py

"""
Chapter16_Lesson2.py
Difference Equations and Discrete-Time State-Space Models
"""

import numpy as np
import matplotlib.pyplot as plt


def simulate_scalar_difference(a, b, u, y0):
    """
    Simulate y[k+1] = a*y[k] + b*u[k]
    Parameters
    ----------
    a, b : float
    u : array_like, length N
    y0 : float
    Returns
    -------
    y : ndarray, length N+1
        y[0] is the initial condition.
    """
    u = np.asarray(u, dtype=float)
    N = len(u)
    y = np.zeros(N + 1, dtype=float)
    y[0] = y0
    for k in range(N):
        y[k + 1] = a * y[k] + b * u[k]
    return y


def simulate_discrete_state_space(A, B, C, D, u_seq, x0):
    """
    Simulate:
        x[k+1] = A x[k] + B u[k]
        y[k]   = C x[k] + D u[k]
    for scalar input u[k].
    """
    A = np.asarray(A, dtype=float)
    B = np.asarray(B, dtype=float).reshape(-1, 1)
    C = np.asarray(C, dtype=float).reshape(1, -1)
    D = float(D)
    x = np.asarray(x0, dtype=float).reshape(-1, 1)

    u_seq = np.asarray(u_seq, dtype=float)
    N = len(u_seq)

    nx = A.shape[0]
    X = np.zeros((N + 1, nx), dtype=float)
    Y = np.zeros(N, dtype=float)

    X[0, :] = x.ravel()
    for k in range(N):
        u = float(u_seq[k])
        yk = (C @ x + D * u).item()
        Y[k] = yk
        x = A @ x + B * u
        X[k + 1, :] = x.ravel()

    return X, Y


def main():
    # Example 1: scalar first-order recursion
    N = 40
    u = np.ones(N)
    a, b = 0.85, 0.2
    y0 = 0.0
    y = simulate_scalar_difference(a, b, u, y0)

    # Example 2: second-order difference equation in companion-form state space
    # y[k+2] - 1.5 y[k+1] + 0.7 y[k] = u[k]
    # Let x[k] = [y[k+1], y[k]]^T
    # Then x[k+1] = [[1.5, -0.7], [1, 0]] x[k] + [[1], [0]] u[k]
    # and y[k] = [0, 1] x[k]
    A = np.array([[1.5, -0.7],
                  [1.0,  0.0]])
    B = np.array([[1.0],
                  [0.0]])
    C = np.array([[0.0, 1.0]])
    D = 0.0

    u2 = np.zeros(N)
    u2[0:10] = 1.0  # finite-duration pulse
    x0 = np.array([0.0, 0.0])

    X, Y = simulate_discrete_state_space(A, B, C, D, u2, x0)

    print("First 10 samples of scalar recursion y[k]:")
    for k in range(10):
        print(f"k={k:2d}, y={y[k]: .6f}")
    print()

    print("First 10 outputs of companion-form model:")
    for k in range(10):
        print(f"k={k:2d}, y={Y[k]: .6f}, x={X[k]}")

    # Plot results
    k1 = np.arange(N + 1)
    plt.figure(figsize=(8, 4))
    plt.stem(k1, y, basefmt=" ", use_line_collection=True)
    plt.xlabel("k")
    plt.ylabel("y[k]")
    plt.title("Scalar Difference Equation: y[k+1] = a y[k] + b u[k]")
    plt.grid(True)
    plt.tight_layout()
    plt.show()

    k2 = np.arange(N)
    plt.figure(figsize=(8, 4))
    plt.stem(k2, Y, basefmt=" ", use_line_collection=True)
    plt.xlabel("k")
    plt.ylabel("y[k]")
    plt.title("Discrete-Time State-Space Output (Companion Form)")
    plt.grid(True)
    plt.tight_layout()
    plt.show()


if __name__ == "__main__":
    main()
      

This implementation explicitly follows the sample-by-sample update equations, so students can directly see the recurrence structure rather than using a black-box solver.

The C++ version uses only the standard library and is suitable for embedded or high-performance contexts where explicit memory and loop control matter.

Chapter16_Lesson2.cpp

// Chapter16_Lesson2.cpp
// Difference Equations and Discrete-Time State-Space Models

#include <iostream>
#include <vector>
#include <array>
#include <iomanip>

std::vector<double> simulateScalarDifference(double a, double b,
                                             const std::vector<double>& u,
                                             double y0) {
    std::vector<double> y(u.size() + 1, 0.0);
    y[0] = y0;
    for (std::size_t k = 0; k < u.size(); ++k) {
        y[k + 1] = a * y[k] + b * u[k];
    }
    return y;
}

struct SimResult {
    std::vector<std::array<double, 2>> X; // states x[k]
    std::vector<double> Y;                // outputs y[k]
};

SimResult simulateDiscreteSS(const std::array<std::array<double, 2>, 2>& A,
                             const std::array<double, 2>& B,
                             const std::array<double, 2>& C,
                             double D,
                             const std::vector<double>& u,
                             std::array<double, 2> x0) {
    SimResult r;
    r.X.resize(u.size() + 1);
    r.Y.resize(u.size());
    r.X[0] = x0;

    auto x = x0;
    for (std::size_t k = 0; k < u.size(); ++k) {
        double uk = u[k];
        double yk = C[0] * x[0] + C[1] * x[1] + D * uk;
        r.Y[k] = yk;

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

        x = xnext;
        r.X[k + 1] = x;
    }
    return r;
}

int main() {
    const int N = 40;

    // Example 1: scalar recursion y[k+1] = a y[k] + b u[k]
    std::vector<double> u1(N, 1.0);
    double a = 0.85, b = 0.2, y0 = 0.0;
    auto y = simulateScalarDifference(a, b, u1, y0);

    std::cout << "Scalar recursion (first 10 samples)\n";
    for (int k = 0; k < 10; ++k) {
        std::cout << "k=" << std::setw(2) << k
                  << "  y=" << std::fixed << std::setprecision(6)
                  << y[k] << "\n";
    }
    std::cout << "\n";

    // Example 2: companion-form state-space for
    // y[k+2] - 1.5 y[k+1] + 0.7 y[k] = u[k]
    std::array<std::array<double, 2>, 2> A{ { {1.5, -0.7}, {1.0, 0.0} } };
    std::array<double, 2> B{ {1.0, 0.0} };
    std::array<double, 2> C{ {0.0, 1.0} };
    double D = 0.0;

    std::vector<double> u2(N, 0.0);
    for (int k = 0; k < 10; ++k) u2[k] = 1.0;

    std::array<double, 2> x0_ss{ {0.0, 0.0} };
    auto result = simulateDiscreteSS(A, B, C, D, u2, x0_ss);

    std::cout << "State-space output (first 10 samples)\n";
    for (int k = 0; k < 10; ++k) {
        std::cout << "k=" << std::setw(2) << k
                  << "  y=" << std::fixed << std::setprecision(6) << result.Y[k]
                  << "  x=[" << result.X[k][0] << ", " << result.X[k][1] << "]\n";
    }

    return 0;
}
      

The code structure is intentionally close to the mathematical recursion, which is useful for control-oriented software development and embedded prototyping.

The Java implementation is also written without external libraries, so students can focus on the recurrence logic and matrix-vector update process.

Chapter16_Lesson2.java

// Chapter16_Lesson2.java
// Difference Equations and Discrete-Time State-Space Models

public class Chapter16_Lesson2 {

    static double[] simulateScalarDifference(double a, double b, double[] u, double y0) {
        double[] y = new double[u.length + 1];
        y[0] = y0;
        for (int k = 0; k < u.length; k++) {
            y[k + 1] = a * y[k] + b * u[k];
        }
        return y;
    }

    static class SimResult {
        double[][] X;
        double[] Y;
        SimResult(int N, int nx) {
            X = new double[N + 1][nx];
            Y = new double[N];
        }
    }

    static SimResult simulateDiscreteSS(double[][] A, double[] B, double[] C, double D,
                                        double[] u, double[] x0) {
        int N = u.length;
        int nx = x0.length;
        SimResult r = new SimResult(N, nx);

        double[] x = new double[nx];
        System.arraycopy(x0, 0, x, 0, nx);
        System.arraycopy(x, 0, r.X[0], 0, nx);

        for (int k = 0; k < N; k++) {
            double uk = u[k];
            double yk = 0.0;
            for (int i = 0; i < nx; i++) {
                yk += C[i] * x[i];
            }
            yk += D * uk;
            r.Y[k] = yk;

            double[] xNext = new double[nx];
            for (int i = 0; i < nx; i++) {
                double sum = 0.0;
                for (int j = 0; j < nx; j++) {
                    sum += A[i][j] * x[j];
                }
                xNext[i] = sum + B[i] * uk;
            }

            x = xNext;
            System.arraycopy(x, 0, r.X[k + 1], 0, nx);
        }

        return r;
    }

    public static void main(String[] args) {
        int N = 40;

        // Example 1: scalar recursion
        double[] u1 = new double[N];
        for (int k = 0; k < N; k++) u1[k] = 1.0;
        double[] y = simulateScalarDifference(0.85, 0.2, u1, 0.0);

        System.out.println("Scalar recursion (first 10 samples)");
        for (int k = 0; k < 10; k++) {
            System.out.printf("k=%2d  y=% .6f%n", k, y[k]);
        }
        System.out.println();

        // Example 2: companion-form discrete state-space
        // y[k+2] - 1.5 y[k+1] + 0.7 y[k] = u[k]
        double[][] A = {
            {1.5, -0.7},
            {1.0,  0.0}
        };
        double[] B = {1.0, 0.0};
        double[] C = {0.0, 1.0};
        double D = 0.0;

        double[] u2 = new double[N];
        for (int k = 0; k < 10; k++) u2[k] = 1.0;
        double[] x0 = {0.0, 0.0};

        SimResult r = simulateDiscreteSS(A, B, C, D, u2, x0);

        System.out.println("State-space output (first 10 samples)");
        for (int k = 0; k < 10; k++) {
            System.out.printf("k=%2d  y=% .6f  x=[% .6f, % .6f]%n",
                              k, r.Y[k], r.X[k][0], r.X[k][1]);
        }
    }
}
      

Java is especially useful for students building educational simulation tools, GUIs, or integration with larger software systems.

MATLAB provides a natural matrix-based environment for discrete-time state-space simulation. The following script implements the same two examples.

Chapter16_Lesson2.m

% Chapter16_Lesson2.m
% Difference Equations and Discrete-Time State-Space Models

clear; clc; close all;

%% Example 1: Scalar first-order difference equation
% y(k+1) = a*y(k) + b*u(k)
N = 40;
u1 = ones(1, N);
a = 0.85;
b = 0.2;
y = zeros(1, N+1);
y(1) = 0.0; % y(0)

for k = 1:N
    y(k+1) = a*y(k) + b*u1(k);
end

disp('Scalar recursion (first 10 samples, y(k))');
for k = 1:10
    fprintf('k=%2d, y=% .6f\n', k-1, y(k));
end

figure;
stem(0:N, y, 'filled');
grid on;
xlabel('k');
ylabel('y[k]');
title('Scalar Difference Equation: y[k+1] = a y[k] + b u[k]');

%% Example 2: Companion-form state-space
% y(k+2) - 1.5 y(k+1) + 0.7 y(k) = u(k)
% Let x(k) = [y(k+1); y(k)]
A = [1.5, -0.7;
     1.0,  0.0];
B = [1.0;
     0.0];
C = [0.0, 1.0];
D = 0.0;

u2 = zeros(1, N);
u2(1:10) = 1.0;  % finite pulse

nx = size(A,1);
X = zeros(nx, N+1);
Y = zeros(1, N);
x = [0.0; 0.0];

X(:,1) = x;
for k = 1:N
    Y(k) = C*x + D*u2(k);
    x = A*x + B*u2(k);
    X(:,k+1) = x;
end

disp('State-space output (first 10 samples)');
for k = 1:10
    fprintf('k=%2d, y=% .6f, x=[% .6f, % .6f]\n', k-1, Y(k), X(1,k), X(2,k));
end

figure;
stem(0:N-1, Y, 'filled');
grid on;
xlabel('k');
ylabel('y[k]');
title('Discrete-Time State-Space Output (Companion Form)');

%% Optional toolbox check (if Control System Toolbox is available)
% sysd = ss(A, B, C, D, 1); % sample time Ts = 1
% [y_toolbox, t_toolbox, x_toolbox] = lsim(sysd, u2, 0:N-1, [0;0]);
      

Simulink setup (conceptual):

To model the same companion-form system in Simulink, use: a Discrete State-Space block with matrices \( \mathbf{A}, \mathbf{B}, \mathbf{C}, \mathbf{D} \) as above, sample time \( T_s \) (e.g., \( T_s=1 \) for normalized samples), and a pulse/step source as input.

This block-level representation is the standard bridge from equations to simulation diagrams in control engineering workflows.

The Mathematica code below implements the same recursions using explicit loops, so the sample-indexed dynamics remain transparent.

Chapter16_Lesson2.nb

(* Chapter16_Lesson2.nb *)
(* Difference Equations and Discrete-Time State-Space Models *)

ClearAll["Global`*"];

(* Example 1: Scalar recursion y[k+1] = a y[k] + b u[k] *)
simulateScalarDifference[a_, b_, u_List, y0_] := Module[
  {n = Length[u], y, k},
  y = ConstantArray[0.0, n + 1];
  y[[1]] = y0; (* y[0] stored at index 1 *)
  For[k = 1, k <= n, k++,
    y[[k + 1]] = a y[[k]] + b u[[k]];
  ];
  y
];

(* Example 2: Discrete state-space simulation
   x[k+1] = A.x[k] + B u[k]
   y[k]   = C.x[k] + D u[k]
*)
simulateDiscreteSS[A_, B_, C_, D_, u_List, x0_List] := Module[
  {n = Length[u], nx = Length[x0], X, Y, x, k},
  X = ConstantArray[0.0, {n + 1, nx}];
  Y = ConstantArray[0.0, n];
  x = N[x0];
  X[[1]] = x;
  For[k = 1, k <= n, k++,
    Y[[k]] = C.x + D u[[k]];
    x = A.x + B u[[k]];
    X[[k + 1]] = x;
  ];
  <|"X" -> X, "Y" -> Y|>
];

Nsteps = 40;

u1 = ConstantArray[1.0, Nsteps];
y1 = simulateScalarDifference[0.85, 0.2, u1, 0.0];

Print["Scalar recursion (first 10 samples)"];
Do[
  Print[Row[{"k=", k - 1, ", y=", NumberForm[y1[[k]], {8, 6}]}]],
  {k, 1, 10}
];

A = { {1.5, -0.7}, {1.0, 0.0} };
B = {1.0, 0.0};
C = {0.0, 1.0};
D = 0.0;

u2 = Join[ConstantArray[1.0, 10], ConstantArray[0.0, Nsteps - 10]];
res = simulateDiscreteSS[A, B, C, D, u2, {0.0, 0.0}];

Print["State-space output (first 10 samples)"];
Do[
  Print[
    Row[{
      "k=", k - 1, 
      ", y=", NumberForm[res["Y"][[k]], {8, 6}],
      ", x=", res["X"][[k]]
    }]
  ],
  {k, 1, 10}
];

ListLinePlot[
  y1,
  PlotMarkers -> Automatic,
  GridLines -> Automatic,
  AxesLabel -> {"k (index+1 in list)", "y"},
  PlotLabel -> "Scalar Difference Equation Response"
]

ListLinePlot[
  res["Y"],
  PlotMarkers -> Automatic,
  GridLines -> Automatic,
  AxesLabel -> {"k", "y[k]"},
  PlotLabel -> "Discrete-Time State-Space Output"
]
      

The code uses explicit recurrence loops (instead of a symbolic transform approach) to reinforce the sample-to-sample dynamics introduced in this lesson.

Problem 1 (Closed-form scalar recursion): Solve the first-order difference equation \( y[k+1] = 0.8\,y[k] + 2\,u[k] \) with \( y[0]=1 \) and constant input \( u[k]=3 \). Also compute the equilibrium value.

Solution:

The equilibrium is

\[ y^\star = \frac{b}{1-a}u_0 = \frac{2}{1-0.8}(3) = 30 \]

The closed-form solution is

\[ y[k] = 0.8^k y[0] + \sum_{i=0}^{k-1}0.8^{k-1-i}(2)(3) \]

\[ y[k] = 0.8^k + 6\sum_{j=0}^{k-1}0.8^j = 0.8^k + 6\frac{1-0.8^k}{1-0.8} = 0.8^k + 30(1-0.8^k) \]

\[ y[k] = 30 - 29(0.8^k) \]

Since \( |0.8| < 1 \), the transient term decays and \( y[k] → 30 \).

Problem 2 (Difference equation to state-space): Convert \( y[k+2] + 0.4\,y[k+1] - 0.12\,y[k] = 2u[k] \) into a discrete-time state-space model using \( \mathbf{x}[k]=[y[k+1]\;\;y[k]]^\top \).

Solution:

Rearranging the difference equation gives

\[ y[k+2] = -0.4\,y[k+1] + 0.12\,y[k] + 2u[k] \]

Therefore

\[ \mathbf{x}[k+1] = \begin{bmatrix} -0.4 & 0.12 \\ 1 & 0 \end{bmatrix}\mathbf{x}[k] + \begin{bmatrix} 2 \\ 0 \end{bmatrix}u[k] \]

\[ y[k] = \begin{bmatrix} 0 & 1 \end{bmatrix}\mathbf{x}[k] \]

This is the companion-form realization associated with the given recurrence.

Problem 3 (k-step state solution): For \( \mathbf{A}=\begin{bmatrix}0.5 & 0\\0 & 0.2\end{bmatrix} \), \( \mathbf{B}=\begin{bmatrix}1\\1\end{bmatrix} \), \( \mathbf{u}[k]=1 \), and \( \mathbf{x}[0]=\begin{bmatrix}0\\0\end{bmatrix} \), compute \( \mathbf{x}[k] \) explicitly.

Solution:

Since \( \mathbf{A} \) is diagonal,

\[ \mathbf{A}^k = \begin{bmatrix} 0.5^k & 0\\ 0 & 0.2^k \end{bmatrix} \]

With zero initial state, the free term vanishes:

\[ \mathbf{x}[k] = \sum_{i=0}^{k-1}\mathbf{A}^{k-1-i}\mathbf{B} \]

Component-wise,

\[ x_1[k] = \sum_{i=0}^{k-1}0.5^{k-1-i} = \sum_{j=0}^{k-1}0.5^j = \frac{1-0.5^k}{1-0.5} = 2(1-0.5^k) \]

\[ x_2[k] = \sum_{i=0}^{k-1}0.2^{k-1-i} = \sum_{j=0}^{k-1}0.2^j = \frac{1-0.2^k}{1-0.2} = 1.25(1-0.2^k) \]

Hence

\[ \mathbf{x}[k] = \begin{bmatrix} 2(1-0.5^k)\\ 1.25(1-0.2^k) \end{bmatrix} \]

Problem 4 (Equilibrium and error dynamics): Consider the scalar recursion \( y[k+1]=1.1y[k]-0.5u[k] \) with constant \( u[k]=4 \). Find the equilibrium and determine whether the system converges to it.

Solution:

The equilibrium is

\[ y^\star = \frac{b}{1-a}u_0 = \frac{-0.5}{1-1.1}(4) = 20 \]

Define \( \tilde{y}[k]=y[k]-20 \). Then

\[ \tilde{y}[k+1]=1.1\tilde{y}[k] \]

Since \( |1.1| > 1 \), the error magnitude grows (unless exactly zero initially). Therefore the equilibrium exists algebraically but is unstable, and trajectories generally do not converge to it.

Problem 5 (Output sequence from state-space): Let \( \mathbf{A}=\begin{bmatrix}0 & 1\\-0.25 & 1.2\end{bmatrix} \), \( \mathbf{B}=\begin{bmatrix}0\\1\end{bmatrix} \), \( \mathbf{C}=\begin{bmatrix}1 & 0\end{bmatrix} \), \( \mathbf{D}=0 \), and \( u[k]=0 \). If \( \mathbf{x}[0]=\begin{bmatrix}2\\0\end{bmatrix} \), compute \( y[0],y[1],y[2] \).

Solution:

Since \( y[k]=\mathbf{C}\mathbf{x}[k] \),

\[ y[0] = \begin{bmatrix}1 & 0\end{bmatrix} \begin{bmatrix}2\\0\end{bmatrix} = 2 \]

Next state:

\[ \mathbf{x}[1]=\mathbf{A}\mathbf{x}[0] = \begin{bmatrix}0 & 1\\-0.25 & 1.2\end{bmatrix} \begin{bmatrix}2\\0\end{bmatrix} = \begin{bmatrix}0\\-0.5\end{bmatrix} \]

Thus \( y[1]=0 \). Then

\[ \mathbf{x}[2]=\mathbf{A}\mathbf{x}[1] = \begin{bmatrix}0 & 1\\-0.25 & 1.2\end{bmatrix} \begin{bmatrix}0\\-0.5\end{bmatrix} = \begin{bmatrix}-0.5\\-0.6\end{bmatrix} \]

Hence \( y[2] = -0.5 \).

In this lesson, we formalized discrete-time dynamics using scalar and higher-order difference equations, derived the k-step state-space solution, and showed how to convert a difference equation into a companion-form state-space model. We also analyzed equilibria and asymptotic stability through characteristic roots and eigenvalues of \( \mathbf{A} \). These results provide the foundation for the next lesson, where the z-transform and discrete-time transfer functions are introduced systematically.