ODE: Systems of First-Order ODEs

Mastering systems of first-order ordinary differential equations is a pivotal skill in engineering, moving you from analyzing single components to modeling entire dynamic systems. Whether you're predicting the vibrations of a vehicle's suspension, tracking heat flow through a composite material, or modeling the population dynamics in an ecosystem, you are dealing with interdependent variables whose rates of change are linked. The matrix-based framework for these systems provides a powerful, unified language for describing, analyzing, and solving such interconnected problems.

From Single Equations to a System of First-Order Equations

A system of first-order linear ODEs is a set of equations where the derivative of each unknown function depends linearly on all the functions themselves. The most compact and powerful way to express such a system is in matrix form. Consider a system with two dependent variables, $x_1(t)$ and $x_2(t)$. Instead of writing: $$x_1^{\prime }=a_{11}{x}_1+a_{12}{x}_2+f_1(t)$$ $$x_2^{\prime }=a_{21}{x}_1+a_{22}{x}_2+f_2(t)$$ you can write it as: $${\mathbf{X}}^{\prime }=A\mathbf{X}+\mathbf{F}(t)$$ where $$\mathbf{X}=\left[\begin{array}{c}{x}_1\\ x_2\end{array}\right],A=\left[\begin{array}{cc}{a}_{11}& a_{12}\\ a_{21}& a_{22}\end{array}\right],\mathbf{F}(t)=\left[\begin{array}{c}{f}_1(t)\\ f_2(t)\end{array}\right].$$ Here, ${\mathbf{X}}^{\prime }$ denotes the derivative of the vector $\mathbf{X}$. This formalism scales seamlessly to $n$ equations, transforming a daunting set of intertwined problems into a single, elegant matrix equation. This is not just notation; it unlocks the entire toolbox of linear algebra for finding solutions.

Converting an nth-Order ODE to a First-Order System

A crucial technique for solving complex, higher-order differential equations is to convert an nth-order ODE to a first-order system. This method allows you to apply matrix system solvers to problems that initially look quite different. The procedure is systematic:

1. Given an nth-order equation: $y^{(n)}+a_{n-1}(t)y^{(n-1)}+...+a_1(t)y^{\prime }+a_0(t)y=g(t)$.
2. Define $n$ new variables, each representing a lower-order derivative of the original $y$:

$$x_1=y,x_2=y^{\prime },x_3=y^{\prime \prime },...,x_n=y^{(n-1)}.$$

1. Take derivatives of these new variables to build your system:

$$\begin{array}{rl}{x}_1^{\prime }& =y^{\prime }=x_2\\ x_2^{\prime }& =y^{\prime \prime }=x_3\\ & ⋮\\ x_{n-1}^{\prime }& =y^{(n-1)}=x_n\\ x_n^{\prime }& =y^{(n)}=-a_0(t)x_1-a_1(t)x_2-...-a_{n-1}(t)x_n+g(t)\end{array}$$

For example, convert the third-order equation $y^{\prime \prime \prime }+2y^{\prime \prime }-y^{\prime }+3y=\sin (t)$. Define $x_1=y$, $x_2=y^{\prime }$, $x_3=y^{\prime \prime }$. The equivalent system is: $$\begin{array}{rl}{x}_1^{\prime }& =x_2\\ x_2^{\prime }& =x_3\\ x_3^{\prime }& =-3x_1+x_2-2x_3+\sin (t)\end{array}$$ In matrix form, this is ${\mathbf{X}}^{\prime }=\left[\begin{array}{ccc}0& 1& 0\\ 0& 0& 1\\ -3& 1& -2\end{array}\right]\mathbf{X}+\left[\begin{array}{c}0\\ 0\\ \sin (t)\end{array}\right]$.

Solution Structure of Linear Homogeneous Systems

For the homogeneous system ${\mathbf{X}}^{\prime }=A\mathbf{X}$ (where $\mathbf{F}(t)=0$), the solution structure mirrors that of single linear equations but in a vector space. The set of all solutions forms an $n$-dimensional vector space. To build a general solution, you need $n$ linearly independent solution vectors.

When the coefficient matrix $A$ is constant, we find solutions of the form $\mathbf{X}(t)=e^{\lambda t}\mathbf{v}$, where $\lambda$ is an eigenvalue of $A$ and $\mathbf{v}$ is its corresponding eigenvector. Substituting this guess into ${\mathbf{X}}^{\prime }=A\mathbf{X}$ leads to the eigenvalue problem $A\mathbf{v}=\lambda \mathbf{v}$.

• Case 1: Real and Distinct Eigenvalues. For each eigenvalue-eigenvector pair $({\lambda }_i,{\mathbf{v}}_i)$, you have a solution ${\mathbf{X}}_i(t)=e^{{\lambda }_{i}t}{\mathbf{v}}_i$. The general solution is a linear combination:

$$\mathbf{X}(t)=c_1{e}^{{\lambda }_{1}t}{\mathbf{v}}_1+c_2{e}^{{\lambda }_{2}t}{\mathbf{v}}_2+...+c_n{e}^{{\lambda }_{n}t}{\mathbf{v}}_n.$$

• Case 2: Complex Eigenvalues. These always come in conjugate pairs. For $\lambda =\alpha \pm i\beta$ with eigenvector $\mathbf{v}=\mathbf{a}\pm i\mathbf{b}$, two real, linearly independent solutions are:

$$e^{\alpha t}[\mathbf{a}\cos (\beta t)-\mathbf{b}\sin (\beta t)]\text{and}{e}^{\alpha t}[\mathbf{a}\sin (\beta t)+\mathbf{b}\cos (\beta t)].$$

• Case 3: Repeated Eigenvalues. The method requires generalized eigenvectors. For a repeated eigenvalue $\lambda$ with deficiency, you look for solutions of the form $e^{\lambda t}{\mathbf{v}}_1$ and $e^{\lambda t}(t{\mathbf{v}}_1+{\mathbf{v}}_2)$, where ${\mathbf{v}}_2$ satisfies $(A-\lambda I){\mathbf{v}}_2={\mathbf{v}}_1$.

The Fundamental Matrix and Solution Framework

A fundamental matrix $\Phi (t)$ is an $n\times n$ matrix whose columns are $n$ linearly independent solution vectors of ${\mathbf{X}}^{\prime }=A\mathbf{X}$. If ${\mathbf{X}}_1(t),{\mathbf{X}}_2(t),...,{\mathbf{X}}_n(t)$ are these solutions, then: $$\Phi (t)=\left[\begin{array}{cccc}|& |& & |\\ {\mathbf{X}}_1(t)& {\mathbf{X}}_2(t)& \cdots & {\mathbf{X}}_n(t)\\ |& |& & |\end{array}\right].$$ This matrix is a powerful conceptual and computational tool. The general solution to the homogeneous system can be written simply as $\mathbf{X}(t)=\Phi (t)\mathbf{c}$, where $\mathbf{c}$ is a column vector of arbitrary constants.

For the nonhomogeneous system ${\mathbf{X}}^{\prime }=A\mathbf{X}+\mathbf{F}(t)$, you find a particular solution using the variation of parameters formula, which generalizes beautifully via the fundamental matrix: $${\mathbf{X}}_p(t)=\Phi (t)\int \Phi (t{)}^{-1}\mathbf{F}(t)dt.$$ The general solution is then $\mathbf{X}(t)=\Phi (t)\mathbf{c}+{\mathbf{X}}_p(t)$. The fundamental matrix organizes the solution space, making this formula straightforward to apply, even for complicated forcing functions $\mathbf{F}(t)$.

Applications to Coupled Oscillators and Multi-Compartment Models

The real power of this theory is revealed in modeling physical and engineering systems. Two classic applications are coupled oscillators and multi-compartment models.

Coupled Oscillators: Consider two masses connected by springs to each other and to walls. Applying Newton's second law to each mass yields two second-order ODEs that are coupled—the equation for each mass's acceleration depends on the other's position. By defining state variables for position and velocity ($x_1,v_1,x_2,v_2$), you convert this to a $4\times 4$ first-order system ${\mathbf{X}}^{\prime }=A\mathbf{X}$. Solving this eigenvalue problem doesn't just give solutions; it reveals the system's normal modes—the specific frequencies and patterns at which the system will naturally oscillate. The eigenvalues ($\lambda =i\omega$) give the frequencies, and the eigenvectors describe the motion pattern of each mass in that mode.

Multi-Compartment Models: These are used to model flow between connected "compartments," such as the exchange of pollutants between lakes, the diffusion of drugs between body organs, or heat transfer between layered materials. If $x_i(t)$ represents the amount of substance (or temperature) in compartment $i$, the rate of change $x_i^{\prime }$ is typically a linear function of the amounts in all connected compartments. This directly creates a first-order linear system ${\mathbf{X}}^{\prime }=A\mathbf{X}$. The eigenvalues of $A$ determine the time scales of the system's approach to equilibrium, and analyzing the solution tells you how long it takes for a drug to clear or how a contaminant spreads.

Common Pitfalls

1. Misapplying Solution Methods: Attempting to use the $e^{\lambda t}\mathbf{v}$ method for a system where $A$ is not constant ($A=A(t)$). This method only works for constant-coefficient matrices. For non-constant matrices, you must generally use numerical methods or seek a fundamental matrix through other means.
2. Incorrect Eigenvector Handling for Complex Eigenvalues: A common mistake is to use the complex solution vectors $\mathbf{X}=e^{(\alpha +i\beta )t}(\mathbf{a}+i\mathbf{b})$ and its conjugate as the two fundamental real solutions. You must extract the real and imaginary parts as described to get two real-valued, linearly independent solutions for your general solution.
3. Forgetting to Convert All Initial Conditions: When converting an nth-order ODE with initial conditions like $y(0)=y_0,y^{\prime }(0)=y_0^{\prime },y^{\prime \prime }(0)=y_0^{\prime \prime }$, you must also convert these conditions for your new vector $\mathbf{X}$. They become $\mathbf{X}(0)=[x_1(0),x_2(0),x_3(0){]}^T=[y_0,y_0^{\prime },y_0^{\prime \prime }{]}^T$. Applying the old conditions to the wrong variables is a frequent error.
4. Matrix Multiplication Errors: When using the variation of parameters formula ${\mathbf{X}}_p=\Phi (t)\int {\Phi }^{-1}(t)\mathbf{F}(t)dt$, the order of multiplication is critical. The product inside the integral is ${\Phi }^{-1}(t)\mathbf{F}(t)$, which results in an $n\times 1$ column vector. You then integrate each component of this vector before multiplying by $\Phi (t)$ on the left.

Summary

• Systems of first-order linear ODEs are expressed compactly as ${\mathbf{X}}^{\prime }=A\mathbf{X}+\mathbf{F}(t)$, which enables the use of linear algebra techniques for solution.
• Any single nth-order linear ODE can be converted into an equivalent system of $n$ first-order ODEs by defining new variables for each derivative of the original function.
• The general solution to a homogeneous constant-coefficient system ${\mathbf{X}}^{\prime }=A\mathbf{X}$ is built from solutions of the form $e^{\lambda t}\mathbf{v}$, found by solving the eigenvalue problem for matrix $A$. The solution structure differs for cases of real distinct, complex conjugate, and repeated eigenvalues.
• A fundamental matrix $\Phi (t)$ is formed from a set of linearly independent solution vectors and provides a framework for writing the general homogeneous solution and for computing a particular solution to nonhomogeneous systems via the variation of parameters formula.
• These methods are directly applicable to essential engineering models like coupled oscillators (revealing normal modes of vibration) and multi-compartment models (analyzing flow and diffusion between connected units).