A linear time-varying (LTV) system becomes slowly time-varying when its matrices change on a time scale much larger than the dominant state response time. A convenient representation is

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

The variable \( s=\varepsilon t \) is called the slow-time variable. Engineering examples include vehicles whose mass changes because of fuel consumption, robotic manipulators whose payload changes slowly, thermal systems whose coefficients drift with temperature, aircraft whose operating point changes along a trajectory, and electrical circuits whose equivalent resistance changes with heating or aging.

In this lesson, "slow" does not mean static. It means that the matrix derivative is small relative to the local modal decay or oscillation rates. For a homogeneous system, a useful dimensionless diagnostic is

\[ \eta(t)= \frac{\|\dot{\mathbf{A} }(t)\|}{\alpha(t)^2},\qquad \alpha(t)= -\max_i \operatorname{Re}\lambda_i(\mathbf{A}(t)). \]

If \( \alpha(t) \) remains positive and \( \eta(t) \) remains small, the dynamics often behave close to a sequence of frozen-time LTI systems. However, as learned in the previous lessons, frozen eigenvalues alone are not a complete stability theorem for arbitrary LTV systems.

flowchart TD
  A["Engineering plant"] --> B["Identify states, inputs, outputs"]
  B --> C["Detect drifting parameters"]
  C --> D["Write A(s), B(s), C(s), D(s) with slow time s"]
  D --> E["Check frozen-time stability margins"]
  E --> F["Check variation rate relative to modal speed"]
  F --> G["Simulate full LTV model"]
  G --> H["Compare with frozen-time approximation"]
        

At a fixed instant \( t_0 \), the frozen-time model is the LTI approximation obtained by holding the matrices constant:

\[ \dot{\mathbf{z} }(t)=\mathbf{A}(t_0)\mathbf{z}(t)+ \mathbf{B}(t_0)\mathbf{u}(t),\qquad t\geq t_0. \]

This model is local in time. It is useful for interpreting modal structure, local damping, local natural frequency, and instantaneous input authority. It is not a replacement for the full state transition matrix \( \boldsymbol{\Phi}(t,t_0) \).

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

The frozen approximation over a short interval \( [t_k,t_{k+1}] \) is

\[ \mathbf{x}(t_{k+1})\approx e^{\mathbf{A}(t_k)h}\mathbf{x}(t_k)+ \int_{0}^{h}e^{\mathbf{A}(t_k)(h-\rho)} \mathbf{B}(t_k)\mathbf{u}(t_k+\rho)\,d\rho,\qquad h=t_{k+1}-t_k. \]

This is the computational basis of frozen-step simulation and gain-scheduled design. Its accuracy improves when \( h \|\dot{\mathbf{A} }\| \) and \( h \|\dot{\mathbf{B} }\| \) are small.

Let \( \mathbf{A}(\cdot) \in C^1([0,\infty),\mathbb{R}^{n\times n}) \). One practical definition of slow variation is the existence of constants \( M_A \) and \( \varepsilon \) such that

\[ \|\dot{\mathbf{A} }(t)\|\leq \varepsilon M_A,\qquad \|\dot{\mathbf{B} }(t)\|\leq \varepsilon M_B,\qquad 0 < \varepsilon \ll 1. \]

Suppose all frozen-time systems have a uniform decay margin:

\[ \max_i\operatorname{Re}\lambda_i(\mathbf{A}(t))\leq -\alpha \qquad\text{for all }t\geq 0,\qquad \alpha > 0. \]

This condition alone does not generally prove uniform exponential stability for an arbitrary LTV system. A stronger and more useful sufficient condition is the existence of a smoothly varying positive definite Lyapunov matrix \( \mathbf{P}(t) \):

\[ m_1\mathbf{I}\preceq\mathbf{P}(t)\preceq m_2\mathbf{I},\qquad \mathbf{A}(t)^{T}\mathbf{P}(t)+\mathbf{P}(t)\mathbf{A}(t)+ \dot{\mathbf{P} }(t)\preceq -q\mathbf{I},\qquad q > 0. \]

Then the homogeneous LTV system is uniformly exponentially stable. The proof follows directly from a time-varying quadratic Lyapunov function.

Proof. Define

\[ V(t,\mathbf{x})=\mathbf{x}^{T}\mathbf{P}(t)\mathbf{x}. \]

Along the trajectories of \( \dot{\mathbf{x} }=\mathbf{A}(t)\mathbf{x} \),

\[ \dot{V}=\mathbf{x}^{T}\left( \mathbf{A}^{T}\mathbf{P}+\mathbf{P}\mathbf{A}+\dot{\mathbf{P} } \right)\mathbf{x}\leq -q\|\mathbf{x}\|^2. \]

Since \( m_1\|\mathbf{x}\|^2\leq V\leq m_2\|\mathbf{x}\|^2 \), we obtain

\[ \dot{V}\leq -\frac{q}{m_2}V,\qquad V(t)\leq V(t_0)e^{-\frac{q}{m_2}(t-t_0)}. \]

Therefore

\[ \|\mathbf{x}(t)\|\leq \sqrt{\frac{m_2}{m_1} }\, e^{-\frac{q}{2m_2}(t-t_0)}\|\mathbf{x}(t_0)\|,\qquad t\geq t_0. \]

Consider a mass-spring-damper oscillator in which the mass, damping, and stiffness vary slowly because of fuel consumption, payload redistribution, thermal drift, or material compliance changes:

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

With states \( x_1=q \) and \( x_2=\dot{q} \),

\[ \dot{\mathbf{x} }= \begin{bmatrix} 0 & 1\\ -\frac{k(t)}{m(t)} & -\frac{c(t)}{m(t)} \end{bmatrix}\mathbf{x}+ \begin{bmatrix} 0\\ \frac{1}{m(t)} \end{bmatrix}u. \]

A typical slowly varying parameterization is

\[ m(t)=m_0(1+\delta_m\sin(\varepsilon t)),\qquad c(t)=c_0(1+\delta_c\cos(\varepsilon t)),\qquad k(t)=k_0(1+\delta_k\sin(\varepsilon t/2)). \]

The frozen-time natural frequency and damping ratio are

\[ \omega_n(t)=\sqrt{\frac{k(t)}{m(t)} },\qquad \zeta(t)=\frac{c(t)}{2\sqrt{k(t)m(t)} }. \]

If \( k(t)>0 \) and \( c(t)>0 \) for all time, each frozen mechanical model is damped. Slow variation allows the actual LTV response to track the changing modal envelope, but rapid parameter changes can inject energy and invalidate frozen-time intuition.

Many vehicle models are nonlinear in the original variables. Around a slowly changing operating trajectory \( \mathbf{x}_0(t),\mathbf{u}_0(t) \), the perturbation model is LTV:

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

\[ \mathbf{A}(t)= \left.\frac{\partial \mathbf{f} }{\partial \mathbf{x} }\right|_{ \mathbf{x}_0(t),\mathbf{u}_0(t)},\qquad \mathbf{B}(t)= \left.\frac{\partial \mathbf{f} }{\partial \mathbf{u} }\right|_{ \mathbf{x}_0(t),\mathbf{u}_0(t)}. \]

For aircraft, \( \mathbf{A}(t) \) may vary with speed, altitude, trim angle, and mass. For marine vehicles, hydrodynamic derivatives may vary with speed, depth, and propulsor operating point. For mobile robots, wheel-ground contact parameters and payload inertia may vary during a mission.

A standard engineering workflow is gain scheduling: design \( \mathbf{K}(s) \) for frozen operating points and implement \( \mathbf{u}(t)=-\mathbf{K}(\varepsilon t)\mathbf{x}(t) \). The closed-loop LTV matrix is

\[ \mathbf{A}_{cl}(t)=\mathbf{A}(t)- \mathbf{B}(t)\mathbf{K}(t). \]

Stability must be verified for the scheduled closed-loop system, not only for the frozen designs.

A robot manipulator linearized around a slowly changing trajectory produces an LTV model. Near a nominal motion \( \mathbf{q}_0(t) \), the perturbation coordinates often satisfy

\[ \delta\dot{\mathbf{x} }= \mathbf{A}(t)\delta\mathbf{x}+ \mathbf{B}(t)\delta\boldsymbol{\tau},\qquad \delta\mathbf{x}= \begin{bmatrix}\delta\mathbf{q}\\ \delta\dot{\mathbf{q} }\end{bmatrix}. \]

If payload mass changes slowly, the inertia matrix \( \mathbf{M}(\mathbf{q},t) \) changes slowly, and so do the linearized matrices. A simplified second-order perturbation model is

\[ \mathbf{M}(t)\delta\ddot{\mathbf{q} }+ \mathbf{D}(t)\delta\dot{\mathbf{q} }+ \mathbf{K}(t)\delta\mathbf{q}=\delta\boldsymbol{\tau}. \]

With state \( \delta\mathbf{x}=[\delta\mathbf{q}^{T}\; \delta\dot{\mathbf{q} }^{T}]^{T} \), the first-order form is

\[ \delta\dot{\mathbf{x} }= \begin{bmatrix} \mathbf{0} & \mathbf{I}\\ -\mathbf{M}(t)^{-1}\mathbf{K}(t) & -\mathbf{M}(t)^{-1}\mathbf{D}(t) \end{bmatrix}\delta\mathbf{x}+ \begin{bmatrix} \mathbf{0}\\ \mathbf{M}(t)^{-1} \end{bmatrix}\delta\boldsymbol{\tau}. \]

This example is important because the frozen modal picture has physical meaning: natural modes and damping ratios move as payload and configuration move.

Thermal and electrical systems often have coefficients that drift slowly because of heating, aging, charge state, or environmental conditions. A simple RC electrical model with temperature-dependent resistance is

\[ C\dot{v}(t)=-\frac{1}{R(t)}v(t)+\frac{1}{R(t)}u(t). \]

In state-space form,

\[ \dot{x}(t)=a(t)x(t)+b(t)u(t),\qquad a(t)=-\frac{1}{R(t)C},\qquad b(t)=\frac{1}{R(t)C}. \]

If \( R(t)=R_0(1+\gamma\theta(t)) \) and temperature \( \theta(t) \) changes slowly, then the electrical state is fast while the coefficient drift is slow. Similar structures appear in heat exchangers and battery equivalent-circuit models.

Consider the homogeneous system \( \dot{\mathbf{x} }=\mathbf{A}(t)\mathbf{x} \). On a short interval \( [t_k,t_k+h] \), write

\[ \mathbf{A}(t_k+\rho)=\mathbf{A}(t_k)+ \int_{0}^{\rho}\dot{\mathbf{A} }(t_k+s)\,ds. \]

If \( \|\dot{\mathbf{A} }(t)\|\leq \varepsilon M_A \), then

\[ \|\mathbf{A}(t_k+\rho)-\mathbf{A}(t_k)\| \leq \varepsilon M_A \rho. \]

Let \( \mathbf{x}_f \) be the frozen solution using \( \mathbf{A}(t_k) \). Variation of constants gives an error equation of the form

\[ \mathbf{e}(t_k+h)= \int_{0}^{h} e^{\mathbf{A}(t_k)(h-\rho)} \left(\mathbf{A}(t_k+\rho)-\mathbf{A}(t_k)\right) \mathbf{x}(t_k+\rho)\,d\rho. \]

If \( \|e^{\mathbf{A}(t_k)r}\|\leq M e^{-\alpha r} \) locally and \( \|\mathbf{x}(t)\|\leq X \) on the interval, then

\[ \|\mathbf{e}(t_k+h)\| \leq M\varepsilon M_A X \int_{0}^{h}e^{-\alpha(h-\rho)}\rho\,d\rho =O(\varepsilon h^2). \]

This explains why frozen-time simulation is accurate for small step lengths and slowly varying matrices.

flowchart TD
  A["Select parameter trajectory"] --> B["Build A(t), B(t)"]
  B --> C["Compute frozen eigenvalues"]
  C --> D["Estimate dA/dt and dB/dt"]
  D --> E["Compute slow-variation index"]
  E --> F["Simulate full LTV system"]
  F --> G["Simulate frozen-step model"]
  G --> H["Compare states, margins, and errors"]
        

Chapter29_Lesson5.py

This program simulates the slowly varying mass-spring-damper system using scipy.integrate.solve_ivp, compares it with a frozen-step matrix-exponential propagation, and computes a numerical slow-variation index.

"""
Chapter29_Lesson5.py

Engineering example for slowly time-varying linear systems.
System: mass-spring-damper oscillator with slowly drifting mass, damping, and stiffness.

State:
    x1 = displacement
    x2 = velocity

Dynamics:
    x_dot = A(t) x + B(t) u(t)

Required libraries:
    numpy
    scipy
    matplotlib
"""

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


EPSILON = 0.03


def mass(t):
    return 1.0 + 0.12 * np.sin(EPSILON * t)


def damping(t):
    return 0.22 + 0.04 * np.cos(EPSILON * t)


def stiffness(t):
    return 2.0 + 0.30 * np.sin(0.5 * EPSILON * t)


def input_force(t):
    return 0.2 * np.sin(0.7 * t)


def A_matrix(t):
    m = mass(t)
    c = damping(t)
    k = stiffness(t)
    return np.array([[0.0, 1.0], [-k / m, -c / m]], dtype=float)


def B_matrix(t):
    return np.array([0.0, 1.0 / mass(t)], dtype=float)


def rhs(t, x):
    return A_matrix(t).dot(x) + B_matrix(t) * input_force(t)


def frozen_step(A, b, u_value, x, dt):
    """
    Exact step for frozen dynamics over [t, t + dt]:
        x_dot = A x + b u_value
    The input is frozen during the step.
    """
    n = A.shape[0]
    aug = np.zeros((n + 1, n + 1))
    aug[:n, :n] = A
    aug[:n, n] = b * u_value
    phi_aug = expm(aug * dt)
    return phi_aug[:n, :n].dot(x) + phi_aug[:n, n]


def simulate_frozen_grid(t_grid, x0):
    x = np.zeros((len(t_grid), 2))
    x[0] = x0
    for i in range(len(t_grid) - 1):
        t = t_grid[i]
        dt = t_grid[i + 1] - t_grid[i]
        x[i + 1] = frozen_step(A_matrix(t), B_matrix(t), input_force(t), x[i], dt)
    return x


def slow_variation_index(t_grid):
    """
    Compute a numerical index:
        eta(t) = ||dA/dt||_F / alpha(t)^2
    where alpha(t) is the frozen-time decay margin:
        alpha(t) = -max Re(lambda_i(A(t)))
    A small eta suggests a slowly varying regime relative to the local decay rate.
    """
    values = []
    for t in t_grid:
        h = 1e-3
        dA = (A_matrix(t + h) - A_matrix(t - h)) / (2.0 * h)
        eigvals = np.linalg.eigvals(A_matrix(t))
        alpha = -np.max(np.real(eigvals))
        eta = np.linalg.norm(dA, ord="fro") / max(alpha * alpha, 1e-12)
        values.append((t, alpha, eta))
    return np.array(values)


def main():
    t0, tf = 0.0, 160.0
    x0 = np.array([1.0, 0.0])
    t_eval = np.linspace(t0, tf, 2001)

    sol = solve_ivp(rhs, (t0, tf), x0, t_eval=t_eval, rtol=1e-9, atol=1e-11)
    frozen = simulate_frozen_grid(t_eval, x0)
    svi = slow_variation_index(t_eval[::100])

    print("Final continuous LTV state:", sol.y[:, -1])
    print("Final frozen-step state:", frozen[-1])
    print("Maximum absolute frozen-step error:", np.max(np.abs(sol.y.T - frozen)))
    print("Maximum slow-variation index eta:", np.max(svi[:, 2]))

    np.savetxt(
        "Chapter29_Lesson5_python_output.csv",
        np.column_stack((t_eval, sol.y.T, frozen)),
        delimiter=",",
        header="t,x1_ltv,x2_ltv,x1_frozen,x2_frozen",
        comments="",
    )

    plt.figure()
    plt.plot(t_eval, sol.y[0], label="LTV displacement")
    plt.plot(t_eval, frozen[:, 0], "--", label="frozen-step displacement")
    plt.xlabel("time")
    plt.ylabel("displacement")
    plt.title("Slowly Time-Varying Mass-Spring-Damper System")
    plt.legend()
    plt.grid(True)
    plt.tight_layout()
    plt.savefig("Chapter29_Lesson5_python_plot.png", dpi=180)
    plt.show()


if __name__ == "__main__":
    main()

Chapter29_Lesson5.cpp

The C++ version uses a fourth-order Runge-Kutta integrator from scratch and writes the simulated state and frozen-time decay margin to a CSV file.

/*
Chapter29_Lesson5.cpp

Slowly time-varying mass-spring-damper example.
No external C++ libraries are required.

Compile:
    g++ -std=c++17 -O2 Chapter29_Lesson5.cpp -o Chapter29_Lesson5

Run:
    ./Chapter29_Lesson5
*/

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

struct State {
    double x1;
    double x2;
};

const double EPSILON = 0.03;

double mass(double t) {
    return 1.0 + 0.12 * std::sin(EPSILON * t);
}

double damping(double t) {
    return 0.22 + 0.04 * std::cos(EPSILON * t);
}

double stiffness(double t) {
    return 2.0 + 0.30 * std::sin(0.5 * EPSILON * t);
}

double inputForce(double t) {
    return 0.2 * std::sin(0.7 * t);
}

State derivative(double t, State x) {
    double m = mass(t);
    double c = damping(t);
    double k = stiffness(t);
    double u = inputForce(t);

    State dx;
    dx.x1 = x.x2;
    dx.x2 = -(k / m) * x.x1 - (c / m) * x.x2 + (1.0 / m) * u;
    return dx;
}

State add(State a, State b, double scale) {
    return State{a.x1 + scale * b.x1, a.x2 + scale * b.x2};
}

State rk4Step(double t, State x, double dt) {
    State k1 = derivative(t, x);
    State k2 = derivative(t + 0.5 * dt, add(x, k1, 0.5 * dt));
    State k3 = derivative(t + 0.5 * dt, add(x, k2, 0.5 * dt));
    State k4 = derivative(t + dt, add(x, k3, dt));

    State next;
    next.x1 = x.x1 + (dt / 6.0) * (k1.x1 + 2.0 * k2.x1 + 2.0 * k3.x1 + k4.x1);
    next.x2 = x.x2 + (dt / 6.0) * (k1.x2 + 2.0 * k2.x2 + 2.0 * k3.x2 + k4.x2);
    return next;
}

double frozenDecayMargin(double t) {
    double m = mass(t);
    double c = damping(t);
    double k = stiffness(t);

    double trace = -c / m;
    double determinant = k / m;
    double discriminant = trace * trace - 4.0 * determinant;

    if (discriminant >= 0.0) {
        double lambda1 = 0.5 * (trace + std::sqrt(discriminant));
        double lambda2 = 0.5 * (trace - std::sqrt(discriminant));
        return -std::max(lambda1, lambda2);
    }

    return -0.5 * trace;
}

int main() {
    const double t0 = 0.0;
    const double tf = 160.0;
    const double dt = 0.01;
    const int steps = static_cast<int>((tf - t0) / dt);

    State x{1.0, 0.0};
    std::ofstream file("Chapter29_Lesson5_cpp_output.csv");
    file << "t,x1,x2,m,c,k,alpha\n";
    file << std::fixed << std::setprecision(10);

    for (int i = 0; i <= steps; ++i) {
        double t = t0 + i * dt;
        file << t << "," << x.x1 << "," << x.x2 << ","
             << mass(t) << "," << damping(t) << "," << stiffness(t) << ","
             << frozenDecayMargin(t) << "\n";

        if (i < steps) {
            x = rk4Step(t, x, dt);
        }
    }

    file.close();

    std::cout << "Final state: x1 = " << x.x1 << ", x2 = " << x.x2 << "\n";
    std::cout << "CSV written to Chapter29_Lesson5_cpp_output.csv\n";
    return 0;
}

Chapter29_Lesson5.java

The Java implementation mirrors the C++ implementation and is useful for students who want a platform-independent object-oriented simulation.

/*
Chapter29_Lesson5.java

Slowly time-varying mass-spring-damper example.
No external Java libraries are required.

Compile:
    javac Chapter29_Lesson5.java

Run:
    java Chapter29_Lesson5
*/

import java.io.FileWriter;
import java.io.IOException;
import java.io.PrintWriter;
import java.util.Locale;

public class Chapter29_Lesson5 {
    private static final double EPSILON = 0.03;

    static class State {
        double x1;
        double x2;

        State(double x1, double x2) {
            this.x1 = x1;
            this.x2 = x2;
        }
    }

    static double mass(double t) {
        return 1.0 + 0.12 * Math.sin(EPSILON * t);
    }

    static double damping(double t) {
        return 0.22 + 0.04 * Math.cos(EPSILON * t);
    }

    static double stiffness(double t) {
        return 2.0 + 0.30 * Math.sin(0.5 * EPSILON * t);
    }

    static double inputForce(double t) {
        return 0.2 * Math.sin(0.7 * t);
    }

    static State derivative(double t, State x) {
        double m = mass(t);
        double c = damping(t);
        double k = stiffness(t);
        double u = inputForce(t);

        double dx1 = x.x2;
        double dx2 = -(k / m) * x.x1 - (c / m) * x.x2 + (1.0 / m) * u;
        return new State(dx1, dx2);
    }

    static State add(State a, State b, double scale) {
        return new State(a.x1 + scale * b.x1, a.x2 + scale * b.x2);
    }

    static State rk4Step(double t, State x, double dt) {
        State k1 = derivative(t, x);
        State k2 = derivative(t + 0.5 * dt, add(x, k1, 0.5 * dt));
        State k3 = derivative(t + 0.5 * dt, add(x, k2, 0.5 * dt));
        State k4 = derivative(t + dt, add(x, k3, dt));

        double nextX1 = x.x1 + (dt / 6.0) * (k1.x1 + 2.0 * k2.x1 + 2.0 * k3.x1 + k4.x1);
        double nextX2 = x.x2 + (dt / 6.0) * (k1.x2 + 2.0 * k2.x2 + 2.0 * k3.x2 + k4.x2);
        return new State(nextX1, nextX2);
    }

    static double frozenDecayMargin(double t) {
        double m = mass(t);
        double c = damping(t);
        double k = stiffness(t);

        double trace = -c / m;
        double determinant = k / m;
        double discriminant = trace * trace - 4.0 * determinant;

        if (discriminant >= 0.0) {
            double lambda1 = 0.5 * (trace + Math.sqrt(discriminant));
            double lambda2 = 0.5 * (trace - Math.sqrt(discriminant));
            return -Math.max(lambda1, lambda2);
        }

        return -0.5 * trace;
    }

    public static void main(String[] args) throws IOException {
        Locale.setDefault(Locale.US);

        double t0 = 0.0;
        double tf = 160.0;
        double dt = 0.01;
        int steps = (int) ((tf - t0) / dt);

        State x = new State(1.0, 0.0);

        try (PrintWriter writer = new PrintWriter(new FileWriter("Chapter29_Lesson5_java_output.csv"))) {
            writer.println("t,x1,x2,m,c,k,alpha");

            for (int i = 0; i <= steps; i++) {
                double t = t0 + i * dt;
                writer.printf(
                    "%.10f,%.10f,%.10f,%.10f,%.10f,%.10f,%.10f%n",
                    t, x.x1, x.x2, mass(t), damping(t), stiffness(t), frozenDecayMargin(t)
                );

                if (i < steps) {
                    x = rk4Step(t, x, dt);
                }
            }
        }

        System.out.printf("Final state: x1 = %.8f, x2 = %.8f%n", x.x1, x.x2);
        System.out.println("CSV written to Chapter29_Lesson5_java_output.csv");
    }
}

Chapter29_Lesson5.m

The MATLAB script uses ode45 and computes frozen-time eigenvalue margins. The comments also indicate the equivalent Simulink block structure using two integrators and time-varying coefficient blocks.

% Chapter29_Lesson5.m
%
% Slowly time-varying mass-spring-damper example for MATLAB/Simulink-oriented study.
%
% Required MATLAB functions:
%   ode45, eig, plot, writematrix
%
% Optional Simulink mapping:
%   x1_dot = x2
%   x2_dot = -(k(t)/m(t))*x1 - (c(t)/m(t))*x2 + u(t)/m(t)
% Use two Integrator blocks, MATLAB Function blocks for m(t), c(t), k(t), and u(t),
% and Sum/Gain/Product blocks from the Simulink Continuous and Math Operations libraries.

clear; clc; close all;

epsilon = 0.03;

mass = @(t) 1.0 + 0.12*sin(epsilon*t);
damping = @(t) 0.22 + 0.04*cos(epsilon*t);
stiffness = @(t) 2.0 + 0.30*sin(0.5*epsilon*t);
inputForce = @(t) 0.2*sin(0.7*t);

A = @(t) [0, 1; -stiffness(t)/mass(t), -damping(t)/mass(t)];
B = @(t) [0; 1/mass(t)];

rhs = @(t, x) A(t)*x + B(t)*inputForce(t);

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

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

alpha = zeros(size(t));
for i = 1:length(t)
    lambda = eig(A(t(i)));
    alpha(i) = -max(real(lambda));
end

output = [t, x(:,1), x(:,2), alpha];
writematrix(output, 'Chapter29_Lesson5_matlab_output.csv');

fprintf('Final state: x1 = %.8f, x2 = %.8f\n', x(end,1), x(end,2));
fprintf('Minimum frozen-time decay margin alpha = %.8f\n', min(alpha));

figure;
plot(t, x(:,1), 'LineWidth', 1.4);
grid on;
xlabel('time');
ylabel('displacement');
title('Slowly Time-Varying Mass-Spring-Damper System');

figure;
plot(t, alpha, 'LineWidth', 1.4);
grid on;
xlabel('time');
ylabel('frozen-time decay margin');
title('Frozen-Time Stability Margin');

Chapter29_Lesson5.nb

The Mathematica source solves the same LTV oscillator using NDSolveValue and plots the state and frozen-time decay margin.

(* Chapter29_Lesson5.nb *)

(* Slowly time-varying mass-spring-damper example for Wolfram Mathematica. *)
(* This text is a Wolfram Language notebook-style source saved with .nb extension. *)

ClearAll["Global`*"];

epsilon = 0.03;

m[t_] := 1.0 + 0.12 Sin[epsilon t];
c[t_] := 0.22 + 0.04 Cos[epsilon t];
k[t_] := 2.0 + 0.30 Sin[0.5 epsilon t];
u[t_] := 0.2 Sin[0.7 t];

A[t_] := { {0, 1}, {-k[t]/m[t], -c[t]/m[t]} };
B[t_] := {0, 1/m[t]};

solution = NDSolveValue[
   {
    x1'[t] == x2[t],
    x2'[t] == -(k[t]/m[t]) x1[t] - (c[t]/m[t]) x2[t] + u[t]/m[t],
    x1[0] == 1,
    x2[0] == 0
    },
   {x1, x2},
   {t, 0, 160},
   MaxStepFraction -> 1/2000
   ];

x1fun = solution[[1]];
x2fun = solution[[2]];

alpha[t_] := -Max[Re[Eigenvalues[A[t]]]];

finalState = {x1fun[160], x2fun[160]};
minDecayMargin = NMinimize[{alpha[s], 0 <= s <= 160}, s];

Print["Final state = ", finalState];
Print["Minimum frozen-time decay margin = ", minDecayMargin];

Plot[
 Evaluate[{x1fun[t], x2fun[t]}],
 {t, 0, 160},
 PlotLegends -> {"x1 displacement", "x2 velocity"},
 GridLines -> Automatic,
 AxesLabel -> {"time", "state"},
 PlotLabel -> "Slowly Time-Varying Mass-Spring-Damper System"
 ]

Plot[
 alpha[t],
 {t, 0, 160},
 GridLines -> Automatic,
 AxesLabel -> {"time", "alpha(t)"},
 PlotLabel -> "Frozen-Time Decay Margin"
 ]

Problem 1 (State-Space Form): Convert \( m(t)\ddot{q}+c(t)\dot{q}+k(t)q=u \) into first-order state-space form using \( x_1=q \) and \( x_2=\dot{q} \).

Solution:

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

Therefore,

\[ \dot{\mathbf{x} }= \begin{bmatrix}0 & 1\\ -k(t)/m(t) & -c(t)/m(t)\end{bmatrix} \mathbf{x}+ \begin{bmatrix}0\\1/m(t)\end{bmatrix}u. \]

Problem 2 (Frozen-Time Eigenvalues): For the mechanical system in Problem 1, derive the frozen-time characteristic polynomial.

Solution: Holding \( t=t_0 \),

\[ \det(\lambda\mathbf{I}-\mathbf{A}(t_0))= \lambda^2+\frac{c(t_0)}{m(t_0)}\lambda+ \frac{k(t_0)}{m(t_0)}. \]

The frozen-time poles are

\[ \lambda_{1,2}(t_0)= -\frac{c(t_0)}{2m(t_0)} \pm \frac{1}{2} \sqrt{\left(\frac{c(t_0)}{m(t_0)}\right)^2 -4\frac{k(t_0)}{m(t_0)} }. \]

Problem 3 (Slow-Variation Bound): Suppose \( \mathbf{A}(t)=\mathbf{A}_0+\varepsilon \sin(\varepsilon t)\mathbf{A}_1 \). Find a bound for \( \|\dot{\mathbf{A} }(t)\| \).

Solution:

\[ \dot{\mathbf{A} }(t)= \varepsilon^2\cos(\varepsilon t)\mathbf{A}_1. \]

Thus, using any submultiplicative matrix norm,

\[ \|\dot{\mathbf{A} }(t)\|\leq \varepsilon^2\|\mathbf{A}_1\|. \]

The coefficient drift is second order in \( \varepsilon \) for this particular parameterization.

Problem 4 (Lyapunov Stability Certificate): Let \( \mathbf{P}(t)=\mathbf{P}^{T}(t)\succ 0 \) satisfy \( m_1\mathbf{I}\preceq\mathbf{P}(t)\preceq m_2\mathbf{I} \) and \( \mathbf{A}^{T}\mathbf{P}+\mathbf{P}\mathbf{A}+\dot{\mathbf{P} } \preceq -q\mathbf{I} \). Prove uniform exponential stability.

Solution: With \( V=\mathbf{x}^{T}\mathbf{P}(t)\mathbf{x} \),

\[ \dot{V}\leq -q\|\mathbf{x}\|^2 \leq -\frac{q}{m_2}V. \]

Integrating gives

\[ V(t)\leq V(t_0)e^{-\frac{q}{m_2}(t-t_0)}. \]

Using the quadratic bounds on \( V \),

\[ \|\mathbf{x}(t)\|\leq \sqrt{\frac{m_2}{m_1} } e^{-\frac{q}{2m_2}(t-t_0)} \|\mathbf{x}(t_0)\|. \]

Problem 5 (Engineering Classification): Explain why a fuel-burning aircraft near a scheduled climb trajectory can be modeled as slowly time-varying.

Solution: The aircraft equations are nonlinear, but linearization about the trajectory gives \( \delta\dot{\mathbf{x} }=\mathbf{A}(t)\delta\mathbf{x}+ \mathbf{B}(t)\delta\mathbf{u} \). During a normal climb, mass, altitude, airspeed, and aerodynamic derivatives change gradually relative to short-period or phugoid modal periods. Hence the local linear model is not constant, but its parameters drift slowly enough that frozen-time analysis and gain scheduling are useful first approximations.

Slowly time-varying systems arise when the physical plant changes on a slow time scale: mass, inertia, stiffness, damping, trim condition, environmental parameters, or electrical/thermal coefficients drift while the state dynamics continue on a faster scale. Frozen-time models are valuable for interpretation and simulation, but rigorous conclusions require uniform stability, Lyapunov inequalities, or direct analysis of the state transition matrix. The practical engineering workflow is to model the coefficient drift, check local stability margins, estimate variation rates, simulate the full LTV dynamics, and compare against frozen-time approximations.