Differential equation

A differential equation is an equation whose unknown is a function and which also involves that function's derivatives (rates of change). Where an ordinary equation such as $x^2=9$ is solved by numbers, a differential equation such as $y'+2y=0$ is solved by functions $y(x)$. Laws of nature state how quantities change, so differential equations describe pendulums, cooling drinks, growing populations, and discharging capacitors. This article covers the standard analytical solution methods, by class of equation, each stated in general and then demonstrated on a concrete numerical example, and the numerical, series, and qualitative routes used when no exact formula exists. It treats ordinary differential equations (one independent variable) and, briefly, partial differential equations (several).

A first example: slopes and a family of solutions

The simplest differential equation prescribes the slope of a function $y(x)$:

$$\frac{dy}{dx}=2x$$

Integration inverts differentiation, so integrating both sides gives

$$\int\frac{dy}{dx}\,dx=\int 2x\,dx\qquad\Longrightarrow\qquad y(x)=x^{2}+C$$

Every $C$ works, since $\frac{d}{dx}\left(x^2+C\right)=2x$; the solutions form the parabola family $y=x^2+C$, the general solution.

If $y(0)=3$, then

$$3=0^2+C\qquad\Longrightarrow\qquad C=3\qquad\Longrightarrow\qquad y=x^2+3$$

A prescribed value such as this is an initial condition.

Classifying differential equations

Three features decide how to solve an equation: its order, its linearity, and how many independent variables it involves.

Order

The order is the order of the highest derivative present. $dy/dx=2x$ is first order; Newton's second law,

m\frac{d^{2}x}{dt^{2}}=F

is second order ($x(t)$ position of mass $m$, $F$ net force). Integration introduces one arbitrary constant per integration, so the general solution of an nth-order equation carries $n$ constants, fixed by $n$ initial conditions. For equations of the special form $y^{(n)}=f(x)$ the constants appear exactly as the integration constants of $n$ successive integrations; the free-fall example in the second-order section below works this out for $n=2$.

Linearity and homogeneity

An equation is linear when the unknown and its derivatives appear only to the first power and never multiplied together. A linear first-order equation can always be written

$$\frac{dy}{dx}+p(x)\,y=q(x)$$

and is homogeneous when $q(x)=0$. The equations $dy/dx=y^2$ and $d^2\theta/dt^2+\sin\theta=0$ are nonlinear (square of $y$; sine of $\theta$).

If $y_1,y_2$ solve a homogeneous linear equation, so does $c_1y_1+c_2y_2$ (the superposition principle): substituting the combination adds the two expressions that already vanish. For a nonlinear equation the combination does not generally solve it: if $y_1'=y_1^2$ and $y_2'=y_2^2$, then

$$(y_1+y_2)'=y_1^2+y_2^2\neq (y_1+y_2)^2$$

so $y_1+y_2$ does not solve $y'=y^2$.

Superposition also joins the homogeneous and non-homogeneous problems of one linear equation. Write the left-hand side as $L(y)$, so the equation reads $L(y)=q(x)$, with $L(y)=0$ its homogeneous form. If $y_p$ is any single solution of $L(y)=q$ (a particular solution) and $y_h$ runs through all solutions of $L(y)=0$, then every solution of the original equation is

$$y=y_p+y_h$$

because $L(y_p+y_h)=L(y_p)+L(y_h)=q+0=q$, and conversely any two solutions of the non-homogeneous equation differ by a solution of the homogeneous one. The constants of integration therefore live entirely in $y_h$: the general solution of a linear equation is one particular solution plus the whole homogeneous family. This is why each linear method below is presented in two parts, the homogeneous case first.

Ordinary and partial

An ordinary differential equation (ODE) has one independent variable. A partial differential equation (PDE) has several, with partial derivatives. For example, the temperature $u(x,t)$ of an insulated metal bar, which depends on position $x$ and time $t$, obeys the heat equation

\frac{\partial u}{\partial t}=\alpha\frac{\partial^{2}u}{\partial x^{2}}

where $\alpha$ is the thermal diffusivity.

Slope fields

A first-order equation can be written

\frac{dy}{dx}=f(x,y)

assigning to each point $(x,y)$ the slope $f(x,y)$ a solution must have there. Drawing short segments of that slope gives a direction field; solution curves run tangent to it.

Direction field of $dy/dx=y$. Credit: jjbeard (public domain).

Numerical methods such as Euler's method follow the field: read the slope, step a short distance along it, repeat.

Solving differential equations

Closed-form solutions are known only for restricted classes of equations; the standard practice is to identify the class by order, linearity, and coefficients, and to apply that class's method. The linear methods below follow the two-step structure of the classification section: solve the homogeneous equation, whose general solution carries all arbitrary constants, then add one particular solution of the non-homogeneous equation.

First-order ODEs

Method 1: separable equations

General case. A first-order equation is separable when it can be brought to the separated form

\frac{dy}{dx}=g(x)\,h(y)\qquad\Longrightarrow\qquad\int\frac{dy}{h(y)}=\int g(x)\,dx

after which both integrals are evaluated directly.

Example: exponential growth and decay. For $\dfrac{dy}{dt}=ky$,

\frac{dy}{dt}=k\,y

$$\int\frac{dy}{y}=\int k\,dt\;\Longrightarrow\;\ln|y|=kt+C_1\;\Longrightarrow\;y=Ce^{kt}$$

The initial condition $y(0)=y_0$ fixes $C=y_0$:

y(t)=y_{0}\,e^{kt}

Numerical case: €1000 at 5% interest compounded continuously, $k=0.05\ \text{yr}^{-1}$:

$$y(t)=1000\,e^{0.05t},\qquad y(10)=1000\,e^{0.5}\approx 1648.7$$

$$t_{\text{double}}=\frac{\ln 2}{k}\approx 13.9\ \text{yr},\qquad t_{1/2}=\frac{\ln 2}{-k}\ (k<0)$$

Example: Newton's law of cooling. For $\dfrac{dT}{dt}=-k(T-T_a)$,

\frac{dT}{dt}=-k\bigl(T-T_{a}\bigr)

$$\int\frac{dT}{T-T_a}=-\int k\,dt\;\Longrightarrow\;\ln|T-T_a|=-kt+C\;\Longrightarrow\;T-T_a=Ce^{-kt}$$

with $C=T_0-T_a$ from $T(0)=T_0$:

$$T(t)=T_a+(T_0-T_a)e^{-kt}$$

Numerical case: a drink at $T_0=80\,^{\circ}\mathrm{C}$ in a room at $T_a=20\,^{\circ}\mathrm{C}$, $k=0.1\ \text{min}^{-1}$:

$$T(t)=20+60e^{-0.1t},\qquad T=40\,^{\circ}\mathrm{C}\text{ at }t=10\ln 3\approx 11\ \text{min}.$$

Method 2: linear first-order equations (integrating factor)

General case. For the linear equation

$$y'+p(x)\,y=q(x)$$

introduce the integrating factor $\mu=e^{\int p\,dx}$, for which $\mu'=p\mu$. The product rule then collapses the left-hand side:

\frac{dy}{dx}+p(x)\,y=q(x),\text{multiply both sides by } \mu(x)=e^{\int p(x)\,dx}\;\Longrightarrow\;\frac{d}{dx}\bigl(\mu(x)\,y\bigr)=\mu(x)\,q(x)

Integrating both sides,

$$\mu\,y=\int\mu\,q\,dx+C\;\Longrightarrow\;y=\frac{1}{\mu}\int\mu\,q\,dx+\frac{C}{\mu}$$

The first term is a particular solution of the non-homogeneous equation; the second, $C/\mu=Ce^{-\int p\,dx}$, is the general solution of the homogeneous equation $y'+py=0$.

Worked demonstration. $y'+y=e^{-x}$: $p=1$, $\mu=e^x$, and $(e^x y)'=e^x(y'+y)=1$, so $e^x y=x+C$:

\frac{dy}{dx}+y=e^{-x}\qquad\Longrightarrow\qquad y=(x+C)e^{-x}

$y(0)=2$ gives $C=2$.

Example: falling with air resistance. Newton's second law with drag $-bv$ gives

$$m\frac{dv}{dt}=mg-bv\;\Longrightarrow\;v'+\frac{b}{m}v=g$$

With $\mu=e^{(b/m)t}$,

$$\frac{d}{dt}\left(e^{(b/m)t}v\right)=g\,e^{(b/m)t}\;\Longrightarrow\;v=\frac{mg}{b}+Ce^{-(b/m)t}$$

$v(0)=0$ fixes $C=-mg/b$:

$$v(t)=\frac{mg}{b}\left(1-e^{-(b/m)t}\right)$$

Numerical case: $m=70\ \mathrm{kg}$, $b=14\ \mathrm{kg\,s^{-1}}$, so $mg/b=49\ \mathrm{m\,s^{-1}}$ (terminal velocity) and $b/m=0.2\ \mathrm{s^{-1}}$:

$$v(t)=49\left(1-e^{-0.2t}\right),\qquad v(5)\approx 31,\quad v(10)\approx 42\ \mathrm{m\,s^{-1}}$$

Method 3: constant-coefficient linear equations (trial solutions)

General case. For $y'+ay=q(x)$, the homogeneous equation is solved by the exponential trial $y=Ce^{bx}$:

$$(b+a)Ce^{bx}=0\;\Longrightarrow\;b=-a\;\Longrightarrow\;y_h=Ce^{-ax}$$

Growth $y'=ky$ is the case $a=-k$. The forced equation then has, by linearity,

$$y=y_h+y_p,\qquad y_p\ \text{any solution of }y'+ay=q$$

When $q$ is constant, exponential, sinusoidal, or polynomial, $y_p$ is guessed in the same family and its coefficient fixed by substitution (method of undetermined coefficients); a guess satisfying the homogeneous equation is multiplied by $x$.

Example: an account with steady withdrawals. $y'=0.1y-100$, $y(0)=5000$:

$$y_h=Ce^{0.1t},\qquad y_p=A:\ 0.1A-100=0\;\Longrightarrow\;A=1000$$

$$y(0)=1000+C=5000\;\Longrightarrow\;C=4000\;\Longrightarrow\;y(t)=1000+4000e^{0.1t}$$

Check: $y'-0.1y=400e^{0.1t}-(100+400e^{0.1t})=-100$. Evaluation:

$$y(10)=1000+4000e\approx 11\,873;\qquad \text{without withdrawals: }5000e\approx 13\,591$$

Further first-order classes, $y'=f(y/x)$, Bernoulli, exact, reduce to these by substitution or by recognising a total differential.

Second-order ODEs

Method 1: direct integration

General case. For $y''=f(x)$,

$$y''=f(x)\;\Longrightarrow\;y'=\int f(x)\,dx+C_1\;\Longrightarrow\;y=\int\!\!\left(\int f(x)\,dx\right)dx+C_1x+C_2$$

and likewise $y^{(n)}=f(x)$ by $n$ integrations.

Example: free fall. $x''=-g$:

$$\frac{dx}{dt}=-gt+v_0\;\Longrightarrow\;x(t)=-\tfrac{g}{2}t^2+v_0t+x_0$$

Dropped from rest at $19.6\ \mathrm{m}$ ($v_0=0$, $x_0=19.6$):

$$0=19.6-4.9t^2\;\Longrightarrow\;t=\sqrt{19.6/4.9}=2\ \text{s}$$

Method 2: linear equations with constant coefficients

General case.

$$y''+a\,y'+b\,y=f(x)$$

Homogeneous case ($f=0$). The exponential trial $y=e^{rx}$,

y''+a\,y'+b\,y=0,\qquad y=e^{rx}\ \Rightarrow\ r^{2}+a\,r+b=0

gives the characteristic equation $r^2+ar+b=0$, whose roots determine $y_h$:

• $r_1\neq r_2$ real: $y_h=C_1e^{r_1x}+C_2e^{r_2x}$;
• $r_1=r_2=r$: $y_h=(C_1+C_2x)e^{rx}$;
• $r=\alpha\pm i\beta$: $y_h=e^{\alpha x}(C_1\cos\beta x+C_2\sin\beta x)$.

Non-homogeneous case ($f\neq 0$). $y=y_h+y_p$, with $y_p$ found by undetermined coefficients as in Method 3.

Worked demonstration (homogeneous). $y''-3y'+2y=0$: $r^2-3r+2=(r-1)(r-2)=0$,

y''-3y'+2y=0\qquad\Longrightarrow\qquad y=C_{1}e^{x}+C_{2}e^{2x}

Check: $e^x$ gives $(1-3+2)e^x=0$.

Worked demonstration (non-homogeneous). $y''-3y'+2y=2e^{3x}$: keep $y_h$ above, try $y_p=Ae^{3x}$:

$$y_p''-3y_p'+2y_p=(9-9+2)Ae^{3x}=2Ae^{3x}\;\Longrightarrow\;A=1$$

$$y=C_1e^x+C_2e^{2x}+e^{3x}$$

Example: the harmonic oscillator (a mass on a spring). Hooke's law $F=-kx$ in Newton's second law,

$$m\frac{d^2x}{dt^2}=-kx\;\Longrightarrow\;x''+\omega_0^2x=0,\qquad \omega_0^2=\frac{k}{m}$$

\frac{d^{2}x}{dt^{2}}+\omega_{0}^{2}x=0

$r^2+\omega_0^2=0$ gives $r=\pm i\omega_0$, hence

$$x(t)=A\cos\omega_0t+B\sin\omega_0t$$

with $A,B$ fixed by initial position and velocity.

A mass on a spring: the harmonic oscillator solution is a sinusoid of fixed amplitude. Credit: Evil saltine (public domain).

Numerical case: $m=2\ \mathrm{kg}$, $k=8\ \mathrm{N\,m^{-1}}$: $\omega_0=2\ \text{rad\,s}^{-1}$; released from rest at $10\ \mathrm{cm}$,

$$x(t)=0.10\cos 2t\ \mathrm{m},\qquad P=\frac{2\pi}{\omega_0}=\pi\approx 3.14\ \text{s},\qquad x(1)\approx -0.042\ \mathrm{m}$$

Partial differential equations

Two elementary classes of linear PDE admit closed-form solutions: diffusion on finite domains, by separation of variables, and first-order transport, by travelling waves.

Method 1: separation of variables (the heat equation)

General case. For a linear, homogeneous PDE on a simple domain, assume $u(x,t)=X(x)T(t)$; substitution splits the PDE into ordinary equations for $X$ and $T$, boundary conditions select the admissible solutions, and their superposition matches the initial profile.

Application. The heat equation on a bar of length $L$ with ends held at $0$,

$$\frac{\partial u}{\partial t}=\alpha\frac{\partial^2u}{\partial x^2},\qquad u(0,t)=u(L,t)=0$$

$$u=X(x)T(t):\qquad XT'=\alpha X''T\;\Longrightarrow\;\frac{X''}{X}=\frac{T'}{\alpha T}$$

u=X(x)\,T(t)\ \Rightarrow\ \frac{X''}{X}=\frac{T'}{\alpha\,T}=-\lambda

Both sides depend on different variables, hence equal a constant $-\lambda$:

$$T'=-\alpha\lambda T\;\Longrightarrow\;T=e^{-\alpha\lambda t}$$

$$X''=-\lambda X\;\Longrightarrow\;X=A\cos(\sqrt\lambda\,x)+B\sin(\sqrt\lambda\,x)$$

The boundary conditions force $X(0)=X(L)=0$: $A=0$ and $\sin(\sqrt\lambda L)=0$, so $\sqrt\lambda\,L=n\pi$, $n=1,2,\dots$; each $\lambda=(n\pi/L)^2$ gives one mode

$$u_n(x,t)=\sin\frac{n\pi x}{L}\,e^{-\alpha(n\pi/L)^2t}$$

and superposition gives the general solution

u(x,t)=\sum_{n=1}^{\infty}b_{n}\sin\Bigl(\frac{n\pi x}{L}\Bigr)\,e^{-\alpha (n\pi/L)^{2}t}

with $b_n$ determined by the Fourier sine series of the initial profile $u(x,0)$; the decay rate $\alpha(n\pi/L)^2$ grows as $n^2$.

Numerical case. A $1\ \mathrm{m}$ iron bar, $\alpha\approx 2.3\times10^{-5}\ \mathrm{m^2s^{-1}}$, initial profile $u(x,0)=100\sin(\pi x/L)$:

$$u(x,t)=100\sin\frac{\pi x}{L}\,e^{-\alpha\pi^2t/L^2},\qquad u\!\left(\tfrac12,t\right)=100\,e^{-2.3\times10^{-4}t}$$

$$u(\tfrac12,1\ \text{h})\approx 44\,^{\circ}\mathrm{C},\qquad 50\,^{\circ}\mathrm{C}\text{ at }t=\frac{\ln 2}{2.3\times10^{-4}}\approx 51\ \text{min}$$

Method 2: travelling waves (method of characteristics)

General case. The transport equation

$$u_t+c\,u_x=0$$

states that $u$ is carried unchanged. The travelling-wave trial $u=f(x-ct)$ gives $u_t=-cf'$, $u_x=f'$, hence $u_t+cu_x=0$ identically:

$$u(x,t)=f(x-ct),\qquad u(x,0)=f(x)$$

so $u$ is constant on the characteristic lines $x-ct=\text{const}$. The non-homogeneous equation $u_t+cu_x=s(x,t)$ accumulates the source along each characteristic:

$$u(x,t)=f(x-ct)+\int_0^t s\bigl(x-c(t-\tau),\tau\bigr)\,d\tau$$

Example: a slug of pollutant in a river. A river at $c=2\ \mathrm{m\,s^{-1}}$ carries a Gaussian release of peak $50\ \mathrm{mg\,L^{-1}}$ and width such that $u$ falls by $e^{-1}$ at $\pm10\ \mathrm{m}$:

$$u(x,0)=50\,e^{-(x/10)^2}\;\Longrightarrow\;u(x,t)=50\,e^{-((x-2t)/10)^2}\ \mathrm{mg\,L^{-1}}$$

After $t=60\ \mathrm{s}$ the peak is at $x=ct=120\ \mathrm{m}$ with the original value $50\ \mathrm{mg\,L^{-1}}$; pure transport does not spread the pulse, which requires the second-order term $\alpha u_{xx}$ of the heat equation.

Wave equation. The second-order wave equation is the two-directional travelling-wave problem:

$$u_{tt}=c^2u_{xx}\;\Longrightarrow\;u(x,t)=f(x-ct)+g(x+ct)$$

(d'Alembert, 1747; see the history section).

When no formula exists

Most equations, especially nonlinear ones, fit none of the classes above and have no solution in terms of familiar functions. They are studied in one of three ways:

• Numerically, when numbers suffice: Euler's method steps along the slope field;
• Qualitative methods: equilibria, stability and long-term behaviour, without formulas;
• Series and transforms, for linear cases: Power series or the Laplace transform.

The exact methods occupy the branches on the left; most equations encountered in research fall through to the routes on the right, each treated in its own article.

A short history

The origins of differential equations coincide with those of the calculus, since the calculus supplies the language in which rates of change are expressed and inverted. Newton's laws of motion and of universal gravitation, published in the Philosophiae Naturalis Principia Mathematica (1687), are differential equations; Newton treated them by the geometrical and infinite-series methods of his fluxional calculus. Although Newton developed a notation for fluxions, the differential notation $dy/dx$ introduced by Leibniz in the 1670s proved the more enduring: it exhibits the structure of the equation directly and is the notation adopted in this article.

Isaac Newton (portrait after Godfrey Kneller, 1689). Newton's laws of motion and of gravitation (Principia, 1687) are differential equations. Credit: James Thronill after Godfrey Kneller (public domain).

The consolidation of these techniques into a systematic theory is due in large measure to Leonhard Euler, whose work in the middle decades of the eighteenth century established the principal exact methods. Euler showed that linear equations with constant coefficients are solved by the substitution $y=e^{rx}$, which reduces the problem to an algebraic equation, and he advanced the theory of series solutions. For equations that admitted no closed-form solution, he introduced the step-by-step numerical procedure, described above as Euler's method, that bears his name. The exact methods presented in this article derive, in large part, from his work.

Leonhard Euler (portrait by Jakob Emanuel Handmann, 1753). Credit: Jakob Emanuel Handmann (public domain).

The theory of partial differential equations arose from the demands of eighteenth-century physics. In 1747, Jean le Rond d'Alembert derived the wave equation for the vibrating string and established that its general solution consists of two waves propagating in opposite directions. The problem of heat conduction proved more demanding, because the initial temperature distribution of a conducting body is arbitrary. In his Théorie analytique de la chaleur (1822), Joseph Fourier derived the heat equation from the physical principles of conduction and solved it by expanding the initial data into a trigonometric series. This work established separation of variables as a standard technique of mathematical physics, and the Fourier series introduced for the purpose has since become fundamental to the analysis of periodic phenomena, from acoustics to signal processing.

The limits of closed-form methods became apparent towards the end of the nineteenth century, and the later history of the subject is concerned principally with equations for which elementary solutions do not exist. In his investigation of the three-body problem of celestial mechanics, Henri Poincaré demonstrated that qualitative properties of the motion, such as its equilibria, stability, and long-term behaviour, can be characterised without solving the equations, thereby founding the qualitative theory of dynamical systems. The subsequent development of electronic computing made numerical approximation, of which Euler's method is the simplest instance, a routine and general technique. The two strands converged in 1963, when Edward Lorenz, studying a simplified system of three ordinary differential equations that models atmospheric convection, established the phenomenon of deterministic chaos: although the equations are deterministic, their solutions are aperiodic and depend so sensitively on initial conditions that long-term weather prediction is not feasible in practice. These later approaches, qualitative analysis, numerical approximation, and series and transform methods, are treated in dedicated articles.

Further reading

• Differential equation, Wikipedia
• Differential Equation, Wolfram MathWorld
• 18.03SC Differential Equations, MIT OpenCourseWare
• Leonhard Euler, MacTutor History of Mathematics