# Linear differential equation

(Redirected from Linear differential equations)

In mathematics, linear differential equations are differential equations having solutions which can be added together in particular linear combinations to form further solutions. They equate 0 to a polynomial that is linear in the value and various derivatives of a variable; its linearity means that each term in the polynomial has degree either 0 or 1.

Linear differential equations can be ordinary (ODEs) or partial (PDEs).

The solutions to (homogeneous) linear differential equations form a vector space (unlike non-linear differential equations).

## Basic features

Linear differential equations are of the form

${\displaystyle Ly=f}$

where the differential operator L is a linear operator, y is the unknown function, and the right hand side f is a given function (called the source term) of the same variable. For a function dependent on time we may write the equation more expressly as

${\displaystyle Ly(t)=f(t)}$

and, even more precisely by bracketing

${\displaystyle L[y(t)]=f(t)}$.

The linear operator L may be considered to be of the form[1]

${\displaystyle L_{n}(y)\equiv {\frac {d^{n}y}{dt^{n}}}+A_{1}(t){\frac {d^{n-1}y}{dt^{n-1}}}+\cdots +A_{n-1}(t){\frac {dy}{dt}}+A_{n}(t)y}$

The linearity condition on L rules out operations such as taking the square of the derivative of y; but permits, for example, taking the second derivative of y. It is convenient to rewrite this equation in an operator form

${\displaystyle L_{n}(y)\equiv \left[\,D^{n}+A_{1}(t)D^{n-1}+\cdots +A_{n-1}(t)D+A_{n}(t)\right]y}$

where D is the differential operator d/dt (i.e. Dy = y' = dy/dt , D2y = y" = d2y/dt2,... ), and the An are given functions.

Such an equation is said to have order n, the index of the highest derivative of y that is involved.

A typical simple example is the linear differential equation used to model radioactive decay.[2] Let N(t) denote the number of radioactive atoms remaining in some sample of material [3] at time t. Then for some constant k > 0, the rate at which the radioactive atoms decay can be modelled by

${\displaystyle {\frac {dN}{dt}}=-kN}$

If y is assumed to be a function of only one variable, one speaks about an ordinary differential equation, else the derivatives and their coefficients must be understood as (contracted) vectors, matrices or tensors of higher rank, and we have a (linear) partial differential equation.

The case where f = 0 is called a homogeneous equation and its solutions are called complementary functions. It is particularly important to the solution of the general case, since any complementary function can be added to a solution of the inhomogeneous equation to give another solution (by a method traditionally called particular integral and complementary function). When the Ai are numbers, the equation is said to have constant coefficients.

## Homogeneous equations with constant coefficients

The first method of solving linear homogeneous ordinary differential equations with constant coefficients is due to Euler, who realized that solutions have the form ezx, for possibly-complex values of z. The exponential function is one of the few functions to keep its shape after differentiation, allowing the sum of its multiple derivatives to cancel out to zero, as required by the equation. Thus, for constant values A1,..., An, to solve:

${\displaystyle y^{(n)}+A_{1}y^{(n-1)}+\cdots +A_{n}y=0\,,}$

we set y = ezx, leading to

${\displaystyle z^{n}e^{zx}+A_{1}z^{n-1}e^{zx}+\cdots +A_{n}e^{zx}=0.}$

Division by ezx gives the nth-order polynomial:

${\displaystyle F(z)=z^{n}+A_{1}z^{n-1}+\cdots +A_{n}=0.\,}$

This algebraic equation F(z) = 0 is the characteristic equation considered later by Gaspard Monge and Augustin-Louis Cauchy.

Formally, the terms ${\displaystyle y^{(k)}(k=1,2,\dots ,n)}$ of the original differential equation are replaced by zk. Solving the polynomial gives n values of z, z1, ..., zn. Substitution of any of those values for z into ezx gives a solution ezix. Since homogeneous linear differential equations obey the superposition principle, any linear combination of these functions also satisfies the differential equation.

When these roots are all distinct, we have n distinct solutions to the differential equation. It can be shown that these are linearly independent, by applying the Vandermonde determinant, and together they form a basis of the space of all solutions of the differential equation.

Examples
${\displaystyle y''''-2y'''+2y''-2y'+y=0}$

has the characteristic equation

${\displaystyle z^{4}-2z^{3}+2z^{2}-2z+1=0.}$

This has zeroes, i, −i, and 1 (multiplicity 2). The solution basis is then

${\displaystyle e^{ix},\,e^{-ix},\,e^{x},\,xe^{x}.}$

This corresponds to the real-valued solution basis

${\displaystyle \cos x,\,\sin x,\,e^{x},\,xe^{x}\,.}$

The preceding gave a solution for the case when all zeros are distinct, that is, each has multiplicity 1. For the general case, if z is a (possibly complex) zero (or root) of F(z) having multiplicity m, then, for ${\displaystyle k\in \{0,1,\dots ,m-1\}\,}$, ${\displaystyle y=x^{k}e^{zx}}$ is a solution of the ordinary differential equation. Applying this to all roots gives a collection of n distinct and linearly independent functions, where n is the degree of F(z). As before, these functions make up a basis of the solution space.

If the coefficients Ai of the differential equation are real, then real-valued solutions are generally preferable. Since non-real roots z then come in conjugate pairs, so do their corresponding basis functions xkezx, and the desired result is obtained by replacing each pair with their real-valued linear combinations Re(y) and Im(y), where y is one of the pair.

A case that involves complex roots can be solved with the aid of Euler's formula.

### Second-order case

In the n=2 case

${\displaystyle y''+ay'+by=0,}$

the characteristic equation is of the form

${\displaystyle z^{2}+az+b=0.\,}$

We then can solve for z. There are three particular cases of interest:

• Case #1: Two distinct roots, z1 and z2
• Case #2: One real repeated root, z
• Case #3: Complex roots, α ± βi

In case #1, the general solution is given by

${\displaystyle y=c_{1}e^{z_{1}x}+c_{2}e^{z_{2}x}.\,}$

In case #2, the general solution is given by

${\displaystyle y=(c_{1}+c_{2}x)e^{zx}.}$

In case #3, the general solution is given, using Euler's equation, by

${\displaystyle y=c_{1}e^{\alpha x}\cos(\beta x)+c_{2}e^{\alpha x}\sin(\beta x),\,}$
${\displaystyle \alpha ={\mathop {\rm {Re}}}(z),\,}$
${\displaystyle \beta ={\mathop {\rm {Im}}}(z).\,}$

In each case, the constants ${\displaystyle c_{1},\,c_{2}}$ are functions of the initial conditions ${\displaystyle y(0),\,y'(0).}$ They can be found by using the values of the initial conditions in the solution equation for y and in the resulting equation for y' , giving two equations in the two unknown parameters.

### Examples

Given ${\displaystyle y''-4y'+5y=0}$. The characteristic equation is ${\displaystyle z^{2}-4z+5=0}$ which has roots "2±i". Thus the solution basis ${\displaystyle \{y_{1},y_{2}\}}$ is ${\displaystyle \{e^{(2+i)x},e^{(2-i)x}\}}$. Now y is a solution if and only if ${\displaystyle y=c_{1}y_{1}+c_{2}y_{2}}$ for ${\displaystyle c_{1},c_{2}\in \mathbf {C} }$.

Because the coefficients are real,

• we are likely not interested in the complex solutions
• our basis elements are mutual conjugates

The linear combinations

${\displaystyle u_{1}={\mbox{Re}}(y_{1})={\tfrac {1}{2}}(y_{1}+y_{2})=e^{2x}\cos(x),}$
${\displaystyle u_{2}={\mbox{Im}}(y_{1})={\tfrac {1}{2i}}(y_{1}-y_{2})=e^{2x}\sin(x),}$

will give us a real basis in ${\displaystyle \{u_{1},u_{2}\}}$.

#### Simple harmonic oscillator

The second order differential equation

${\displaystyle D^{2}y=-k^{2}y,}$

which represents a simple harmonic oscillator, can be restated as

${\displaystyle (D^{2}+k^{2})y=0.}$

The expression in parenthesis can be factored out, yielding

${\displaystyle ((D+ik)(D-ik))y=0,}$

which has a pair of linearly independent solutions:

${\displaystyle (D-ik)y=0}$
${\displaystyle (D+ik)y=0.}$

The solutions are, respectively,

${\displaystyle y_{0}=A_{0}e^{ikx}}$

and

${\displaystyle y_{1}=A_{1}e^{-ikx}.}$

These solutions provide a basis for the two-dimensional solution space of the second order differential equation: meaning that linear combinations of these solutions will also be solutions. In particular, the following solutions can be constructed

${\displaystyle y_{0'}={C_{0}e^{ikx}+C_{0}e^{-ikx} \over 2}=C_{0}\cos(kx)}$

and

${\displaystyle y_{1'}={C_{1}e^{ikx}-C_{1}e^{-ikx} \over 2i}=C_{1}\sin(kx).}$

These last two trigonometric solutions are linearly independent, so they can serve as another basis for the solution space, yielding the following general solution:

${\displaystyle y_{H}=C_{0}\cos(kx)+C_{1}\sin(kx).}$

#### Damped harmonic oscillator

Given the equation for the damped harmonic oscillator:

${\displaystyle \left(D^{2}+{\frac {b}{m}}D+\omega _{0}^{2}\right)y=0,}$

the expression in parentheses can be factored out: first obtain the characteristic equation by replacing D with z. This equation must be satisfied for all y, thus:

${\displaystyle z^{2}+{\frac {b}{m}}z+\omega _{0}^{2}=0.}$

${\displaystyle z={\tfrac {1}{2}}\left(-{\frac {b}{m}}\pm {\sqrt {{\frac {b^{2}}{m^{2}}}-4\omega _{0}^{2}}}\right).}$

Use these characteristic roots to factor the left side of the original differential equation:

${\displaystyle \left(D+{\frac {b}{2m}}-{\sqrt {{\frac {b^{2}}{4m^{2}}}-\omega _{0}^{2}}}\right)\left(D+{\frac {b}{2m}}+{\sqrt {{\frac {b^{2}}{4m^{2}}}-\omega _{0}^{2}}}\right)y=0.}$

This implies a pair of solutions, one corresponding to

${\displaystyle \left(D+{\frac {b}{2m}}-{\sqrt {{\frac {b^{2}}{4m^{2}}}-\omega _{0}^{2}}}\right)y=0}$
${\displaystyle \left(D+{\frac {b}{2m}}+{\sqrt {{\frac {b^{2}}{4m^{2}}}-\omega _{0}^{2}}}\right)y=0}$

The solutions are, respectively,

${\displaystyle y_{0}=A_{0}e^{-\omega x+{\sqrt {\omega ^{2}-\omega _{0}^{2}}}x}=A_{0}e^{-\omega x}e^{{\sqrt {\omega ^{2}-\omega _{0}^{2}}}x}}$
${\displaystyle y_{1}=A_{1}e^{-\omega x-{\sqrt {\omega ^{2}-\omega _{0}^{2}}}x}=A_{1}e^{-\omega x}e^{-{\sqrt {\omega ^{2}-\omega _{0}^{2}}}x}}$

where ω = b/2m. From this linearly independent pair of solutions can be constructed another linearly independent pair which thus serve as a basis for the two-dimensional solution space:

${\displaystyle y_{H}(A_{0},A_{1})(x)=\left(A_{0}\sinh \left({\sqrt {\omega ^{2}-\omega _{0}^{2}}}x\right)+A_{1}\cosh \left({\sqrt {\omega ^{2}-\omega _{0}^{2}}}x\right)\right)e^{-\omega x}.}$

However, if |ω| < |ω0| then it is preferable to get rid of the consequential imaginaries, expressing the general solution as

${\displaystyle y_{H}(A_{0},A_{1})(x)=\left(A_{0}\sin \left({\sqrt {\omega _{0}^{2}-\omega ^{2}}}x\right)+A_{1}\cos \left({\sqrt {\omega _{0}^{2}-\omega ^{2}}}x\right)\right)e^{-\omega x}.}$

This latter solution corresponds to the underdamped case, whereas the former one corresponds to the overdamped case: the solutions for the underdamped case oscillate whereas the solutions for the overdamped case do not.

## Nonhomogeneous equation with constant coefficients

To obtain the solution to the nonhomogeneous equation (sometimes called inhomogeneous equation), find a particular integral yP(x) by the method of undetermined coefficients, method of variation of parameters or the use of the exponential response formula (below); the general solution to the linear differential equation is the sum of the general solution of the related homogeneous equation and the particular integral. Or, when the initial conditions are set, use Laplace transform to obtain the particular solution directly.

Suppose we face

${\displaystyle {\frac {d^{n}y(x)}{dx^{n}}}+A_{1}{\frac {d^{n-1}y(x)}{dx^{n-1}}}+\cdots +A_{n}y(x)=f(x).}$

For later convenience, define the characteristic polynomial

${\displaystyle P(z)=z^{n}+A_{1}z^{n-1}+\cdots +A_{n}.}$

We find a solution basis ${\displaystyle \{y_{1}(x),y_{2}(x),\ldots ,y_{n}(x)\}}$ for the homogeneous (f(x) = 0) case. We now seek a particular integral yp(x) by the variation of parameters method. Let the coefficients of the linear combination be functions of x:

${\displaystyle y_{p}(x)=u_{1}(x)y_{1}(x)+u_{2}(x)y_{2}(x)+\cdots +u_{n}(x)y_{n}(x).}$

For ease of notation we will drop the dependency on x (i.e. the various (x)). Using the operator notation D = d/dx, the ODE in question is P(D)y = f; so

${\displaystyle f=P(D)y_{p}=P(D)(u_{1}y_{1})+P(D)(u_{2}y_{2})+\cdots +P(D)(u_{n}y_{n}).}$

With the constraints

${\displaystyle 0=u'_{1}y_{1}+u'_{2}y_{2}+\cdots +u'_{n}y_{n}}$
${\displaystyle 0=u'_{1}y'_{1}+u'_{2}y'_{2}+\cdots +u'_{n}y'_{n}}$
${\displaystyle \cdots }$
${\displaystyle 0=u'_{1}y_{1}^{(n-2)}+u'_{2}y_{2}^{(n-2)}+\cdots +u'_{n}y_{n}^{(n-2)}}$

the parameters commute out,

${\displaystyle f=u_{1}P(D)y_{1}+u_{2}P(D)y_{2}+\cdots +u_{n}P(D)y_{n}+u'_{1}y_{1}^{(n-1)}+u'_{2}y_{2}^{(n-1)}+\cdots +u'_{n}y_{n}^{(n-1)}.}$

But P(D)yj = 0, therefore

${\displaystyle f=u'_{1}y_{1}^{(n-1)}+u'_{2}y_{2}^{(n-1)}+\cdots +u'_{n}y_{n}^{(n-1)}.}$

This, with the constraints, gives a linear system in the u′j. This much can always be solved; in fact, combining Cramer's rule with the Wronskian,

${\displaystyle u'_{j}=(-1)^{n+j}{\frac {W(y_{1},\ldots ,y_{j-1},y_{j+1}\ldots ,y_{n})_{0 \choose f}}{W(y_{1},y_{2},\ldots ,y_{n})}}.}$

In the very non-standard notation used above, one should take the i,n-minor of W and multiply it by f. That's why we get a minus-sign. Alternatively, forget about the minus sign and just compute the determinant of the matrix obtained by substituting the j-th W column with (0, 0, ..., f).

The rest is a matter of integrating u′j.

The particular integral is not unique; ${\displaystyle y_{p}+c_{1}y_{1}+\cdots +c_{n}y_{n}}$ also satisfies the ODE for any set of constants cj.

### Exponential response formula

The particular solution of

${\displaystyle P(D)y=\sum _{i}\left({a_{i}e^{r_{i}x}}\right)}$

can be found as

${\displaystyle y_{p}=\sum _{i}\left({{\frac {a_{i}}{P\left(r_{i}\right)}}e^{r_{i}x}}\right).}$

### Example

Suppose ${\displaystyle y''-4y'+5y=\sin(kx)}$. We take the solution basis found above ${\displaystyle \{e^{(2+i)x}=y_{1}(x),e^{(2-i)x}=y_{2}(x)\}}$.

{\displaystyle {\begin{aligned}W&={\begin{vmatrix}e^{(2+i)x}&e^{(2-i)x}\\(2+i)e^{(2+i)x}&(2-i)e^{(2-i)x}\end{vmatrix}}=e^{4x}{\begin{vmatrix}1&1\\2+i&2-i\end{vmatrix}}=-2ie^{4x}\\u'_{1}&={\frac {1}{W}}{\begin{vmatrix}0&e^{(2-i)x}\\\sin(kx)&(2-i)e^{(2-i)x}\end{vmatrix}}=-{\tfrac {i}{2}}\sin(kx)e^{(-2-i)x}\\u'_{2}&={\frac {1}{W}}{\begin{vmatrix}e^{(2+i)x}&0\\(2+i)e^{(2+i)x}&\sin(kx)\end{vmatrix}}={\tfrac {i}{2}}\sin(kx)e^{(-2+i)x}.\end{aligned}}}
${\displaystyle u_{1}=-{\tfrac {i}{2}}\int \sin(kx)e^{(-2-i)x}\,dx={\frac {ie^{(-2-i)x}}{2(3+4i+k^{2})}}\left((2+i)\sin(kx)+k\cos(kx)\right)}$
${\displaystyle u_{2}={\tfrac {i}{2}}\int \sin(kx)e^{(-2+i)x}\,dx={\frac {ie^{(i-2)x}}{2(3-4i+k^{2})}}\left((i-2)\sin(kx)-k\cos(kx)\right).}$

And so

{\displaystyle {\begin{aligned}y_{p}&=u_{1}(x)y_{1}(x)+u_{2}(x)y_{2}(x)={\frac {i}{2(3+4i+k^{2})}}\left((2+i)\sin(kx)+k\cos(kx)\right)+{\frac {i}{2(3-4i+k^{2})}}\left((i-2)\sin(kx)-k\cos(kx)\right)\\&={\frac {(5-k^{2})\sin(kx)+4k\cos(kx)}{(3+k^{2})^{2}+16}}.\end{aligned}}}

(Notice that u1 and u2 had factors that canceled y1 and y2; that is typical.)

For interest's sake, this ODE has a physical interpretation as a driven damped harmonic oscillator; yp represents the steady state, and ${\displaystyle c_{1}y_{1}+c_{2}y_{2}}$ is the transient.

As ${\displaystyle \sin(kx)={\frac {e^{ikx}-e^{-ikx}}{2i}},}$ the method of the exponential response formula produces

{\displaystyle {\begin{aligned}y_{p}&={\frac {{\frac {e^{ikx}}{5-k^{2}-4ik}}-{\frac {e^{-ikx}}{5-k^{2}+4ik}}}{2i}}\\&={\frac {(5-k^{2})\sin(kx)+4k\cos(kx)}{(5-k^{2})^{2}+(4k)^{2}}},\end{aligned}}}

## Equation with variable coefficients

A linear ODE of order n with variable coefficients has the general form

${\displaystyle p_{n}(x)y^{(n)}(x)+p_{n-1}(x)y^{(n-1)}(x)+\cdots +p_{0}(x)y(x)=r(x).}$

### Examples

A simple example is the Cauchy–Euler equation often used in engineering

${\displaystyle x^{n}y^{(n)}(x)+a_{n-1}x^{n-1}y^{(n-1)}(x)+\cdots +a_{0}y(x)=0.}$

## First-order equation with variable coefficients

Examples
Solve the equation
${\displaystyle y'(x)+3y(x)=2}$

with the initial condition

${\displaystyle y(0)=2.}$

Using the general solution method:

${\displaystyle y=e^{-3x}\left(\int 2e^{3x}\,dx+\kappa \right).\,}$

The indefinite integral is solved to give:

${\displaystyle y=e^{-3x}\left(\left(2/3\right)e^{3x}+\kappa \right).\,}$

Then we can reduce to:

${\displaystyle y=2/3+\kappa e^{-3x}.\,}$

where κ = 4/3 from the initial condition.

A linear ODE of order 1 with variable coefficients has the general form

${\displaystyle Dy(x)+f(x)y(x)=g(x).}$

Where D is the differential operator. Equations of this form can be solved by multiplying the integrating factor

${\displaystyle e^{\int f(x)\,dx}}$

throughout to obtain

${\displaystyle Dy(x)e^{\int f(x)\,dx}+f(x)y(x)e^{\int f(x)\,dx}=g(x)e^{\int f(x)\,dx},}$

which simplifies due to the product rule (applied backwards) to

${\displaystyle D\left(y(x)e^{\int f(x)\,dx}\right)=g(x)e^{\int f(x)\,dx}}$

which, on integrating both sides and solving for y(x) gives:

${\displaystyle y(x)=e^{-\int f(x)\,dx}\left(\int g(x)e^{\int f(x)\,dx}\,dx+\kappa \right).}$

In other words: The solution of a first-order linear ODE

${\displaystyle y'(x)+f(x)y(x)=g(x),}$

with coefficients that may or may not vary with x, is:

${\displaystyle y=e^{-a(x)}\left(\int g(x)e^{a(x)}\,dx+\kappa \right)}$

where κ is the constant of integration, and

${\displaystyle a(x)=\int {f(x)\,dx}.}$

A compact form of the general solution based on a Green's function is[4]

${\displaystyle y(x)=\int _{a}^{x}\!{[y(a)\delta (t-a)+g(t)]e^{-\int _{t}^{x}\!f(u)du}\,dt}\,.}$

where δ(x) is the generalized Dirac delta function.

### Examples

Consider a first order differential equation with constant coefficients:

${\displaystyle {\frac {dy}{dx}}+by=1.}$

This equation is particularly relevant to first order systems such as RC circuits and mass-damper systems.

In this case, f(x) = b, g(x) = 1.

Hence its solution is

${\displaystyle y(x)=e^{-bx}\left({\frac {e^{bx}}{b}}+C\right)={\frac {1}{b}}+Ce^{-bx}.}$

## Systems of linear differential equations

An arbitrary linear ordinary differential equation or even a system of such equations can be converted into a first order system of linear differential equations by adding variables for all but the highest order derivatives. A linear system can be viewed as a single equation with a vector-valued variable. The general treatment is analogous to the treatment above of ordinary first order linear differential equations, but with complications stemming from noncommutativity of matrix multiplication.

To solve

${\displaystyle \left\{{\begin{array}{rl}\mathbf {y} '(x)&=A(x)\mathbf {y} (x)+\mathbf {b} (x)\\\mathbf {y} (x_{0})&=\mathbf {y} _{0}\end{array}}\right.}$

(here ${\displaystyle \mathbf {y} (x)}$ is a vector or matrix, and ${\displaystyle A(x)}$ is a matrix), let ${\displaystyle U(x)}$ be the solution of ${\displaystyle \mathbf {y} '(x)=A(x)\mathbf {y} (x)}$ with ${\displaystyle U(x_{0})=I}$ (the identity matrix). ${\displaystyle U}$ is a fundamental matrix for the equation — the columns of ${\displaystyle U}$ form a complete linearly independent set of solutions for the homogeneous equation. After substituting ${\displaystyle \mathbf {y} (x)=U(x)\mathbf {z} (x)}$, the equation ${\displaystyle \mathbf {y} '(x)=A(x)\mathbf {y} (x)+\mathbf {b} (x)}$ simplifies to ${\displaystyle U(x)\mathbf {z} '(x)=\mathbf {b} (x).}$ Thus,

${\displaystyle \mathbf {y} (x)=U(x)\mathbf {y_{0}} +U(x)\int _{x_{0}}^{x}U^{-1}(t)\mathbf {b} (t)\,dt}$

If ${\displaystyle A(x_{1})}$ commutes with ${\displaystyle A(x_{2})}$ for all ${\displaystyle x_{1}}$ and ${\displaystyle x_{2}}$, then[citation needed]

${\displaystyle U(x)=e^{\int _{x_{0}}^{x}A(x)\,dx}}$

and thus

${\displaystyle U^{-1}(x)=e^{-\int _{x_{0}}^{x}A(x)\,dx},}$

but in the general case there is no closed form solution, and an approximation method such as Magnus expansion may have to be used. Note that the exponentials are matrix exponentials.