Consider the (possibly time-varying) linear homogeneous system from Lesson 2: \( \dot{\mathbf{x}}(t)=\mathbf{A}(t)\mathbf{x}(t) \), where \( \mathbf{x}(t)\in\mathbb{R}^n \) and \( \mathbf{A}(t)\in\mathbb{R}^{n\times n} \). A central idea is to collect n linearly independent solutions as columns of a matrix \( \mathbf{\Phi}(t)\in\mathbb{R}^{n\times n} \). This matrix is called a fundamental matrix.

flowchart TD
  S["Start: x_dot = A(t) x, given t0 and x0"] --> F["Build a fundamental matrix Phi(t)"]
  F --> C["Compute constant vector c = inv(Phi(t0)) * x0"]
  C --> X["Solution: x(t) = Phi(t) * c"]
  X --> V["Validate: substitute into x_dot = A(t) x"]
        

When \( \mathbf{A}(t)=\mathbf{A} \) is constant, Lesson 3 introduced the matrix exponential \( e^{\mathbf{A}t} \). In that special case, a canonical fundamental matrix is \( \mathbf{\Phi}(t)=e^{\mathbf{A}(t-t_0)} \) (normalized so that \( \mathbf{\Phi}(t_0)=\mathbf{I} \)).

Definition (Fundamental Matrix). Let \( \mathbf{x}_1(t),\dots,\mathbf{x}_n(t) \) be solutions of \( \dot{\mathbf{x}}=\mathbf{A}(t)\mathbf{x} \). Define \( \mathbf{\Phi}(t)=[\mathbf{x}_1(t)\ \cdots\ \mathbf{x}_n(t)] \). If the columns are linearly independent for some (hence all) times in an interval, then \( \mathbf{\Phi}(t) \) is a fundamental matrix on that interval.

Key matrix differential equation. Because differentiation acts columnwise, \( \mathbf{\Phi}(t) \) satisfies

\[ \frac{d}{dt}\mathbf{\Phi}(t) = \mathbf{A}(t)\mathbf{\Phi}(t). \]

Conversely, any matrix solution \( \mathbf{\Phi}(t) \) of \( \dot{\mathbf{\Phi}}=\mathbf{A}(t)\mathbf{\Phi} \) with \( \det(\mathbf{\Phi}(t_0))\neq 0 \) is a fundamental matrix (its columns are independent).

We now prove the most important structural result: every solution is a linear combination of the columns of a fundamental matrix, with constant coefficients.

Theorem (General Solution via Fundamental Matrix). If \( \mathbf{\Phi}(t) \) is a fundamental matrix for \( \dot{\mathbf{x}}=\mathbf{A}(t)\mathbf{x} \), then any solution \( \mathbf{x}(t) \) can be written uniquely as \( \mathbf{x}(t)=\mathbf{\Phi}(t)\mathbf{c} \) for some constant vector \( \mathbf{c}\in\mathbb{R}^n \).

Proof. Fix any time \( t_\star \) in the interval. Since the columns of \( \mathbf{\Phi}(t_\star) \) are linearly independent, the matrix \( \mathbf{\Phi}(t_\star) \) is invertible, hence any vector \( \mathbf{x}(t_\star) \) can be expressed uniquely as \( \mathbf{x}(t_\star)=\mathbf{\Phi}(t_\star)\mathbf{c} \) with \( \mathbf{c}=\mathbf{\Phi}(t_\star)^{-1}\mathbf{x}(t_\star) \).

Define \( \tilde{\mathbf{x}}(t)=\mathbf{\Phi}(t)\mathbf{c} \). Then

\[ \dot{\tilde{\mathbf{x}}}(t) = \dot{\mathbf{\Phi}}(t)\mathbf{c} = \mathbf{A}(t)\mathbf{\Phi}(t)\mathbf{c} = \mathbf{A}(t)\tilde{\mathbf{x}}(t), \]

so \( \tilde{\mathbf{x}}(t) \) is also a solution and matches \( \mathbf{x}(t) \) at \( t_\star \). By the uniqueness principle for linear IVPs (proved formally in Lesson 5), they coincide on the interval, hence \( \mathbf{x}(t)=\mathbf{\Phi}(t)\mathbf{c} \). Uniqueness of \( \mathbf{c} \) follows from invertibility of \( \mathbf{\Phi}(t_\star) \). □

Corollary (IVP Solution). For the IVP \( \dot{\mathbf{x}}(t)=\mathbf{A}(t)\mathbf{x}(t),\ \mathbf{x}(t_0)=\mathbf{x}_0 \), if \( \mathbf{\Phi}(t) \) is any fundamental matrix, then

\[ \mathbf{x}(t) = \mathbf{\Phi}(t)\mathbf{\Phi}(t_0)^{-1}\mathbf{x}_0. \]

The product \( \mathbf{\Phi}(t)\mathbf{\Phi}(t_0)^{-1} \) is the normalized mapping that takes \( \mathbf{x}_0 \) at time \( t_0 \) to the state at time \( t \). In later chapters this object will be given a dedicated name and additional properties; here we focus on the fundamental-matrix logic that produces it.

Theorem. If \( \mathbf{\Phi}_1(t) \) and \( \mathbf{\Phi}_2(t) \) are two fundamental matrices for the same system \( \dot{\mathbf{x}}=\mathbf{A}(t)\mathbf{x} \) on an interval, then there exists a constant invertible matrix \( \mathbf{C} \) such that

\[ \mathbf{\Phi}_2(t)=\mathbf{\Phi}_1(t)\mathbf{C}\quad \text{for all } t. \]

Proof. Define \( \mathbf{C}(t)=\mathbf{\Phi}_1(t)^{-1}\mathbf{\Phi}_2(t) \). Differentiate using the product rule and the identity \( \frac{d}{dt}\mathbf{\Phi}_1^{-1}(t)=-\mathbf{\Phi}_1^{-1}(t)\dot{\mathbf{\Phi}}_1(t)\mathbf{\Phi}_1^{-1}(t) \):

\[ \dot{\mathbf{C}}(t) = \frac{d}{dt}\!\big(\mathbf{\Phi}_1^{-1}(t)\big)\mathbf{\Phi}_2(t) + \mathbf{\Phi}_1^{-1}(t)\dot{\mathbf{\Phi}}_2(t) = -\mathbf{\Phi}_1^{-1}\mathbf{A}\mathbf{\Phi}_1\mathbf{\Phi}_1^{-1}\mathbf{\Phi}_2 + \mathbf{\Phi}_1^{-1}\mathbf{A}\mathbf{\Phi}_2 = \mathbf{0}. \]

Hence \( \mathbf{C}(t) \) is constant: \( \mathbf{C}(t)\equiv \mathbf{C} \). Since both fundamental matrices are invertible, \( \mathbf{C} \) is invertible. □

Practical meaning. Choosing a fundamental matrix is like choosing a basis of solutions; changing the basis multiplies by a constant invertible matrix.

The determinant \( \det(\mathbf{\Phi}(t)) \) plays the role of a matrix Wronskian. A fundamental matrix must remain invertible; Liouville’s formula explains how its determinant evolves.

Lemma (Jacobi’s Formula). If \( \mathbf{M}(t) \) is differentiable and invertible, then

\[ \frac{d}{dt}\det(\mathbf{M}(t)) = \det(\mathbf{M}(t))\,\operatorname{tr}\!\big(\mathbf{M}(t)^{-1}\dot{\mathbf{M}}(t)\big). \]

Sketch of proof. Consider \( g(t)=\ln\det(\mathbf{M}(t)) \). Using the first-order expansion \( \det(\mathbf{I}+\epsilon\mathbf{K})=1+\epsilon\,\operatorname{tr}(\mathbf{K})+o(\epsilon) \), one obtains \( \frac{d}{dt}\ln\det(\mathbf{M})=\operatorname{tr}(\mathbf{M}^{-1}\dot{\mathbf{M}}) \), and then multiply by \( \det(\mathbf{M}) \). □

Theorem (Liouville). If \( \mathbf{\Phi}(t) \) satisfies \( \dot{\mathbf{\Phi}}(t)=\mathbf{A}(t)\mathbf{\Phi}(t) \), then

\[ \frac{d}{dt}\det(\mathbf{\Phi}(t)) = \operatorname{tr}\!\big(\mathbf{A}(t)\big)\,\det(\mathbf{\Phi}(t)). \]

Proof. Apply Jacobi’s formula with \( \mathbf{M}=\mathbf{\Phi} \): \( \det(\mathbf{\Phi})'=\det(\mathbf{\Phi})\operatorname{tr}(\mathbf{\Phi}^{-1}\dot{\mathbf{\Phi}}) \). Substitute \( \dot{\mathbf{\Phi}}=\mathbf{A}\mathbf{\Phi} \): \( \mathbf{\Phi}^{-1}\dot{\mathbf{\Phi}}=\mathbf{\Phi}^{-1}\mathbf{A}\mathbf{\Phi} \). Use cyclic invariance of trace: \( \operatorname{tr}(\mathbf{\Phi}^{-1}\mathbf{A}\mathbf{\Phi})=\operatorname{tr}(\mathbf{A}\mathbf{\Phi}\mathbf{\Phi}^{-1})=\operatorname{tr}(\mathbf{A}) \). □

Integrating the scalar ODE for \( w(t)=\det(\mathbf{\Phi}(t)) \) yields

\[ \det(\mathbf{\Phi}(t)) = \det(\mathbf{\Phi}(t_0))\, \exp\!\left(\int_{t_0}^{t}\operatorname{tr}(\mathbf{A}(\sigma))\,d\sigma\right). \]

Therefore, if \( \det(\mathbf{\Phi}(t_0))\neq 0 \), then \( \det(\mathbf{\Phi}(t))\neq 0 \) for all times in the interval, establishing global invertibility of the fundamental matrix on that interval.

Assume \( \mathbf{A}(t)=\mathbf{A} \) is constant. Lesson 3 defined \( e^{\mathbf{A}t}=\sum_{k=0}^{\infty}\frac{(\mathbf{A}t)^k}{k!} \). We now show it solves the matrix IVP for the fundamental matrix.

Theorem. The matrix \( \mathbf{\Phi}(t)=e^{\mathbf{A}(t-t_0)} \) satisfies

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

Proof. Differentiate the series term-by-term (justified by uniform convergence on bounded time intervals for finite matrices):

\[ \frac{d}{dt}e^{\mathbf{A}(t-t_0)} = \frac{d}{dt}\sum_{k=0}^{\infty}\frac{\mathbf{A}^k(t-t_0)^k}{k!} = \sum_{k=1}^{\infty}\frac{\mathbf{A}^k k (t-t_0)^{k-1}}{k!} = \mathbf{A}\sum_{k=1}^{\infty}\frac{\mathbf{A}^{k-1}(t-t_0)^{k-1}}{(k-1)!} = \mathbf{A}e^{\mathbf{A}(t-t_0)}. \]

Also \( e^{\mathbf{A}(t_0-t_0)}=e^{\mathbf{0}}=\mathbf{I} \). □

Consequently, the IVP solution is \( \mathbf{x}(t)=e^{\mathbf{A}(t-t_0)}\mathbf{x}_0 \). This is the most computationally convenient case (via robust expm-type routines).

Computationally, the fundamental matrix is obtained by solving the matrix ODE \( \dot{\mathbf{\Phi}}=\mathbf{A}(t)\mathbf{\Phi} \) with \( \mathbf{\Phi}(t_0)=\mathbf{I} \). For constant \( \mathbf{A} \), use a matrix exponential routine; for time-varying \( \mathbf{A}(t) \), integrate the matrix ODE (often by vectorizing \( \mathbf{\Phi} \) into \( n^2 \) states).

flowchart TD
  A0["Input: A(t), t0, t_grid"] --> Q["Is A(t) constant?"]
  Q -->|yes| E["Compute Phi(t) = expm(A*(t-t0))"]
  Q -->|no| O["Integrate Phi_dot = A(t)*Phi, Phi(t0)=I"]
  E --> X["Solve IVP: x(t)=Phi(t)*inv(Phi(t0))*x0"]
  O --> X
        

In this lesson, we focus strictly on the homogeneous system. (Non-homogeneous forcing terms will be handled later when we study forced responses.)

Libraries. Use SciPy for expm and for ODE integration. For modern control workflows later, the python-control ecosystem is common, but for this lesson we stay with the linear-algebra/ODE core.

import numpy as np
from scipy.linalg import expm
from scipy.integrate import solve_ivp

def fundamental_matrix_LTI(A, t0, t):
    """
    Phi(t) for constant A with normalization Phi(t0) = I:
      Phi(t) = expm(A*(t-t0))
    """
    return expm(A * (t - t0))

def fundamental_matrix_LTV(Afun, t0, t_eval):
    """
    Compute Phi(t) for time-varying A(t) by integrating:
      Phi_dot = A(t) Phi,  Phi(t0)=I
    Vectorize Phi into a length n*n vector to use solve_ivp.
    """
    t_eval = np.asarray(t_eval)
    n = Afun(t0).shape[0]
    Phi0 = np.eye(n).reshape(-1)

    def ode(t, phi_vec):
        Phi = phi_vec.reshape(n, n)
        dPhi = Afun(t) @ Phi
        return dPhi.reshape(-1)

    sol = solve_ivp(ode, (t0, float(t_eval[-1])), Phi0, t_eval=t_eval, rtol=1e-9, atol=1e-12)
    Phi_list = [sol.y[:, k].reshape(n, n) for k in range(sol.y.shape[1])]
    return Phi_list

# Example 1 (LTI): x_dot = A x
A = np.array([[0.0, 1.0],
              [-2.0, -3.0]])
t0 = 0.0
x0 = np.array([1.0, 0.0])

t = 1.25
Phi = fundamental_matrix_LTI(A, t0, t)
x_t = Phi @ x0
print("Phi(t):\n", Phi)
print("x(t): ", x_t)

# Example 2 (LTV): A(t) = [[0, 1], [-(2+0.1 t), -3]]
def Afun(t):
    return np.array([[0.0, 1.0],
                     [-(2.0 + 0.1*t), -3.0]])

t_grid = np.linspace(0.0, 2.0, 21)
Phi_grid = fundamental_matrix_LTV(Afun, 0.0, t_grid)

# Solve IVP x(t) = Phi(t) x0 when Phi(0)=I
x_grid = np.array([Phi_grid[k] @ x0 for k in range(len(t_grid))])
print("x(2.0) approx:", x_grid[-1])

The LTV routine computes \( \mathbf{\Phi}(t) \) directly; since we enforce \( \mathbf{\Phi}(t_0)=\mathbf{I} \), the IVP solution reduces to \( \mathbf{x}(t)=\mathbf{\Phi}(t)\mathbf{x}_0 \).

Libraries. Eigen is a standard C++ linear algebra library. For matrix exponential, Eigen provides unsupported/Eigen/MatrixFunctions. For LTV systems, we integrate \( \dot{\mathbf{\Phi}}=\mathbf{A}(t)\mathbf{\Phi} \) using a basic RK4 integrator (from scratch).

#include <iostream>
#include <vector>
#include <Eigen/Dense>
#include <unsupported/Eigen/MatrixFunctions>

using Eigen::MatrixXd;
using Eigen::VectorXd;

MatrixXd fundamentalMatrixLTI(const MatrixXd& A, double t0, double t) {
  return (A * (t - t0)).exp(); // Phi(t) = expm(A*(t-t0))
}

// Example time-varying A(t)
MatrixXd Afun(double t) {
  MatrixXd A(2,2);
  A << 0.0, 1.0,
       -(2.0 + 0.1*t), -3.0;
  return A;
}

// One RK4 step for Phi_dot = A(t)*Phi
MatrixXd rk4StepPhi(double t, double h, const MatrixXd& Phi) {
  MatrixXd k1 = Afun(t) * Phi;
  MatrixXd k2 = Afun(t + 0.5*h) * (Phi + 0.5*h*k1);
  MatrixXd k3 = Afun(t + 0.5*h) * (Phi + 0.5*h*k2);
  MatrixXd k4 = Afun(t + h) * (Phi + h*k3);
  return Phi + (h/6.0) * (k1 + 2.0*k2 + 2.0*k3 + k4);
}

std::vector<MatrixXd> fundamentalMatrixLTV(double t0, const std::vector<double>& tGrid) {
  int n = 2;
  MatrixXd Phi = MatrixXd::Identity(n,n); // Phi(t0)=I
  std::vector<MatrixXd> out;
  out.reserve(tGrid.size());
  out.push_back(Phi);

  double t = t0;
  for (size_t k = 1; k < tGrid.size(); ++k) {
    double tNext = tGrid[k];
    int steps = 50; // simple fixed sub-stepping for stability/accuracy
    double h = (tNext - t) / static_cast<double>(steps);
    for (int s = 0; s < steps; ++s) {
      Phi = rk4StepPhi(t + s*h, h, Phi);
    }
    t = tNext;
    out.push_back(Phi);
  }
  return out;
}

int main() {
  // LTI example
  MatrixXd A(2,2);
  A << 0.0, 1.0,
       -2.0, -3.0;

  double t0 = 0.0;
  double t = 1.25;
  MatrixXd PhiLTI = fundamentalMatrixLTI(A, t0, t);

  VectorXd x0(2);
  x0 << 1.0, 0.0;
  VectorXd x_t = PhiLTI * x0;

  std::cout << "Phi_LTI(t):\n" << PhiLTI << "\n";
  std::cout << "x(t):\n" << x_t << "\n";

  // LTV example
  std::vector<double> tGrid;
  for (int k = 0; k <= 20; ++k) tGrid.push_back(0.1 * k);

  auto PhiGrid = fundamentalMatrixLTV(0.0, tGrid);
  VectorXd x_end = PhiGrid.back() * x0;
  std::cout << "x(2.0) approx:\n" << x_end << "\n";
  return 0;
}

The RK4 code is intentionally minimal for transparency. In production numerical work, prefer adaptive solvers and error control; however, the matrix ODE structure is exactly as derived in this lesson.

Libraries. Apache Commons Math provides ODE integrators and linear algebra utilities. We integrate the vectorized matrix equation by stacking columns of \( \mathbf{\Phi} \) into one vector of size \( n^2 \).

import java.util.Arrays;
import org.apache.commons.math3.ode.FirstOrderDifferentialEquations;
import org.apache.commons.math3.ode.nonstiff.DormandPrince853Integrator;
import org.apache.commons.math3.ode.sampling.StepHandler;
import org.apache.commons.math3.ode.sampling.StepInterpolator;

public class FundamentalMatrixLTV {

  // Example A(t) for n=2
  static double[][] A(double t) {
    return new double[][]{
      {0.0, 1.0},
      {-(2.0 + 0.1*t), -3.0}
    };
  }

  // Multiply 2x2 by 2x2
  static double[][] matMul(double[][] M, double[][] N) {
    double[][] R = new double[2][2];
    for (int i = 0; i < 2; i++) {
      for (int j = 0; j < 2; j++) {
        R[i][j] = 0.0;
        for (int k = 0; k < 2; k++) R[i][j] += M[i][k] * N[k][j];
      }
    }
    return R;
  }

  // Stack columns of Phi into a length-4 vector: [Phi11, Phi21, Phi12, Phi22]
  static double[] vecPhi(double[][] Phi) {
    return new double[]{Phi[0][0], Phi[1][0], Phi[0][1], Phi[1][1]};
  }

  static double[][] unvecPhi(double[] v) {
    return new double[][]{
      {v[0], v[2]},
      {v[1], v[3]}
    };
  }

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

    // Phi(t0)=I
    double[] y0 = new double[]{1.0, 0.0, 0.0, 1.0};

    FirstOrderDifferentialEquations ode = new FirstOrderDifferentialEquations() {
      public int getDimension() { return 4; }

      public void computeDerivatives(double t, double[] y, double[] yDot) {
        double[][] Phi = unvecPhi(y);
        double[][] dPhi = matMul(A(t), Phi);
        double[] v = vecPhi(dPhi);
        System.arraycopy(v, 0, yDot, 0, 4);
      }
    };

    DormandPrince853Integrator integrator =
        new DormandPrince853Integrator(1.0e-6, 0.1, 1.0e-10, 1.0e-12);

    integrator.addStepHandler(new StepHandler() {
      public void init(double t0, double[] y0, double t) {}

      public void handleStep(StepInterpolator interpolator, boolean isLast) {
        double t = interpolator.getCurrentTime();
        double[] y = interpolator.getInterpolatedState();
        if (isLast) {
          double[][] Phi = unvecPhi(y);
          System.out.println("Phi(tf) = " + Arrays.deepToString(Phi));

          // Example IVP x(tf)=Phi(tf)*x0 (since Phi(t0)=I)
          double[] x0 = new double[]{1.0, 0.0};
          double[] x = new double[]{
            Phi[0][0]*x0[0] + Phi[0][1]*x0[1],
            Phi[1][0]*x0[0] + Phi[1][1]*x0[1]
          };
          System.out.println("x(tf) approx = " + Arrays.toString(x));
        }
      }
    });

    double[] y = y0.clone();
    integrator.integrate(ode, t0, y, tf, y);
  }
}

For constant \( \mathbf{A} \), you would typically use a dedicated matrix exponential routine (often via a numerical library). The ODE approach above works uniformly for both LTI and LTV, at the cost of more computation.

MATLAB provides robust matrix exponentials via expm. For time-varying systems, integrate the matrix ODE for \( \mathbf{\Phi}(t) \) using ode45 by vectorizing the matrix state.

% LTI example: Phi(t) = expm(A*(t-t0))
A = [0 1; -2 -3];
t0 = 0;
t  = 1.25;
Phi = expm(A*(t-t0));

x0 = [1; 0];
x_t = Phi*x0;

disp('Phi(t) ='); disp(Phi);
disp('x(t) ='); disp(x_t);

% LTV example: integrate Phi_dot = A(t)*Phi, Phi(t0)=I using ode45
Afun = @(tt) [0 1; -(2 + 0.1*tt) -3];
n = 2;

phi0 = reshape(eye(n), [], 1); % vec(Phi(t0))

ode = @(tt, phi_vec) reshape( Afun(tt) * reshape(phi_vec, n, n), [], 1 );

tspan = [0 2];
opts = odeset('RelTol',1e-9,'AbsTol',1e-12);
[tt, phi_sol] = ode45(ode, tspan, phi0, opts);

Phi_tf = reshape(phi_sol(end,:).', n, n);
x_tf = Phi_tf * x0;

disp('Phi(2) approx ='); disp(Phi_tf);
disp('x(2) approx ='); disp(x_tf);

Simulink approach (conceptual, no images). A convenient way to approximate a fundamental matrix in Simulink is to simulate the same state-space model multiple times with different initial conditions: \( \mathbf{x}(t_0)=\mathbf{e}_i \), where \( \mathbf{e}_i \) is the i-th standard basis vector. The resulting trajectories become the columns of \( \mathbf{\Phi}(t) \).

Suggested block setup: Use a State-Space block (continuous) with matrix \( \mathbf{A} \), and set \( \mathbf{B}=\mathbf{0} \), \( \mathbf{u}=0 \). Then:

• Run \( n \) simulations with initial condition \( \mathbf{x}(t_0)=\mathbf{e}_1,\dots,\mathbf{e}_n \).
• Log output \( \mathbf{x}(t) \) each run.
• Assemble \( \mathbf{\Phi}(t)=[\mathbf{x}^{(1)}(t)\ \cdots\ \mathbf{x}^{(n)}(t)] \) offline in MATLAB.

% Offline assembly idea after n runs (pseudo-structure)
% Suppose logs are stored as x1(t), x2(t), ..., xn(t) each size n-by-N
% Then Phi(t_k) = [x1(:,k) x2(:,k) ... xn(:,k)] for each sample k.

Mathematica supports symbolic/numeric matrix exponentials and can solve matrix ODEs directly.

(* LTI example *)
A = { {0, 1}, {-2, -3} };
t0 = 0;
t  = 5/4;

Phi[t_] := MatrixExp[A (t - t0)];
x0 = {1, 0};
x[t_] := Phi[t].x0;

PhiVal = Phi[t];
xVal = x[t];

Print["Phi(t) = ", PhiVal];
Print["x(t) = ", xVal];

(* LTV example: Phi'(t)=A(t) Phi(t), Phi(t0)=I *)
Afun[t_] := { {0, 1}, {-(2 + 0.1 t), -3} };
PhiSym = Array[Subscript[\[CapitalPhi], ##] &, {2, 2}];

eqs = Flatten@Table[
   D[PhiSym[[i, j]], t] == Sum[Afun[t][[i, k]] PhiSym[[k, j]], {k, 1, 2}],
   {i, 1, 2}, {j, 1, 2}
];

ics = Flatten@Table[
   PhiSym[[i, j]] /. t -> 0 == If[i == j, 1, 0],
   {i, 1, 2}, {j, 1, 2}
];

sol = NDSolve[Join[eqs, ics], PhiSym, {t, 0, 2}][[1]];

PhiNum[tt_] := Evaluate[PhiSym /. sol /. t -> tt];
xNum[tt_] := PhiNum[tt].x0;

Print["Phi(2) approx = ", PhiNum[2]];
Print["x(2) approx = ", xNum[2]];

Problem 1 (Fundamental Matrix Representation): Let \( \mathbf{\Phi}(t) \) be a fundamental matrix of \( \dot{\mathbf{x}}=\mathbf{A}(t)\mathbf{x} \). Show that if \( \mathbf{x}(t)=\mathbf{\Phi}(t)\mathbf{c}(t) \) for a differentiable vector \( \mathbf{c}(t) \), then \( \dot{\mathbf{c}}(t)=\mathbf{0} \).

Solution: Differentiate and substitute:

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

But \( \dot{\mathbf{x}}(t)=\mathbf{A}(t)\mathbf{x}(t)=\mathbf{A}(t)\mathbf{\Phi}(t)\mathbf{c}(t) \). Subtracting yields \( \mathbf{\Phi}(t)\dot{\mathbf{c}}(t)=\mathbf{0} \). Since \( \mathbf{\Phi}(t) \) is invertible, \( \dot{\mathbf{c}}(t)=\mathbf{0} \). Hence \( \mathbf{c}(t) \) is constant. □

Problem 2 (Relation Between Two Fundamental Matrices): Suppose \( \mathbf{\Phi}_1(t) \) and \( \mathbf{\Phi}_2(t) \) are fundamental matrices for the same system. Prove that \( \mathbf{\Phi}_2(t)=\mathbf{\Phi}_1(t)\mathbf{C} \) for a constant invertible matrix \( \mathbf{C} \).

Solution: Define \( \mathbf{C}(t)=\mathbf{\Phi}_1(t)^{-1}\mathbf{\Phi}_2(t) \). As shown in Section 4, \( \dot{\mathbf{C}}(t)=\mathbf{0} \), so \( \mathbf{C}(t)\equiv \mathbf{C} \) constant. Invertibility follows because both fundamental matrices are invertible. □

Problem 3 (Compute a Fundamental Matrix for an LTI System): Let \( \mathbf{A}=\begin{bmatrix}0 & 1\\ -2 & -3\end{bmatrix} \). Find a normalized fundamental matrix \( \mathbf{\Phi}(t) \) satisfying \( \mathbf{\Phi}(0)=\mathbf{I} \).

Solution: For constant \( \mathbf{A} \), take \( \mathbf{\Phi}(t)=e^{\mathbf{A}t} \). Compute eigenvalues of \( \mathbf{A} \): characteristic polynomial \( \lambda^2+3\lambda+2=0 \) gives \( \lambda_1=-1,\ \lambda_2=-2 \).

With distinct eigenvalues, \( \mathbf{A} \) is diagonalizable: \( \mathbf{A}=\mathbf{V}\mathbf{\Lambda}\mathbf{V}^{-1} \) with \( \mathbf{\Lambda}=\operatorname{diag}(-1,-2) \). Then

\[ e^{\mathbf{A}t} = \mathbf{V}e^{\mathbf{\Lambda}t}\mathbf{V}^{-1} = \mathbf{V}\operatorname{diag}(e^{-t},e^{-2t})\mathbf{V}^{-1}. \]

Any correct \( \mathbf{V} \) from eigenvectors yields a valid closed form. In practice, the most reliable method is to compute \( e^{\mathbf{A}t} \) numerically via expm (MATLAB/Python) or library routines (C++/Java), which exactly implements \( \mathbf{\Phi}(t) \). □

Problem 4 (Liouville Determinant Formula): Assume \( \mathbf{\Phi}(t) \) satisfies \( \dot{\mathbf{\Phi}}=\mathbf{A}(t)\mathbf{\Phi} \). Show that \( \det(\mathbf{\Phi}(t))=\det(\mathbf{\Phi}(t_0))\exp\!\left(\int_{t_0}^{t}\operatorname{tr}(\mathbf{A}(\sigma))\,d\sigma\right) \).

Solution: By Section 5, \( \frac{d}{dt}\det(\mathbf{\Phi}(t))=\operatorname{tr}(\mathbf{A}(t))\det(\mathbf{\Phi}(t)) \). This is a scalar first-order linear ODE in \( w(t)=\det(\mathbf{\Phi}(t)) \):

\[ \dot{w}(t)=\operatorname{tr}(\mathbf{A}(t))w(t), \qquad w(t_0)=\det(\mathbf{\Phi}(t_0)). \]

Separation and integration give \( \ln w(t)-\ln w(t_0)=\int_{t_0}^{t}\operatorname{tr}(\mathbf{A}(\sigma))\,d\sigma \), hence the claimed formula. □

Problem 5 (Diagonal Time-Varying System): Let \( \mathbf{A}(t)=\operatorname{diag}(a_1(t),\dots,a_n(t)) \). Find a normalized fundamental matrix \( \mathbf{\Phi}(t) \) with \( \mathbf{\Phi}(t_0)=\mathbf{I} \).

Solution: The system decouples into scalar ODEs: \( \dot{x}_i(t)=a_i(t)x_i(t) \). The normalized scalar fundamental solution for component i is

\[ \phi_i(t) = \exp\!\left(\int_{t_0}^{t} a_i(\sigma)\,d\sigma\right), \qquad \phi_i(t_0)=1. \]

Therefore,

\[ \mathbf{\Phi}(t) = \operatorname{diag}\!\left( \exp\!\left(\int_{t_0}^{t} a_1(\sigma)\,d\sigma\right), \dots, \exp\!\left(\int_{t_0}^{t} a_n(\sigma)\,d\sigma\right) \right). \]

This also matches Liouville’s formula since \( \operatorname{tr}(\mathbf{A}(t))=\sum_{i=1}^n a_i(t) \) and \( \det(\mathbf{\Phi}(t))=\prod_{i=1}^n \phi_i(t) \). □

We introduced the fundamental matrix \( \mathbf{\Phi}(t) \) as a matrix whose columns form a basis of solutions to \( \dot{\mathbf{x}}=\mathbf{A}(t)\mathbf{x} \). We proved the representation \( \mathbf{x}(t)=\mathbf{\Phi}(t)\mathbf{\Phi}(t_0)^{-1}\mathbf{x}_0 \) for IVPs, established that any two fundamental matrices differ by a constant invertible factor, and derived Liouville’s formula governing \( \det(\mathbf{\Phi}(t)) \). Finally, we implemented fundamental-matrix computation in major technical languages using matrix exponentials (LTI) and matrix ODE integration (LTV).