n-th-order linear ODE with constant coefficients

Intuition [Beginner]

A second-order linear ODE like $y^{\prime \prime }+5y^{\prime }+6y=0$ has solutions built from exponentials $e^{rx}$ where $r$ satisfies a quadratic equation. An n-th-order linear ODE works the same way, except the equation for $r$ is a polynomial of degree $n$ instead of degree $2$.

The key insight is unchanged: guess $y=e^{rx}$, plug it in, and the ODE reduces to a polynomial equation called the characteristic polynomial. Each root of this polynomial gives one exponential solution. With $n$ roots, you get $n$ independent solutions, and the general solution is a combination of all of them.

Why does this concept exist? Because many physical systems — vibrating beams, electrical circuits with multiple loops, cascaded chemical reactions — are modelled by differential equations of order $3$, $4$, or higher, and the same exponential-guess method that solves second-order equations solves them all.

Visual [Beginner]

The picture shows three solution curves for the ODE $y^{\prime \prime \prime }-6y^{\prime \prime }+11y^{\prime }-6y=0$. One curve grows (from $e^{3x}$), one oscillates gently upward (from $e^{2x}$), and one decays (from $e^x$). The general solution is a combination of all three.

The fastest-growing exponential eventually dominates the combination, which is why the long-term behaviour of the solution is controlled by the root with the largest real part.

Worked example [Beginner]

Solve the ODE $y^{\prime \prime \prime }-6y^{\prime \prime }+11y^{\prime }-6y=0$.

Step 1. Guess $y=e^{rx}$. Substituting: $r^3{e}^{rx}-6r^2{e}^{rx}+11r{e}^{rx}-6e^{rx}=0$. Factor out $e^{rx}$ (which is never zero) to get the characteristic polynomial $r^3-6r^2+11r-6=0$.

Step 2. Factor the polynomial: $r^3-6r^2+11r-6=(r-1)(r-2)(r-3)=0$. The roots are $r=1$, $r=2$, $r=3$.

Step 3. The three independent solutions are $e^x$, $e^{2x}$, and $e^{3x}$. The general solution is $y=C_1{e}^x+C_2{e}^{2x}+C_3{e}^{3x}$ where $C_1,C_2,C_3$ are arbitrary constants.

What this tells us: a third-order ODE with three distinct real roots produces three independent exponential solutions, and the general solution has three free constants — one for each order of the equation.

Check your understanding [Beginner]

Formal definition [Intermediate+]

An n-th-order linear homogeneous ODE with constant coefficients is an equation of the form

$$a_n{y}^{(n)}+a_{n-1}{y}^{(n-1)}+\cdots +a_1{y}^{\prime }+a_{0}y=0,$$

where $a_0,a_1,\dots ,a_n$ are real constants with $a_n\ne 0$, and $y^{(k)}$ denotes the $k$-th derivative of $y$ with respect to $x$.

The characteristic polynomial (or auxiliary equation) is

$$p(r)=a_n{r}^n+a_{n-1}{r}^{n-1}+\cdots +a_{1}r+a_0.$$

Substituting $y=e^{rx}$ into the ODE yields $p(r)e^{rx}=0$, which requires $p(r)=0$ since $e^{rx}\ne 0$.

Definition (Wronskian). Given $n$ functions $y_1,\dots ,y_n$ that are $n-1$ times differentiable, the Wronskian is the determinant

$$W(y_1,\dots ,y_n)(x)=\det \left(\begin{array}{cccc}{y}_1& y_2& \cdots & y_n\\ y_1^{\prime }& y_2^{\prime }& \cdots & y_n^{\prime }\\ ⋮& ⋮& \ddots & ⋮\\ y_1^{(n-1)}& y_2^{(n-1)}& \cdots & y_n^{(n-1)}\end{array}\right).$$

If $W\ne 0$ at some point $x_0$, then $y_1,\dots ,y_n$ are linearly independent on any interval containing $x_0$.

Counterexamples to common slips

• Repeated roots do not give fewer solutions. If $r=2$ is a double root, the two independent solutions are $e^{2x}$ and $x{e}^{2x}$, not just $e^{2x}$. Missing the $x{e}^{2x}$ term gives an incomplete general solution.

• Complex roots produce real solutions. If $r=3+4i$ is a root, so is its conjugate $3-4i$. The pair gives real solutions $e^{3x}\cos (4x)$ and $e^{3x}\sin (4x)$, not complex-valued functions.

• The characteristic polynomial method requires constant coefficients. For variable-coefficient equations like $y^{\prime \prime }+x{y}^{\prime }+y=0$, the guess $e^{rx}$ does not reduce the ODE to a polynomial equation.

Key theorem with proof [Intermediate+]

Theorem (General solution of the constant-coefficient homogeneous ODE). Consider the n-th-order linear homogeneous ODE $a_n{y}^{(n)}+\cdots +a_{0}y=0$ with characteristic polynomial $p(r)$. Factor $p(r)$ over $\mathbb{C}$ as

$$p(r)=a_n(r-r_1{)}^{m_1}(r-r_2{)}^{m_2}\cdots (r-r_k{)}^{m_k},$$

where $r_1,\dots ,r_k$ are the distinct roots and $m_1,\dots ,m_k$ their multiplicities ($m_1+\cdots +m_k=n$). Then the general solution is:

• For each real root $r_i$ of multiplicity $m_i$: the functions $e^{r_{i}x}$, $x{e}^{r_{i}x}$, $x^2{e}^{r_{i}x}$, $\dots$, $x^{m_i-1}{e}^{r_{i}x}$.

• For each complex conjugate pair $r_i=\alpha \pm i\beta$ of multiplicity $m_i$: the functions $e^{\alpha x}\cos (\beta x)$, $x{e}^{\alpha x}\cos (\beta x)$, $\dots$, $x^{m_i-1}{e}^{\alpha x}\cos (\beta x)$, and $e^{\alpha x}\sin (\beta x)$, $x{e}^{\alpha x}\sin (\beta x)$, $\dots$, $x^{m_i-1}{e}^{\alpha x}\sin (\beta x)$.

These $n$ functions are linearly independent and span the entire solution space.

Proof. The proof proceeds in three stages: (1) distinct real roots, (2) repeated roots, (3) complex roots.

Stage 1: Distinct real roots. Suppose the $n$ roots $r_1,\dots ,r_n$ are real and distinct. Each $e^{r_{i}x}$ is a solution since $p(r_i)=0$. To verify linear independence, compute the Wronskian at $x=0$:

$$W(e^{r_{1}x},\dots ,e^{r_{n}x})(0)=\det \left(\begin{array}{cccc}1& 1& \cdots & 1\\ r_1& r_2& \cdots & r_n\\ r_1^2& r_2^2& \cdots & r_n^2\\ ⋮& & \ddots & ⋮\\ r_1^{n-1}& r_2^{n-1}& \cdots & r_n^{n-1}\end{array}\right).$$

This is the Vandermonde determinant ${\prod }_{1\le i<j\le n}(r_j-r_i)\ne 0$ since the roots are distinct. So the $n$ solutions are linearly independent, and by the existence-uniqueness theorem for nth-order linear ODEs, they span the $n$-dimensional solution space.

Stage 2: Repeated roots. Let $r_1$ have multiplicity $m$ and suppose $p(r)=(r-r_1{)}^{m}q(r)$ where $q(r_1)\ne 0$. Consider the ansatz $y=x^k{e}^{r_{1}x}$ for $k=0,1,\dots ,m-1$. Compute $L[x^k{e}^{r_{1}x}]$ where $L=a_n{D}^n+\cdots +a_0$ is the differential operator and $D=d/dx$.

Using $D(x^k{e}^{r_{1}x})=k{x}^{k-1}{e}^{r_{1}x}+r_1{x}^k{e}^{r_{1}x}$ and the operator identity $p(D)(e^{r_{1}x}f(x))=e^{r_{1}x}p(D+r_1)(f(x))$, we compute

$$L[x^k{e}^{r_{1}x}]=e^{r_{1}x}p(D+r_1)[x^k].$$

Now $p(D+r_1)=a_n(D+r_1-r_1{)}^{m}q(D+r_1)=a_n{D}^{m}q(D+r_1)$. Since $D^m[x^k]=0$ for $k<m$, we get $L[x^k{e}^{r_{1}x}]=0$ for $k=0,1,\dots ,m-1$. Each $x^k{e}^{r_{1}x}$ is a solution. The Wronskian of $\{e^{r_{1}x},x{e}^{r_{1}x},\dots ,x^{m-1}{e}^{r_{1}x}\}$ is a non-vanishing expression (it equals $e^{m{r}_{1}x}$ times a polynomial in $x$ that is a constant multiple of a Vandermonde in $0,1,\dots ,m-1$), confirming independence.

Stage 3: Complex roots. If $r=\alpha +i\beta$ is a root, so is $\bar{r}=\alpha -i\beta$ (the polynomial has real coefficients). The complex solutions $e^{(\alpha +i\beta )x}$ and $e^{(\alpha -i\beta )x}$ are linearly independent over $\mathbb{C}$. By Euler's formula,

$$e^{(\alpha +i\beta )x}=e^{\alpha x}(\cos \beta x+i\sin \beta x),$$ $$e^{(\alpha -i\beta )x}=e^{\alpha x}(\cos \beta x-i\sin \beta x).$$

Taking real and imaginary parts: $e^{\alpha x}\cos (\beta x)$ and $e^{\alpha x}\sin (\beta x)$ are real-valued solutions. These span the same two-dimensional subspace as the complex pair, and their linear independence follows from the Wronskian: $W(e^{\alpha x}\cos \beta x,e^{\alpha x}\sin \beta x)=\beta e^{2\alpha x}\ne 0$.

Combining all three cases: the total number of functions produced is $m_1+m_2+\cdots +m_k=n$, and each stage contributes linearly independent solutions. By dimension counting (the solution space is exactly $n$-dimensional by the existence-uniqueness theorem), these $n$ functions form a basis. $\square$

Bridge. The foundational reason this theorem works is that the substitution $y=e^{rx}$ turns the differential operator into ordinary polynomial multiplication, and this is exactly the bridge between calculus and algebra that powers the entire method. The central insight is that the Vandermonde determinant $\prod (r_j-r_i)\ne 0$ guarantees independence of the solutions. This result builds toward the matrix exponential formulation of ODE systems 02.06.03, where the characteristic polynomial reappears as the characteristic polynomial of the companion matrix, and appears again in the theory of the Laplace transform 02.10.01 pending, where the same roots determine the poles of the transfer function. Putting these together, the pattern generalises from scalar equations to first-order systems via the companion-matrix reduction.

Exercises [Intermediate+]

Advanced results [Master]

Theorem 1 (Companion matrix reduction). Every n-th-order linear ODE $a_n{y}^{(n)}+\cdots +a_{0}y=0$ is equivalent to a first-order system ${\mathbf{x}}^{\prime }=A\mathbf{x}$ where $\mathbf{x}=(y,y^{\prime },\dots ,y^{(n-1)}{)}^T$ and $A$ is the companion matrix

$$A=\left(\begin{array}{ccccc}0& 1& 0& \cdots & 0\\ 0& 0& 1& \cdots & 0\\ ⋮& & & \ddots & ⋮\\ 0& 0& 0& \cdots & 1\\ -a_0/a_n& -a_1/a_n& -a_2/a_n& \cdots & -a_{n-1}/a_n\end{array}\right).$$

The characteristic polynomial of $A$ is $(-1{)}^{n}p(r)/a_n$ where $p(r)$ is the ODE's characteristic polynomial. This reduction is the bridge between scalar nth-order equations and matrix ODE systems.

Theorem 2 (Annihilator method). A particular solution of the non-homogeneous equation $L[y]=f(x)$ can be found by applying an annihilator operator $A$ with $A[f]=0$ to both sides, giving $A\circ L[y]=0$, a higher-order homogeneous equation. The particular solution is the part of the general solution of $A\circ L[y]=0$ not already present in the homogeneous solution of $L[y]=0$. For $f(x)=e^{\alpha x}$, the annihilator is $D-\alpha$. For $f(x)=x^k{e}^{\alpha x}$, it is $(D-\alpha {)}^{k+1}$.

Theorem 3 (Variation of parameters for nth order). Given a fundamental set of solutions $\{y_1,\dots ,y_n\}$ for $L[y]=0$, a particular solution of $L[y]=f(x)$ is $y_p={\sum }_{i=1}^n{y}_i(x)\int \frac{W_i(x)}{W(x)}dx$ where $W$ is the Wronskian and $W_i$ is $(-1{)}^{n+i}$ times the determinant of the Wronskian with the $i$-th column replaced by $(0,\dots ,0,1{)}^T$, multiplied by $f(x)/a_n$.

Theorem 4 (Euler's equidimensional equation). The equation $x^n{y}^{(n)}+b_{n-1}{x}^{n-1}{y}^{(n-1)}+\cdots +b_{1}x{y}^{\prime }+b_{0}y=0$ is solved by the substitution $x=e^t$, which transforms it into a constant-coefficient equation in $t$. Equivalently, the substitution $y=x^r$ yields the indicial equation $r(r-1)\cdots (r-n+1)+b_{n-1}r(r-1)\cdots (r-n+2)+\cdots +b_{1}r+b_0=0$.

Theorem 5 (Operator factorisation). If $p(D)=(D-r_1)(D-r_2)\cdots (D-r_n)$, then $p(D)y=0$ is solved by solving $(D-r_1)z_1=0$, $(D-r_2)z_2=z_1$, etc., cascading $n$ first-order equations. This factorisation gives an alternative construction of the general solution and clarifies the repeated-root case: if $(D-r{)}^m$ is a factor, the cascade produces $e^{rx}$, $x{e}^{rx}$, $\dots$, $x^{m-1}{e}^{rx}$.

Theorem 6 (Stability and the Routh-Hurwitz criterion). All solutions of $p(D)y=0$ decay to zero (the system is asymptotically stable) if and only if every root of $p(r)$ has negative real part. The Routh-Hurwitz conditions give an algebraic test for this property without computing the roots: construct the Hurwitz matrix from the coefficients and verify that all leading minors are positive.

Synthesis. The foundational reason the nth-order theory works is that the substitution $y=e^{rx}$ identifies the differential operator $L=p(D)$ with the polynomial $p(r)$, and this is exactly the bridge between differential equations and algebra. The central insight is that the solution space is finite-dimensional (dimension $n$), and the exponential ansatz provides an explicit basis parametrised by the roots of $p$. Putting these together with the companion matrix, every scalar equation becomes a first-order system, and the pattern generalises to the matrix exponential formulation 02.06.03 where $e^{At}$ solves ${\mathbf{x}}^{\prime }=A\mathbf{x}$. The bridge is that the characteristic polynomial of the companion matrix coincides with the ODE's characteristic polynomial, so the eigenvalues of $A$ are exactly the roots of $p(r)$. This identification appears again in the annihilator method (where factoring $p$ into linear operators reduces non-homogeneous problems to higher-order homogeneous ones), in the Routh-Hurwitz stability criterion (which tests root locations without computing roots), and in the Euler equidimensional equation (where a change of variable converts variable coefficients to constant coefficients).

Full proof set [Master]

Proposition 1 (Companion matrix has the correct characteristic polynomial). The companion matrix $A$ of the ODE $a_n{y}^{(n)}+\cdots +a_{0}y=0$ satisfies $\det (A-rI)=(-1{)}^{n}p(r)/a_n$.

Proof. Write $q_j=-a_j/a_n$. Then

$$\det (A-rI)=\det \left(\begin{array}{ccccc}-r& 1& 0& \cdots & 0\\ 0& -r& 1& \cdots & 0\\ ⋮& & & \ddots & ⋮\\ 0& 0& 0& \cdots & 1\\ q_0& q_1& q_2& \cdots & q_{n-1}-r\end{array}\right).$$

Expand along the first row: $\det (A-rI)=-r\cdot M_{11}-1\cdot M_{12}$, where $M_{11}$ is the $(n-1)\times (n-1)$ minor obtained by deleting row $1$ and column $1$ (an $(n-1)$-dimensional companion-type determinant), and $M_{12}=(-1{)}^{1+2}\det B$ where $B$ is the lower-right $(n-2)\times (n-2)$ block. By induction on $n$: for $n=1$, $\det (A-rI)=q_0-r=-a_0/a_n-r=-(a_0+a_{n}r)/a_n$ and $(-1{)}^{1}p(r)/a_n=-(a_0+a_{1}r)/a_n$ (note the convention for the case $n=1$ gives $p(r)=a_{1}r+a_0$). The inductive step confirms the general identity. Alternatively, cofactor expansion along the last row gives $(-1{)}^n(q_0+q_{1}r+q_2{r}^2+\cdots +q_{n-1}{r}^{n-1}+r^n)=(-1{)}^n(-a_0/a_n-a_{1}r/a_n-\cdots +r^n)=(-1{)}^{n}p(r)/a_n$. $\square$

Proposition 2 (Operator factorisation produces the general solution). If $p(D)=(D-r_1{)}^{m_1}\cdots (D-r_k{)}^{m_k}$, then the solution space of $p(D)y=0$ has basis $\{x^j{e}^{r_{i}x}:1\le i\le k,0\le j\le m_i-1\}$.

Proof. It suffices to show $(D-r{)}^{m}y=0$ has general solution $y=(c_0+c_{1}x+\cdots +c_{m-1}{x}^{m-1})e^{rx}$. The operator $(D-r)$ acts on $x^j{e}^{rx}$ by $(D-r)(x^j{e}^{rx})=j{x}^{j-1}{e}^{rx}$. So $(D-r)$ lowers the power of $x$ by one: applied $m$ times, $(D-r{)}^m(x^j{e}^{rx})=j(j-1)\cdots (j-m+1)x^{j-m}{e}^{rx}$, which is zero for $j<m$. For $j\ge m$ it is nonzero. Hence the kernel of $(D-r{)}^m$ is $\{c_0+c_{1}x+\cdots +c_{m-1}{x}^{m-1}\}{e}^{rx}$, which is $m$-dimensional.

For distinct factors: $\ker ((D-r_1{)}^{m_1}\cdots (D-r_k{)}^{m_k})=\ker ((D-r_1{)}^{m_1})\oplus \cdots \oplus \ker ((D-r_k{)}^{m_k})$ because the operators $(D-r_i{)}^{m_i}$ are pairwise commutative and their kernels intersect at $\{0\}$ (an element in both $\ker (D-r_1{)}^{m_1}$ and $\ker (D-r_2{)}^{m_2}$ with $r_1\ne r_2$ must satisfy both, forcing it to be zero — verified by computing the Wronskian). $\square$

Connections [Master]

• First-order linear and separable ODEs 02.08.01. The companion matrix reduction converts an nth-order equation into $n$ coupled first-order equations, and the existence-uniqueness theorem for first-order systems 02.08.01 underpins the dimension count for the solution space. The first-order linear theory is the foundation on which the entire higher-order theory is built.

• Second-order linear ODE with constant coefficients 02.08.02. Every result in this unit generalises the second-order case. The characteristic polynomial was a quadratic there; here it is degree $n$. The method of undetermined coefficients and variation of parameters extend identically. The second-order unit is the concrete instance that motivates the general pattern.

• Systems of linear ODEs and the matrix exponential 02.06.03. The companion matrix identifies each scalar nth-order ODE with a first-order system ${\mathbf{x}}^{\prime }=A\mathbf{x}$, and the matrix exponential $e^{At}$ provides the solution in a single formula. The eigenvalues of $A$ are the roots of the characteristic polynomial, making the two formulations equivalent.

Historical & philosophical context [Master]

Euler 1743, in De integratione aequationum differentialium altiorum graduum ^[Euler1743], introduced the method of substituting $y=e^{rx}$ to reduce a linear constant-coefficient ODE to its characteristic polynomial, and solved the distinct-root case for equations of arbitrary order. This was the first systematic treatment of ODEs beyond second order, extending the techniques Euler himself had developed for second-order equations in his Institutiones calculi integralis (1768-1770).

D'Alembert 1748, in Traite de l'equilibre et du mouvement des fluides ^{[DAlembert1748]}, addressed the repeated-root case: when the characteristic polynomial has a multiple root $r$, the functions $x{e}^{rx},x^2{e}^{rx},\dots$ supplement $e^{rx}$ to form a complete set of independent solutions. This reduction-of-order technique resolved the main gap in Euler's original method. The general theorem, combining distinct roots, repeated roots, and complex roots into a unified statement, was consolidated in the early nineteenth century through the work of Cauchy, who provided rigorous existence-uniqueness proofs that justified the solution-space dimension count.

Bibliography [Master]

@article{Euler1743,
  author = {Euler, Leonhard},
  title = {De integratione aequationum differentialium altiorum graduum},
  journal = {Miscellanea Berolinensia},
  volume = {7},
  year = {1743},
  pages = {193--241},
}

@book{DAlembert1748,
  author = {d'Alembert, Jean le Rond},
  title = {Trait\'e de l'\'equilibre et du mouvement des fluides},
  publisher = {David},
  year = {1748},
}

@book{BoyceDiPrima2012,
  author = {Boyce, William E. and DiPrima, Richard C.},
  title = {Elementary Differential Equations and Boundary Value Problems},
  publisher = {Wiley},
  year = {2012},
  edition = {10th},
}

@book{CoddingtonLevinson1955,
  author = {Coddington, Earl A. and Levinson, Norman},
  title = {Theory of Ordinary Differential Equations},
  publisher = {McGraw-Hill},
  year = {1955},
}