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. Each method below is stated in general and then applied to a concrete scenario: the scenario is first translated into a differential equation, which the method then solves step by step. The linear methods 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. A quantity whose rate of change is proportional to its own size, such as an unchecked population or a radioactive sample, obeys

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

Separating variables and integrating,

$$\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$, giving

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

If $k>0$ the quantity grows, with doubling time $t_{\text{double}}=(\ln 2)/k$; if $k<0$ it decays, with half-life $t_{1/2}=(\ln 2)/(-k)$. Radioactive decay is the standard decay instance: writing $k=-\lambda$ with $\lambda>0$ gives $y(t)=y_0e^{-\lambda t}$.

Numeric scenario: carbon-14, the isotope used in radiocarbon dating, has $\lambda=1.21\times10^{-4}\ \mathrm{yr}^{-1}$, i.e. a half-life of 5730 years. A sample that retains 60% of its initial ${}^{14}\mathrm{C}$ satisfies $e^{-\lambda t}=0.6$, so its age is

$$t=\frac{1}{\lambda}\ln\frac{1}{0.6}\approx 4.2\times10^{3}\ \mathrm{yr}$$

Example: Newton's law of cooling. A hot object in a cooler room loses heat through its surface, and the larger the temperature gap, the faster it cools: the gap shrinks at a rate proportional to itself. The temperature is modelled by

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

Separating variables and integrating,

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

so with $T(0)=T_0$, hence $C=T_0-T_a$,

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

Numeric scenario: a drink at $T_0=80\,^{\circ}\mathrm{C}$ cools in a room at $T_a=20\,^{\circ}\mathrm{C}$ with rate constant $k=0.1\ \mathrm{min}^{-1}$. It obeys $dT/dt=-0.1(T-20)$, $T(0)=80$, and the solution above becomes

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

The drink reaches $40\,^{\circ}\mathrm{C}$ when $20+60e^{-0.1t}=40$, i.e. at $t=10\ln 3\approx 11$ 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}$, chosen so that $\mu'=p\mu$; multiplying by $\mu$ collapses the left-hand side 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)

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 one particular solution of the non-homogeneous equation, and the second term, $C/\mu=Ce^{-\int p\,dx}$, is the general solution of its homogeneous part, so the formula realises the structure $y=y_p+y_h$.

Worked demonstration. Solve $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}

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

Example: falling with air resistance. A falling body is pulled down by gravity $mg$ and slowed by air drag $-bv$ proportional to its speed; as the speed grows the drag grows too, until it balances gravity and the speed stops increasing. Newton's second law models the fall,

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

a linear first-order equation with constant $p=b/m$ and $q=g$, so $\mu=e^{(b/m)t}$ and

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

Dropped from rest, $v(0)=0$ fixes $C=-mg/b$:

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

Numeric scenario: a skydiver of mass $m=70\ \mathrm{kg}$ and drag coefficient $b=14\ \mathrm{kg\,s^{-1}}$ jumps from rest. The equation of motion is $dv/dt=9.8-0.2v$, $v(0)=0$ (since $mg/b=49\ \mathrm{m\,s^{-1}}$ and $b/m=0.2\ \mathrm{s^{-1}}$), and the solution reads

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

so $v(5)\approx 31\ \mathrm{m\,s^{-1}}$ and $v(10)\approx 42\ \mathrm{m\,s^{-1}}$, approaching but never reaching the terminal velocity $mg/b=49\ \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}$$

Exponential growth $y'=ky$ is the case $a=-k$. By linearity the forced equation has general solution

$$y=y_h+y_p$$

where $y_p$ is any single solution of $y'+ay=q$. When $q$ is constant, exponential, sinusoidal, or polynomial, $y_p$ is tried in the same family and its coefficient fixed by substitution (method of undetermined coefficients); a trial that already solves the homogeneous equation is multiplied by $x$.

Example: a fish stock with constant harvesting. A fish population growing at 10% per year is harvested at the constant rate of 100 fish per year; the population $y(t)$ changes at the rate of natural growth minus the catch:

$$y'=0.1y-100,\qquad y(0)=5000$$

The homogeneous part is solved by the exponential trial $y=Ce^{bt}$: substituting into $y'=0.1y$ gives $b=0.1$, hence $y_h=Ce^{0.1t}$. Since the forcing is constant, try the constant particular solution $y_p=A$, the population at which growth exactly balances the catch; substituting gives $0.1A-100=0$, so $A=1000$, and

$$y(t)=1000+Ce^{0.1t}$$

The initial stock, $y(0)=1000+C=5000$, fixes $C=4000$:

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

Check: $y'-0.1y=400e^{0.1t}-(100+400e^{0.1t})=-100$, so the equation is satisfied. After ten years

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

whereas without harvesting the same initial stock would have grown to $5000e\approx 13\,591$ fish.

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. A ball released above the ground falls under gravity alone, which accelerates it downward at the constant rate $g\approx 9.8\ \mathrm{m\,s^{-2}}$; Newton's second law gives the second-order equation $x''=-g$, of the form above with $f(x)=-g$. Integrating twice,

$$x''=-9.8\;\Longrightarrow\;\frac{dx}{dt}=-9.8t+v_0\;\Longrightarrow\;x(t)=-4.9t^2+v_0t+x_0$$

Numeric scenario: the ball is dropped from rest, $v_0=0$, at height $x_0=19.6\ \mathrm{m}$. It reaches the ground, $x=0$, when

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

so the two initial conditions have pinned down the whole trajectory.

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). A mass attached to a spring is pulled back towards its rest position by a force $-kx$ proportional to the displacement (Hooke's law); once released it oscillates. Newton's second law models the motion,

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

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

The trial $x=e^{rt}$ gives the characteristic equation $r^2+\omega_0^2=0$ with roots $r=\pm i\omega_0$, the complex-pair case with $\alpha=0$, hence

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

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

Numeric scenario: a mass $m=2\ \mathrm{kg}$ hangs on a spring with $k=8\ \mathrm{N\,m^{-1}}$, so $\omega_0=\sqrt{8/2}=2\ \mathrm{rad\,s^{-1}}$ and the displacement obeys $x''+4x=0$. Pulled $0.10\ \mathrm{m}$ from rest and released, the conditions $x(0)=0.10$, $x'(0)=0$ give $A=0.10$, $B=0$:

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

After one second $x(1)=0.10\cos 2\approx -0.042\ \mathrm{m}$, and the motion is simple harmonic motion.

Partial differential equations

Two elementary classes of linear PDE admit closed-form solutions: separable problems on bounded domains, by separation of variables, and transport and wave problems, by travelling waves.

Method 1: separation of variables

General form. For a linear homogeneous PDE in two variables on a bounded domain with homogeneous boundary conditions, seek a solution of the separated form

$$u(x,t)=X(x)\,T(t)$$

Substituting into the PDE and dividing by $XT$ separates the variables into one ordinary differential equation in $x$ and one in $t$. Since the two sides are functions of different variables, they can be identically equal only if each equals the same constant, the separation constant $-\lambda$. The $x$-equation together with the boundary conditions is an eigenvalue problem: only a discrete sequence of constants $\lambda_n$, with eigenfunctions $X_n(x)$, is admissible. The $t$-equation then has a solution $T_n(t)$ for each $n$, and every product $X_nT_n$ solves the PDE.

General algebraic solution. The PDE is linear and homogeneous, so the separated modes superimpose:

$$u(x,t)=\sum_n c_n\,X_n(x)\,T_n(t)$$

with the coefficients $c_n$ chosen so that the series equals the initial profile $u(x,0)$; orthogonality of the eigenfunctions $X_n$ determines them.

Example: the heat equation. The temperature of a bar of length $L$ with insulated sides and both ends held at $0$ obeys

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

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

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

The $t$-equation $T'=-\alpha\lambda T$ has solution $T=e^{-\alpha\lambda t}$, and the $x$-equation

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

together with the boundary conditions forces $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 $\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 the general algebraic solution above becomes

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$.

Numeric scenario: a $1\ \mathrm{m}$ iron bar, heated so that its centre is at $100\,^{\circ}\mathrm{C}$ while both ends are held at $0\,^{\circ}\mathrm{C}$, cools by conduction with iron's diffusivity $\alpha\approx 2.3\times10^{-5}\ \mathrm{m^2s^{-1}}$. The initial profile $u(x,0)=100\sin(\pi x/L)$ is exactly the first mode, so only $n=1$ contributes and

$$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 (method of characteristics)

General form. The transport equation for a quantity $u(x,t)$ carried by a flow at constant speed $c$ is the first-order linear PDE

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

General algebraic solution. For any differentiable profile $f$, the travelling wave $u(x,t)=f(x-ct)$ satisfies the equation identically, since $u_t=-cf'(x-ct)$ and $u_x=f'(x-ct)$; conversely, every solution is of this form, with $f$ fixed by the initial data:

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

so $u$ is constant on each characteristic line $x-ct=\text{const}$, and the initial profile is transported unchanged at speed $c$. For the forced equation $u_t+cu_x=s(x,t)$, the source is accumulated along the characteristic:

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

Example: transport of a pollutant. A river flows steadily at speed $c=2\ \mathrm{m\,s^{-1}}$, and a factory releases a concentrated slug of pollutant at one point; as long as mixing and diffusion are negligible, the current simply carries the whole slug downstream without changing it. The concentration obeys the transport equation $u_t+2u_x=0$ with the Gaussian initial profile

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

(peak $50\ \mathrm{mg\,L^{-1}}$ at the release point, falling by $e^{-1}$ ten metres away). The solution above gives

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

After one minute the peak has moved from $x=0$ to $x=ct=120\ \mathrm{m}$, still reading $50\ \mathrm{mg\,L^{-1}}$; pure transport does not spread the slug, which would require the second-order term $\alpha u_{xx}$ of the heat equation.

Example: the wave equation. The second-order wave equation is the two-directional version of the same idea, and its general algebraic solution is a pair of waves travelling left and right (d'Alembert, 1747):

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

The term $f(x-ct)$ obeys the transport equation with speed $c$ and $g(x+ct)$ with speed $-c$; each therefore satisfies the wave equation, and linearity superposes them into the general solution. The functions $f$ and $g$ are fixed by the initial displacement and velocity: for a string released from rest with initial displacement $\phi(x)$, the conditions $u(x,0)=\phi(x)$ and $u_t(x,0)=0$ give $f=g=\phi/2$, so

$$u(x,t)=\frac{\phi(x-ct)+\phi(x+ct)}{2}$$

and the initial hump separates into two half-size copies travelling apart at speed $c$.

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