Differential Equations

Ordinary differential equations (ODEs)

Scalar ODEs

Ordinary differential equation (ODE): equation of the form $F(x,y,y^{\prime },\dots ,y^{(n)})=0$, where $y$ is a function of $x$

Order of an ODE: order of the highest derivative appearing in the equation

Autonomous ODE: $F$ does not (explicitly) depend on $x$

Linear ODE: $F$ is a linear combination of the derivatives of $y$, plus a function of $x$ called the source term/forcing function; that is, an equation of the form $$\sum _{i=0}^n{a}_i(x)y^{(i)}=b(x)$$ ​ Homogeneous linear ODE: $b(x)\equiv 0$ (in which case $y=0$ is a trivial solution)

​ Nonhomogeneous linear ODE: $b(x)\ne 0$

The general solution to a linear ODE is $y=y_c+y_p$, where $y_c$ is the general solution of the homogeneous equation (the complementary solution) and $y_p$ is a solution of the nonhomogeneous equation (a particular solution).

Superposition principle: a linear combination of solutions to a homogeneous linear ODE is also a solution to the ODE ¹

A basis for the space of solutions to a homogeneous linear ODE is called a fundamental system of solutions.

First-order scalar ODEs

Picard-Lindelöf theorem

Given the initial value problem $y^{\prime }=f(x,y),y(x_0)=y_0$, suppose that $f$ is continuous in $x$ and Lipschitz continuous in $y$ (uniformly in $x$). Then there exists a unique solution to the IVP on $[x_0-ϵ,x_0+ϵ]$ for some $ϵ>0$.

Separable equations

To solve a separable equation $y^{\prime }=f(x)g(y)$:

• Integrate: $$\int \frac{dy}{g(y)}=\int f(x)dx$$

First-order linear ODEs

To solve a first-order linear ODE of the form $y^{\prime }+p(x)y=f(x)$:

• Multiply both sides by an integrating factor $r(x)=e^{\int p(x)dx}$, which yields $(r(x)\cdot y{)}^{\prime }=r(x)f(x)$

• Integrate: $$y=\frac{1}{r(x)}\int r(x)f(x)dx=e^{-\int p(x)dx}\int e^{\int p(x)dx}f(x)dx$$

Bernoulli equations

To solve a Bernoulli equation $y^{\prime }+p(x)y=q(x)y^n$:

• Substitute $u=y^{1-n}$, which yields the first-order linear equation $\frac{1}{1-n}{u}^{\prime }+p(x)u=q(x)$

Given a first-order autonomous equation $y^{\prime }=f(y)$, the points on the $y$-axis where $f(y)=0$ are called critical points. If $y_0$ is a critical point, the constant solution $y=y_0$ is called an equilibrium solution.

The phase line of such an equation consists of the $y$-axis along with arrows indicating the sign of $f$ between each of the critical points – ▲ if $f>0$; ▼ if $f<0$. If both arrows point toward a critical point, it is stable; if both point away, it is unstable; otherwise, it is semi-stable (sometimes just called unstable).

A general first-order scalar ODE $y^{\prime }=f(x,y)$ may be visualized using a slope field, a plot of the $xy$-plane with small line segments of slope $f(x,y)$ at each point $(x,y)$.

Second-order scalar ODEs

We will consider only linear second-order ODEs.

Existence and uniqueness of solutions to second-order linear ODEs

Given the equation $y^{″}+p(x)y^{\prime }+q(x)y=f(x)$, suppose that $p,q,f$ are continuous on some interval $I$. Then for any fixed $x_0\in I$ and $y_0,y_0^{\prime }$, there exists a unique solution to the ODE on $I$ satisfying $y(x_0)=y_0,y^{\prime }(x_0)=y_0^{\prime }$.

Superposition of solutions to second-order linear homogeneous ODEs

Given the equation $y^{″}+p(x)y^{\prime }+q(x)y=0$, suppose that $p,q$ are continuous. Then the solution space of the ODE is two-dimensional; viz., if $y_1,y_2$ are linearly independent solutions to the ODE, the general solution is $y=c_1{y}_1+c_2{y}_2$.

Now we consider the constant-coefficient second-order linear homogeneous ODE $a{y}^{″}+b{y}^{\prime }+cy=0$, where $a,b,c\in \mathbb{R}$. Assuming the solution is of the form $y=e^{\lambda x}$ (with $\lambda$ possibly non-real), we obtain the characteristic equation $a{\lambda }^2+b\lambda +c=0$.

Second-order linear homogeneous ODEs with constant coefficients

To solve an ODE of the form $a{y}^{″}+b{y}^{\prime }+cy=0$:

• Compute the roots ${\lambda }_1,{\lambda }_2$ of its characteristic equation

• If ${\lambda }_1,{\lambda }_2$ are real and distinct ($b^2-4ac>0$), the general solution is $y=c_1{e}^{{\lambda }_{1}x}+c_2{e}^{{\lambda }_{2}x}$ ²

• If ${\lambda }_1,{\lambda }_2$ are real and equal ($b^2-4ac=0$), the general solution is $y=c_1{e}^{{\lambda }_{1}x}+c_{2}x{e}^{{\lambda }_{1}x}$

• If ${\lambda }_1,{\lambda }_2$ are complex ($b^2-4ac<0$) with ${\lambda }_{1,2}=\mu \pm \omega i$, the general solution is $y=c_1{e}^{\mu x}\cos (\omega x)+c_2{e}^{\mu x}\sin (\omega x)$

Reduction of order of a second-order linear ODE

Given the equation $y^{″}+p(x)y^{\prime }+q(x)y=f(x)$:

• Let $y_1$ be a solution to the homogeneous equation $y^{″}+p(x)y^{\prime }+q(x)y=0$
• Substitute the ansatz $y_2=v(x)y_1$, which yields a first-order linear ODE in $y^{\prime }$
• Solving for $y_2$ yields a second, linearly independent solution to the second-order ODE

In particular, reduction of order applied to $a{y}^{″}+b{y}^{\prime }+cy=0$ in the case ${\lambda }_1={\lambda }_2=-b/2a$ with $y_1=e^{{\lambda }_{1}x}$ yields $y_2=x{e}^{{\lambda }_{1}x}$.

The (translational mechanical) harmonic oscillator is modelled by the constant-coefficient second-order linear ODE $$m{x}^{″}(t)+c{x}^{\prime }(t)+kx(t)=F(t),$$ where $x(t)$ is the position of the object as a function of time, $m>0$ is the mass of the object, $c\ge 0$ is the damping coefficient, $k>0$ is the spring constant, and $F$ is a driving/forcing function.

Several physical systems are harmonic oscillators, as exemplified by the following table.

Damped; undamped harmonic oscillator: $c>0$; $c=0$

Forced; unforced/free harmonic oscillator: $F\not\equiv 0$; $F\equiv 0$

Undamped/natural angular frequency: ${\omega }_0=\sqrt{k/m}$ (i.e., the angular frequency of the oscillator if it were undamped) ³

Damping ratio: $\zeta =c/(2\sqrt{mk})$

We first consider the behaviour of unforced oscillators.

Behaviour of an unforced undamped harmonic oscillator

If $c=0$ (i.e., $\zeta =0$), the oscillator is undamped. The solution $$x=\sqrt{c_1^2+c_2^2}\cos ({\omega }_{0}t-{\tan }^{-1}\left(\frac{c_2}{c_1}\right))$$ oscillates with amplitude $\sqrt{c_1^2+c_2^2}$, angular frequency ${\omega }_0$, and phase shift ${\tan }^{-1}(c_2/c_1)$.

If the oscillator is damped, we distinguish three cases, depending as above on the sign of the discriminant of the characteristic polynomial.

Behaviour of an unforced damped harmonic oscillator

• If $c^2-4mk>0$ (i.e., $\zeta >1$), the oscillator is overdamped

  • The solution $x=c_1{e}^{{\lambda }_{1}t}+c_2{e}^{{\lambda }_{2}t}$ decays as $t\to \mathrm{\infty }$ and does not oscillate
  • Its decay constants are $-{\lambda }_{1,2}={\omega }_0\zeta \mp {\omega }_0\sqrt{{\zeta }^2-1}$

• If $c^2-4mk=0$ (i.e., $\zeta =1$), the oscillator is critically damped

  • The solution $x=c_1{e}^{{\lambda }_{1}t}+c_{2}t{e}^{{\lambda }_{1}t}$ decays as $t\to \mathrm{\infty }$ and does not oscillate
  • Its decay constant is $-{\lambda }_1={\omega }_0\zeta$

• If $c^2-4mk<0$ (i.e., $\zeta <1$), the oscillator is underdamped

  • The solution $x=e^{\mu x}(c_1\cos ({\omega }_{1}t)+c_2\sin ({\omega }_{1}t))$ decays as $t\to \mathrm{\infty }$ while oscillating ⁴

    • Its decay constant is $-\mu ={\omega }_0\zeta$ and the angular pseudo-frequency of its oscillations is ${\omega }_1={\omega }_0\sqrt{1-{\zeta }^2}$

Now suppose the oscillator is subject to sinusoidal forcing $F(t)=F_0\cos (\omega t)$.

Behaviour of a forced undamped harmonic oscillator

• If $\omega \ne {\omega }_0$, there are beats
  • The particular solution is $x_p=F_0\cos (\omega t)/(m({\omega }_0^2-{\omega }^2))$
  • This combines with $x_c$ to produce beats; the closer $\omega$ is to ${\omega }_0$, the greater the amplitude of the beats
• If $\omega ={\omega }_0$, there is resonance
  • The particular solution is $x_p=F_{0}t\sin (\omega t)/(2m\omega )$
  • This dominates $x_c$ as $t\to \mathrm{\infty }$, producing increasingly large oscillations at the resonant frequency

Behaviour of a forced damped harmonic oscillator

The particular solution $x_p$ is a sinusoid with amplitude $$\frac{F_0}{\sqrt{[m({\omega }_0^2-{\omega }^2){]}^2+(c\omega {)}^2}}$$ and angular frequency $\omega$, and is called the steady-state solution. The complementary solution decays exponentially and is therefore called the transient solution.

The amplitude of $x_p$ is maximized at the practical resonance frequency $\omega ={\omega }_0\sqrt{1-2{\zeta }^2}$, provided that $\zeta <1/\sqrt{2}$.

If $\zeta \ge 1/\sqrt{2}$, there is no maximum (for $\omega >0$), but the amplitude increases as $\omega \to 0^{+}$.

Higher-order scalar ODEs

$n$^th-order linear homogeneous ODEs with constant coefficients

To solve an ODE of the form $\sum _{i=0}^n{a}_i{y}^{(i)}=0$:

• Compute the roots of its characteristic polynomial $\sum _{i=0}^n{a}_i{\lambda }^i$
• Each real root $\lambda$ of multiplicity $k$ contributes $(\sum _{j=0}^{k-1}{c}_j{x}^j)e^{\lambda x}$ to the general solution
• Each pair of complex roots $\mu \pm \omega i$ of (individual) multiplicity $k$ contributes $(\sum _{j=0}^{k-1}{c}_j{x}^j)e^{\mu x}\cos (\omega x)+(\sum _{j=0}^{k-1}{d}_j{x}^j)e^{\mu x}\sin (\omega x)$ to the general solution

There are two commonly used methods for finding a particular solution to a nonhomogeneous linear ODE with source term $f(x)$.

The method of undetermined coefficients

Suppose that $f(x)=P(x)e^{\lambda x}$ for some polynomial $P$ of degree $m$ and some $\lambda \in \mathbb{C}$. Given the equation $\sum _{i=0}^n{a}_i{y}^{(i)}=f(x)$:

• Let $k\ge 0$ be the multiplicity of $\lambda$ as a root of its characteristic polynomial

• Substitute the ansatz $y_p=x^k\left(\sum _{j=0}^m{A}_j{x}^j\right)e^{\lambda x}$ and match coefficients to determine the $A_j$

  • If $\lambda =\mu \pm \omega i$, the ansatz $y_p=x^k\left(\sum _{j=0}^m{A}_j{x}^j\right)e^{\mu x}\cos (\omega x)+x^k\left(\sum _{j=0}^m{B}_j{x}^j\right)e^{\mu x}\sin (\omega x)$ can be used instead (which is useful when $f(x)=P(x)\cos (\omega x)$ or $f(x)=P(x)\sin (\omega x)$, for example)

• If $f$ is a sum of such terms, apply the principle of superposition; i.e., solve the equation with each term as the source term separately, then add the solutions

To verify that a solution set of an ODE is linearly independent, it is occasionally useful to compute the Wronskian determinant.

Wronskian (determinant) of $f_1,\dots ,f_n$: $$W(f_1,\dots ,f_n)(x)=det\left[\begin{array}{ccc}{f}_1(x)& \cdots & f_n(x)\\ f_1^{\prime }(x)& \cdots & f_n^{\prime }(x)\\ ⋮& \ddots & ⋮\\ f_1^{(n-1)}(x)& \cdots & f_n^{(n-1)}(x)\end{array}\right]$$

Linear dependence and the Wronskian

Given a set of $(n-1)$-times differentiable functions, if the functions are linearly dependent on an interval, then their Wronskian vanishes identically thereon. ⁵

Of course, the contrapositive of this result is used in practice.

The method of variation of parameters

Given the equation $y^{(n)}+\sum _{i=0}^{n-1}{a}_i(x)y^{(i)}=f(x)$:

• Let $y_1,\dots ,y_n$ be a fundamental system of solutions to the homogeneous equation

• Substitute the ansatz $y_p=\sum _{i=1}^n{c}_i(x)y_i$ and impose the constraints $\sum _{i=1}^n{c}_i^{\prime }(x)y_i^{(j)}=0$ for $j=0,\dots ,n-2$

• Solve the resulting linear system $$\left[\begin{array}{ccc}{y}_1& \cdots & y_n\\ ⋮& \ddots & ⋮\\ y_1^{(n-2)}& \cdots & y_n^{(n-2)}\\ y_1^{(n-1)}& \cdots & y_n^{(n-1)}\end{array}\right]\left[\begin{array}{c}{c}_1^{\prime }(x)\\ ⋮\\ c_n^{\prime }(x)\end{array}\right]=\left[\begin{array}{c}0\\ ⋮\\ 0\\ f(x)\end{array}\right]$$

  • Note that the matrix is the Wronskian matrix of the fundamental system, so by Cramer’s rule, $c_i(x)=\int W_i(y_1,\dots ,y_n)(x)/W(y_1,\dots ,y_n)(x)dx$, where $W_i(y_1,\dots ,y_n)$ denotes the Wronskian determinant with the $i$^th column of the matrix replaced by the right-hand side

The Laplace transform

Laplace transform of $f:[0,\mathrm{\infty })\to \mathbb{C}$: $$\mathcal{L}\{f\}(s)=F(s)={\int }_0^{\mathrm{\infty }}f(t)e^{-st}dt$$

• If $f\in L_{\mathrm{loc}}^1([0,\mathrm{\infty }))$ and $f$ is of exponential type (that is, $f\in \mathcal{O}(e^{at})$ as $t\to \mathrm{\infty }$ for some $a\in \mathbb{R}$), then $\mathcal{L}\{f\}$ is defined for all $\mathrm{\Re }(s)>a$

  • Moreover, $\underset{s\to \mathrm{\infty }}{lim}\mathcal{L}\{f\}(s)=0$

• The Laplace transform is linear

• The Laplace transform is injective in the sense that if $\mathcal{L}\{f\}=\mathcal{L}\{g\}$, then $f=g$ a.e. (Lerch’s theorem)

  • In particular, if $f$ and $g$ are continuous, then $f=g$ for all $t\ge 0$

• The variables $t$ and $s$ are typically thought of as “time” and “frequency”, respectively

Inverse Laplace transform of $F:\mathbb{C}\to \mathbb{C}$: $${\mathcal{L}}^{-1}\{F\}(t)=f(t)=\frac{1}{2\pi i}{\int }_{a-\mathrm{\infty }i}^{a+\mathrm{\infty }i}F(s)e^{st}ds,$$ where $a\in \mathbb{R}$ is such that the contour of integration (i.e., the line $\mathrm{\Re }(s)=a$) lies in the region of convergence of $F$

• The formula above is called Mellin’s integral formula and is derived from the Fourier inversion theorem
• In practice, the inverse Laplace transform is computed by inspection (using tables of known Laplace transforms; see below)

Dirac delta “function”: Borel measure defined by $\delta (E)=1_E(0)$; $\mathcal{L}\{\delta \}$ is interpreted as $\mathcal{L}\{\delta \}(s)={\int }_0^{\mathrm{\infty }}{e}^{-st}d\delta =1$

Heaviside step function: $H(t)=1_{[0,\mathrm{\infty })}(t)$ (sometimes denoted $\theta$ or $u$)

Unit step functions: $u_a(t)=H(t-a)$

We observe that $1_{[a,b]}=u_a-u_b$, which is useful for expressing piecewise functions.

When applicable, the Laplace transform can be used to solve ODEs by transforming both sides of the equation, solving for the transform of the independent variable, and computing its inverse transform.

In the tables below, we assume where necessary that $f(t)=0$ for all $t<0$. Thus, for instance, $t^n$ represents $t^{n}H(t)$ and $f(t-a)$ represents $f(t-a)u_a(t)$.

These “elementary” transforms can be combined with the general properties in the table below.

For example, $\mathcal{L}\{u_a\}(s)=e^{-as}/s$ by the time shift property.

Initial value theorem

If $f:[0,\mathrm{\infty })\to \mathbb{C}$ is of exponential type and $\underset{t\to 0^{+}}{lim}f(t)$ exists (and is finite), then $\underset{t\to 0^{+}}{lim}f(t)=\underset{s\to \mathrm{\infty }}{lim}sF(s)$.

Final value theorem

If $f:[0,\mathrm{\infty })\to \mathbb{C}$ is bounded and $\underset{t\to \mathrm{\infty }}{lim}f(t)$ exists (and is finite), then $\underset{t\to \mathrm{\infty }}{lim}f(t)=\underset{s\to 0^{+}}{lim}sF(s)$.

If $L$ is a linear differential operator with constant coefficients (that is, a polynomial differential operator), then the solution to $Lx=\delta (t)$ is called the impulse response of the system.

Linear ODEs with constant coefficients and the impulse response

If $L$ is a linear differential operator and $x_{\delta }$ is the impulse response of $Lx=\delta (t)$, then the solution to $Lx=f(t)$ is $x_{\delta }\ast f$.

(To see this, we can convolve both sides of the equation $L{x}_{\delta }=\delta$ with $f$ and use the fact that $L{x}_{\delta }\ast f=L(x_{\delta }\ast f)$ and $\delta \ast f=f$. We assume that $x_{\delta }$ is sufficiently smooth to justify the manipulation on the left-hand side.)

Power series methods

Consider the second-order linear homogeneous ODE $$p(x)y^{″}+q(x)y^{\prime }+r(x)y=0,$$ where $p,q,r$ are polynomials. ⁶

Ordinary; singular point $x_0$: $p(x_0)\ne 0$; $p(x_0)=0$ ⁷

Regular singular point $x_0$: $x_0$ is a singular point but $(x-x_0)q(x)/p(x)$ and $(x-x_0{)}^{2}r(x)/p(x)$ tend to finite limits as $x\to x_0$ ⁸

At an ordinary point $x_0$, we can substitute the ansatz $y=\sum _{k=0}^{\mathrm{\infty }}{a}_k(x-x_0{)}^k$ and derive a recurrence relation for the coefficients $a_k$.

At a regular singular point, we substitute the ansatz $y=\sum _{k=0}^{\mathrm{\infty }}{a}_k(x-x_0{)}^{k+r}$. By equating the trailing coefficient (usually that of $x^r$) of the resulting series to zero, we obtain the indicial equation.

Method of Frobenius

Suppose that $x_0$ is a regular singular point of $p(x)y^{″}+q(x)y^{\prime }+r(x)y=0$ and that $r_1,r_2$ are the roots of the indicial equation. Then the ODE has two linearly independent solutions $y_1,y_2$ as given below.

• If $r_1-r_2\notin \mathbb{Z}$, then $$y_1=\sum _{k=0}^{\mathrm{\infty }}{a}_k(x-x_0{)}^{k+r_1},y_2=\sum _{k=0}^{\mathrm{\infty }}{b}_k(x-x_0{)}^{k+r_2}$$

• If $r_1=r_2$, then $$y_1=\sum _{k=0}^{\mathrm{\infty }}{a}_k(x-x_0{)}^{k+r_1},y_2=y_1\ln (x-x_0)+\sum _{k=0}^{\mathrm{\infty }}{b}_k(x-x_0{)}^{k+r_2}$$

• If $r_1-r_2\in \mathbb{Z}\setminus \{0\}$, then $$y_1=\sum _{k=0}^{\mathrm{\infty }}{a}_k(x-x_0{)}^{k+r_1},y_2=C{y}_1\ln (x-x_0)+\sum _{k=0}^{\mathrm{\infty }}{b}_k(x-x_0{)}^{k+r_2}$$

If $r_1,r_2\in \mathbb{C}$ with $r_1=\overline{r_2}$, we can take $\mathrm{\Re }(y_1)$ and $\mathrm{\Im }(y_1)$ to obtain linearly independent real solutions.

Vector ODEs

A vector ODE is simply an equation of the form $\mathbf{F}(x,\mathbf{y},{\mathbf{y}}^{\prime },\dots ,{\mathbf{y}}^{(n)})=\mathbf{0}$, where $\mathbf{y}$ is a function of $x$. Most definitions and properties given in the section on Scalar ODEs generalize to vector ODEs in the obvious manner.

We note that a linear vector ODE takes the form $$\sum _{i=0}^n{A}_i(x){\mathbf{y}}^{(i)}=\mathbf{b}(x),$$ where the $A_i$ are matrix-valued functions and the source term $\mathbf{b}$ is vector-valued.

When the (vector) functions of a fundamental system of solutions to a homogeneous linear vector ODE are concatenated horizontally, the resulting matrix is called a fundamental matrix (solution).

The analogue of the phase line for first-order autonomous vector ODEs ${\mathbf{y}}^{\prime }=\mathbf{F}(\mathbf{y})$ is the direction/vector field. A plot of the vector field $\mathbf{F}$ is made in phase space (i.e., the domain of $\mathbf{y}$), in contradistinction to slope fields, which are plotted in solution space. Solution curves (also called trajectories) parametrized by $x$ may also be drawn in phase space, even for non-autonomous ODEs. Such curves remain tangent to the vector field. A plot of several trajectories in phase space is called a phase portrait.

Vector ODEs are useful for expressing (systems of) scalar ODEs. An $n$^th-order (scalar) ODE $F(x,y,y^{\prime },\dots ,y^{(n)})=0$ can equivalently be viewed as a system of $n$ 1^st-order ODEs $$\{\begin{array}{l}{y}_0^{\prime }=y_1\\ y_1^{\prime }=y_2\\ ⋮\\ y_{n-2}^{\prime }=y_{n-1}\\ F(x,y_0,y_1,\dots ,y_{n-1},y_{n-1}^{\prime })=0\end{array},$$ where $y_i=y^{(i)}$ for $i=0,\dots ,n-1$. Similarly, a system of $k$ ODEs of order $n$ (with $k$ dependent variables) can be converted to a system of $kn$ ODEs of order 1 (with $kn$ dependent variables).

We will therefore restrict our attention to systems of 1^st-order ODEs, which can themselves be regarded as individual 1^st-order vector ODEs. In the example above, we can define $\mathbf{y}=(y_0,\dots ,y_{n-1})$ and $\mathbf{F}(x,\mathbf{y},{\mathbf{y}}^{\prime })=(0,0,\dots ,F(x,{\mathbf{y}}_1,{\mathbf{y}}_2,\dots ,{\mathbf{y}}_n,{\mathbf{y}}_n^{\prime }))$, where ${\mathbf{v}}_i$ denotes the $i$^th component of $\mathbf{v}$.

Linear vector ODEs

A first-order linear vector ODE (which we shall refer to interchangeably as a “linear system of ODEs”) can be written in the form ${\mathbf{x}}^{\prime }=A(t)\mathbf{x}+\mathbf{f}(t)$ (cf. first-order linear scalar ODEs). We will primarily be interested in an even more specific type of first-order linear vector ODEs: those with constant coefficients (i.e., $A(t)$ is constant).

As with second-order linear scalar ODEs, we begin with the homogeneous case.

First-order linear homogeneous vector ODEs with constant coefficients

The general solution of ${\mathbf{x}}^{\prime }=A\mathbf{x}$ is $\mathbf{x}=e^{tA}\mathbf{c}$, where $\mathbf{c}=\mathbf{x}(0)$. In other words, $X(t)=e^{tA}$ is a fundamental matrix solution.

(See Appendix B for the definition of $e^{tA}$ and further remarks.)

Corollary:

If $X$ is a fundamental matrix solution to ${\mathbf{x}}^{\prime }=A\mathbf{x}$, then $e^{tA}=X(t)X^{-1}(0)$. ⁹

We note that ${\mathbf{x}}^{\prime }=A\mathbf{x}$ is an autonomous system whose sole critical point is the origin.

In general, if ${\mathbf{x}}_0$ is an isolated critical point of a first-order autonomous system (not necessarily linear), we can classify the critical point according to the behaviour of trajectories near it. We can also describe the stability of the point itself as follows.

Stable critical point: given a distance $ϵ>0$ and any initial condition within a possibly smaller distance from ${\mathbf{x}}_0$, the trajectory of the system therefrom will never travel further than $ϵ$ from ${\mathbf{x}}_0$

Asymptotically stable critical point: ${\mathbf{x}}_0$ is stable and moreover every trajectory beginning sufficiently close to ${\mathbf{x}}_0$ converges to it

Unstable critical point: ${\mathbf{x}}_0$ is not stable

A stable critical point that is not asymptotically so is sometimes called marginally/neutrally stable.

When the system is two-dimensional and $A$ is invertible ¹⁰ and has distinct eigenvalues, the signs of the eigenvalues of $A$ determine the behaviour near and stability of the origin:

(The “centre” is so called because trajectories are ellipses centred at the origin.)

If $A$ has a repeated eigenvalue, we further distinguish between proper and improper nodes according as the eigenvalue is nondefective or defective: ¹¹

We now consider the general (i.e., non-homogeneous) case (cf. First-order scalar ODEs).

First-order linear vector ODEs with constant coefficients

To solve a first-order linear vector ODE of the form ${\mathbf{x}}^{\prime }+P\mathbf{x}=\mathbf{f}(t)$:

• Multiply both sides by the integrating factor $R(t)=e^{tP}$, which yields $(R(t)\mathbf{x}{)}^{\prime }=R(t)\mathbf{f}(t)$

• Integrate: $$\mathbf{x}={R(t)}^{-1}\int R(t)\mathbf{f}(t)dt=e^{-tP}\int e^{tP}\mathbf{f}(t)dt$$

However, the integrating factor method is not applicable when $P$ is a function of $t$, since in general $(\exp (\int P(t)dt){)}^{\prime }\ne P(t)\exp (\int P(t)dt)$ as matrix multiplication is noncommutative.

(For constant-coefficient linear systems, the method of undetermined coefficients may also be used, where the ‘coefficients’ are vectors. One difference is that ansatz terms augmented by powers of $t$ must be included in addition to, rather than in place of other terms. For instance, if $(0,1)e^t$ is part of the complementary solution, then both $\mathbf{a}{e}^t$ and $\mathbf{b}t{e}^t$ must be included in the ansatz. In general, this method is unwieldy for vector ODEs.)

To handle variable-coefficient systems, we can use variation of parameters:

First-order linear vector ODEs with variable coefficients

To solve a first-order linear vector ODE of the form ${\mathbf{x}}^{\prime }=A(t)\mathbf{x}+\mathbf{f}(t)$:

• Let $X(t)$ be a fundamental matrix solution to the homogeneous equation

• Substitute the ansatz ${\mathbf{x}}_p=X(t)\mathbf{c}(t)$, which yields $X(t){\mathbf{c}}^{\prime }(t)=\mathbf{f}(t)$

• Integrate: $${\mathbf{x}}_p=X(t)\int X(t{)}^{-1}\mathbf{f}(t)dt$$

(Note that this reduces to the integrating factor method when $A(t)\equiv -P$ and $X(t)=e^{-tP}$.)

When $A$ is diagonalizable, we can solve first- and second-order linear vector ODEs of the form ${\mathbf{x}}^{″}=A\mathbf{x}+\mathbf{f}(t)$ by exploiting the eigendecomposition of $A$.

Linear vector ODEs with constant coefficients

To solve a first- or second-order linear vector ODE of the form ${\mathbf{x}}^{(\nu )}=A\mathbf{x}+\mathbf{f}(t)$, where $A$ is diagonalizable:

• Let $({\lambda }_1,{\mathbf{v}}_1),\dots ,({\lambda }_n,{\mathbf{v}}_n)$ be the eigenpairs of $A$

• Write $\mathbf{x}(t)=\sum _{i=1}^n{\mathbf{v}}_i{\xi }_i(t)$ and $\mathbf{f}(t)=\sum _{i=1}^n{\mathbf{v}}_i{\varphi }_i(t)$

• Solve for the ${\varphi }_i$ (in the linear system $\left[\begin{array}{ccc}{\mathbf{v}}_1& \cdots & {\mathbf{v}}_n\end{array}\right]{\left[\begin{array}{ccc}{\varphi }_1& \cdots & {\varphi }_n\end{array}\right]}^T=\mathbf{f}$)

• Substitute these decompositions into the equation and equate the coefficients of the eigenvectors, which yields $n$ decoupled first- or second-order linear scalar ODEs $$\begin{array}{rl}{\xi }_1^{(\nu )}& ={\lambda }_1{\xi }_1+{\varphi }_1,\\ & ⋮\\ {\xi }_n^{(\nu )}& ={\lambda }_n{\xi }_n+{\varphi }_n,\end{array}$$

  • Solve the scalar ODEs using the methods described in First-order scalar ODEs or Second-order scalar ODEs

Nonlinear vector ODEs

Suppose that ${\mathbf{x}}_0$ is an isolated critical point of the first-order autonomous vector ODE ${\mathbf{x}}^{\prime }=\mathbf{F}(\mathbf{x})$ (i.e., $\mathbf{F}({\mathbf{x}}_0)=\mathbf{0}$). By the Taylor expansion of $\mathbf{F}$ about ${\mathbf{x}}_0$, we have $${\mathbf{x}}^{\prime }\approx \mathbf{F}({\mathbf{x}}_0)+d{\mathbf{F}}_{{\mathbf{x}}_0}(\mathbf{x}-{\mathbf{x}}_0)=d{\mathbf{F}}_{{\mathbf{x}}_0}(\mathbf{x}-{\mathbf{x}}_0).$$ In other words, near ${\mathbf{x}}_0$, the nonlinear ODE behaves like the constant-coefficient linear ODE ${\mathbf{v}}^{\prime }=d{\mathbf{F}}_{{\mathbf{x}}_0}\mathbf{v}$ – its linearization at ${\mathbf{x}}_0$ – provided that the linearized system itself has an isolated critical point (see Linear vector ODEs). In this case, the nonlinear ODE is called almost linear at ${\mathbf{x}}_0$.

However, we note that a centre of the linearized system tends to correspond to a spiral point in the nonlinear ODE, since it is unlikely that the real parts of the eigenvalues of $d\mathbf{F}$ all vanish in a whole neighbourhood of ${\mathbf{x}}_0$. (Similarly, a node arising from a repeated real eigenvalue could also correspond to a spiral point.)

Conservative equation: equation of the form $x^{″}(t)+f(x)=0$

A conservative equation can be written as the first-order autonomous vector ODE $${\underset{\mathbf{x}}{\underbrace{\left[\begin{array}{c}x\\ y\end{array}\right]}}}^{\prime }=\underset{\mathbf{F}(\mathbf{x})}{\underbrace{\left[\begin{array}{c}y\\ -f(x)\end{array}\right]}},$$ where $y=x^{\prime }$.

By integrating both sides of the scalar equation with respect to $x$, we find that $$\begin{array}{rl}C=\int x^{″}dx+\int f(x)dx& =\int y\frac{dy}{dx}dx+\int f(x)dx\\ & =\int ydy+\int f(x)dx\\ & =\frac{1}{2}{y}^2+\int f(x)dx\end{array}$$ since $x^{″}=y^{\prime }=dy/dx\cdot dx/dt$. Hence the quantity $\frac{1}{2}{y}^2+\int f(x)dx$ is conserved for any solution, which can be used to find trajectories in phase (i.e., $xy$) space. ¹²

The critical points are clearly the points on the $x$-axis where $f(x)=0$. Since $$d\mathbf{F}=\left[\begin{array}{cc}0& 1\\ -f^{\prime }(x)& 0\end{array}\right],$$ a critical point is almost linear when $f^{\prime }(x)\ne 0$. The eigenvalues of $d\mathbf{F}$ at such a point are $\pm \sqrt{-f^{\prime }(x)}$, so the point is a saddle point if $f^{\prime }(x)<0$ and a centre if $f^{\prime }(x)>0$. ^{13 14}

Partial differential equations (PDEs)

Separation of variables

Boundary-value problems

Let $A:\mathcal{D}(A)\subseteq L^2([a,b])\to L^2([a,b])$ be the linear operator given by $A=-\frac{d^2}{d{t}^2}$. We will take $\mathcal{D}(A)$ to be some subspace of $C^2([a,b])$ of functions satisfying certain boundary conditions (BCs). The three main types are:

• Dirichlet boundary conditions: $x(a)=x(b)=0$
• Neumann boundary conditions: $x^{\prime }(a)=x^{\prime }(b)=0$
• Periodic boundary conditions: $x(a)=x(b),x^{\prime }(a)=x^{\prime }(b)$

Dirichlet and Neumann boundary conditions can also be mixed; e.g., $x^{\prime }(a)=0,x(b)=0$ are “Neumann-Dirichlet” boundary conditions.

All such conditions ensure that $A$ is a symmetric operator (by integrating by parts, we see that $A$ is symmetric if and only if $x_1^{\prime }{x}_2-x_1{x}_2^{\prime }{|}_a^b=0$ for all $x_1,x_2\in \mathcal{D}(A)$).

The eigenvectors of $A$ are the nontrivial solutions to the boundary value problem (BVP) $x^{″}+\lambda x=0$, where $x$ satisfies one of the boundary conditions above. In this context, they are called eigenfunctions. As $A$ is symmetric, its eigenvalues are real, and eigenfunctions of distinct eigenvalues are orthogonal.

Fredholm alternative

For any given $\lambda$, either:

• $\lambda$ is an eigenvalue of $A$
• $A-\lambda I=f$ has a unique solution for every $f\in C([a,b])$

Trigonometric series

Fourier series of $f\in L^2([-L,L])$: $$\frac{1}{2}{a}_0+\sum _{n=1}^{\mathrm{\infty }}{a}_n\cos \left(\frac{n\pi }{L}t\right)+b_n\sin \left(\frac{n\pi }{L}t\right),$$ where $$\begin{array}{rlrl}{a}_0& =2\frac{⟨f,1⟩}{\Vert 1{\Vert }^2}& & =\frac{1}{L}{\int }_{-L}^{L}f(t)dt,\\ a_n& =\frac{⟨f,\cos (n\pi t/L)⟩}{\Vert \cos (n\pi t/L){\Vert }^2}& & =\frac{1}{L}{\int }_{-L}^{L}f(t)\cos \left(\frac{n\pi }{L}t\right)dt,\\ b_n& =\frac{⟨f,\sin (n\pi t/L)⟩}{\Vert \sin (n\pi t/L){\Vert }^2}& & =\frac{1}{L}{\int }_{-L}^{L}f(t)\sin \left(\frac{n\pi }{L}t\right)dt,\end{array}$$ and the inner product is $⟨f,g⟩={\int }_{-L}^{L}f(t)g(t)dt$

Periodic extension of $f:(-L,L]\to \mathbb{R}$: the function $F:\mathbb{R}\to \mathbb{R}$ given by $F(t)=f(t-2kL)$ for $t\in (-L,L]+2kL,k\in \mathbb{Z}$

Convergence of Fourier series

If the periodic extension $F$ of $f$ is piecewise smooth, then the Fourier series of $f$ converges pointwise to $[F(t^{-})+F(t^{+})]/2$ for every $t$.

Differentiation and antidifferentiation of Fourier series

Suppose that the Fourier series of $f$ is piecewise smooth.

If $f$ is continuous and $f^{\prime }$ is piecewise smooth, then the Fourier series of $f^{\prime }$ may be computed by differentiating that of $f$ term-wise.

Similarly, an antiderivative of $f$ may be computed by antidifferentiating the Fourier series of $f$ term-wise. (N.B.: In general, the resulting series will not be a Fourier series.)

Decay of Fourier series coefficients

If $f\in C^k([a,b])$, then $a_n,b_n\in \mathcal{O}(n^{-(k+2)})$ as $n\to \mathrm{\infty }$.

Parseval’s identity $$\Vert f{\Vert }^2=\frac{a_0^2}{4}+\sum _{n=1}^{\mathrm{\infty }}\frac{a_n^2+b_n^2}{2}.$$

Even extension of $f:[0,L]\to \mathbb{R}$: the function $F_{\mathrm{even}}:(-L,L]\to \mathbb{R}$ given by $F_{\mathrm{even}}(t)=\{\begin{array}{ll}f(t),& t\in [0,L]\\ f(-t),& t\in (-L,0)\end{array}$

Odd extension of $f:[0,L]\to \mathbb{R}$: the function $F_{\mathrm{odd}}:(-L,L]\to \mathbb{R}$ given by $F_{\mathrm{odd}}(t)=\{\begin{array}{ll}f(t),& t\in [0,L]\\ -f(-t),& t\in (-L,0)\end{array}$

The periodic extensions of $F_{\mathrm{even}}$ and $F_{\mathrm{odd}}$ are called the even and odd periodic extensions of $f$, respectively.

Fourier cosine series of $f\in L^2([0,L])$: the Fourier series of its even extension, viz., $$\frac{1}{2}{a}_0+\sum _{n=1}^{\mathrm{\infty }}{a}_n\cos \left(\frac{n\pi }{L}t\right),a_n=\frac{2}{L}{\int }_0^{L}f(t)\cos \left(\frac{n\pi }{L}t\right)dt$$ Fourier sine series of $f\in L^2([0,L])$: the Fourier series of its odd extension, viz., $$\sum _{n=1}^{\mathrm{\infty }}{b}_n\sin \left(\frac{n\pi }{L}t\right),b_n=\frac{2}{L}{\int }_0^{L}f(t)\sin \left(\frac{n\pi }{L}t\right)dt$$ Thus, by definition, the Fourier series of an even (resp., odd) function coincides with its Fourier cosine (resp., sine) series. As above, if the even (resp., odd) periodic extension of $f$ is piecewise smooth, then the Fourier cosine (resp., sine) series of $f$ converges pointwise a.e. to $F_{\mathrm{even}}$ (resp., $F_{\mathrm{odd}}$).

To solve a BVP of the form $x^{″}+\lambda x=f(t)$ (where $\lambda$ is not an eigenvalue of the homogeneous problem), we can expand $x$ and $f$ as Fourier sine, Fourier cosine, or Fourier series according as the BCs are Dirichlet, Neumann, or periodic. We then equate coefficients and solve for those of $x$ as with power series.

Fourier series are also effective in solving periodically forced harmonic oscillator equations. If necessary, resonant terms in the Fourier series are multiplied by $t$.

Second-order linear PDEs

Every second-order linear PDE in two independent variables $x,y$ can be written in the form $A{u}_{xx}+2B{u}_{xy}+C{u}_{yy}+\cdots +G=0$, where $A,B,C,\dots ,G$ are functions of $x$ and $y$.