In the previous chapters, most state-space systems were linear time-invariant: \( \dot{\mathbf{x} } = A\mathbf{x}+B\mathbf{u} \), \( \mathbf{y}=C\mathbf{x}+D\mathbf{u} \). In many engineering systems, however, the local dynamics change with time. Examples include aircraft with changing mass and speed, robotic systems linearized along a trajectory, switched operating regimes approximated by smooth parameter variation, and linearization of nonlinear systems about time-dependent nominal motions.

The continuous-time linear time-varying model is

\[ \dot{\mathbf{x} }(t)=A(t)\mathbf{x}(t)+B(t)\mathbf{u}(t),\qquad \mathbf{y}(t)=C(t)\mathbf{x}(t)+D(t)\mathbf{u}(t). \]

Here \( \mathbf{x}(t)\in\mathbb{R}^{n} \) is the state, \( \mathbf{u}(t)\in\mathbb{R}^{m} \) is the input, and \( \mathbf{y}(t)\in\mathbb{R}^{p} \) is the output. The matrices have dimensions

\[ A(t)\in\mathbb{R}^{n\times n},\quad B(t)\in\mathbb{R}^{n\times m},\quad C(t)\in\mathbb{R}^{p\times n},\quad D(t)\in\mathbb{R}^{p\times m}. \]

flowchart TD
  S["Engineering system with changing operating condition"] --> L["Linearize or model at each time"]
  L --> M["Time-varying matrices A(t), B(t), C(t), D(t)"]
  M --> DYN["State equation: xdot = A(t)x + B(t)u"]
  DYN --> OUT["Output equation: y = C(t)x + D(t)u"]
  OUT --> SIM["Analyze, simulate, and later study stability and controllability"]
        

Let the time interval of interest be \( \mathcal{I}=[t_0,t_f] \). A continuous-time LTV state-space system is specified by four matrix-valued functions \( A(\cdot),B(\cdot),C(\cdot),D(\cdot) \) and the equations

\[ \Sigma_{\mathrm{LTV} }: \begin{cases} \dot{\mathbf{x} }(t)=A(t)\mathbf{x}(t)+B(t)\mathbf{u}(t),\\ \mathbf{y}(t)=C(t)\mathbf{x}(t)+D(t)\mathbf{u}(t),\\ \mathbf{x}(t_0)=\mathbf{x}_0. \end{cases} \]

The system is linear because the map from \( (\mathbf{x}_0,\mathbf{u}) \) to \( (\mathbf{x}(t),\mathbf{y}(t)) \) satisfies superposition. It is time-varying because the matrices may depend explicitly on \( t \).

A typical regularity assumption in this introductory setting is piecewise continuity:

\[ A(t),B(t),C(t),D(t)\ \text{are piecewise continuous on}\ \mathcal{I}, \qquad \mathbf{u}(t)\ \text{is piecewise continuous on}\ \mathcal{I}. \]

Under these assumptions, the state equation has a unique absolutely continuous solution for every initial condition \( \mathbf{x}_0 \).

The interpretation of the four matrices is similar to the LTI case, but now the interpretation is local in time:

\( A(t) \): instantaneous internal coupling matrix. It describes how the current state affects the state derivative at time \( t \).

\[ A_{ij}(t)= \frac{\partial \dot{x}_i(t)}{\partial x_j(t)} \quad \text{for a linearization about a nominal trajectory.} \]

\( B(t) \): instantaneous input distribution matrix. It describes actuator authority at time \( t \).

\[ B_{ij}(t)= \frac{\partial \dot{x}_i(t)}{\partial u_j(t)}. \]

\( C(t) \): instantaneous output or sensor matrix. It maps the state to the measured output at time \( t \).

\[ C_{ij}(t)= \frac{\partial y_i(t)}{\partial x_j(t)}. \]

\( D(t) \): instantaneous feedthrough matrix. It maps the input directly to the output at the same time.

\[ D_{ij}(t)= \frac{\partial y_i(t)}{\partial u_j(t)}. \]

In mechanical systems, \( A(t) \) may change because the operating point changes; \( B(t) \) may change because actuator effectiveness varies; \( C(t) \) may change because sensors measure a moving coordinate frame; and \( D(t) \) may change because direct transmission effects depend on configuration.

First consider the homogeneous state equation:

\[ \dot{\mathbf{x} }(t)=A(t)\mathbf{x}(t),\qquad \mathbf{x}(t_0)=\mathbf{x}_0. \]

A matrix function \( X(t) \) is called a fundamental matrix if its columns are linearly independent solutions of the homogeneous system:

\[ \dot{X}(t)=A(t)X(t),\qquad \det X(t)\neq 0. \]

The state-transition matrix is the particular fundamental matrix normalized at the initial time:

\[ \Phi(t,t_0)=X(t)X^{-1}(t_0),\qquad \Phi(t_0,t_0)=I. \]

Therefore the homogeneous solution is

\[ \mathbf{x}(t)=\Phi(t,t_0)\mathbf{x}_0. \]

For LTI systems, \( \Phi(t,t_0)=e^{A(t-t_0)} \). For LTV systems, this simplification generally fails because \( A(t_1)A(t_2) \) need not equal \( A(t_2)A(t_1) \).

\[ A(t_1)A(t_2)\neq A(t_2)A(t_1) \quad \Longrightarrow \quad \Phi(t,t_0)\neq \exp\!\left(\int_{t_0}^{t}A(\tau)\,d\tau\right) \ \text{in general.} \]

For the forced LTV system

\[ \dot{\mathbf{x} }(t)=A(t)\mathbf{x}(t)+B(t)\mathbf{u}(t), \qquad \mathbf{x}(t_0)=\mathbf{x}_0, \]

the solution is

\[ \mathbf{x}(t)=\Phi(t,t_0)\mathbf{x}_0+ \int_{t_0}^{t}\Phi(t,\tau)B(\tau)\mathbf{u}(\tau)\,d\tau. \]

The output is then obtained from

\[ \mathbf{y}(t)=C(t)\Phi(t,t_0)\mathbf{x}_0+ C(t)\int_{t_0}^{t}\Phi(t,\tau)B(\tau)\mathbf{u}(\tau)\,d\tau+ D(t)\mathbf{u}(t). \]

Proof. Define \( \mathbf{x}(t)=\Phi(t,t_0)\mathbf{z}(t) \). Since \( \dot{\Phi}(t,t_0)=A(t)\Phi(t,t_0) \),

\[ \dot{\mathbf{x} }(t)=A(t)\Phi(t,t_0)\mathbf{z}(t)+ \Phi(t,t_0)\dot{\mathbf{z} }(t). \]

Matching this with \( A(t)\mathbf{x}(t)+B(t)\mathbf{u}(t) \) gives

\[ \Phi(t,t_0)\dot{\mathbf{z} }(t)=B(t)\mathbf{u}(t), \qquad \dot{\mathbf{z} }(t)=\Phi^{-1}(t,t_0)B(t)\mathbf{u}(t). \]

Integrating and using \( \Phi(t,\tau)=\Phi(t,t_0)\Phi^{-1}(\tau,t_0) \) yields the stated forced-response formula.

It is tempting to treat an LTV system as an LTI system whose matrices are slowly updated. This intuition is useful in simulation, but it is mathematically incomplete. In LTI systems, eigenvalues of \( A \) determine modal behavior. In LTV systems, instantaneous eigenvalues of \( A(t) \) do not by themselves determine stability or response.

The central object is not the spectrum of \( A(t) \) at one time instant. The central object is the mapping \( \Phi(t,t_0) \), which describes how the state at \( t_0 \) propagates to time \( t \).

\[ \mathbf{x}(t_0)\mapsto \mathbf{x}(t) \quad \text{is represented by} \quad \Phi(t,t_0). \]

This observation is the bridge to the next lessons: stability, controllability, and observability for LTV systems must be stated using time intervals and state-transition operators rather than only constant matrices.

flowchart TD
  LTI["LTI model"] --> EA["Constant A"]
  EA --> EXP["State transition: exp(A times time difference)"]
  EXP --> SPEC["Eigenvalues of A guide stability"]

  LTV["LTV model"] --> AT["A(t) changes with time"]
  AT --> STM["State transition: Phi(t,t0)"]
  STM --> INT["Behavior depends on evolution over an interval"]
        

Consider a second-order system whose natural frequency and damping vary with time:

\[ \ddot{q}(t)+2\zeta(t)\omega(t)\dot{q}(t)+\omega^2(t)q(t)=b(t)u(t). \]

Define the state \( \mathbf{x}(t)=\begin{bmatrix}q(t)&\dot{q}(t)\end{bmatrix}^{T} \). Then

\[ \dot{\mathbf{x} }(t)= \begin{bmatrix} 0 & 1\\ -\omega^2(t) & -2\zeta(t)\omega(t) \end{bmatrix}\mathbf{x}(t)+ \begin{bmatrix} 0\\ b(t) \end{bmatrix}u(t). \]

If the measured output is a time-varying sensor gain on position,

\[ y(t)= \begin{bmatrix} c(t) & 0 \end{bmatrix}\mathbf{x}(t), \]

then

\[ A(t)= \begin{bmatrix} 0 & 1\\ -\omega^2(t) & -2\zeta(t)\omega(t) \end{bmatrix}, \quad B(t)= \begin{bmatrix} 0\\ b(t) \end{bmatrix}, \quad C(t)= \begin{bmatrix} c(t) & 0 \end{bmatrix}, \quad D(t)=0. \]

To simulate an LTV state equation numerically, one evaluates \( A(t),B(t),C(t),D(t) \) at the current integration time and advances the differential equation using an ODE solver.

\[ \mathbf{x}_{k+1}\approx \mathbf{x}_{k}+ h\left[A(t_k)\mathbf{x}_k+B(t_k)\mathbf{u}(t_k)\right] \]

Higher-order solvers, such as Runge-Kutta methods, improve the local approximation. When the state-transition matrix is needed, integrate the matrix differential equation

\[ \dot{\Phi}(t,t_0)=A(t)\Phi(t,t_0), \qquad \Phi(t_0,t_0)=I. \]

This matrix equation contains \( n^2 \) scalar differential equations, so computing \( \Phi(t,t_0) \) is usually more expensive than simulating a single trajectory.

The Python implementation uses NumPy, SciPy solve_ivp, and Matplotlib. Related modern-control libraries include python-control, Slycot, SciPy, and NumPy.

Chapter29_Lesson1.py

# Chapter29_Lesson1.py
# Linear time-varying (LTV) state-space simulation:
#     x_dot = A(t)x + B(t)u(t)
#     y     = C(t)x + D(t)u(t)
#
# Requirements:
#     pip install numpy scipy matplotlib

import numpy as np
from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt


def matrices(t: float):
    """Return A(t), B(t), C(t), D(t) for a second-order LTV oscillator."""
    omega = 1.0 + 0.30 * np.sin(0.70 * t)
    damping = 0.15 + 0.05 * np.cos(0.50 * t)

    A = np.array([
        [0.0, 1.0],
        [-(omega ** 2), -2.0 * damping * omega]
    ], dtype=float)

    B = np.array([
        [0.0],
        [1.0 + 0.20 * np.sin(0.30 * t)]
    ], dtype=float)

    C = np.array([
        [1.0 + 0.10 * np.sin(0.40 * t), 0.0]
    ], dtype=float)

    D = np.array([[0.0]], dtype=float)
    return A, B, C, D


def input_signal(t: float) -> np.ndarray:
    """Scalar input written as a length-1 vector."""
    return np.array([np.sin(1.2 * t)], dtype=float)


def rhs_state(t: float, x: np.ndarray) -> np.ndarray:
    A, B, _, _ = matrices(t)
    u = input_signal(t)
    return A @ x + B @ u


def rhs_augmented_state_transition(t: float, z: np.ndarray) -> np.ndarray:
    """
    Integrate x(t) and Phi(t,t0) together.

    z = [x; vec(Phi)] where Phi is flattened column-wise.
    Phi_dot = A(t) Phi, Phi(t0,t0) = I.
    """
    n = 2
    x = z[:n]
    Phi = z[n:].reshape((n, n), order="F")
    A, B, _, _ = matrices(t)
    u = input_signal(t)

    x_dot = A @ x + B @ u
    Phi_dot = A @ Phi

    return np.concatenate([x_dot, Phi_dot.reshape(n * n, order="F")])


def main():
    t0, tf = 0.0, 20.0
    x0 = np.array([1.0, 0.0], dtype=float)
    n = len(x0)

    sol_x = solve_ivp(rhs_state, (t0, tf), x0, dense_output=True, rtol=1e-9, atol=1e-11)

    Phi0 = np.eye(n)
    z0 = np.concatenate([x0, Phi0.reshape(n * n, order="F")])
    sol_z = solve_ivp(
        rhs_augmented_state_transition,
        (t0, tf),
        z0,
        dense_output=True,
        rtol=1e-9,
        atol=1e-11,
    )

    t_grid = np.linspace(t0, tf, 600)
    x_grid = sol_z.sol(t_grid)[:n, :].T

    y_grid = []
    for t, x in zip(t_grid, x_grid):
        _, _, C, D = matrices(float(t))
        u = input_signal(float(t))
        y_grid.append(float((C @ x + D @ u)[0]))
    y_grid = np.array(y_grid)

    Phi_tf_t0 = sol_z.y[n:, -1].reshape((n, n), order="F")
    print("Phi(tf,t0) =")
    print(Phi_tf_t0)
    print("x(tf) =", sol_z.y[:n, -1])
    print("y(tf) =", y_grid[-1])

    plt.figure()
    plt.plot(t_grid, x_grid[:, 0], label="x1(t)")
    plt.plot(t_grid, x_grid[:, 1], label="x2(t)")
    plt.plot(t_grid, y_grid, label="y(t)", linestyle="--")
    plt.xlabel("time")
    plt.ylabel("response")
    plt.title("Chapter 29 Lesson 1: LTV state-space response")
    plt.grid(True)
    plt.legend()
    plt.show()


if __name__ == "__main__":
    main()
      

The C++ implementation uses fixed-step RK4 from scratch. For larger modern-control projects, useful libraries include Eigen, Armadillo, Boost.Odeint, and CasADi.

Chapter29_Lesson1.cpp

// Chapter29_Lesson1.cpp
// Linear time-varying (LTV) state-space simulation using fixed-step RK4.
// Compile:
//     g++ -std=c++17 -O2 Chapter29_Lesson1.cpp -o Chapter29_Lesson1

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

using Vec2 = std::array<double, 2>;
using Mat2 = std::array<std::array<double, 2>, 2>;

struct Matrices {
    Mat2 A;
    Vec2 B;
    std::array<double, 2> C;
    double D;
};

Matrices matrices(double t) {
    const double omega = 1.0 + 0.30 * std::sin(0.70 * t);
    const double damping = 0.15 + 0.05 * std::cos(0.50 * t);

    Matrices m{};
    m.A = { {
        { {0.0, 1.0} },
        { {-(omega * omega), -2.0 * damping * omega} }
    } };
    m.B = { {0.0, 1.0 + 0.20 * std::sin(0.30 * t)} };
    m.C = { {1.0 + 0.10 * std::sin(0.40 * t), 0.0} };
    m.D = 0.0;
    return m;
}

double inputSignal(double t) {
    return std::sin(1.2 * t);
}

Vec2 rhs(double t, const Vec2& x) {
    Matrices m = matrices(t);
    const double u = inputSignal(t);

    return { {
        m.A[0][0] * x[0] + m.A[0][1] * x[1] + m.B[0] * u,
        m.A[1][0] * x[0] + m.A[1][1] * x[1] + m.B[1] * u
    } };
}

Vec2 addScaled(const Vec2& x, const Vec2& k, double scale) {
    return { {x[0] + scale * k[0], x[1] + scale * k[1]} };
}

Vec2 rk4Step(double t, const Vec2& x, double h) {
    Vec2 k1 = rhs(t, x);
    Vec2 k2 = rhs(t + 0.5 * h, addScaled(x, k1, 0.5 * h));
    Vec2 k3 = rhs(t + 0.5 * h, addScaled(x, k2, 0.5 * h));
    Vec2 k4 = rhs(t + h, addScaled(x, k3, h));

    return { {
        x[0] + (h / 6.0) * (k1[0] + 2.0 * k2[0] + 2.0 * k3[0] + k4[0]),
        x[1] + (h / 6.0) * (k1[1] + 2.0 * k2[1] + 2.0 * k3[1] + k4[1])
    } };
}

int main() {
    const double t0 = 0.0;
    const double tf = 20.0;
    const double h = 0.01;

    Vec2 x = { {1.0, 0.0} };

    std::cout << std::fixed << std::setprecision(8);
    std::cout << "t,x1,x2,y\n";

    for (double t = t0; t <= tf + 1e-12; t += h) {
        Matrices m = matrices(t);
        double u = inputSignal(t);
        double y = m.C[0] * x[0] + m.C[1] * x[1] + m.D * u;

        if (static_cast<int>(std::round(t / h)) % 100 == 0) {
            std::cout << t << "," << x[0] << "," << x[1] << "," << y << "\n";
        }

        if (t + h <= tf + 1e-12) {
            x = rk4Step(t, x, h);
        }
    }

    return 0;
}
      

The Java implementation also uses RK4 from scratch. For larger projects, useful Java numerical libraries include Apache Commons Math, EJML, ojAlgo, and Hipparchus.

Chapter29_Lesson1.java

// Chapter29_Lesson1.java
// Linear time-varying (LTV) state-space simulation using fixed-step RK4.
// Compile and run:
//     javac Chapter29_Lesson1.java
//     java Chapter29_Lesson1

public class Chapter29_Lesson1 {
    static class Matrices {
        double[][] A = new double[2][2];
        double[] B = new double[2];
        double[] C = new double[2];
        double D;
    }

    static Matrices matrices(double t) {
        double omega = 1.0 + 0.30 * Math.sin(0.70 * t);
        double damping = 0.15 + 0.05 * Math.cos(0.50 * t);

        Matrices m = new Matrices();
        m.A[0][0] = 0.0;
        m.A[0][1] = 1.0;
        m.A[1][0] = -(omega * omega);
        m.A[1][1] = -2.0 * damping * omega;
        m.B[0] = 0.0;
        m.B[1] = 1.0 + 0.20 * Math.sin(0.30 * t);
        m.C[0] = 1.0 + 0.10 * Math.sin(0.40 * t);
        m.C[1] = 0.0;
        m.D = 0.0;
        return m;
    }

    static double inputSignal(double t) {
        return Math.sin(1.2 * t);
    }

    static double[] rhs(double t, double[] x) {
        Matrices m = matrices(t);
        double u = inputSignal(t);
        return new double[] {
            m.A[0][0] * x[0] + m.A[0][1] * x[1] + m.B[0] * u,
            m.A[1][0] * x[0] + m.A[1][1] * x[1] + m.B[1] * u
        };
    }

    static double[] addScaled(double[] x, double[] k, double scale) {
        return new double[] {x[0] + scale * k[0], x[1] + scale * k[1]};
    }

    static double[] rk4Step(double t, double[] x, double h) {
        double[] k1 = rhs(t, x);
        double[] k2 = rhs(t + 0.5 * h, addScaled(x, k1, 0.5 * h));
        double[] k3 = rhs(t + 0.5 * h, addScaled(x, k2, 0.5 * h));
        double[] k4 = rhs(t + h, addScaled(x, k3, h));

        return new double[] {
            x[0] + (h / 6.0) * (k1[0] + 2.0 * k2[0] + 2.0 * k3[0] + k4[0]),
            x[1] + (h / 6.0) * (k1[1] + 2.0 * k2[1] + 2.0 * k3[1] + k4[1])
        };
    }

    public static void main(String[] args) {
        double t0 = 0.0;
        double tf = 20.0;
        double h = 0.01;
        double[] x = {1.0, 0.0};

        System.out.println("t,x1,x2,y");
        for (double t = t0; t <= tf + 1e-12; t += h) {
            Matrices m = matrices(t);
            double u = inputSignal(t);
            double y = m.C[0] * x[0] + m.C[1] * x[1] + m.D * u;

            int step = (int) Math.round(t / h);
            if (step % 100 == 0) {
                System.out.printf("%.8f,%.8f,%.8f,%.8f%n", t, x[0], x[1], y);
            }

            if (t + h <= tf + 1e-12) {
                x = rk4Step(t, x, h);
            }
        }
    }
}
      

MATLAB directly supports numerical integration using ode45. The Control System Toolbox is primarily oriented toward LTI and some LTV workflows through gain scheduling and simulation scripting. In Simulink, the standard approach is to compute \( A(t),B(t),C(t),D(t) \) inside MATLAB Function blocks and connect them to an Integrator block.

Chapter29_Lesson1.m

% Chapter29_Lesson1.m
% Linear time-varying (LTV) state-space simulation:
%     x_dot = A(t)x + B(t)u(t)
%     y     = C(t)x + D(t)u(t)

clear; clc; close all;

tspan = [0, 20];
x0 = [1; 0];

opts = odeset('RelTol', 1e-9, 'AbsTol', 1e-11);
[t, x] = ode45(@rhs_state, tspan, x0, opts);

y = zeros(size(t));
for k = 1:length(t)
    [~, ~, C, D] = ltv_matrices(t(k));
    u = input_signal(t(k));
    y(k) = C * x(k, :)' + D * u;
end

figure;
plot(t, x(:, 1), 'LineWidth', 1.5); hold on;
plot(t, x(:, 2), 'LineWidth', 1.5);
plot(t, y, '--', 'LineWidth', 1.5);
grid on;
xlabel('time');
ylabel('response');
title('Chapter 29 Lesson 1: LTV state-space response');
legend('x_1(t)', 'x_2(t)', 'y(t)');

% Optional Simulink guidance:
% Implement A(t), B(t), C(t), D(t) as MATLAB Function blocks, use an
% Integrator block for x_dot, and feed x and u into y = C(t)x + D(t)u.

function dx = rhs_state(t, x)
    [A, B, ~, ~] = ltv_matrices(t);
    u = input_signal(t);
    dx = A * x + B * u;
end

function [A, B, C, D] = ltv_matrices(t)
    omega = 1.0 + 0.30 * sin(0.70 * t);
    damping = 0.15 + 0.05 * cos(0.50 * t);

    A = [0, 1;
        -(omega^2), -2.0 * damping * omega];

    B = [0;
        1.0 + 0.20 * sin(0.30 * t)];

    C = [1.0 + 0.10 * sin(0.40 * t), 0];
    D = 0;
end

function u = input_signal(t)
    u = sin(1.2 * t);
end
      

Mathematica can simulate the LTV differential equations using NDSolve. Although StateSpaceModel is most convenient for LTI models, direct differential-equation representation is natural for LTV models.

Chapter29_Lesson1.nb

(* Chapter29_Lesson1.nb *)
(* Wolfram Mathematica code for a linear time-varying state-space model. *)

ClearAll["Global`*"];

omega[t_] := 1.0 + 0.30 Sin[0.70 t];
damping[t_] := 0.15 + 0.05 Cos[0.50 t];

A[t_] := { {0, 1}, {-omega[t]^2, -2 damping[t] omega[t]} };
B[t_] := { {0}, {1.0 + 0.20 Sin[0.30 t]} };
Cmat[t_] := { {1.0 + 0.10 Sin[0.40 t], 0} };
Dmat[t_] := { {0} };
u[t_] := Sin[1.2 t];

t0 = 0;
tf = 20;
x10 = 1;
x20 = 0;

sol = NDSolve[
   {
    x1'[t] == x2[t],
    x2'[t] == -omega[t]^2 x1[t] - 2 damping[t] omega[t] x2[t] +
      (1.0 + 0.20 Sin[0.30 t]) u[t],
    x1[t0] == x10,
    x2[t0] == x20
    },
   {x1, x2},
   {t, t0, tf},
   WorkingPrecision -> MachinePrecision
   ];

y[t_] := (Cmat[t].{ {x1[t]}, {x2[t]} } + Dmat[t].{ {u[t]} })[[1, 1]] /. sol[[1]];

Plot[
 Evaluate[{x1[t], x2[t], y[t]} /. sol[[1]]],
 {t, t0, tf},
 PlotLegends -> {"x1(t)", "x2(t)", "y(t)"},
 AxesLabel -> {"time", "response"},
 PlotLabel -> "Chapter 29 Lesson 1: LTV state-space response"
 ]

(* State-transition matrix Phi(t,t0): solve Phi_dot = A(t) Phi, Phi(t0,t0) = I. *)
phiSol = NDSolve[
   Flatten[
    Join[
     Table[Derivative[1, 0][phi][i, j][t] ==
       Sum[A[t][[i, k]] phi[k, j][t], {k, 1, 2}], {i, 1, 2}, {j, 1, 2}],
     Table[phi[i, j][t0] == KroneckerDelta[i, j], {i, 1, 2}, {j, 1, 2}]
     ]
    ],
   Flatten[Table[phi[i, j], {i, 1, 2}, {j, 1, 2}]],
   {t, t0, tf}
   ];

Phi[t_] := Table[phi[i, j][t], {i, 1, 2}, {j, 1, 2}] /. phiSol[[1]];
MatrixForm[Phi[tf]]
      

Problem 1: Consider \( \dot{\mathbf{x} }(t)=A(t)\mathbf{x}(t) \) with \( \Phi(t,t_0) \) satisfying \( \dot{\Phi}(t,t_0)=A(t)\Phi(t,t_0) \) and \( \Phi(t_0,t_0)=I \). Show that \( \mathbf{x}(t)=\Phi(t,t_0)\mathbf{x}_0 \) solves the homogeneous initial-value problem.

Solution: Differentiate the candidate solution:

\[ \dot{\mathbf{x} }(t)=\dot{\Phi}(t,t_0)\mathbf{x}_0 =A(t)\Phi(t,t_0)\mathbf{x}_0=A(t)\mathbf{x}(t). \]

At \( t=t_0 \),

\[ \mathbf{x}(t_0)=\Phi(t_0,t_0)\mathbf{x}_0=I\mathbf{x}_0=\mathbf{x}_0. \]

Therefore the formula satisfies both the differential equation and the initial condition.

Problem 2: For the scalar LTV system \( \dot{x}(t)=a(t)x(t)+b(t)u(t) \), derive the state-transition function and the forced response.

Solution: Since the system is scalar, multiplication commutes. The transition function is

\[ \phi(t,t_0)= \exp\!\left(\int_{t_0}^{t}a(\sigma)\,d\sigma\right). \]

Hence

\[ x(t)=\phi(t,t_0)x_0+ \int_{t_0}^{t}\phi(t,\tau)b(\tau)u(\tau)\,d\tau. \]

Problem 3: Let \( A(t)=\begin{bmatrix}0 & 1\\ -k(t) & -c(t)\end{bmatrix} \). Write the corresponding second-order differential equation when \( B(t)=\begin{bmatrix}0\\ 1\end{bmatrix} \) and \( \mathbf{x}=\begin{bmatrix}q & \dot{q}\end{bmatrix}^{T} \).

Solution: The state equations are

\[ \dot{x}_1=x_2,\qquad \dot{x}_2=-k(t)x_1-c(t)x_2+u(t). \]

Since \( x_1=q \) and \( x_2=\dot{q} \), we obtain

\[ \ddot{q}(t)+c(t)\dot{q}(t)+k(t)q(t)=u(t). \]

Problem 4: Explain why the formula \( \Phi(t,t_0)=\exp\!\left(\int_{t_0}^{t}A(\tau)d\tau\right) \) is not generally valid for matrix-valued LTV systems.

Solution: If \( A(t_1) \) and \( A(t_2) \) do not commute, then the order in which the infinitesimal dynamics occur matters. A single exponential of an integrated matrix loses this order information. Therefore the formula is valid only in special cases, for example when \( A(t_1)A(t_2)=A(t_2)A(t_1) \) for all relevant times.

Problem 5: Suppose \( C(t)=\begin{bmatrix}1+\epsilon\sin t & 0\end{bmatrix} \) and \( D(t)=0 \). If \( \mathbf{x}(t)=\begin{bmatrix}x_1(t)&x_2(t)\end{bmatrix}^{T} \), find \( y(t) \).

Solution:

\[ y(t)=C(t)\mathbf{x}(t)= \begin{bmatrix} 1+\epsilon\sin t & 0 \end{bmatrix} \begin{bmatrix} x_1(t)\\ x_2(t) \end{bmatrix} = \left(1+\epsilon\sin t\right)x_1(t). \]

Linear time-varying systems preserve linearity but lose time invariance. Their state-space form is \( \dot{\mathbf{x} }(t)=A(t)\mathbf{x}(t)+B(t)\mathbf{u}(t) \) and \( \mathbf{y}(t)=C(t)\mathbf{x}(t)+D(t)\mathbf{u}(t) \). The main mathematical object is the state-transition matrix \( \Phi(t,t_0) \), not a constant matrix exponential. The forced response is obtained through a time-varying convolution-like integral involving \( \Phi(t,\tau)B(\tau)\mathbf{u}(\tau) \). These ideas prepare the ground for the next lesson, where the time-varying state-transition matrix is studied in more detail.