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

There is no formula that solves every differential equation. The practical route is to recognise the class of the equation and apply that class's method. For first-order ODEs there are two standard classes, separable and linear, with the constant-coefficient case of the latter solved fastest by an exponential trial (Method 3), the same idea that grows into the characteristic equation of order two; for second-order ODEs there are direct integration and constant-coefficient linear equations; for linear PDEs there are separation of variables and the method of characteristics (travelling waves). Each linear method below follows the structure derived above: solve the homogeneous equation (its solutions carry the constants), then add one particular solution of the non-homogeneous one. Each method is stated for its general case and then applied to a concrete numerical example. When an equation fits none of these classes, it is treated by the numerical, series, or qualitative methods summarised at the end of this section.

First-order ODEs

Method 1: separable equations

General case. A first-order equation is separable when the right-hand side factors into a function of $x$ alone times a function of $y$ alone:

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

Divide both sides by $h(y)$ and integrate: all the $y$'s land on one side and all the $x$'s on the other, and if the two integrals can be evaluated the resulting relation between $y$ and $x$ is the general solution.

Example: exponential growth and decay. When a quantity changes at a rate proportional to its own size — a bank balance earning interest, a population with unlimited food, a radioactive sample — the equation is separable with $g(t)=k$, $h(y)=y$:

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

Divide by $y$ and integrate:

$$\int\frac{dy}{y}=\int k\,dt\qquad\Longrightarrow\qquad \ln|y|=kt+C$$

Exponentiating, $|y|=e^C e^{kt}$. The sign of $y$ never changes, so absorbing it into the constant and writing $y(0)=y_0$,

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

Numbers: €1000 at 5% interest compounded continuously ($k=0.05\ \text{yr}^{-1}$) gives $y(t)=1000\,e^{0.05t}$, and

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

Doubling time: $1000\,e^{0.05t}=2000\Rightarrow t=\ln 2/0.05\approx 13.9$ years. For $k<0$ the same solution describes decay; the half-life $y=y_0/2$ is $t_{1/2}=(\ln 2)/(-k)$.

Example: Newton's law of cooling. A body hotter than its surroundings cools at a rate proportional to the temperature gap, which makes the equation separable:

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

Separate and integrate:

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

Exponentiating and folding the (constant-sign) factor $T-T_a$ into the constant, with $T(0)=T_0$:

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

The gap $T-T_a$ decays exponentially; $T$ itself does not. Example: a drink at $80\,^{\circ}\mathrm{C}$ in a $20\,^{\circ}\mathrm{C}$ room, $k=0.1\ \text{min}^{-1}$:

$$T(t)=20+60\,e^{-0.1t}$$

Reaches $40\,^{\circ}\mathrm{C}$ when $20+60e^{-0.1t}=40$, i.e. $t=10\ln 3\approx 11$ min.

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

General case. Many first-order equations are linear but not separable. A linear first-order equation has the form

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

with $q(x)$ the forcing. The homogeneous case $q=0$ is the separable equation of Method 1 with $h(y)=y$, and its solutions are

$$y_h=C\,e^{-\int p(x)\,dx}$$

For the non-homogeneous case multiply both sides by the integrating factor

$$\mu(x)=e^{\int p(x)\,dx}$$

so that the left-hand side collapses into a single derivative:

\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)

Both sides then integrate directly, giving the general solution

$$y(x)=\frac{1}{\mu(x)}\left(\int \mu(x)\,q(x)\,dx + C\right)$$

which is the homogeneous solution $y_h=C/\mu$ plus the particular solution $y_p=(1/\mu)\int\mu q\,dx$, the structure $y=y_p+y_h$ of the classification section.

Worked demonstration: solve $y'+y=e^{-x}$. Here $p=1$, $\mu=e^x$, and

$$e^x y'+e^x y=1\qquad\Longrightarrow\qquad (e^x y)'=1\qquad\Longrightarrow\qquad e^x y=x+C\qquad\Longrightarrow\qquad y=(x+C)e^{-x}$$

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

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

Example: falling with air resistance. A falling body of mass $m$ is pulled down by gravity $mg$ and slowed by air drag proportional to its speed, $-bv$. Newton's second law gives the linear first-order equation

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

Here $p=b/m$ and $q=g$, both constant, so $\mu=e^{(b/m)t}$ and

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

hence, with $v(0)=0$,

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

As $t$ grows the exponential fades and the speed approaches the constant terminal velocity $mg/b$; in the language of the classification section the terminal velocity is the particular solution and the fading exponential the homogeneous part. Numbers: a skydiver of $m=70\ \mathrm{kg}$ with $b=14\ \mathrm{kg\,s^{-1}}$ has $mg/b = 70\times9.8/14 = 49\ \mathrm{m\,s^{-1}}$ and $b/m = 0.2\ \mathrm{s^{-1}}$, so

$$v(t)=49\left(1-e^{-0.2t}\right)\ \mathrm{m\,s^{-1}}$$

After 5 seconds $v=49(1-e^{-1})\approx 31\ \mathrm{m\,s^{-1}}$; after 10 seconds $v=49(1-e^{-2})\approx 42\ \mathrm{m\,s^{-1}}$; the terminal $49\ \mathrm{m\,s^{-1}}$ is approached but never quite reached.

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

General case. When the coefficient $p(x)$ is the constant $a$, the linear first-order equation has constant coefficients, and the homogeneous and non-homogeneous cases are solved separately and then added, as in the classification section.

Homogeneous case. The equation $y'+a\,y=0$ is solved by an exponential trial. Try $y=C e^{bx}$: substituting gives $(b+a)\,C e^{bx}=0$, and since the exponential never vanishes the only possible exponent is $b=-a$, so the homogeneous solution is

$$y_h=C e^{-ax}$$

Exponential growth $y'=k\,y$ is the special case $a=-k$, giving $y_h=C e^{kt}$ with no integration at all.

Non-homogeneous case. For the forced equation $y'+a\,y=q(x)$, the general solution is the particular-plus-homogeneous sum

$$y=y_p+y_h=C e^{-ax}+y_p$$

because substituting $y_p+y_h$ leaves $q(x)+0=q(x)$: the homogeneous part already vanishes on its own, so it can be added to any particular solution without spoiling it. When $q(x)$ is a constant, an exponential, a sine or cosine, or a polynomial, a particular solution of the same form can be guessed and its coefficient fixed by substitution — this is the method of undetermined coefficients. A trial that duplicates the homogeneous solution is multiplied by $x$ instead. The same exponential trial reappears, as the characteristic equation, for the constant-coefficient equations of order two below.

Example: an account with steady withdrawals. An account earning 10% interest compounded continuously, from which €100 is withdrawn each year, is forced exponential growth:

$$y'=0.1\,y-100$$

The homogeneous part, found by the trial $y=C e^{bt}$, is $y_h=C e^{0.1t}$. The withdrawal term is constant, so guess the constant particular solution $y_p=A$; substituting gives $0.1A-100=0$, hence $A=1000$. Starting with €5000, the condition $y(0)=1000+C=5000$ gives $C=4000$ (the constant lives, as always, in the homogeneous part), and

$$y(t)=1000+4000\,e^{0.1t}$$

Check: $y'-0.1y=400e^{0.1t}-(100+400e^{0.1t})=-100$, as required. The €1000 is the balance whose annual interest (10% of €1000) exactly offsets the withdrawals: a balance above €1000 grows, and one below it shrinks. After ten years

$$y(10)=1000+4000\,e\approx 11\,873$$

whereas without the withdrawals the €5000 would have grown to $5000e\approx 13\,591$.

Textbooks add further first-order classes — equations homogeneous in the sense $y'=f(y/x)$, Bernoulli equations, and exact equations — each solvable by an extra change of variables or, for exact equations, by recognising a total differential.

Second-order ODEs

Method 1: direct integration

General case. When the equation has the form $y''=f(x)$, with the right-hand side depending only on the independent variable, each derivative is undone by one integration:

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

The two constants are fixed by two conditions, typically the initial value and the initial derivative (see the discussion of order above). The same pattern applies to $y^{(n)}=f(x)$ with $n$ integrations and $n$ constants.

Example: free fall. With only gravity acting, the height $x(t)$ of a falling object obeys $x''=-g$. Integrating twice,

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

where $v_0$ and $x_0$ are the speed and height at $t=0$. A ball dropped from rest ($v_0=0$) at height $19.6\ \mathrm{m}$ reaches the ground, $x=0$, when

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

The two initial conditions have fixed the whole trajectory.

Method 2: linear equations with constant coefficients

General case. Linear equations with constant coefficients,

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

model a mass on a spring, a small-angle pendulum, and an RLC circuit. As with every linear equation, the homogeneous and non-homogeneous cases are solved separately and joined by superposition.

Homogeneous case ($f=0$). Try the exponential $y=e^{rx}$ (the trial of Method 3), since $y'=re^{rx}$ and $y''=r^2e^{rx}$, giving

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

The factor $e^{rx}$ is never zero, so the exponential solves the equation exactly when $r$ solves the algebraic characteristic equation

$$r^2+ar+b=0$$

whose roots determine the homogeneous solution:

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

The two arbitrary constants, fixed by the initial conditions, live here and only here.

Non-homogeneous case ($f\neq 0$). By the classification result, the general solution is

$$y=y_p+y_h$$

with $y_h$ from above and $y_p$ any single solution of the full equation: substituting $y_p+y_h$ gives $f(x)+0$, since the homogeneous part vanishes on its own. When $f(x)$ is a constant, an exponential, a sine or cosine, or a polynomial, $y_p$ is found by undetermined coefficients, exactly as in Method 3, and a trial duplicating $y_h$ is multiplied by $x$.

Worked demonstration (homogeneous): solve $y''-3y'+2y=0$. The characteristic equation $r^2-3r+2=(r-1)(r-2)=0$ has the two distinct real roots $1$ and $2$, so

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

Check: for $y=e^x$, $y''-3y'+2y=(1-3+2)e^x=0$.

Worked demonstration (non-homogeneous): the same operator forced by $2e^{3x}$, i.e. $y''-3y'+2y=2e^{3x}$, keeps the homogeneous solution found above, $y_h=C_1e^x+C_2e^{2x}$, and adds a trial $y_p=Ae^{3x}$. Substituting,

$$y_p''-3y_p'+2y_p=(9-9+2)Ae^{3x}=2Ae^{3x}=2e^{3x}$$

so $A=1$ and

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

Check: substituting $y=e^{3x}$ gives $(9-9+2)e^{3x}=2e^{3x}$; the homogeneous terms $e^x,e^{2x}$ vanish by the earlier check, so the sum solves the forced equation for any $C_1,C_2$.

Example: the harmonic oscillator (a mass on a spring). A mass displaced $x$ from rest is pulled back by $-kx$ (Hooke's law), so Newton's second law gives

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

Writing $\omega_0^2=k/m$ gives the harmonic oscillator equation

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

whose characteristic equation $r^2+\omega_0^2=0$ has the purely imaginary roots $r=\pm i\omega_0$, the complex-pair case with $\alpha=0$. Since $\frac{d^2}{dt^2}\cos\omega_0 t=-\omega_0^2\cos\omega_0 t$, and likewise for sine, superposition gives

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

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

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

Numbers: $m=2\ \mathrm{kg}$, $k=8\ \mathrm{N/m}$, so $\omega_0=\sqrt{8/2}=2\ \text{rad/s}$. Pulled $10\ \mathrm{cm}$ out and released from rest, $B=0$ and $x(t)=0.10\cos 2t$ metres. The period is

$$P=\frac{2\pi}{\omega_0}=\pi\approx 3.14\ \text{s}$$

and after one second

$$x(1)=0.10\cos 2\approx 0.10(-0.416)\approx -0.042\ \text{m}$$

Such fixed-amplitude sinusoidal motion is simple harmonic motion.

Partial differential equations

Two elementary exact methods cover the standard introductory cases: separation of variables, for diffusion problems on finite domains, and travelling waves (the method of characteristics), for transport. Each is stated in general and then applied to a concrete numerical example.

Method 1: separation of variables (the heat equation)

General case. For a linear, homogeneous PDE on a simple domain, look for a solution that is a product of functions of the separate independent variables, $u(x,t)=X(x)T(t)$. Substituting splits the PDE into two linked ordinary equations; boundary conditions pick out which solutions survive; and superposition of those basic solutions then gives the general solution.

Application. Solve the heat equation on a bar of length $L$ with insulated sides and both ends held at $0$:

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

Seek $u(x,t)=X(x)T(t)$. Substituting gives $XT'=\alpha X''T$; dividing by $\alpha XT$,

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

The left side depends only on $t$ and the right only on $x$, so both must equal one and the same constant, written $-\lambda$. Each side is now an ODE:

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

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

The end conditions $u(0,t)=u(L,t)=0$ force $X(0)=X(L)=0$: hence $A=0$ and $\sin(\sqrt\lambda\,L)=0$, so $\sqrt\lambda\,L=n\pi$, $n=1,2,\dots$ Each allowed value $\lambda=(n\pi/L)^2$ gives one basic solution, a mode

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

The equation is linear and homogeneous, so superposition applies and the general solution is

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 the $b_n$ fixed by the initial profile $u(x,0)$ (a Fourier sine series). The decay rate $\alpha(n\pi/L)^2$ grows as $n^2$, so higher modes die out first.

Numbers: a $1\ \mathrm{m}$ iron bar, $\alpha\approx 2.3\times10^{-5}\ \mathrm{m^2s^{-1}}$, heated so $u(x,0)=100\sin(\pi x/L)$ (ends at $0\,^{\circ}\mathrm{C}$, centre $100\,^{\circ}\mathrm{C}$). Only the $n=1$ mode is present:

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

At the centre, with $L=1$ and $\alpha\pi^2\approx 2.3\times10^{-4}\ \text{s}^{-1}$,

$$u\!\left(\tfrac12,t\right)=100\,e^{-2.3\times10^{-4}t}$$

so after one hour $u\approx 100e^{-0.82}\approx 44\,^{\circ}\mathrm{C}$, and $50\,^{\circ}\mathrm{C}$ is reached at $t=\ln 2/(2.3\times10^{-4})\approx 3050\ \text{s}\approx 51$ min.

Method 2: travelling waves (the method of characteristics)

General case. Some first-order PDEs are solved not by separating variables but by noticing that their solutions travel. The transport equation with constant speed $c$,

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

states that $u$ is carried along unchanged. Trying a travelling wave $u(x,t)=f(x-ct)$, with $f$ arbitrary, the chain rule gives $u_t=-c\,f'$ and $u_x=f'$, so

$$u_t+c\,u_x=-c\,f'+c\,f'=0$$

whatever $f$ is. Hence the general solution is any profile sliding rigidly to the right at speed $c$ (to the left if $c<0$):

$$u(x,t)=f(x-ct)$$

with $f$ fixed by the initial profile, $u(x,0)=f(x)$. Equivalently, $u$ is constant on each of the straight characteristic lines $x-ct=\text{constant}$. The equation is linear and homogeneous, so sums of travelling waves are again solutions. A source term makes it non-homogeneous, $u_t+cu_x=s(x,t)$, and adds, along each characteristic, the accumulated contribution of $s$, the same particular-plus-homogeneous structure as for ODEs:

$$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 flows at $2\ \mathrm{m\,s^{-1}}$, and at $t=0$ a release at $x=0$ gives the Gaussian concentration profile $u(x,0)=50\,e^{-(x/10)^2}\ \mathrm{mg\,L^{-1}}$: $50\ \mathrm{mg\,L^{-1}}$ at the release point, falling by a factor $e^{-1}\approx 0.37$ ten metres away. The solution slides the whole profile downstream without changing it:

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

After one minute the peak, which started at $x=0$, has reached $x=2\times 60=120$ m downstream and still reads $50\ \mathrm{mg\,L^{-1}}$; a monitor $120$ m downstream sees the peak arrive after $60$ s. The transport equation describes pure advection, with no spreading: that is why the heat example above needed the extra second-order term $\alpha u_{xx}$ to spread its pulse out.

Remark: the wave equation. The second-order wave equation $u_{tt}=c^2u_{xx}$ (a plucked string) is solved by the same travelling-wave idea applied in both directions at once; its general solution is

$$u(x,t)=f(x-ct)+g(x+ct)$$

two arbitrary shapes moving left and right (d'Alembert, 1747, described in 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

Differential equations are as old as the calculus itself, because the calculus is the mathematics of change: the laws of physics say how quantities change, and predicting the future means undoing those changes. When Newton published his laws of motion and of gravitation in the Principia in 1687, the equations he needed were differential equations, and he solved them by geometry and by infinite series. Newton had his own notation for rates of change, but it was Leibniz's $dy/dx$, first written in the 1670s, that survived: it displays the whole equation on the page, and it is the notation used throughout this article.

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

The scattered tricks of the early calculus became a subject when Leonhard Euler took them up in the middle decades of the 18th century. He recognised that a linear equation with constant coefficients is solved by substituting $y=e^{rx}$, turning calculus into algebra; he developed series solutions; and, for equations that resisted formulas, he invented the step-by-step numerical scheme, described in this article, that still bears his name. Most of the exact methods above descend from his work.

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

Meanwhile physics began asking for equations with more than one independent variable. A plucked string takes a shape that depends on position along the string and on time, and in 1747 Jean le Rond d'Alembert wrote down the wave equation for it and solved it, showing that its solutions are two waves travelling in opposite directions. Heat conduction posed a subtler problem, because the initial temperature of a bar can have any shape at all. Joseph Fourier derived the heat equation from the physics of conduction and, to solve it, had to express an arbitrary initial profile as a sum of sine ripples. His Théorie analytique de la chaleur of 1822 turned separation of variables into a cornerstone of applied mathematics, and the Fourier series invented for the purpose now appears wherever signals are analysed, from acoustics to image compression.

The exact formulas, however, have their limits, and the history of the subject since the late 19th century is largely the story of what to do when no formula exists. Studying the three-body problem of celestial mechanics, Henri Poincaré realised that the shape of the motion can be understood without solving the equations, founding the qualitative theory of dynamical systems. The electronic computer then made the numerical route routine: approximate the solution step by step, as Euler's method does, refining the steps until the error is acceptable. The two strands met in 1963, when the meteorologist Edward Lorenz found that a simple system of three differential equations, meant to model atmospheric convection, behaved chaotically: the equations were deterministic, yet their solutions were aperiodic and so sensitive to initial conditions that long-term weather prediction is impossible in practice. Each of these later routes, qualitative study, numerical stepping, series, and transforms, is the subject of its own article.

Further reading

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