First-order ordinary differential equations

1.1 Separable ODEs

A separable first-order ODE can be written in the form

$\begin{equation}\label{eqn:first:separable} y'=g(x)h(y).\tag{1.1}\end{equation}$

We will assume that $g(x)$ is continuous over some range of values of $x$ , e.g. $x\in (a,b)$ , possibly the whole real line $(-\infty ,\infty )$ . Continuity of $g(x)$ ensures that we can take the necessary integrals.

1.1.1 A simple solution

Any (constant) solution of the equation $h(y)=0$ is also a possible solution of (1.1).

Example 1.1

Show that $y=n\pi$ ( $n\in \mathbb {Z}$ ) is a possible solution of the first-order ODE,

$y’=x\sin (y).$

Solution

Here $g(x)=x$ and $h(y)=\sin (y)$ . Now, $h(y)=\sin (y)=0$ if $y=n\pi$ (where $n\in \mathbb {Z}$ ), so the right-hand side of the ODE is zero. Clearly $y=n\pi$ (constant) $\forall x$ $\Longrightarrow$ $y'=0$ (i.e. the left-hand side is also zero). Therefore $y=n\pi$ is a possible solution.

1.1.2 Separation of variables

To derive a more general solution, we must assume that $h(y)\neq 0$ , so that we can divide each side by $h(y)$ ,

$\frac{y'}{h(y)}=g(x).$

Integrating each side with respect to $x$ ,

$\int \frac{y'}{h(y)}~ {\rm d} x=\int g(x)~ {\rm d}x.$

since $\mathrm{d}y=y'\, \mathrm{d}x$ , this implies that

$\begin{equation}\label{eqn:first:sepsolution} \int \frac{{\rm d}y}{h(y)}=\int g(x)~ {\rm d}x.\tag{1.2}\end{equation}$

Evaluating these integrals leads to a solution for $y(x)$ . Sometimes the general solution will contain solutions of $h(y)=0$ , but this is not always the case.

Example 1.2

Find the general solution to

$y’=x(y+1)^2.$

Solution

Clearly $y=-1$ is a possible constant solution. If $y\neq -1$ , we separate variables,

$\begin{eqnarray*} \int \frac{{\rm d}y}{(y+1)^2}& =& \int x~ {\rm d}x\\ \Longrightarrow -\frac{1}{y+1} & =& \frac{x^2}{2} + C, \end{eqnarray*}$

where C is an arbitrary constant of integration. Rearranging this expression,

$\begin{eqnarray*} -\frac{2}{y+1} & =& x^2+2C\\ \Longrightarrow y+1 & =& -\frac{2}{x^2+2C}\\ \Longrightarrow y & =& -1 - \left(\frac{2}{x^2+D}\right), \end{eqnarray*}$

where $D=2C$ is a rescaled constant. Alongside the special solution $y=-1$ , this is the most general solution to this first order ODE. Different values of $D$ lead to a family of solution curves (or integral curves) in the $xy$ -plane.

Check (not lectured): Differentiating this function, we see that

$y’ = \frac{4x}{\left(x^2+D\right)^2} = x(y+1)^2,$

as required.

Example 1.3

Find the general solution of the ODE $y'=(y+xy)/x^2$ and then the particular solution which passes through the point $(1,1)$ , i.e. $y(1)=1$ .

Solution

We have

$\begin{eqnarray*} y’& =& \frac{(y+xy)}{x^2} \\ & =& \frac{(1+x)y}{x^2}. \end{eqnarray*}$

Clearly $y=0$ is a possible solution, but this is incompatible with $y(1)=1$ so can be ignored. Assuming $y\neq 0$ , separation of variables leads to

$\int \frac{{\rm d}y}{y}=\int \left(\frac{1+x}{x^2}\right){\rm d}x = \int \left(\frac{1}{x^2}+\frac{1}{x}\right)\, \mathrm{d}x \Longrightarrow \ln |y|=-\frac{1}{x}+\ln |x| + C$

(where $C$ is the constant of integration). Taking the exponential of both sides,

$|y|=\exp \left[-\frac{1}{x}+\ln |x| + C\right]=|x|\exp \left[-\frac{1}{x}+C\right]= |x|\exp [C]\exp \left[-\frac{1}{x}\right],$

where $\displaystyle {\exp (x)\equiv \mathrm{e}^ x}$ . Dropping the modulus operators from both sides,

$y = \pm \, x \exp [C]\exp \left[-\frac{1}{x}\right].$

Absorbing the plus/minus sign into a new constant $D$ (where $|D|=\exp [C]$ ),

$y=D x \exp \left[-\frac{1}{x}\right].$

If $y(1)=1$ , this implies that

$1=D\cdot 1 \cdot \exp [-1] \equiv D\mathrm{e}^{-1} \Rightarrow D=\mathrm{e}.$

So the required solution is,

$y=\mathrm{e}x\exp \left[-\frac{1}{x}\right] = x \exp \left[1-\frac{1}{x}\right] \left(\equiv x \mathrm{e}^{\left[1-\left(1/x\right)\right]}\right).$

Example 1.4

If $N(t)$ is the size of a radioactive sample at time $t$ , then

$\frac{\mathrm{d}N}{\mathrm{d}t} = -\alpha N,$

where $\alpha$ is a positive constant. How long does it take for half the sample to decay?

Solution

Assuming that $N\neq 0$ , we again proceed by separating variables (note the change of notation!):

$\int \frac{{\rm d}N}{N}=-\alpha \int \, \mathrm{d}t \Longrightarrow \ln |N| = -\alpha t + C,$

where $C$ is a constant of integration. Hence,

$|N(t)| = \mathrm{e}^ C \mathrm{e}^{-\alpha t} \Longrightarrow N(t) =\pm \, \mathrm{e}^ C \mathrm{e}^{-\alpha t}\Longrightarrow N(t) = A \mathrm{e}^{-\alpha t},$

where $A$ is a new constant ( $|A|=\mathrm{e}^ C$ ) into which the plus/minus sign has been absorbed. If $N(0)=N_0$ (initial sample size) then,

$N_0 = A \mathrm{e}^0 = A \Longrightarrow N(t) = N_0 \mathrm{e}^{-\alpha t}.$

Let $t_ h$ be the length of time it takes for half the sample to decay, i.e. $N(t_ h) = N_0/2$ . This implies that,

$\frac{N_0}{2} = N(t_ h) = N_0 \mathrm{e}^{-\alpha t_ h}.$

Hence

$\begin{eqnarray*} \frac{1}{2} & =& \mathrm{e}^{-\alpha t_ h}\\ \mbox{(taking logs)}\; \Longrightarrow \ln \left[\frac{1}{2}\right] & =& -\alpha t_ h\\ \Longrightarrow \ln 1 - \ln 2 & =& -\alpha t_ h \quad \mbox{(recall that $\ln (x/y)=\ln x - \ln y)$}\\ \Longrightarrow 0 - \ln 2 & =& -\alpha t_ h\\ \Longrightarrow t_ h & =& \frac{\ln 2}{\alpha }. \end{eqnarray*}$

This is the time it takes for half the sample to decay (known as the sample’s half-life).

1.2 Homogeneous ODEs

A homogeneous first-order ODE can be written in the form

$\begin{equation}y'=f(x,y)=F\left(\frac{y}{x} \right). \label{eqn:first:homoneode}\tag{1.3}\end{equation}$

We will assume that $F(y/x)$ is a continuous function of its argument.

Idea: Define a new variable, $v(x)$ , such that

$\begin{equation}\label{eqn:first:homoneodeidea} y(x)=xv(x).\tag{1.4}\end{equation}$

Rewriting (1.3) in terms of $v(x)$ and $x$ ,

$y’=F\left(\frac{y}{x}\right) \Longrightarrow \left(xv\right)’=F(v) \Longrightarrow xv’ + v = F(v).$

Hence,

$\begin{equation}\label{eqn:first:homoneodesolution} xv' = F(v) - v,\tag{1.5}\end{equation}$

which can be solved by separating variables.

Example 1.5

Find the general solution to $y'=(y^2+2xy)/x^2$ . Then find the specific solution satisfying $y(1)=0.5$ .

Solution

Rearranging this equation,

$y’=\left(\frac{y}{x}\right)^2+2\left(\frac{y}{x}\right) \equiv F\left(\frac{y}{x}\right).$

This is a homogeneous, first-order ODE, so we set $y(x)=xv(x)$ . Following the procedure leading to (1.5),

$\begin{eqnarray*} \left(xv\right)’ & =& v^2+2v \\ \Longrightarrow xv’ + v & =& v^2 +2v\\ \Longrightarrow xv’ & =& v^2+v\\ \Longrightarrow xv’ & =& v(1+v). \end{eqnarray*}$

Separating variables and integrating both sides,

$\begin{equation*} \int \frac{{\rm d}v}{v(1+v)}=\int \frac{{\rm d}x}{x}. \end{equation*}$

Note that we have implicitly assumed that $v=y/x\neq 0$ and $v=y/x\neq -1$ . Neither of these special solutions are consistent with $y(1)=0.5$ , so this is valid.

To evaluate the integral on the left-hand side, use a partial fractions decomposition,

$\frac{1}{v(1+v)}\equiv \frac{A}{v}+\frac{B}{1+v},$

where $A$ and $B$ are constants. Hence,

$\frac{1}{v(1+v)} = \frac{A(1+v)+Bv}{v(1+v)} = \frac{A + (A+B)v}{v(1+v)}.$

For this to be true for any $v$ , we require that $A=1$ and $A+B=0 \Longrightarrow B=-A=-1$ . Therefore

$\frac{1}{v(1+v)}=\frac{1}{v}-\frac{1}{1+v}.$

Substituting this into the separated ODE,

$\begin{eqnarray*} \int \left[\frac{1}{v}-\frac{1}{1+v}\right]\, \mathrm{d}v & =& \int \frac{{\rm d}x}{x}\\ \Longrightarrow \ln |v|-\ln |1+v|& =& \ln |x|+C, \end{eqnarray*}$

where $C$ is a constant. Hence

$\ln \left|\frac{v}{1+v}\right| = \ln |x| + C \Longrightarrow \frac{v}{1+v} = \pm \exp \left[\ln |x|+C\right] = \pm \mathrm{e}^ C |x| = Dx,$

where $D$ is a new constant. Rearranging,

$\frac{v}{1+v} = Dx \Longrightarrow v = Dx(1+v) \Longrightarrow v(1 - Dx) = Dx \Longrightarrow v = \frac{Dx}{1-Dx}.$

Recalling that $v=y/x$ , this leads to the general solution

$\frac{y}{x} = \frac{Dx}{1-Dx} \Longrightarrow y = \frac{Dx^2}{1-Dx}.$

Now, given that $y(1)=0.5$ ,

$\frac{1}{2} = \frac{D}{1-D} \Longrightarrow 1-D = 2D \Longrightarrow 3D = 1 \Longrightarrow D= \frac{1}{3}.$

Hence,

$y = \frac{\frac{1}{3}x^2}{1- \frac{1}{3}x} = \frac{x^2}{3-x},$

which is the required answer.

1.3 Linear first-order ODEs

1.3.1 Terminology

Recalling our definition of a linear operator,

$L(y) = c_ n(x)y^{(n)} + c_{n-1}(x)y^{(n-1)}+\ldots + c_1(x)y’ + c_0(x)y$

(where the $c_ i(x)$ are specified functions of $x$ ), any linear first-order ODE can be written in the form

$\begin{equation}\label{eqn:first:pregenerallinear} L(y)\equiv c_1(x)y'+c_0(x)y = f(x),\tag{1.6}\end{equation}$

where $f(x)$ is another given function of $x$ . Without loss of generality, we shall assume $c_1(x)\neq 0$ .

Note

The case of $c_1(x)=0$ is trivial, i.e. $\displaystyle {c_0(x)y = f(x) \Longrightarrow y = \frac{f(x)}{c_0(x)}}$ .

Dividing (1.6) by $c_1(x)$ ,

$\begin{equation}\label{eqn:first:generallinear} y'+p(x)y = q(x),\tag{1.7}\end{equation}$

where $p(x)=c_0(x)/c_1(x)$ and $q(x)=f(x)/c_1(x)$ . This is the most general form for a (non-trivial) first-order linear ODE.

1.3.2 The integrating factor method

The challenge is to solve $y'+p(x)y = q(x)$ . We shall assume that $p(x)$ and $q(x)$ are continuous (and hence integrable) over some interval $(a,b)$ , possibly the whole real line $(-\infty ,\infty )$ . Let

$P(x) = \int p(x)\, \mathrm{d}x.$

For any $x\in (a,b)$ , the fundamental theorem of calculus tells us that

$P’(x) = p(x).$

This implies that

$\begin{eqnarray*} \frac{\mathrm{d}}{\mathrm{d}x}\left[y\mathrm{e}^{P(x)}\right] & =& y’\mathrm{e}^{P(x)} + P’(x)y \mathrm{e}^{P(x)}\\ & =& y’\mathrm{e}^{P(x)} + p(x)y \mathrm{e}^{P(x)}\\ & =& \left[y’+p(x)y\right] \mathrm{e}^{P(x)}. \end{eqnarray*}$

The quantity $\displaystyle {\mathrm{e}^{P(x)}\equiv \exp \left[P(x)\right]}$ is known as the integrating factor. As we have just shown, if we multiply our original ODE by this integrating factor, i.e.

$\left[y’+p(x)y\right] \mathrm{e}^{P(x)}=q(x) \mathrm{e}^{P(x)},$

this is equivalent to

$\begin{equation}\label{eqn:first:preintfactor} \frac{\mathrm{d}}{\mathrm{d}x}\left[y\mathrm{e}^{P(x)}\right] = q(x) \mathrm{e}^{P(x)}.\tag{1.8}\end{equation}$

Integrating this equation,

$y\mathrm{e}^{P(x)} = C + \int \left[q(x) \mathrm{e}^{P(x)}\right]\, \mathrm{d}x,$

where $C$ is a constant. Hence, recalling that $\displaystyle {P(x)=\int p(x)\, \mathrm{d}x}$ ,

$\begin{equation}\label{eqn:first:intfactor} y = C \mathrm{e}^{-P(x)} + \mathrm{e}^{-P(x)}\int \left[q(x) \mathrm{e}^{P(x)}\right]\, \mathrm{d}x.\tag{1.9}\end{equation}$

This is the general solution. The term proportional to $C$ is often referred to as the complementary function, whilst the rest of the right-hand side is known as the particular integral.

Example 1.6

Find the general solutions of the following first-order ODEs:

i)	$\displaystyle {y’+\frac{y}{x} = x}$ ;
ii)	$\displaystyle {5xy’+6y = 1+x^2}$ .

Solution

i) This equation is already in the form $y' + p(x)y = q(x)$ . Here, $p(x)=1/x$ and $q(x)=x$ . Clearly,

$P(x) = \int p(x)\, \mathrm{d}x= \int \frac{\mathrm{d}x}{x}= \ln |x|$

(note that the integrating factor is simply a computational tool, so we can ignore any constants of integration here). Hence,

$\mathrm{e}^{P(x)} = \mathrm{e}^{\ln |x|} = |x|.$

Without loss of generality, we can ignore the modulus operator here, setting $\displaystyle {\mathrm{e}^{P(x)}=x}$ (this is because a negative integrating factor would also work in the same way). Multiplying our original equation by this integrating factor, we obtain the following ODE:

$x\left(y’+\frac{y}{x}\right) = x^2 \Longrightarrow xy’ + y = x^2.$

Referrring to (1.8), we have already shown that this is equivalent to

$\frac{\mathrm{d}}{\mathrm{d}x}\left[y\mathrm{e}^{P(x)}\right] = q(x) \mathrm{e}^{P(x)} \Longrightarrow \frac{\mathrm{d}}{\mathrm{d}x}\left[xy\right] = x^2.$

Integrating this equation,

$\begin{eqnarray*} xy & =& C+ \frac{x^3}{3}\quad \mbox{(where $C$ is constant)}\\ \Longrightarrow y & =& \frac{C}{x} +\frac{x^2}{3}, \end{eqnarray*}$

which is the required result.

[Check: $\displaystyle {y’+\frac{y}{x}=-\frac{C}{x^2}+\frac{2x}{3}+\frac{C}{x^2}+\frac{x}{3}=x}$ , as required.]

ii) This equation ( $\displaystyle {5xy’+6y = 1+x^2}$ ) is not in the correct form to apply the integrating factor directly. Dividing through by $5x$ ,

$y’ + \frac{6y}{5x} = \frac{1}{5x} + \frac{x}{5},$

so $p(x)=6/5x$ and $q(x)=(1/5x)+(x/5)$ . Here,

$P(x) = \int p(x)\, \mathrm{d}x= \int \frac{\mathrm{6d}x}{5x}= \frac{6}{5}\ln |x|= \ln \left(|x|^{6/5}\right).$

Hence,

$\mathrm{e}^{P(x)} = |x|^{6/5}.$

Again we can drop the modulus operator, so that $\displaystyle {\mathrm{e}^{P(x)}=x^{6/5}}$ . Applying this integrating factor,

$x^{6/5}y’ + \frac{6yx^{1/5}}{5} = x^{6/5} \left(\frac{1}{5x} + \frac{x}{5}\right),$

which, by (1.8), is equivalent to

$\frac{\mathrm{d}}{\mathrm{d}x} \left[x^{6/5}y\right]= \frac{x^{1/5}}{5} + \frac{x^{11/5}}{5}.$

Integrating each side of this equation,

$\begin{eqnarray*} x^{6/5}y & =& C + \frac{x^{6/5}}{6}+\frac{x^{16/5}}{16}\quad \mbox{(where $C$ is constant)}\\ \Longrightarrow y & =& Cx^{-6/5} + \frac{1}{6} + \frac{x^2}{16}, \end{eqnarray*}$

which is the required solution.

[Check: $\displaystyle {5xy’+6y=-6Cx^{-6/5}+\frac{5x^2}{8}+6Cx^{-6/5}+1+\frac{3x^2}{8}=1+x^2}$ , as required.]

1.3.3 Existence and uniqueness

Theorem: If the functions $p(x)$ and $q(x)$ are continuous $\forall x \in (a,b)$ , then there exists a unique function $y(x)$ , defined on the open interval $a < x < b$ , satisfying the linear first-order ODE

$y’+p(x)y=q(x), \quad y\left(x_0\right)=y_0,$

for some given $x_0\in (a,b)$ and $y_0\in \mathbb {R}$ .

Proof: Existence has already been demonstrated by the construction of an explicit solution (1.9). For uniqueness, we assume that there are two solutions $y_1(x)$ and $y_2(x)$ . Consider the function $Y(x)=y_1(x)-y_2(x)$ . Clearly

$\begin{eqnarray*} Y’ + p(x)Y & =& \left(y_1-y_2\right)’ + p(x)\left(y_1-y_2\right)\\ & =& \left[y_1’+p(x)y_1\right] - \left[y_2’+p(x)y_2\right]\\ & =& q(x)-q(x)\\ \Longrightarrow Y’ + p(x)Y & =& 0, \end{eqnarray*}$

whilst $\displaystyle {Y\left(x_0\right) = y_1\left(x_0\right)-y_2\left(x_0\right) = y_0-y_0 = 0}$ . Obviously $Y=0$ is a possible solution of this ODE, in which case $y_1(x)=y_2(x)$ (i.e. the solution is unique). For $Y\neq 0$ , separation of variables leads to

$Y’ = -p(x)Y \iff \int \frac{\mathrm{d}Y}{Y} = -\int p(x)\, \mathrm{d}x \iff \ln |Y| = C-P(x),$

(where, as before. $\displaystyle {P(x) = \int p(x)\, \mathrm{d}x}$ and $C$ is a constant). Rearranging this expression, we see that

$Y = A\mathrm{e}^{-P(x)},$

where $A$ is a redefined constant. The only way that this can be compatible with $Y(0)=0$ is if $A=0 \Longrightarrow Y(x)=0 \; \forall x$ . This again implies that $y_1(x)=y_2(x)$ , i.e. unique solution.

Note

Linearity is crucial to this result. There are many nonlinear ODEs for which solutions are not unique. For example, consider the following differentiable function,

$y(x)=\left\{ \begin{array}{l} 0 \hspace{0.83in} \mbox{if} \hspace{0.1in} x\le C \\ (x-C)^2 \hspace{0.3in} \mbox{if} \hspace{0.1in} x > C, \end{array} \right\}$

where $C$ is any positive constant (i.e. $C\ge 0$ ). For any valid choice of $C$ (of which there are infinitely many), this function satisfies the nonlinear ODE

$y’=2\sqrt{y},$

subject to $y(0)=0$ .

1.4 Isoclines

Consider an equation of the form $y'=f(x,y)$ . The function $f(x,y)$ defines the slope of the solution curve at each point in the $xy$ -plane. Lines of constant slope, known as isoclines, satisfy

$\begin{equation}\label{eqn:first:isocline} f(x,y)=m,\tag{1.10}\end{equation}$

where $m$ is a constant ( $m\in \mathbb {R}$ ). Isoclines can be used to sketch the solution curves.

Example 1.7

Consider the following ODE:

$y’ = y+2x.$

By finding the isoclines, sketch some integral curves (i.e. solution curves) for this ODE (Exercise: show that the general solution is $\displaystyle {y=A\mathrm{e}^ x-2-2x}$ , where $A$ is a constant).

Solution

The isoclines correspond to

$y+2x =m \Longrightarrow y=m-2x.$

These are straight lines on which the slope of the solution curve is $m$ . The short line elements (or flow vectors) on the sketch indicate the corresponding direction field, i.e. the slope of the solution curves at each point in space. By connecting these, we can sketch the solution.

Example 1.8

Find the isoclines and sketch some solution curves for the following ODE,

$y’ =-\frac{x}{y}.$

Check your solution curves by finding the exact solution.

Solution

The isoclines are given by

$-\frac{x}{y} = m \Longrightarrow y = -\frac{x}{m},$

which are straight lines, of gradient $-1/m$ , passing through the origin.

Connecting up the flow vectors (see sketch) we obtain the solution curves. These look like concentric circles, centred at the origin. To find the exact solution, we can separate variables:

$\int y\, \mathrm{d}y = -\int x \, \mathrm{d}x \Longrightarrow \frac{y^2}{2} = - \frac{x^2}{2}+C,$

where $C$ is a constant of integration. Rearranging,

$y^2+x^2 = 2C,$

which describes a circle, centred at the origin in the $xy$ -plane, of radius $\sqrt{2C}$ . This is consistent with the solution curves obtained from the isoclines.

Terminology: Consider an ODE of the form

$y’ =\frac{g(x,y)}{h(x,y)},$

with isoclines defined by

$\frac{g(x,y)}{h(x,y)} = m.$

Lines of zero slope (i.e. $m=0$ ) are defined by $g(x,y)=0$ , whilst lines of infinite slope (i.e. $1/m=0$ ) are defined by $h(x,y)=0$ . If

$g(x,y)=h(x,y)=0$

at any point in the $xy$ -plane, then the point is said to be singular. At this point, the slope is not single-valued: it lies on both the $m=0$ and $1/m=0$ isoclines. In Example 1.8, all of the isoclines intersect at the origin, which is the only singular point in this case. See handout for a more complicated example.