In the previous lesson, we introduced the homogeneous LTV state equation

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

For an LTI system, the solution is \( \mathbf{x}(t)=e^{\mathbf{A}(t-t_0)}\mathbf{x}_0 \). For an LTV system, a tempting but generally false guess is \( \mathbf{x}(t)= \exp\!\left(\int_{t_0}^t \mathbf{A}(\tau)\,d\tau\right) \mathbf{x}_0 \) . The obstruction is noncommutativity:

\[ \mathbf{A}(\tau_1)\mathbf{A}(\tau_2) \neq \mathbf{A}(\tau_2)\mathbf{A}(\tau_1) \quad \text{for some } \tau_1,\tau_2 . \]

Therefore, the solution operator must preserve the chronological order in which the matrices act on the state. This operator is the state transition matrix.

flowchart TD
  A["LTV homogeneous equation: x_dot = A(t) x"] --> B["Need a map from x(t0) to x(t)"]
  B --> C["Define Phi(t,t0)"]
  C --> D["Phi solves matrix IVP"]
  D --> E{"Do A(t1) and A(t2) \ncommute?"}
  E -->|"yes"| F["Exponential of integral may work"]
  E -->|"no"| G["Use time-ordered Peano-Baker series"]
  G --> H["Compute or approximate Phi(t,t0)"]
        

The state transition matrix \( \boldsymbol{\Phi}(t,t_0) \in \mathbb{R}^{n\times n} \) is the unique matrix satisfying

\[ \frac{\partial}{\partial t}\boldsymbol{\Phi}(t,t_0) =\mathbf{A}(t)\boldsymbol{\Phi}(t,t_0),\qquad \boldsymbol{\Phi}(t_0,t_0)=\mathbf{I} . \]

Once this matrix is known, the homogeneous solution is

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

The columns of \( \boldsymbol{\Phi}(t,t_0) \) are simply the solutions obtained by using the standard basis vectors as initial states:

\[ \boldsymbol{\Phi}(t,t_0) =\begin{bmatrix} \mathbf{x}_1(t) & \mathbf{x}_2(t) & \cdots & \mathbf{x}_n(t) \end{bmatrix}, \qquad \mathbf{x}_i(t_0)=\mathbf{e}_i . \]

Thus, computing the transition matrix is equivalent to solving \( n \) homogeneous initial-value problems at once.

The LTV transition matrix retains the most important structural properties of the LTI matrix exponential, but the order of the two time arguments matters.

Identity property:

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

Composition property: for any intermediate time \( t_1 \),

\[ \boxed{ \boldsymbol{\Phi}(t_2,t_0) = \boldsymbol{\Phi}(t_2,t_1) \boldsymbol{\Phi}(t_1,t_0) } . \]

Inverse property:

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

To prove the composition property, define \( \mathbf{Z}(t)= \boldsymbol{\Phi}(t,t_1)\boldsymbol{\Phi}(t_1,t_0) \) . Then

\[ \dot{\mathbf{Z}}(t) = \mathbf{A}(t)\boldsymbol{\Phi}(t,t_1) \boldsymbol{\Phi}(t_1,t_0) = \mathbf{A}(t)\mathbf{Z}(t), \qquad \mathbf{Z}(t_1)=\boldsymbol{\Phi}(t_1,t_0). \]

The matrix \( \boldsymbol{\Phi}(t,t_0) \) satisfies the same matrix IVP at \( t=t_1 \). Uniqueness of solutions gives the result.

Consider the nonhomogeneous LTV state equation

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

The solution is

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

This is the LTV counterpart of the LTI convolution integral. The kernel \( \boldsymbol{\Phi}(t,\tau) \) depends on both the observation time \( t \) and the input time \( \tau \); therefore, the response is generally not a time-invariant convolution.

Proof sketch. Let

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

Since \( \frac{\partial}{\partial t}\boldsymbol{\Phi}(t_0,t) =-\boldsymbol{\Phi}(t_0,t)\mathbf{A}(t) \) , differentiation gives

\[ \dot{\mathbf{z}}(t) = \boldsymbol{\Phi}(t_0,t)\mathbf{B}(t)\mathbf{u}(t). \]

Integrating from \( t_0 \) to \( t \) and multiplying by \( \boldsymbol{\Phi}(t,t_0) \) gives the stated formula.

Integrating the homogeneous matrix IVP gives the Volterra integral equation

\[ \boldsymbol{\Phi}(t,t_0) = \mathbf{I}+ \int_{t_0}^t \mathbf{A}(\tau_1) \boldsymbol{\Phi}(\tau_1,t_0) \,d\tau_1 . \]

Substitute this equation recursively into itself. The first terms are

\[ \begin{aligned} \boldsymbol{\Phi}(t,t_0) &= \mathbf{I} + \int_{t_0}^t \mathbf{A}(\tau_1)\,d\tau_1 \\ &\quad+ \int_{t_0}^t \mathbf{A}(\tau_1) \int_{t_0}^{\tau_1} \mathbf{A}(\tau_2)\,d\tau_2\,d\tau_1 \\ &\quad+ \int_{t_0}^t \mathbf{A}(\tau_1) \int_{t_0}^{\tau_1} \mathbf{A}(\tau_2) \int_{t_0}^{\tau_2} \mathbf{A}(\tau_3)\,d\tau_3\,d\tau_2\,d\tau_1 +\cdots . \end{aligned} \]

This infinite expansion is the Peano–Baker series:

\[ \boxed{ \boldsymbol{\Phi}(t,t_0) = \mathbf{I}+ \sum_{k=1}^\infty \int_{t_0}^t \int_{t_0}^{\tau_1} \cdots \int_{t_0}^{\tau_{k-1}} \mathbf{A}(\tau_1)\mathbf{A}(\tau_2)\cdots \mathbf{A}(\tau_k) \,d\tau_k\cdots d\tau_2\,d\tau_1 } . \]

The product order is essential: the latest-time matrix appears at the left and the earliest-time matrix appears at the right.

Suppose \( \mathbf{A}(t) \) is continuous on \( [t_0,T] \). Let

\[ M=\max_{t_0\le t\le T}\|\mathbf{A}(t)\| . \]

The \( k \)-th Peano–Baker term is bounded by

\[ \left\| \int_{t_0}^t \int_{t_0}^{\tau_1} \cdots \int_{t_0}^{\tau_{k-1}} \mathbf{A}(\tau_1)\cdots\mathbf{A}(\tau_k) \,d\tau_k\cdots d\tau_1 \right\| \le \frac{\left(M|t-t_0|\right)^k}{k!} . \]

Therefore, the series is absolutely and uniformly convergent on compact time intervals. If \( \boldsymbol{\Phi}_q(t,t_0) \) denotes the truncation through order \( q \), then a conservative remainder bound is

\[ \left\|\boldsymbol{\Phi}(t,t_0)- \boldsymbol{\Phi}_q(t,t_0)\right\| \le \sum_{k=q+1}^\infty \frac{\left(M|t-t_0|\right)^k}{k!} . \]

This bound is often pessimistic, but it proves convergence and explains why truncation works best on short time intervals or when \( \|\mathbf{A}(t)\| \) is modest.

If the matrix family commutes at all pairs of times,

\[ \mathbf{A}(\tau_1)\mathbf{A}(\tau_2) = \mathbf{A}(\tau_2)\mathbf{A}(\tau_1) \qquad \text{for all } \tau_1,\tau_2 , \]

then the Peano–Baker series collapses to

\[ \boldsymbol{\Phi}(t,t_0)= \exp\!\left( \int_{t_0}^t \mathbf{A}(\tau)\,d\tau \right). \]

This includes the LTI case \( \mathbf{A}(t)=\mathbf{A} \). However, if commutation fails, the above formula is generally wrong. The second-order Peano–Baker term contains the ordered product \( \mathbf{A}(\tau_1)\mathbf{A}(\tau_2) \), while the ordinary exponential would symmetrize the product order.

flowchart TD
  S["Given A(t) on interval"] --> Q{"Check pairwise \ncommutation"}
  Q -->|"commutes"| C["Use exp of integral"]
  Q -->|"does not commute"| P["Use Peano-Baker or numerical matrix IVP"]
  P --> R["Short interval: \ntruncate series"]
  P --> N["Long interval: \nintegrate Phi_dot = A(t) Phi"]
  R --> V["Validate by composition property"]
  N --> V
        

For numerical work, the most direct method is to integrate

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

This is an \( n^2 \)-dimensional linear differential equation. In vectorized form,

\[ \frac{d}{dt}\operatorname{vec} \boldsymbol{\Phi}(t,t_0) = \left(\mathbf{I}\otimes\mathbf{A}(t)\right) \operatorname{vec}\boldsymbol{\Phi}(t,t_0), \]

where column-wise vectorization is assumed. The Peano–Baker recursion is also useful pedagogically and for short-interval approximations:

\[ \mathbf{T}_0(t)=\mathbf{I},\qquad \mathbf{T}_k(t)= \int_{t_0}^t \mathbf{A}(s)\mathbf{T}_{k-1}(s)\,ds,\qquad \boldsymbol{\Phi}_q(t,t_0)=\sum_{k=0}^q \mathbf{T}_k(t). \]

This implementation compares direct matrix-IVP integration with a truncated Peano–Baker approximation.

Chapter29_Lesson2.py

# Chapter29_Lesson2.py
# Time-varying state transition matrix for x_dot = A(t) x
# Methods:
#   1) solve dPhi/dt = A(t) Phi, Phi(t0,t0)=I
#   2) approximate the Peano-Baker series by recursive quadrature

import numpy as np
from scipy.integrate import solve_ivp


def A(t: float) -> np.ndarray:
    """A smooth 2x2 time-varying state matrix."""
    return np.array(
        [
            [0.0, 1.0],
            [-2.0 - 0.5 * np.sin(t), -0.4 + 0.2 * np.cos(2.0 * t)],
        ],
        dtype=float,
    )


def transition_matrix_ode(t0: float, tf: float, n: int = 2) -> np.ndarray:
    """Compute Phi(tf,t0) by integrating the matrix IVP."""
    def rhs(t, y):
        Phi = y.reshape(n, n)
        dPhi = A(t) @ Phi
        return dPhi.reshape(-1)

    y0 = np.eye(n).reshape(-1)
    sol = solve_ivp(rhs, (t0, tf), y0, rtol=1e-10, atol=1e-12)
    if not sol.success:
        raise RuntimeError(sol.message)
    return sol.y[:, -1].reshape(n, n)


def peano_baker(t0: float, tf: float, order: int = 6, steps: int = 4000) -> np.ndarray:
    """
    Approximate Phi(tf,t0) = I + T1(tf) + T2(tf) + ...
    using the recursion
        T_0(t) = I,
        T_k(t) = integral_{t0}^{t} A(s) T_{k-1}(s) ds.

    The nested time ordering is embedded in the recursion.
    """
    if order < 0:
        raise ValueError("order must be nonnegative")
    if steps < 1:
        raise ValueError("steps must be positive")

    n = 2
    grid = np.linspace(t0, tf, steps + 1)
    h = (tf - t0) / steps

    terms = [[np.zeros((n, n)) for _ in range(steps + 1)] for _ in range(order + 1)]
    for j in range(steps + 1):
        terms[0][j] = np.eye(n)

    for k in range(1, order + 1):
        terms[k][0] = np.zeros((n, n))
        for j in range(1, steps + 1):
            smid = 0.5 * (grid[j - 1] + grid[j])
            # Left endpoint for previous term; midpoint for A(t).
            integrand = A(smid) @ terms[k - 1][j - 1]
            terms[k][j] = terms[k][j - 1] + h * integrand

    Phi = sum(terms[k][-1] for k in range(order + 1))
    return Phi


def check_composition(t0: float, tm: float, tf: float) -> None:
    """Check Phi(tf,t0) = Phi(tf,tm) Phi(tm,t0)."""
    Phi_f0 = transition_matrix_ode(t0, tf)
    Phi_fm = transition_matrix_ode(tm, tf)
    Phi_m0 = transition_matrix_ode(t0, tm)
    err = np.linalg.norm(Phi_f0 - Phi_fm @ Phi_m0, ord="fro")
    print("composition error:", err)


if __name__ == "__main__":
    np.set_printoptions(precision=8, suppress=True)

    t0, tf = 0.0, 1.5
    Phi_ode = transition_matrix_ode(t0, tf)
    print("Phi(tf,t0) from matrix ODE:")
    print(Phi_ode)

    for q in [1, 2, 3, 4, 6, 8]:
        Phi_pb = peano_baker(t0, tf, order=q, steps=5000)
        err = np.linalg.norm(Phi_pb - Phi_ode, ord="fro")
        print(f"Peano-Baker order {q:2d} error: {err:.6e}")

    check_composition(0.0, 0.7, 1.5)
      

This C++ version uses no external linear-algebra library; it implements only the required 2-by-2 matrix operations and RK4 integration.

Chapter29_Lesson2.cpp

// Chapter29_Lesson2.cpp
// From-scratch 2x2 implementation for the LTV transition matrix.
// Compile:
//   g++ -std=c++17 -O2 Chapter29_Lesson2.cpp -o Chapter29_Lesson2

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

struct Mat2 {
    std::array<double, 4> a{};
};

Mat2 zero() {
    return Mat2{{0.0, 0.0, 0.0, 0.0}};
}

Mat2 eye() {
    return Mat2{{1.0, 0.0, 0.0, 1.0}};
}

Mat2 add(const Mat2& X, const Mat2& Y) {
    return Mat2{{X.a[0] + Y.a[0], X.a[1] + Y.a[1],
                 X.a[2] + Y.a[2], X.a[3] + Y.a[3]}};
}

Mat2 sub(const Mat2& X, const Mat2& Y) {
    return Mat2{{X.a[0] - Y.a[0], X.a[1] - Y.a[1],
                 X.a[2] - Y.a[2], X.a[3] - Y.a[3]}};
}

Mat2 scale(double c, const Mat2& X) {
    return Mat2{{c * X.a[0], c * X.a[1], c * X.a[2], c * X.a[3]}};
}

Mat2 mul(const Mat2& X, const Mat2& Y) {
    return Mat2{{
        X.a[0] * Y.a[0] + X.a[1] * Y.a[2],
        X.a[0] * Y.a[1] + X.a[1] * Y.a[3],
        X.a[2] * Y.a[0] + X.a[3] * Y.a[2],
        X.a[2] * Y.a[1] + X.a[3] * Y.a[3]
    }};
}

double normF(const Mat2& X) {
    double s = 0.0;
    for (double v : X.a) s += v * v;
    return std::sqrt(s);
}

Mat2 A(double t) {
    return Mat2{{
        0.0, 1.0,
        -2.0 - 0.5 * std::sin(t), -0.4 + 0.2 * std::cos(2.0 * t)
    }};
}

Mat2 rhs(double t, const Mat2& Phi) {
    return mul(A(t), Phi);
}

Mat2 rk4Phi(double t0, double tf, int steps) {
    double h = (tf - t0) / static_cast<double>(steps);
    double t = t0;
    Mat2 Phi = eye();

    for (int i = 0; i < steps; ++i) {
        Mat2 k1 = rhs(t, Phi);
        Mat2 k2 = rhs(t + 0.5 * h, add(Phi, scale(0.5 * h, k1)));
        Mat2 k3 = rhs(t + 0.5 * h, add(Phi, scale(0.5 * h, k2)));
        Mat2 k4 = rhs(t + h, add(Phi, scale(h, k3)));

        Mat2 incr = scale(h / 6.0, add(add(k1, scale(2.0, k2)), add(scale(2.0, k3), k4)));
        Phi = add(Phi, incr);
        t += h;
    }
    return Phi;
}

Mat2 peanoBaker(double t0, double tf, int order, int steps) {
    double h = (tf - t0) / static_cast<double>(steps);
    std::vector<std::vector<Mat2>> terms(order + 1, std::vector<Mat2>(steps + 1, zero()));

    for (int j = 0; j <= steps; ++j) terms[0][j] = eye();

    for (int k = 1; k <= order; ++k) {
        terms[k][0] = zero();
        for (int j = 1; j <= steps; ++j) {
            double smid = t0 + (j - 0.5) * h;
            Mat2 integrand = mul(A(smid), terms[k - 1][j - 1]);
            terms[k][j] = add(terms[k][j - 1], scale(h, integrand));
        }
    }

    Mat2 Phi = zero();
    for (int k = 0; k <= order; ++k) Phi = add(Phi, terms[k][steps]);
    return Phi;
}

void printMat(const Mat2& X) {
    std::cout << "[" << std::setw(12) << X.a[0] << " " << std::setw(12) << X.a[1] << "]\n";
    std::cout << "[" << std::setw(12) << X.a[2] << " " << std::setw(12) << X.a[3] << "]\n";
}

int main() {
    std::cout << std::setprecision(8) << std::fixed;

    double t0 = 0.0;
    double tf = 1.5;
    Mat2 PhiRK = rk4Phi(t0, tf, 20000);

    std::cout << "Phi(tf,t0) by RK4 matrix IVP:\n";
    printMat(PhiRK);

    for (int order : {1, 2, 3, 4, 6, 8}) {
        Mat2 PhiPB = peanoBaker(t0, tf, order, 6000);
        std::cout << "Peano-Baker order " << order
                  << " Frobenius error = " << normF(sub(PhiPB, PhiRK)) << "\n";
    }

    Mat2 Phi60 = rk4Phi(0.0, 1.5, 20000);
    Mat2 Phi62 = rk4Phi(0.7, 1.5, 12000);
    Mat2 Phi20 = rk4Phi(0.0, 0.7, 8000);
    std::cout << "Composition error = " << normF(sub(Phi60, mul(Phi62, Phi20))) << "\n";
    return 0;
}
      

The Java implementation mirrors the C++ version and is useful for checking the algorithm independent of scientific-computing libraries.

Chapter29_Lesson2.java

// Chapter29_Lesson2.java
// From-scratch 2x2 implementation for the LTV transition matrix.
// Compile and run:
//   javac Chapter29_Lesson2.java
//   java Chapter29_Lesson2

public class Chapter29_Lesson2 {
    static final class Mat2 {
        final double a00, a01, a10, a11;

        Mat2(double a00, double a01, double a10, double a11) {
            this.a00 = a00;
            this.a01 = a01;
            this.a10 = a10;
            this.a11 = a11;
        }

        static Mat2 zero() {
            return new Mat2(0.0, 0.0, 0.0, 0.0);
        }

        static Mat2 eye() {
            return new Mat2(1.0, 0.0, 0.0, 1.0);
        }

        Mat2 add(Mat2 b) {
            return new Mat2(a00 + b.a00, a01 + b.a01, a10 + b.a10, a11 + b.a11);
        }

        Mat2 sub(Mat2 b) {
            return new Mat2(a00 - b.a00, a01 - b.a01, a10 - b.a10, a11 - b.a11);
        }

        Mat2 scale(double c) {
            return new Mat2(c * a00, c * a01, c * a10, c * a11);
        }

        Mat2 multiply(Mat2 b) {
            return new Mat2(
                a00 * b.a00 + a01 * b.a10,
                a00 * b.a01 + a01 * b.a11,
                a10 * b.a00 + a11 * b.a10,
                a10 * b.a01 + a11 * b.a11
            );
        }

        double normF() {
            return Math.sqrt(a00 * a00 + a01 * a01 + a10 * a10 + a11 * a11);
        }

        void print() {
            System.out.printf("[%12.8f %12.8f]%n", a00, a01);
            System.out.printf("[%12.8f %12.8f]%n", a10, a11);
        }
    }

    static Mat2 A(double t) {
        return new Mat2(
            0.0, 1.0,
            -2.0 - 0.5 * Math.sin(t), -0.4 + 0.2 * Math.cos(2.0 * t)
        );
    }

    static Mat2 rhs(double t, Mat2 phi) {
        return A(t).multiply(phi);
    }

    static Mat2 rk4Phi(double t0, double tf, int steps) {
        double h = (tf - t0) / steps;
        double t = t0;
        Mat2 phi = Mat2.eye();

        for (int i = 0; i < steps; i++) {
            Mat2 k1 = rhs(t, phi);
            Mat2 k2 = rhs(t + 0.5 * h, phi.add(k1.scale(0.5 * h)));
            Mat2 k3 = rhs(t + 0.5 * h, phi.add(k2.scale(0.5 * h)));
            Mat2 k4 = rhs(t + h, phi.add(k3.scale(h)));

            Mat2 incr = k1.add(k2.scale(2.0)).add(k3.scale(2.0)).add(k4).scale(h / 6.0);
            phi = phi.add(incr);
            t += h;
        }
        return phi;
    }

    static Mat2 peanoBaker(double t0, double tf, int order, int steps) {
        double h = (tf - t0) / steps;
        Mat2[][] terms = new Mat2[order + 1][steps + 1];

        for (int k = 0; k <= order; k++) {
            for (int j = 0; j <= steps; j++) {
                terms[k][j] = Mat2.zero();
            }
        }
        for (int j = 0; j <= steps; j++) terms[0][j] = Mat2.eye();

        for (int k = 1; k <= order; k++) {
            terms[k][0] = Mat2.zero();
            for (int j = 1; j <= steps; j++) {
                double smid = t0 + (j - 0.5) * h;
                Mat2 integrand = A(smid).multiply(terms[k - 1][j - 1]);
                terms[k][j] = terms[k][j - 1].add(integrand.scale(h));
            }
        }

        Mat2 phi = Mat2.zero();
        for (int k = 0; k <= order; k++) phi = phi.add(terms[k][steps]);
        return phi;
    }

    public static void main(String[] args) {
        double t0 = 0.0;
        double tf = 1.5;

        Mat2 phiRK = rk4Phi(t0, tf, 20000);
        System.out.println("Phi(tf,t0) by RK4 matrix IVP:");
        phiRK.print();

        int[] orders = {1, 2, 3, 4, 6, 8};
        for (int order : orders) {
            Mat2 phiPB = peanoBaker(t0, tf, order, 6000);
            System.out.printf("Peano-Baker order %d Frobenius error = %.8e%n",
                              order, phiPB.sub(phiRK).normF());
        }

        Mat2 phi60 = rk4Phi(0.0, 1.5, 20000);
        Mat2 phi62 = rk4Phi(0.7, 1.5, 12000);
        Mat2 phi20 = rk4Phi(0.0, 0.7, 8000);
        System.out.printf("Composition error = %.8e%n", phi60.sub(phi62.multiply(phi20)).normF());
    }
}
      

The MATLAB script integrates the vectorized transition-matrix equation. The same right-hand side can be placed in a Simulink MATLAB Function block with an Integrator initialized by \( \operatorname{vec}(\mathbf{I}) \).

Chapter29_Lesson2.m

% Chapter29_Lesson2.m
% Time-varying state transition matrix for x_dot = A(t) x.
% Includes a MATLAB computation and a Simulink-ready vectorized RHS.

clear; clc;

t0 = 0.0;
tf = 1.5;
n = 2;

Phi0 = eye(n);
y0 = Phi0(:);

opts = odeset('RelTol', 1e-10, 'AbsTol', 1e-12);
[t, y] = ode45(@ltv_phi_rhs, [t0 tf], y0, opts);

Phi_ode = reshape(y(end, :), n, n);
disp('Phi(tf,t0) from MATLAB ode45 matrix IVP:');
disp(Phi_ode);

orders = [1 2 3 4 6 8];
for q = orders
    Phi_pb = peano_baker(t0, tf, q, 5000);
    err = norm(Phi_pb - Phi_ode, 'fro');
    fprintf('Peano-Baker order %2d error: %.6e\n', q, err);
end

Phi60 = transition_matrix_ode(0.0, 1.5);
Phi62 = transition_matrix_ode(0.7, 1.5);
Phi20 = transition_matrix_ode(0.0, 0.7);
fprintf('Composition error: %.6e\n', norm(Phi60 - Phi62 * Phi20, 'fro'));

% Simulink use:
% Put the ltv_phi_rhs equations inside a MATLAB Function block.
% Feed Clock time t and the 4-state vector PhiVec into the block.
% The output dPhiVec integrates through an Integrator block initialized by eye(2)(:).

function dy = ltv_phi_rhs(t, y)
    Phi = reshape(y, 2, 2);
    dPhi = A_matrix(t) * Phi;
    dy = dPhi(:);
end

function A = A_matrix(t)
    A = [0, 1;
        -2 - 0.5 * sin(t), -0.4 + 0.2 * cos(2 * t)];
end

function Phi = transition_matrix_ode(t0, tf)
    y0 = eye(2);
    y0 = y0(:);
    opts = odeset('RelTol', 1e-10, 'AbsTol', 1e-12);
    [~, y] = ode45(@ltv_phi_rhs, [t0 tf], y0, opts);
    Phi = reshape(y(end, :), 2, 2);
end

function Phi = peano_baker(t0, tf, order, steps)
    h = (tf - t0) / steps;
    grid = linspace(t0, tf, steps + 1);

    terms = cell(order + 1, steps + 1);
    for j = 1:(steps + 1)
        terms{1, j} = eye(2);
    end

    for k = 2:(order + 1)
        terms{k, 1} = zeros(2);
        for j = 2:(steps + 1)
            smid = 0.5 * (grid(j - 1) + grid(j));
            integrand = A_matrix(smid) * terms{k - 1, j - 1};
            terms{k, j} = terms{k, j - 1} + h * integrand;
        end
    end

    Phi = zeros(2);
    for k = 1:(order + 1)
        Phi = Phi + terms{k, steps + 1};
    end
end
      

Mathematica can compute the same transition matrix using NDSolveValue and then compare it against the truncated Peano–Baker approximation.

Chapter29_Lesson2.nb

(* Chapter29_Lesson2.nb *)
(* Time-varying state transition matrix and Peano-Baker approximation. *)

ClearAll["Global`*"];

A[t_] := {{0, 1}, {-2 - 1/2 Sin[t], -2/5 + 1/5 Cos[2 t]}};

t0 = 0;
tf = 1.5;
n = 2;

phiSol = NDSolveValue[
   {Phi'[t] == A[t].Phi[t], Phi[t0] == IdentityMatrix[n]},
   Phi,
   {t, t0, tf},
   WorkingPrecision -> 30,
   AccuracyGoal -> 18,
   PrecisionGoal -> 18
];

PhiODE = N[phiSol[tf], 12];

PeanoBaker[t0_, tf_, order_Integer, steps_Integer] := Module[
   {h, grid, terms, smid, integrand, k, j},
   h = N[(tf - t0)/steps];
   grid = N[Subdivide[t0, tf, steps]];
   terms = Table[ConstantArray[0, {2, 2}], {order + 1}, {steps + 1}];

   Do[terms[[1, j]] = IdentityMatrix[2], {j, 1, steps + 1}];

   Do[
    terms[[k, 1]] = ConstantArray[0, {2, 2}];
    Do[
     smid = 0.5 (grid[[j - 1]] + grid[[j]]);
     integrand = A[smid].terms[[k - 1, j - 1]];
     terms[[k, j]] = terms[[k, j - 1]] + h integrand;
     ,
     {j, 2, steps + 1}
     ];
    ,
    {k, 2, order + 1}
    ];

   Total[Table[terms[[k, steps + 1]], {k, 1, order + 1}]]
];

Print["Phi(tf,t0) from NDSolve:"];
Print[MatrixForm[PhiODE]];

Do[
  PhiPB = N[PeanoBaker[t0, tf, q, 5000], 12];
  err = Norm[PhiPB - PhiODE, "Frobenius"];
  Print["Peano-Baker order ", q, " error: ", ScientificForm[err, 6]];
  ,
  {q, {1, 2, 3, 4, 6, 8}}
];

Phi60 = NDSolveValue[
   {P'[t] == A[t].P[t], P[0] == IdentityMatrix[2]},
   P, {t, 0, 1.5}
][1.5];

Phi62 = NDSolveValue[
   {P'[t] == A[t].P[t], P[0.7] == IdentityMatrix[2]},
   P, {t, 0.7, 1.5}
][1.5];

Phi20 = NDSolveValue[
   {P'[t] == A[t].P[t], P[0] == IdentityMatrix[2]},
   P, {t, 0, 0.7}
][0.7];

Print["Composition error: ", Norm[Phi60 - Phi62.Phi20, "Frobenius"]];
      

Problem 1 (Composition property): Prove that \( \boldsymbol{\Phi}(t_2,t_0)= \boldsymbol{\Phi}(t_2,t_1)\boldsymbol{\Phi}(t_1,t_0) \) .

Solution: Define \( \mathbf{Z}(t)= \boldsymbol{\Phi}(t,t_1)\boldsymbol{\Phi}(t_1,t_0) \) . Then \( \dot{\mathbf{Z}}(t)=\mathbf{A}(t)\mathbf{Z}(t) \) and \( \mathbf{Z}(t_1)=\boldsymbol{\Phi}(t_1,t_0) \). The matrix \( \boldsymbol{\Phi}(t,t_0) \) satisfies the same matrix IVP at \( t=t_1 \). By uniqueness, \( \mathbf{Z}(t)=\boldsymbol{\Phi}(t,t_0) \), and setting \( t=t_2 \) proves the claim.

Problem 2 (Scalar LTV system): For \( \dot{x}(t)=a(t)x(t) \), show that the transition matrix is \( \Phi(t,t_0)=\exp\!\left(\int_{t_0}^t a(\tau)d\tau\right) \) .

Solution: Since scalar multiplication commutes at all times, the Peano–Baker series becomes

\[ \Phi(t,t_0) = 1+ \int_{t_0}^t a(\tau_1)d\tau_1+ \frac{1}{2!} \left(\int_{t_0}^t a(\tau)d\tau\right)^2+\cdots = \exp\!\left(\int_{t_0}^t a(\tau)d\tau\right). \]

Problem 3 (Commuting matrix family): Suppose \( \mathbf{A}(\tau_1)\mathbf{A}(\tau_2)= \mathbf{A}(\tau_2)\mathbf{A}(\tau_1) \) for every pair of times. Prove that \( \boldsymbol{\Phi}(t,t_0)= \exp\!\left(\int_{t_0}^t\mathbf{A}(\tau)d\tau\right) \) .

Solution: Let \( \mathbf{G}(t)=\int_{t_0}^t\mathbf{A}(\tau)d\tau \). The commutation hypothesis implies \( \mathbf{G}(t)\mathbf{A}(t)=\mathbf{A}(t)\mathbf{G}(t) \). Therefore,

\[ \frac{d}{dt}e^{\mathbf{G}(t)} = \dot{\mathbf{G}}(t)e^{\mathbf{G}(t)} = \mathbf{A}(t)e^{\mathbf{G}(t)} . \]

Also \( e^{\mathbf{G}(t_0)}=\mathbf{I} \). Thus \( e^{\mathbf{G}(t)} \) satisfies the defining matrix IVP for \( \boldsymbol{\Phi}(t,t_0) \).

Problem 4 (First two Peano–Baker terms): Compute the Peano–Baker approximation through second order for \( \mathbf{A}(t)=\mathbf{A}_0+t\mathbf{A}_1 \) .

Solution: Let \( h=t-t_0 \). The first-order term is

\[ \mathbf{T}_1 = \int_{t_0}^t (\mathbf{A}_0+\tau_1\mathbf{A}_1)d\tau_1 = h\mathbf{A}_0+ \frac{t^2-t_0^2}{2}\mathbf{A}_1 . \]

The second-order term is

\[ \mathbf{T}_2 = \int_{t_0}^t (\mathbf{A}_0+\tau_1\mathbf{A}_1) \int_{t_0}^{\tau_1} (\mathbf{A}_0+\tau_2\mathbf{A}_1) d\tau_2d\tau_1 . \]

If \( \mathbf{A}_0 \) and \( \mathbf{A}_1 \) commute, this equals \( \frac{1}{2}\mathbf{T}_1^2 \) . If they do not commute, the ordered integral must be kept in the displayed form or evaluated by quadrature.

Problem 5 (Forced response): For \( \dot{\mathbf{x}}(t)=\mathbf{A}(t)\mathbf{x}(t)+ \mathbf{B}(t)\mathbf{u}(t) \) , verify the variation-of-constants formula by direct differentiation.

Solution: Let

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

Differentiating the first term gives \( \mathbf{A}(t)\boldsymbol{\Phi}(t,t_0)\mathbf{x}_0 \). For the integral, Leibniz's rule gives a boundary term \( \mathbf{B}(t)\mathbf{u}(t) \) and an integral term \( \int_{t_0}^t \mathbf{A}(t)\boldsymbol{\Phi}(t,\tau) \mathbf{B}(\tau)\mathbf{u}(\tau)d\tau \) . Combining terms yields \( \dot{\mathbf{x}}(t)= \mathbf{A}(t)\mathbf{x}(t)+\mathbf{B}(t)\mathbf{u}(t) \) .

The LTV state transition matrix \( \boldsymbol{\Phi}(t,t_0) \) maps initial states at \( t_0 \) to states at \( t \). It satisfies a matrix differential equation, has identity, inverse, and composition properties, and gives the variation-of-constants formula for forced systems. The Peano–Baker series generalizes the LTI matrix exponential by preserving the time order of matrix products. When \( \mathbf{A}(t) \) commutes with itself at all times, the series reduces to the exponential of the integral; otherwise, direct matrix-IVP integration or time-ordered approximations are needed.