Differential equation

A differential equation is an equation that relates an unknown function to its own derivatives.

Two features are common to all differential equations:

• the unknown is a function, not a number;
• a solution is rarely unique, but instead often a family of functions containing arbitrary constants. A particular member of the solution family can be singled out by extra conditions, such as the state of the system at some time;

Types of differential equations

Ordinary and partial differential equations

Ordinary differential equations

An ordinary differential equation (ODE) describes an unknown function which depends on a single independent variable.

Writing the unknown as y(x), an ODE can be expressed in its implicit form:

F\bigl(x,y,y',y'',\ldots,y^{(n)}\bigr)=0

where F is a given function.

If the equation can be solved for the highest derivative, it takes an explicit form.

For a first-order equation:

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

The right-hand side f(x, y) prescribes the slope that any solution function must have at the point (x, y).

For an equation of order n that is solved for its highest derivative:

y^{(n)}=f\bigl(x,y,y',y'',\ldots,y^{(n-1)}\bigr)

Integrating the equation once removes one derivative and introduces one arbitrary constant. Therefore, a problem of order n generally needs n conditions to single out one particular solution of the solution family.

Partial differential equations

A partial differential equation describes an unknown function that depends on two or more independent variables. Its general form is:

F\Bigl(x_1,x_2,\ldots,x_n,\ u,\ \frac{\partial u}{\partial x_1},\ldots,\ \frac{\partial^{2}u}{\partial x_i\,\partial x_j},\ldots\Bigr)=0

Order and linearity

The order of a differential equation is the order of the highest derivative that appears in it.

A differential equation is linear if the unknown function and its derivatives appear only to the first power and are never multiplied together (they may be multiplied by functions of the independent variables).

An n-th order linear ODE therefore has the form

a_n(x)\,y^{(n)}+a_{n-1}(x)\,y^{(n-1)}+\cdots+a_1(x)\,y'+a_0(x)\,y=g(x)

We consider this equation homogeneous when g(x) = 0.

Linear homogeneous equations obey the superposition principle: any linear combination of solutions is again a solution.

Differential equations as mathematical models

Mechanics

Force and acceleration

Newton's second law of motion states that the acceleration of a body is proportional to the net force acting on it:

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

where x(t) is the position of the body, m its mass, and the force F may itself depend on time, position and velocity.

Harmonic oscillations

For a mass attached to a spring whose restoring force is linear in the displacement, this becomes the harmonic oscillator equation:

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

whose solutions are sinusoidal oscillations at the natural frequency ω₀; such motion is called simple harmonic motion.

A mass on a spring executing simple harmonic motion, the solution of the harmonic oscillator equation. Credit: Evil saltine (public domain).

Exponential growth and decay

A quantity that changes at a rate proportional to its current size satisfies the exponential growth or decay equation:

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

Its solutions are exponentials:

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

When k > 0, it can describe a population with unlimited resources that grows exponentially, whereas k < 0 describes a radioactive substance that decays exponentially towards zero.

Heat transfer

Newton's law of cooling describes a body whose temperature T differs from a constant ambient temperature Tₐ:

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

The temperature difference decays exponentially, so the body approaches the ambient temperature but never quite reaches it.

Geometry of first-order equations

Direction fields

As described above, a first-order equation can be written explicitly

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

which assigns to every point (x, y) the slope that a solution curve must have there.

Drawing a short line segment of that slope at many points produces a direction field (or slope field). Every solution curve must stay tangent to the field directions everywhere.

Direction field of the equation dy/dx = y, with several solution curves. Each segment shows the local slope; the solutions are exponential curves. Credit: jjbeard (public domain).

As such, direction fields can visualise the qualitative behaviour of solutions.

Solving differential equations

Solutions to differential equations can be classified into 3 categories:

• Closed-form solution: a finite combination of elementary functions (polynomials, rational functions, exponentials, logarithms, trigonometric functions and their inverses) built with the usual algebraic operations, which can be evaluated directly at a given point. Only a small fraction of differential equations admit a closed-form solution.

• Open form solution: a combination of elementary functions, and/or infinite series (a power series or a Fourier series) or special functions defined by such series, and integrals.

• numerical estimation: approximation of the real solution to a desired accuracy.

Analytical methods

Classical methods convert the equation into simpler, solvable pieces, giving a closed form when possible and a series or an integral otherwise.

linear first order equations

Separation of variables: Reorganise by variable and reintegrate:

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

For the equation dy/dx = −x/y, this gives the family of circles:

\frac{dy}{dx}=-\frac{x}{y}\qquad\Longrightarrow\qquad x^{2}+y^{2}=C

Integrating factors

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

For instance:

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

homogeneous equations up to 2nd order

Characteristic equations.

For a homogeneous linear equation with constant coefficients, solutions of the exponential form y = e^rx are sought. Substituting the trial solution together with its derivatives into the equation, and using that e^rx never vanishes, reduces the differential equation to an algebraic characteristic equation. For the second-order equation y″ + ay′ + by = 0:

y=e^{rx},\quad y'=r\,e^{rx},\quad y''=r^{2}\,e^{rx};\qquad y''+a\,y'+b\,y=\bigl(r^{2}+a\,r+b\bigr)e^{rx}=0\ \Rightarrow\ r^{2}+a\,r+b=0,\quad e^{rx}\neq 0

The roots r of the characteristic equation build the general solution. For the example

y=e^{rx}\ \Rightarrow\ (r^{2}-3r+2)\,e^{rx}=0\ \Rightarrow\ r^{2}-3r+2=(r-1)(r-2)=0

whose roots r = 1 and r = 2 are real and distinct, each root contributes one exponential, and the general solution is the sum

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

When the roots form a complex-conjugate pair r = α ± iβ, Euler's formula merges the two complex exponentials into real oscillations, and the general solution is

r=\alpha\pm i\beta\ \Rightarrow\ y=e^{\alpha x}\bigl(C_{1}\cos\beta x+C_{2}\sin\beta x\bigr)

For the harmonic oscillator y″ + y = 0, the characteristic roots are r = ±i, so the general solution is sinusoidal:

y''+y=0\ \Rightarrow\ r^{2}+1=0\ \Rightarrow\ r=\pm i\ \Rightarrow\ y=C_{1}\cos x+C_{2}\sin x

For an equation of order n with constant coefficients, the same trial solution leads to a characteristic polynomial of degree n.

non-homogeneous equations up to 2nd order

A forced (or non-homogeneous) linear equation is one whose right-hand side is a non-zero function of x, called the forcing term g(x), which drives the system from outside, as a driving force drives an oscillator.

The structure of the solutions is clearest in the first-order case. For y′ − y = e^x, the complementary solution y_c = Ce^x solves the homogeneous equation; the forcing term e^x is itself a homogeneous solution, so the trial particular solution must be multiplied by x:

y'-y=e^{x}:\qquad y_{c}=C\,e^{x};\qquad y_{p}=A\,x\,e^{x}\ \Rightarrow\ y_{p}'-y_{p}=A\,e^{x}=e^{x}\ \Rightarrow\ A=1;\qquad y=C\,e^{x}+x\,e^{x}

In general, the solution of a forced linear equation is the sum of the complementary solution y_c, the general solution of the homogeneous equation that carries the arbitrary constants, and one particular solution y_p of the full equation:

y''+p(x)\,y'+q(x)\,y=g(x),\qquad y=y_{c}+y_{p},\qquad y_{c}''+p(x)\,y_{c}'+q(x)\,y_{c}=0

Variation of parameters is a general way to construct y_p from the homogeneous solutions y₁ and y₂: replace the constants in y_c = C₁y₁ + C₂y₂ by unknown functions, y_p = u₁(x)y₁ + u₂(x)y₂, and choose these functions so that y_p solves the forced equation. Differentiating the ansatz,

y_{p}=u_{1}y_{1}+u_{2}y_{2}\quad\Rightarrow\quad y_{p}'=u_{1}'y_{1}+u_{2}'y_{2}+u_{1}y_{1}'+u_{2}y_{2}';\qquad u_{1}'y_{1}+u_{2}'y_{2}=0\ \Rightarrow\ y_{p}'=u_{1}y_{1}'+u_{2}y_{2}'

The first condition, u₁′y₁ + u₂′y₂ = 0, is imposed so that y_p′ takes the same form as if the uᵢ were constants. With that, substituting y_p into the full left-hand side leaves only the terms in the uᵢ′, because the brackets vanish as y₁ and y₂ each solve the homogeneous equation:

y_{p}''+p\,y_{p}'+q\,y_{p}=u_{1}\bigl[y_{1}''+p\,y_{1}'+q\,y_{1}\bigr]+u_{2}\bigl[y_{2}''+p\,y_{2}'+q\,y_{2}\bigr]+u_{1}'y_{1}'+u_{2}'y_{2}'=u_{1}'y_{1}'+u_{2}'y_{2}'

What remains of y_p″ + py_p′ + qy_p is u₁′y₁′ + u₂′y₂′. Imposing the second condition, u₁′y₁′ + u₂′y₂′ = g(x), therefore makes y_p solve the forced equation. The two conditions form a 2×2 linear system for the unknown derivatives u₁′ and u₂′,

\begin{pmatrix}y_{1}&y_{2}\\ y_{1}'&y_{2}'\end{pmatrix}\binom{u_{1}'}{u_{2}'}=\binom{0}{g},\qquad W=y_{1}y_{2}'-y_{2}y_{1}';\qquad u_{1}'=\frac{\begin{vmatrix}0&y_{2}\\ g&y_{2}'\end{vmatrix}}{W}=-\frac{y_{2}\,g}{W},\qquad u_{2}'=\frac{\begin{vmatrix}y_{1}&0\\ y_{1}'&g\end{vmatrix}}{W}=\frac{y_{1}\,g}{W}

whose determinant is the Wronskian W = y₁y₂′ − y₂y₁′. Because y₁ and y₂ are linearly independent solutions of a second-order linear equation, W(x) is never zero; division by W is therefore always legitimate, and Cramer's rule yields u₁′ = −y₂g/W and u₂′ = y₁g/W. Integrating the two derivatives gives the functions u₁ and u₂.

For the forced example y″ − 3y′ + 2y = e^x, whose homogeneous solutions are y₁ = e^x and y₂ = e^2x:

y_{1}=e^{x},\quad y_{2}=e^{2x},\quad W=y_{1}y_{2}'-y_{2}y_{1}'=e^{3x};\qquad u_{1}'=-\frac{y_{2}\,g}{W}=-1,\quad u_{2}'=\frac{y_{1}\,g}{W}=e^{-x}

Thus u₁′ = −1 and u₂′ = e^−x, so u₁ = −x and u₂ = −e^−x. The function produced, −xe^x − e^x, contains the homogeneous solution e^x; it is absorbed into C₁, leaving the particular solution y_p = −xe^x and the general solution y = C₁e^x + C₂e^2x − xe^x.

Power series and the Frobenius method

When the coefficients of a linear equation vary with x, elementary closed forms are rare, but a solution can still be sought as a power series y = Σaₙxⁿ: substituting the series into the equation and equating the coefficients of like powers turns the differential equation into recurrence relations that fix the coefficients one after another. Near a regular singular point, the Frobenius method extends the same idea to series with a shifted leading power x^ρ. Such series solve equations exactly (in open form) even when no elementary solution exists, as for the Airy equation y″ − xy = 0, whose solutions are the Airy functions. The full step-by-step demonstrations, including the worked Airy equation and the Bessel equation, are given in the article Power series.

Laplace transforms

For initial value problems, the Laplace transform turns differentiation into multiplication by the transform variable, so the problem becomes an algebraic equation for the transformed function; solving it and inverting with a table of standard transforms yields the solution, with the initial conditions entering automatically. The derivation of the transform rules and a complete worked example are given in the article Laplace transform.

Separation of variables for PDEs

With the product ansatz u(x, t) = X(x)T(t), a linear PDE splits into ordinary equations tied together by a separation constant λ:

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

Superposition of the product solutions then fits the initial data. For the heat equation on a rod of length L whose ends are held at temperature zero, this gives the Fourier sine series

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

whose terms decay at rates set by the thermal diffusivity α.

The theory of these methods is treated in standard textbooks on ordinary and partial differential equations. Still, only a minority of differential equations, almost all of them linear, can be solved explicitly; most equations that arise in applications are treated by numerical or qualitative methods.

Numerical methods

When no closed-form solution is available, solutions are approximated by numerical integration. The simplest scheme, Euler's method, advances an approximate solution along the tangent of the exact solution: with step size h, each step

y_{n+1}=y_{n}+h\,f(x_{n},y_{n})

takes the approximation from xₙ to xₙ₊₁ = xₙ + h. The error of one step is of order h², so over a fixed interval the accumulated error is of order h; halving h halves the error, and the polygonal approximation converges to the exact solution as h → 0. More accurate and more stable methods, such as the Runge–Kutta family, evaluate the slope several times per step and are implemented in essentially every numerical computing environment. The derivation, error analysis, worked examples and code are given in the article Euler's method.

Qualitative methods

For nonlinear equations, and especially for systems of two or more equations, exact solutions are usually impossible, and the geometry of the solutions matters as much as their formulas. Qualitative analysis studies the equilibria, the points at which the system does not change; their stability under small perturbations; the periodic orbits of the system; the bifurcations at which this behaviour changes as a parameter varies; and, for systems in three or more dimensions, chaotic motion. These ideas originated in the geometric methods that Henri Poincaré developed for celestial mechanics at the end of the 19th century and grew, with the help of numerical computing, into the modern theory of dynamical systems. The methods and their full worked demonstrations, the pendulum in the phase plane, the collapse of an overharvested fishery and the Lorenz attractor, are treated in the article Qualitative methods.

History

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

Differential equations arose with the calculus that Isaac Newton and Gottfried Wilhelm Leibniz developed independently in the late 17th century. Newton's laws of motion and of universal gravitation, published in his Principia of 1687, are differential equations. However, it was Leibniz who introduced the notation dy/dx still used today.

In the 18th century Leonhard Euler turned a collection of ad hoc tricks into a systematic theory, developing the general treatment of linear equations with constant coefficients, series methods, and the first numerical integration scheme, which still bears his name.

Leonhard Euler (portrait by Jakob Emanuel Handmann, 1753). Euler created much of the systematic theory of differential equations and gave his name to Euler's method. Credit: Jakob Emanuel Handmann (public domain).

Alongside the theory of ordinary equations, the 18th and 19th centuries saw the development of partial differential equations of physics. Jean le Rond d'Alembert wrote down and solved the one-dimensional wave equation for the vibrating string in the 1740s, and Joseph Fourier derived the heat equation from the physics of conduction and solved it with trigonometric series in his Théorie analytique de la chaleur of 1822, founding Fourier analysis. Augustin-Louis Cauchy and, later, Rudolf Lipschitz and Émile Picard placed the existence and uniqueness of solutions on a rigorous footing. Henri Poincaré developed qualitative, geometric methods to study differential equations at the end of the 19th century to tackle the three-body problem of celestial mechanics. In 1892, Aleksandr Lyapunov gave a general definition of stability.

The growing power of numerical computing in the 20th century enabled the development of the theory of dynamical systems that now pervades the sciences. In 1963, meteorologist Edward Lorenz found chaotic behaviour in a simple system of three differential equations modelling atmospheric convection, a landmark in the study of nonlinear dynamics. Although the system was deterministic, its solutions were aperiodic and so sensitive to initial conditions that nearby states rapidly diverged. He popularised the phenomenon as the butterfly effect.

Further reading

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