Systems of Differential Equations

Coupled linear differential equation systems, state-space representations, and dynamic modal solutions.

What it is

A first-order system of ordinary differential equations is a coupled set of $n$ dynamic equations governing $n$ state variables $\mathbf{x}(t) = \left(x_1(t), x_2(t), \dots, x_n(t)\right)^T$: $\frac{d\mathbf{x}}{dt} = \mathbf{f}(t, \mathbf{x})$ Any $n$-th order scalar ODE $y^{(n)} = F(t, y, y', \dots, y^{(n-1)})$ canonically reduces to an $n$-dimensional first-order system by defining state variables $x_1 = y, x_2 = y', \dots, x_n = y^{(n-1)}$.

For homogeneous linear systems with constant coefficients: $\frac{d\mathbf{x}}{dt} = \mathbf{A} \mathbf{x}, \quad \mathbf{A} \in \mathbb{R}^{n \times n}$ The general solution is built from the eigenvalues $\lambda_i$ and eigenvectors $\mathbf{v}_i$ satisfying $\mathbf{A} \mathbf{v}_i = \lambda_i \mathbf{v}_i$:

• When $\mathbf{A}$ possesses $n$ linearly independent eigenvectors: $\mathbf{x}(t) = \sum_{i=1}^n c_i e^{\lambda_i t} \mathbf{v}_i$
• Subject to Cauchy initial data $\mathbf{x}(0) = \mathbf{x}_0$, the unique trajectory is generated via the matrix exponential: $\mathbf{x}(t) = e^{\mathbf{A} t} \mathbf{x}_0, \quad e^{\mathbf{A} t} = \sum_{k=0}^\infty \frac{(\mathbf{A} t)^k}{k!} = \mathbf{P} e^{\mathbf{D} t} \mathbf{P}^{-1}$ where $\mathbf{P}$ is the modal transformation matrix whose columns are eigenvectors, and $\mathbf{D} = \operatorname{diag}(\lambda_1, \dots, \lambda_n)$.

Why it appears in physics

Nature operates through coupled multi-degree-of-freedom interactions:

• Coupled harmonic oscillators: Masses joined by elastic springs exchange vibrational energy ($m_1 \ddot{x}_1 = -k_1 x_1 + k_{12}(x_2 - x_1)$). Diagonalization decouples the system into independent normal modes, each oscillating at a pure eigenfrequency.
• Multiloop electrical networks: Multi-mesh circuits governed by Kirchhoff's current and voltage laws.
• Bateman radioactive decay cascades: Sequential decay series ($A \to B \to C$), where daughter populations track parent decay rates.
• Hamiltonian mechanics: Classical physics expresses Newton's law as a canonical first-order symplectic system of $2n$ equations: $\dot{q}_i = \frac{\partial H}{\partial p_i}, \quad \dot{p}_i = -\frac{\partial H}{\partial q_i}$

Visual or geometric intuition

The state vector $\mathbf{x}(t)$ traces a continuous path in $n$-dimensional state space (phase space). At every coordinate point $\mathbf{x}$, matrix $\mathbf{A}$ assigns a local tangent velocity vector $\mathbf{v} = \mathbf{A}\mathbf{x}$. Eigenvectors define invariant rays: along an eigenspace direction $\mathbf{v}_i$, the vector field points along the line itself. A state initiated on an eigenvector line never leaves it, expanding ($e^{\lambda t}, \lambda > 0$) or decaying ($e^{\lambda t}, \lambda < 0$) rectilinearly. All other trajectories are curved superpositions of these modal directions.

Minimal symbolic language

• Fundamental matrix $\mathbf{\Phi}(t)$: An $n \times n$ matrix whose columns are $n$ linearly independent solution vectors $\mathbf{x}^{(1)}(t), \dots, \mathbf{x}^{(n)}(t)$: $\frac{d\mathbf{\Phi}}{dt} = \mathbf{A} \mathbf{\Phi}(t), \quad W(t) = \det\mathbf{\Phi}(t) \neq 0$
• Liouville-Jacobi formula for matrix Wronskians: $W(t) = W(0) \exp(\operatorname{tr}(\mathbf{A}) t)$ where $\operatorname{tr}(\mathbf{A}) = \sum \lambda_i$ measures phase volume expansion rate (Liouville's theorem).
• Duhamel convolution formula for forced systems: $\mathbf{x}' = \mathbf{A} \mathbf{x} + \mathbf{g}(t) \implies \mathbf{x}(t) = e^{\mathbf{A} t} \mathbf{x}_0 + \int_0^t e^{\mathbf{A}(t - s)} \mathbf{g}(s)\,ds$

What the result means

1. Modal decoupling: Transforming to modal coordinates $\mathbf{y} = \mathbf{P}^{-1} \mathbf{x}$ splits the coupled dynamics into $n$ uncoupled scalar ODEs $\dot{y}_i = \lambda_i y_i$.
2. Volume preservation and dissipation: If $\operatorname{tr}(\mathbf{A}) = 0$ (Hamiltonian conservative systems), the phase flow strictly preserves volume. If $\operatorname{tr}(\mathbf{A}) < 0$, phase volume contracts toward attractors.
3. Complex eigenvalues: A pair $\lambda = \alpha \pm i\omega$ generates rotational rotation with frequency $\omega$ modulated by radial envelope $e^{\alpha t}$.

When not to use it

• Defective matrices (non-diagonalizable): When algebraic multiplicity exceeds geometric multiplicity, generalized eigenvectors and Jordan forms are mandatory.
• Far-from-equilibrium non-linear systems: Local Jacobian linearization $\mathbf{x}' = \mathbf{J} \mathbf{x}$ fails beyond small perturbations, missing limit cycles and chaos.
• Time-varying non-commuting matrices $\mathbf{A}(t)$: If $\mathbf{A}(t) \mathbf{A}(s) \neq \mathbf{A}(s) \mathbf{A}(t)$, the exponential solution fails; Magnus expansions or Dyson series are required.

Operational summary

1. Express the system in canonical form $\mathbf{x}' = \mathbf{A} \mathbf{x}$.
2. Solve $\det(\mathbf{A} - \lambda \mathbf{I}) = 0$ for eigenvalues $\lambda_i$.
3. Compute eigenvectors from $(\mathbf{A} - \lambda_i \mathbf{I}) \mathbf{v}_i = \mathbf{0}$.
4. For complex roots $\lambda = \alpha \pm i\omega$, extract real and imaginary components of $e^{(\alpha + i\omega) t} \mathbf{v}$.
5. Assemble the fundamental matrix $\mathbf{\Phi}(t)$ and apply initial data $\mathbf{x}_0$ via $\mathbf{\Phi}(0) \mathbf{c} = \mathbf{x}_0$.

Procedure: Solving Linear Systems of ODEs

Before calculating

1. Define the state vector: Explicitly write $\mathbf{x}(t) = (x_1, x_2, \dots)^T$ with physical dimensions.
2. Build matrix $\mathbf{A}$: Check that equations are arranged consistently.
3. Compute trace and determinant: $\operatorname{tr}(\mathbf{A}) = a_{11} + a_{22}$ and $\det(\mathbf{A}) = a_{11} a_{22} - a_{12} a_{21}$. Characteristic polynomial: $\lambda^2 - \operatorname{tr}(\mathbf{A}) \lambda + \det(\mathbf{A}) = 0$.

Step 1: Identify the mathematical object

Write the system matrix: $\mathbf{A} = \begin{pmatrix} a_{11} & a_{12} \\ a_{21} & a_{22} \end{pmatrix}$ Confirm whether the system is homogeneous.

Step 2: Choose the operation

1. Solve the characteristic equation: $\det(\mathbf{A} - \lambda \mathbf{I}) = 0 \implies \lambda^2 - \operatorname{tr}(\mathbf{A})\lambda + \det(\mathbf{A}) = 0$
2. Determine eigenvectors for each root: $(\mathbf{A} - \lambda_k \mathbf{I}) \mathbf{v}_k = \mathbf{0}$
3. Assemble based on eigenvalue category:
   • Distinct real roots: $\mathbf{x}(t) = c_1 e^{\lambda_1 t} \mathbf{v}_1 + c_2 e^{\lambda_2 t} \mathbf{v}_2$
   • Complex conjugate pair ($\lambda = \alpha \pm i\beta$, eigenvector $\mathbf{v} = \mathbf{u} + i\mathbf{w}$): $\mathbf{x}_1(t) = e^{\alpha t}\left(\mathbf{u} \cos(\beta t) - \mathbf{w} \sin(\beta t)\right)$ $\mathbf{x}_2(t) = e^{\alpha t}\left(\mathbf{u} \sin(\beta t) + \mathbf{w} \cos(\beta t)\right)$ $\mathbf{x}(t) = c_1 \mathbf{x}_1(t) + c_2 \mathbf{x}_2(t)$
   • Repeated defective root: Find generalized eigenvector $(\mathbf{A} - \lambda \mathbf{I}) \mathbf{w} = \mathbf{v}$: $\mathbf{x}(t) = c_1 e^{\lambda t} \mathbf{v} + c_2 e^{\lambda t} (t \mathbf{v} + \mathbf{w})$

Step 3: Conserve units and signs

• Eigenvalues $\lambda_i$ carry dimensions of inverse time $\text{s}^{-1}$.
• Eigenvectors fix relative modal scaling.

Step 4: Calculate

1. Build the fundamental matrix $\mathbf{\Phi}(t) = [\mathbf{x}_1(t) \quad \mathbf{x}_2(t)]$.
2. Evaluate at $t = 0$: $\mathbf{\Phi}(0) \mathbf{c} = \mathbf{x}_0$.
3. Solve for constants $\mathbf{c} = \mathbf{\Phi}(0)^{-1} \mathbf{x}_0$.

Step 5: Interpret the result

1. $\operatorname{Re}(\lambda_i) < 0 \implies$ asymptotically stable sink (attractor).
2. Pure imaginary $\alpha = 0 \implies$ closed elliptic orbits (conservative center).
3. Opposite real signs $\implies$ saddle point (unstable).

Self-check checklist

• Does $\lambda_1 + \lambda_2 = \operatorname{tr}(\mathbf{A})$?
• Does $\lambda_1 \lambda_2 = \det(\mathbf{A})$?
• Did I verify $(\mathbf{A} - \lambda \mathbf{I}) \mathbf{v} = \mathbf{0}$?
• Did I eliminate all imaginary units $i$ from the physical solution vector?
• Does $\mathbf{x}(0)$ match initial condition vector $\mathbf{x}_0$?

Worked Examples: Systems of Differential Equations

Step-by-step guided example

Problem

Two interconnected brine tanks contain $V_1 = 100\,\text{L}$ and $V_2 = 100\,\text{L}$ respectively. Fluid is pumped from Tank 1 to Tank 2 at rate $r = 4\,\text{L}/\text{min}$, and from Tank 2 to Tank 1 at identical rate $r = 4\,\text{L}/\text{min}$, maintaining constant liquid volumes. At $t = 0$, Tank 1 holds $x_1(0) = 20\,\text{kg}$ of dissolved salt while Tank 2 holds pure fresh water $x_2(0) = 0\,\text{kg}$.

1. Formulate the matrix system of ODEs for salt masses $x_1(t)$ and $x_2(t)$.
2. Calculate eigenvalues and eigenvectors.
3. Solve analytically for $x_1(t)$ and $x_2(t)$.
4. Determine the asymptotic equilibrium state as $t \to \infty$.

Detailed solution

Step 1: Matrix system formulation Salt rate of change in each tank equals inflow rate minus outflow rate: $\frac{dx_1}{dt} = -\frac{4}{100} x_1 + \frac{4}{100} x_2 = -0.04 x_1 + 0.04 x_2$ $\frac{dx_2}{dt} = \frac{4}{100} x_1 - \frac{4}{100} x_2 = 0.04 x_1 - 0.04 x_2$ Matrix format: $\frac{d}{dt} \begin{pmatrix} x_1 \\ x_2 \end{pmatrix} = \begin{pmatrix} -0.04 & 0.04 \\ 0.04 & -0.04 \end{pmatrix} \begin{pmatrix} x_1 \\ x_2 \end{pmatrix}$

Step 2: Eigenvalue and eigenvector computation Characteristic polynomial: $\det(\mathbf{A} - \lambda \mathbf{I}) = \lambda^2 + 0.08 \lambda = 0 \implies \lambda_1 = 0, \quad \lambda_2 = -0.08\,\text{min}^{-1}$

• For $\lambda_1 = 0$: $\begin{pmatrix} -0.04 & 0.04 \\ 0.04 & -0.04 \end{pmatrix} \begin{pmatrix} v_{11} \\ v_{12} \end{pmatrix} = \mathbf{0} \implies v_{11} = v_{12} \implies \mathbf{v}_1 = \begin{pmatrix} 1 \\ 1 \end{pmatrix}$
• For $\lambda_2 = -0.08$: $\begin{pmatrix} 0.04 & 0.04 \\ 0.04 & 0.04 \end{pmatrix} \begin{pmatrix} v_{21} \\ v_{22} \end{pmatrix} = \mathbf{0} \implies v_{21} = -v_{22} \implies \mathbf{v}_2 = \begin{pmatrix} 1 \\ -1 \end{pmatrix}$

Step 3: General solution and initial condition fitting $\mathbf{x}(t) = c_1 \begin{pmatrix} 1 \\ 1 \end{pmatrix} + c_2 e^{-0.08 t} \begin{pmatrix} 1 \\ -1 \end{pmatrix}$ At $t = 0$: $\mathbf{x}(0) = \begin{pmatrix} 20 \\ 0 \end{pmatrix}$: $\begin{cases} c_1 + c_2 = 20 \\ c_1 - c_2 = 0 \end{cases} \implies c_1 = 10, \quad c_2 = 10$ Yielding: $x_1(t) = 10 + 10 e^{-0.08 t}\,\text{kg}$ $x_2(t) = 10 - 10 e^{-0.08 t}\,\text{kg}$

Step 4: Asymptotic equilibrium

• Total dissolved salt: $x_1(t) + x_2(t) = 20\,\text{kg}$ (conserved).
• As $t \to \infty$, the transient mode vanishes: $x_1(\infty) = 10\,\text{kg}, \quad x_2(\infty) = 10\,\text{kg}$ Both reach a uniform concentration of $0.1\,\text{kg}/\text{L}$.

Inverse example

Problem

Two neural signals recorded by microelectrodes evolve as: $x_1(t) = 2 e^{-t} \cos(3 t) - e^{-t} \sin(3 t)$ $x_2(t) = e^{-t} \cos(3 t) + 2 e^{-t} \sin(3 t)$ Reconstruct the governing matrix $\mathbf{A} \in \mathbb{R}^{2 \times 2}$ such that $\dot{\mathbf{x}} = \mathbf{A} \mathbf{x}$.

Detailed solution

Step 1: Differentiating state variables $\dot{x}_1 = e^{-t}(-5 \cos 3t - 5 \sin 3t)$ $\dot{x}_2 = e^{-t}(5 \cos 3t - 5 \sin 3t)$

Step 2: Sampling at $t = 0$ $\mathbf{x}(0) = \begin{pmatrix} 2 \\ 1 \end{pmatrix}, \quad \dot{\mathbf{x}}(0) = \begin{pmatrix} -5 \\ 5 \end{pmatrix} \implies \mathbf{A} \begin{pmatrix} 2 \\ 1 \end{pmatrix} = \begin{pmatrix} -5 \\ 5 \end{pmatrix}$

Step 3: Sampling at $t = \pi/6$ $\mathbf{x}(\pi/6) = e^{-\pi/6} \begin{pmatrix} -1 \\ 2 \end{pmatrix}, \quad \dot{\mathbf{x}}(\pi/6) = e^{-\pi/6} \begin{pmatrix} -5 \\ -5 \end{pmatrix} \implies \mathbf{A} \begin{pmatrix} -1 \\ 2 \end{pmatrix} = \begin{pmatrix} -5 \\ -5 \end{pmatrix}$

Step 4: Matrix inversion $\mathbf{A} \begin{pmatrix} 2 & -1 \\ 1 & 2 \end{pmatrix} = \begin{pmatrix} -5 & -5 \\ 5 & -5 \end{pmatrix}$ $\mathbf{A} = \begin{pmatrix} -5 & -5 \\ 5 & -5 \end{pmatrix} \frac{1}{5} \begin{pmatrix} 2 & 1 \\ -1 & 2 \end{pmatrix} = \begin{pmatrix} -1 & -3 \\ 3 & -1 \end{pmatrix}$ Eigenvalues check: $\lambda^2 - \operatorname{tr}(\mathbf{A}) \lambda + \det(\mathbf{A}) = \lambda^2 + 2 \lambda + 10 = 0 \implies \lambda = -1 \pm 3 i$ Matches observed damping $\alpha = -1$ and oscillation $\beta = 3$.

Common Errors: Systems of Differential Equations

Error 1: Mismatching eigenvalues with foreign eigenvectors

In $\mathbf{x}(t) = c_1 e^{\lambda_1 t} \mathbf{v}_1 + c_2 e^{\lambda_2 t} \mathbf{v}_2$, multiplying $e^{\lambda_1 t}$ by $\mathbf{v}_2$.

• Correction: Each temporal modal carrier $e^{\lambda_k t}$ is coupled strictly to its own directional eigenvector $\mathbf{v}_k$.

Error 2: Leaving imaginary $i$ in real physical state solutions

Expressing complex-conjugate eigensolutions as formal complex linear combinations without extracting the real and imaginary components.

• Correction: Physical observables require real phase-space coordinates; decompose via Euler's identity into damped real sines and cosines.

Error 3: Assuming defective repeated roots have full eigenvector bases

Assuming any repeated eigenvalue admits two independent eigenvectors without testing the geometric nullspace dimension.

• Consequence: Defective matrices produce a missing basis vector, demanding generalized eigenvector chains $(t \mathbf{v} + \mathbf{w})$.

Error 4: Computing matrix exponentials element-wise

Erroneously writing $e^{\mathbf{A} t} = [e^{a_{ij} t}]$.

• Correction: Matrix exponentials require power series or modal similarity transformations $\mathbf{P} e^{\mathbf{D} t} \mathbf{P}^{-1}$.

Error 5: Dropping cross-coupling flux terms

In conservation balance models, omitting reciprocal exchange rates between nodes.

Quick control rule

Always verify the trace and determinant identities: $\lambda_1 + \lambda_2 = \operatorname{tr}(\mathbf{A})$ and $\lambda_1 \lambda_2 = \det(\mathbf{A})$. If violated, an algebraic error occurred.

Physical Applications: Systems of Differential Equations

1. Molecular normal modes and vibrational spectroscopy

In polyatomic molecules (such as $\text{H}_2\text{O}$ or $\text{CO}_2$), nuclear displacements satisfy: $\mathbf{M} \ddot{\mathbf{x}} + \mathbf{K} \mathbf{x} = \mathbf{0}$ Diagonalizing mass tensor $\mathbf{M}$ alongside stiffness matrix $\mathbf{K}$ yields molecular normal modes (symmetric/asymmetric stretching and bending). These eigenfrequencies match absorption lines in infrared (FTIR) and Raman spectroscopy, enabling chemical identification in astrophysics and pharmacology.

2. Power grid transient stability and smart grids

Transient angular stability across high-voltage transmission networks couples synchronous generators through swing equations: $M_i \ddot{\delta}_i + D_i \dot{\delta}_i = P_{m,i} - \sum_{j=1}^N V_i V_j Y_{ij} \sin(\delta_i - \delta_j)$ Linearizing around steady operation reveals system eigenvalues, preventing cascading blackouts under fault conditions.

3. Ecological population dynamics (Lotka-Volterra)

Two-species competition and predator-prey dynamics obey: $\frac{dN_1}{dt} = N_1 (r_1 - a_{11} N_1 - a_{12} N_2)$ $\frac{dN_2}{dt} = N_2 (r_2 - a_{21} N_1 - a_{22} N_2)$ Matrix linearization at coexistence fixed points establishes competitive exclusion vs. stable ecological equilibrium.

4. Multi-stream heat exchangers

Counterflow shell-and-tube heat exchangers track hot ($T_h$) and cold ($T_c$) stream profiles via: $\dot{m}_h c_{p,h} \frac{dT_h}{dx} = -U P (T_h - T_c), \quad -\dot{m}_c c_{p,c} \frac{dT_c}{dx} = U P (T_h - T_c)$ Solving this linear system derives the Logarithmic Mean Temperature Difference (LMTD) critical for boiler and nuclear steam generator engineering.

5. Orbital rendezvous and Clohessy-Wiltshire relative dynamics

Relative proximity operations between two spacecraft approaching the International Space Station are governed by the Clohessy-Wiltshire (Hill's) equations: $\ddot{x} - 2 \omega \dot{y} - 3 \omega^2 x = F_x / m$ $\ddot{y} + 2 \omega \dot{x} = F_y / m$ $\ddot{z} + \omega^2 z = F_z / m$ Formulating this as a 6D first-order state system enables real-time Kalman filtering and automated docking thruster guidance.

Interdisciplinary connection

Systems of ODEs underpin Dynamic Stochastic General Equilibrium (DSGE) macroeconomic models, biochemical signal cascade networks, and multi-link robotic manipulator inverse kinematics.

Guided Practice: Systems of Differential Equations

Exercise 1: Normal modes of undamped coupled oscillators

Two identical unit masses $m = 1\,\text{kg}$ slide along a frictionless horizontal track. Mass 1 is anchored to a rigid wall by a spring of stiffness $k = 3\,\text{N}/\text{m}$, the two masses are coupled by an intermediate spring $k_c = 2\,\text{N}/\text{m}$, and mass 2 has a free open end.

1. Write Newton's equations and assemble the matrix system for displacement vector $\mathbf{x} = (x_1, x_2)^T$.
2. Determine squared eigenfrequencies $\omega^2$ and normal mode eigenvectors.
3. Interpret the physical motion of each normal mode.

Solution

1. Dynamic equations: $\ddot{x}_1 = -5 x_1 + 2 x_2, \quad \ddot{x}_2 = 2 x_1 - 2 x_2$ Matrix form $\ddot{\mathbf{x}} = \mathbf{K} \mathbf{x}$: $\begin{pmatrix} \ddot{x}_1 \\ \ddot{x}_2 \end{pmatrix} = \begin{pmatrix} -5 & 2 \\ 2 & -2 \end{pmatrix} \begin{pmatrix} x_1 \\ x_2 \end{pmatrix}$
2. Trial $\mathbf{x} = \mathbf{v} e^{i \omega t} \implies \lambda = -\omega^2$: $\det(\mathbf{K} - \lambda \mathbf{I}) = \lambda^2 + 7 \lambda + 6 = 0 \implies \lambda_1 = -1, \quad \lambda_2 = -6$
   • Mode 1: $\omega_1 = 1\,\text{rad}/\text{s}$: $(\mathbf{K} + \mathbf{I})\mathbf{v}_1 = \begin{pmatrix} -4 & 2 \\ 2 & -1 \end{pmatrix} \mathbf{v}_1 = \mathbf{0} \implies \mathbf{v}_1 = \begin{pmatrix} 1 \\ 2 \end{pmatrix}$
   • Mode 2: $\omega_2 = \sqrt{6} \approx 2.45\,\text{rad}/\text{s}$: $(\mathbf{K} + 6\mathbf{I})\mathbf{v}_2 = \begin{pmatrix} 1 & 2 \\ 2 & 4 \end{pmatrix} \mathbf{v}_2 = \mathbf{0} \implies \mathbf{v}_2 = \begin{pmatrix} -2 \\ 1 \end{pmatrix}$
3. Physical description:
   • Mode 1: Masses oscillate in phase in the same direction, with mass 2 moving twice as far as mass 1.
   • Mode 2: Masses move in anti-phase (opposite directions), heavily compressing the coupling spring.

Exercise 2: Bateman radioactive decay chain

Parent nuclide $A$ decays with rate constant $\lambda_1 = 0.5\,\text{h}^{-1}$ into unstable daughter $B$, which decays with $\lambda_2 = 0.1\,\text{h}^{-1}$ into stable isotope $C$: $\dot{N}_A = -\lambda_1 N_A, \quad \dot{N}_B = \lambda_1 N_A - \lambda_2 N_B$ Initial inventory: $N_A(0) = N_0, N_B(0) = 0$.

1. State the lower triangular system matrix $\mathbf{A}$.
2. Solve analytically for daughter inventory $N_B(t)$.
3. Determine peak intermediate daughter activity time.

Solution

1. Matrix equation: $\frac{d}{dt} \begin{pmatrix} N_A \\ N_B \end{pmatrix} = \begin{pmatrix} -0.5 & 0 \\ 0.5 & -0.1 \end{pmatrix} \begin{pmatrix} N_A \\ N_B \end{pmatrix}$
2. Integrating recursively: $N_A(t) = N_0 e^{-0.5 t}$ $\dot{N}_B + 0.1 N_B = 0.5 N_0 e^{-0.5 t}$ Using integrating factor $\mu = e^{0.1 t}$: $N_B(t) = 1.25 N_0 \left(e^{-0.1 t} - e^{-0.5 t}\right)$
3. Peak activity: $\dot{N}_B = 0 \implies -0.1 e^{-0.1 t} + 0.5 e^{-0.5 t} = 0$: $e^{-0.4 t} = 0.2 \implies t_{\text{max}} = \frac{\ln 5}{0.4} \approx 4.024\,\text{hours}$

Exercise 3: Matrix exponential via Taylor expansion

For system $\mathbf{x}' = \mathbf{A} \mathbf{x}$ with $\mathbf{A} = \begin{pmatrix} 0 & 1 \\ -4 & 0 \end{pmatrix}$:

1. Compute $\mathbf{A}^2, \mathbf{A}^3, \mathbf{A}^4$.
2. Form matrix exponential $e^{\mathbf{A} t}$.
3. Match to simple harmonic oscillator motion.

Solution

1. Matrix powers: $\mathbf{A}^2 = -4 \mathbf{I} = -(2)^2 \mathbf{I}, \quad \mathbf{A}^3 = -4 \mathbf{A}, \quad \mathbf{A}^4 = 16 \mathbf{I} = (2)^4 \mathbf{I}$
2. Taylor summation: $e^{\mathbf{A} t} = \cos(2 t) \mathbf{I} + \frac{1}{2} \sin(2 t) \mathbf{A} = \begin{pmatrix} \cos(2 t) & \frac{1}{2} \sin(2 t) \\ -2 \sin(2 t) & \cos(2 t) \end{pmatrix}$
3. Applying $\mathbf{x}(0) = (x_0, v_0)^T$: $x_1(t) = x_0 \cos(2 t) + \frac{v_0}{2} \sin(2 t)$ Matches $\ddot{x} + 4 x = 0$ exactly.