Ordinary Differential Equations

3.1 An introduction

Basically, an ODE is any equation with derivatives in it, and the order is just the highest derivative you can find. The equations are called 'ordinary' because they involve ordinary derivatives; this distinguishes them from partial differential equations, which involve partial derivatives (second year material). For example,

$$\frac{d^{3}y}{d{x}^3}+\sin x\frac{dy}{dx}=3x^{2}y$$

is an ODE of order 3 (independent variable $x$, with $y$ a function of $x$), while

$${\left(\frac{d^{2}x}{d{t}^2}\right)}^{3/2}+\frac{dx}{dt}-tx=0$$

is an ODE of order 2. Notice how the $3/2$ power does not affect the order; only the highest derivative present matters.

Also note that the same ODE can be dressed up in several notations; the equations

$$\frac{d^{2}y}{d{x}^2}+4x\frac{dy}{dx}=e^x,f^{\prime \prime }(x)+4x{f}^{\prime }(x)=e^x,y^{\prime \prime }+4x{y}^{\prime }=e^x$$

all represent the same ODE.

Example. Consider the ODE

$$\frac{dy}{dx}=x^2+5.$$

Then $y(x)=\frac{x^3}{3}+5x$ is a solution (differentiate it and you get back $x^2+5$). But so is $y(x)=\frac{x^3}{3}+5x+6$, and $y(x)=\frac{x^3}{3}+5x-45$; each of these is called a particular solution. Every solution to this ODE can be written in the form

$$y(x)=\frac{x^3}{3}+5x+C,$$

where $C\in \mathbb{R}$; this family of solutions is called the general solution. Basically, a particular solution is one specific function, whilst the general solution is the whole family with the arbitrary constant(s) still floating around.

In the example above, $y$ could be written explicitly as a function of $x$, giving an explicit solution. This cannot always be done; sometimes we must settle for an implicit solution, where $y$ is only defined implicitly by some equation.

Example. Show that $y$, given implicitly by the equation

$$y^2=\cos (x^2+y^2),$$

is a particular solution to the ODE

$$2x\sin (x^2+y^2)+(2y\sin (x^2+y^2)+2y)\frac{dy}{dx}=0.$$

To verify that $y$ solves the ODE, we first need $\frac{dy}{dx}$. Implicitly differentiating both sides with respect to $x$ (using the chain rule on the right),

$$\begin{array}{rl}2y\frac{dy}{dx}& =-\sin (x^2+y^2)\times \left(2x+2y\frac{dy}{dx}\right)\\ 2y\frac{dy}{dx}+2y\sin (x^2+y^2)\frac{dy}{dx}& =-2x\sin (x^2+y^2)\\ 2x\sin (x^2+y^2)+(2y\sin (x^2+y^2)+2y)\frac{dy}{dx}& =0,\end{array}$$

which is exactly the ODE. Hence the ODE is satisfied and $y$ is a (particular) solution. Note we only verified a given solution here; in Section 3.5 we will actually find the general solution to this ODE from scratch.

3.2 Initial value problems

In most practical applications, we also know the value of the unknown function (and possibly its derivatives) at some particular point. This extra information, together with the ODE, forms an initial value problem.

'Initial value problem' is usually abbreviated as IVP. The strategy for solving an IVP is: find the general solution to the ODE (with its unspecified constants), then impose the initial conditions to pin down the constants. Note that an $n$​th order ODE needs $n$ conditions since the general solution has $n$ arbitrary constants.

Example. Solve the IVP

$$\frac{d^{2}y}{d{x}^2}=6x,y^{\prime }(0)=2,y(0)=-1.$$

Integrating the ODE once,

$$\frac{dy}{dx}=3x^2+C,$$

where $C\in \mathbb{R}$. Imposing $y^{\prime }(0)=2$ gives $C=2$. Integrating again,

$$y=x^3+2x+D,$$

where $D\in \mathbb{R}$. Imposing $y(0)=-1$ gives $D=-1$. Therefore the solution to the IVP is

$$y=x^3+2x-1,$$

which is valid for all $x\in \mathbb{R}$.

Not every IVP is this friendly. In general, three questions arise:

1. Does the IVP have a solution at all?
2. Is the solution unique?
3. If the initial values are given at a point $a$, how far on either side of $a$ does the solution extend?

The next two examples show that even 'simple' looking IVPs can misbehave.

Example. Solve the IVP

$$\frac{dy}{dx}=\sqrt{y},y(0)=0.$$

Assuming $y\ne 0$, we can separate:

$$\begin{array}{rl}\frac{1}{\sqrt{y}}dy& =dx\\ 2\sqrt{y}& =x+C.\end{array}$$

The initial condition gives $C=0$, and rearranging,

$$y(x)=\frac{x^2}{4}.$$

However, notice that the constant function $y(x)=0$ is also a solution to the IVP (both sides are $0$). Hence this IVP does not have a unique solution.

Example. Solve the IVP

$$\frac{dy}{dx}=\frac{1}{x},y(1)=2.$$

How far does the solution extend on either side of the point $1$?
The ODE says $\frac{dy}{dx}$ does not exist at $0$, so the solution cannot cross $0$. On the interval $(0,\infty )$ the general solution is $y(x)=\ln x+C$, and $y(1)=2$ gives $C=2$. Hence

$$y(x)=\ln x+2,x>0.$$

So the solution only extends over $(0,\infty )$; you could patch on $\ln |x|+D$ for $x<0$, but for practical applications the break in the domain at $0$ means we would not use it there.

3.3 Separable ODEs

A separable ODE is one where the two variables (say $x$ and $y$) can be separated, so that all the $y$'s end up on one side of the equation and all the $x$'s on the other. In general, a separable ODE is one that can be written in the form

$$\boxed{\frac{dy}{dx}=\frac{g(x)}{h(y)}.}$$

To solve it, write

$$h(y)dy=g(x)dx,$$

integrate both sides to obtain an implicit solution $H(y)=G(x)+C$, and then (whenever possible) isolate $y$ to get the explicit solution. If initial conditions are given, use them to determine $C$.

The fact that this symbolic manipulation with $dy$ and $dx$ 'works' traces back to the chain rule: $h(y)\frac{dy}{dx}=g(x)$, and integrating both sides with respect to $x$ turns $\int h(y)\frac{dy}{dx}dx$ into $\int h(y)dy$. Do not assume you can manipulate the symbols $dy$ and $dx$ in other creative ways and still get valid results.

Example. Solve the initial value problem

$$\frac{dy}{dx}=y^2(1+x^2),y(0)=1.$$

Separating the variables and integrating,

$$\begin{array}{rl}\frac{1}{y^2}dy& =(1+x^2)dx\\ \int \frac{1}{y^2}dy& =\int (1+x^2)dx\\ -\frac{1}{y}& =x+\frac{x^3}{3}+C,\end{array}$$

where $C\in \mathbb{R}$; this is the implicit solution. Imposing $y(0)=1$,

$$-\frac{1}{1}=0+0+C,$$

so $C=-1$. Making $y$ the subject,

$$\begin{array}{rl}-\frac{1}{y}& =x+\frac{x^3}{3}-1\\ y& =\frac{-3}{3x+x^3-3}.\end{array}$$

Therefore the explicit solution is $y=\frac{3}{3-3x-x^3}$.

Example. Solve the equation

$$\sinh y{\cos }^{2}x\frac{dy}{dx}=\tan x+4.$$

Separating and simplifying,

$$\begin{array}{rl}\sinh ydy& =\frac{\tan x+4}{{\cos }^{2}x}dx\\ \sinh ydy& =(\tan x{\sec }^{2}x+4{\sec }^{2}x)dx\\ \int \sinh ydy& =\int (\tan x{\sec }^{2}x+4{\sec }^{2}x)dx\\ \cosh y& =\frac{1}{2}{\tan }^{2}x+4\tan x+C,\end{array}$$

where $C\in \mathbb{R}$. In this case it is best to leave the solution in implicit form, since $\cosh$ is not one-to-one; 'inverting' it would silently throw away solutions.

Example. (Newton's law of cooling)
(a) Newton's law of cooling states that the rate of heat loss of a body is proportional to the difference in temperature between the body and its surroundings. Set up an ODE for this and solve it.
(b) A hot object is placed into a room of temperature $20^{\circ }$​C. The object is initially too hot for the thermometer to measure, but after 6 minutes its temperature is $80^{\circ }$​C, and after 8 minutes it is $50^{\circ }$​C. What was the original temperature?

(a) Suppose $T$ is the temperature of the object at time $t$, $A$ is the ambient (room) temperature and $k$ is the constant of proportionality. Then

$$\frac{dT}{dt}=k(T-A).$$

This is separable:

$$\begin{array}{rl}\frac{1}{T-A}dT& =kdt\\ \int \frac{1}{T-A}dT& =\int kdt\\ \ln (T-A)& =kt+C\\ T-A& =e^{kt+C}\\ T& =A+K{e}^{kt},\end{array}$$

where $K=e^C>0$.

(b) Here $A=20$, $T(6)=80$ and $T(8)=50$. Substituting into the solution,

$$\{\begin{array}{l}80=20+K{e}^{6k}\\ 50=20+K{e}^{8k}.\end{array}$$

So $K{e}^{6k}=60$ and $K{e}^{8k}=30$; dividing the second equation by the first,

$$\begin{array}{rl}{e}^{2k}& =\frac{1}{2}\\ k& =-\frac{\ln 2}{2}.\end{array}$$

Then, from $K{e}^{6k}=60$,

$$\begin{array}{rl}K& =60e^{-6k}\\ & =60e^{3\ln 2}\\ & =60\times 8\\ & =480.\end{array}$$

Finally,

$$T(0)=20+480e^0=500.$$

Therefore the initial temperature of the object was $500^{\circ }$​C.

3.4 First order linear ODEs

A first order linear ODE can be written in the form

$$\boxed{\frac{dy}{dx}+f(x)y=g(x),}$$

where $f$ and $g$ are given functions of the single variable $x$. The ODE is called linear because there are no non-linear terms (like $y^2$, $\sin y$ or $\sqrt{y^{\prime }}$) involving $y$ or $y^{\prime }$; note that $f$ and $g$ are allowed to be as ugly as they like in $x$.

The method for solving these is very slick:

1. Write the ODE in the standard form above.
2. Calculate $h(x)=e^{\int f(x)dx}$ (ignoring the constant of integration). This is called the integrating factor.
3. Multiply the ODE through by $h(x)$:

$$h(x)\frac{dy}{dx}+h(x)f(x)y=g(x)h(x).$$

By the product rule, the left-hand side always collapses to

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

4. Integrate both sides and rearrange for $y$. Do not forget the constant of integration in this step; if it is accidentally omitted then you lose half the solution.

Basically, the integrating factor is a magic multiplier engineered so that the messy left-hand side becomes the derivative of a single product $h(x)y$; this works because $h^{\prime }(x)=f(x)h(x)$ by the chain rule, which is exactly the coefficient needed for the product rule to run backwards.

Example. Solve $\frac{dy}{dx}+3y=e^{-x}$.
The ODE is already in standard form with $f(x)=3$, so the integrating factor is

$$h(x)=e^{\int 3dx}=e^{3x}.$$

Multiplying through by $e^{3x}$ and contracting the left-hand side,

$$\begin{array}{rl}{e}^{3x}\frac{dy}{dx}+3e^{3x}y& =e^{2x}\\ \frac{d}{dx}\left(e^{3x}y\right)& =e^{2x}\\ e^{3x}y& =\frac{1}{2}{e}^{2x}+C\\ y& =\frac{1}{2}{e}^{-x}+C{e}^{-3x},\end{array}$$

where $C\in \mathbb{R}$. Therefore the general solution is $y=\frac{1}{2}{e}^{-x}+C{e}^{-3x}$; notice how forgetting the $+C$ would have lost the entire $C{e}^{-3x}$ family.

Example. Solve the IVP

$$(x-1{)}^3\frac{dy}{dx}+4(x-1{)}^{2}y=x+1,y(0)=2.$$

First rewrite in standard form by dividing by $(x-1{)}^3$:

$$\frac{dy}{dx}+\frac{4}{x-1}y=\frac{x+1}{(x-1{)}^3}.$$

The integrating factor is

$$h(x)=e^{\int \frac{4}{x-1}dx}=e^{4\ln (x-1)}=e^{\ln (x-1{)}^4}=(x-1{)}^4.$$

It is important to simplify $h(x)$ like this before continuing; leaving it as $e^{4\ln (x-1)}$ makes the next steps unusable.
Multiplying the standard form by $(x-1{)}^4$,

$$\begin{array}{rl}(x-1{)}^4\frac{dy}{dx}+4(x-1{)}^{3}y& =x^2-1\\ \frac{d}{dx}((x-1{)}^{4}y)& =x^2-1\\ (x-1{)}^{4}y& =\frac{x^3}{3}-x+C.\end{array}$$

Imposing $y(0)=2$,

$$(0-1{)}^4\times 2=0-0+C,$$

so $C=2$. Therefore

$$y=\frac{\frac{x^3}{3}-x+2}{(x-1{)}^4}=\frac{x^3-3x+6}{3(x-1{)}^4}.$$

Note the solution is valid on $(-\infty ,1)$ or $(1,\infty )$ but never at $x=1$; in fact no real-valued solution of the original ODE can have $1$ in its domain.

Example. An investor has a salary of $\$60,000$ per year, expected to increase at $\$1000$ per annum. An initial deposit of $\$1000$ is invested in a program paying $8\%$ per annum, and the investor deposits $5\%$ of their salary each year. Find the amount invested after $t$ years (assuming interest and deposits are calculated continuously).
Let $y(t)$ denote the dollars invested after $t$ years. Then

$$\frac{dy}{dt}=0.08y+0.05(60000+1000t),$$

i.e. (rate of increase of investment) $=$ ($8\%$ of investment) $+$ ($5\%$ of salary). In standard form,

$$\frac{dy}{dt}-0.08y=3000+50t,$$

which is first order linear with integrating factor $h(t)=e^{\int -0.08dt}=e^{-0.08t}$. Hence

$$\frac{d}{dt}\left(e^{-0.08t}y\right)=(3000+50t)e^{-0.08t}.$$

Integrating the right-hand side by parts, with $u=3000+50t$ and $v^{\prime }=e^{-0.08t}$ (so $v=-12.5e^{-0.08t}$),

$$\begin{array}{rl}\int (3000+50t)e^{-0.08t}dt& =-12.5(3000+50t)e^{-0.08t}+\int 625e^{-0.08t}dt\\ & =-(37500+625t)e^{-0.08t}-7812.5e^{-0.08t}+C\\ & =-(625t+45312.5)e^{-0.08t}+C.\end{array}$$

So $e^{-0.08t}y=-(625t+45312.5)e^{-0.08t}+C$, and dividing through,

$$y(t)=-625t-45312.5+C{e}^{0.08t}.$$

Imposing $y(0)=1000$ gives $C=46312.5$, hence

$$y(t)=46312.5e^{0.08t}-625t-\mathrm{45312.5.}$$

Therefore after 10 years, the investment totals $y(10)\approx \$51507.86$.

Example. Some ODEs are both separable and linear; solve $\frac{dy}{dx}+2y=4$ both ways and check the answers agree.
Method 1 (separable). Write $\frac{dy}{dx}=4-2y$ and separate:

$$\begin{array}{rl}\frac{1}{4-2y}dy& =dx\\ -\frac{1}{2}\ln |4-2y|& =x+C\\ \ln |4-2y|& =-2x+C_1\\ 4-2y& =K{e}^{-2x}\\ y& =2+C_2{e}^{-2x}.\end{array}$$

Method 2 (linear). Here $f(x)=2$, so $h(x)=e^{2x}$ and

$$\begin{array}{rl}\frac{d}{dx}\left(e^{2x}y\right)& =4e^{2x}\\ e^{2x}y& =2e^{2x}+C\\ y& =2+C{e}^{-2x}.\end{array}$$

Both methods produce the same family $y=2+C{e}^{-2x}$; the arbitrary constants just get relabelled along the way. So when an ODE is both types, use whichever method you spot first (separation is usually less writing).

3.5 Exact ODEs

Here is another approach to solving (some) first order ODEs. Suppose $H$ is a function of two variables $x$ and $y$ satisfying

$$H(x,y)=C,$$

where $C$ is a real constant. If we consider $y$ as a function of $x$ and differentiate both sides with respect to $x$, then the chain rule from Chapter 1 gives

$$\frac{\partial H}{\partial x}+\frac{\partial H}{\partial y}\frac{dy}{dx}=0.$$

Denoting $\frac{\partial H}{\partial x}$ and $\frac{\partial H}{\partial y}$ by $F$ and $G$ respectively, we obtain the differential equation

$$F(x,y)+G(x,y)\frac{dy}{dx}=0.$$

Conversely, if we want to solve an ODE of this form and there happens to exist an $H$ with $F=\frac{\partial H}{\partial x}$ and $G=\frac{\partial H}{\partial y}$, then the solution is simply $H(x,y)=C$. The difficulty is that this condition on $F$ and $G$ is not easy to verify directly. Fortunately, the mixed derivative theorem says that for a 'nice' $H$,

$$\frac{{\partial }^{2}H}{\partial y\partial x}=\frac{{\partial }^{2}H}{\partial x\partial y},\text{i.e.}\frac{\partial F}{\partial y}=\frac{\partial G}{\partial x},$$

and this second condition is easy to check.

Basically: an exact ODE is one whose left-hand side is secretly the total derivative of some two-variable function $H$, and the exactness test (equality of the mixed partials) is how you detect this without knowing $H$. Note the ODE can equivalently be written as $\frac{dy}{dx}=-\frac{F(x,y)}{G(x,y)}$ or as $F(x,y)dx+G(x,y)dy=0$.

Example. Show that the differential equation

$$\frac{dy}{dx}=-\frac{2x+y+1}{2y+x+1}$$

is exact, and hence find its solution.
First rewrite the ODE as

$$(2x+y+1)+(2y+x+1)\frac{dy}{dx}=0,$$

and write $F=2x+y+1$ and $G=2y+x+1$. Then

$$\frac{\partial F}{\partial y}=1=\frac{\partial G}{\partial x},$$

so the ODE is exact. Hence there is a function $H$ with $\frac{\partial H}{\partial x}=F$ and $\frac{\partial H}{\partial y}=G$. Integrating $F$ with respect to $x$ (treating $y$ as a constant), and integrating $G$ with respect to $y$ (treating $x$ as a constant),

$$\begin{array}{rl}H(x,y)& =x^2+xy+x+C_1(y),\\ H(x,y)& =y^2+xy+y+C_2(x),\end{array}$$

where the 'constants of integration' $C_1(y)$ and $C_2(x)$ are functions of $y$ and $x$ respectively. Comparing the two expressions,

$$H(x,y)=x^2+xy+y^2+x+y,$$

and so the solution to the ODE is

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

where $C\in \mathbb{R}$. (Technically $H$ should carry its own arbitrary constant $+K$, but it just merges into $C$ on the right-hand side, so it is customary to ignore it.) Since this is a quadratic in $y$ you could use the quadratic formula to write $y$ explicitly, but the answer is much cleaner left in implicit form.

Example. Solve the differential equation

$$2x\sin (x^2+y^2)+(2y\sin (x^2+y^2)+2y)\frac{dy}{dx}=0.$$

Write $F(x,y)=2x\sin (x^2+y^2)$ and $G(x,y)=2y\sin (x^2+y^2)+2y$. Then

$$\frac{\partial F}{\partial y}=4xy\cos (x^2+y^2)=\frac{\partial G}{\partial x},$$

so the equation is exact. This time we use a slightly different approach to find $H$. Integrating $F$ with respect to $x$,

$$H(x,y)=-\cos (x^2+y^2)+C_1(y).$$

To determine $C_1(y)$, differentiate this with respect to $y$ and compare with $G$:

$$\begin{array}{rl}\frac{\partial H}{\partial y}& =2y\sin (x^2+y^2)+C_1^{\prime }(y)\\ & =2y\sin (x^2+y^2)+2y,\end{array}$$

so $C_1^{\prime }(y)=2y$, whence $C_1(y)=y^2$. Hence $H(x,y)=-\cos (x^2+y^2)+y^2$ and the solution is

$$y^2-\cos (x^2+y^2)=C,$$

where $C\in \mathbb{R}$. Notice how taking $C=0$ recovers $y^2=\cos (x^2+y^2)$, which is exactly the particular solution we verified back in Section 3.1; and here it is not possible to make $y$ explicit at all.

Example. Solve the differential equation

$$(e^x-\sin y)dx+\cos ydy=0.$$

Write $F(x,y)=e^x-\sin y$ and $G(x,y)=\cos y$. Since

$$\frac{\partial F}{\partial y}=-\cos y\text{but}\frac{\partial G}{\partial x}=0,$$

the ODE is not exact. What happens if we barge ahead with the method regardless? Integrating $F$ with respect to $x$ gives $H(x,y)=e^x-x\sin y+C_1(y)$, and then

$$\frac{\partial H}{\partial y}=-x\cos y+C_1^{\prime }(y)=\cos y$$

forces $C_1^{\prime }(y)=(1+x)\cos y$, which contradicts $C_1(y)$ being independent of $x$; no such $H$ exists.
Fortunately, not all is lost. Multiplying the whole ODE by $e^{-x}$ gives

$$(1-e^{-x}\sin y)dx+e^{-x}\cos ydy=0,$$

and since $e^{-x}$ is never zero, this has the same solutions as the original. Checking the exactness test again with $F=1-e^{-x}\sin y$ and $G=e^{-x}\cos y$,

$$\frac{\partial F}{\partial y}=-e^{-x}\cos y=\frac{\partial G}{\partial x},$$

so the new equation is exact. Integrating $F$ with respect to $x$,

$$H(x,y)=x+e^{-x}\sin y+C_1(y),$$

and comparing $\frac{\partial H}{\partial y}=e^{-x}\cos y+C_1^{\prime }(y)$ with $G$ shows $C_1^{\prime }(y)=0$. Therefore the solution is

$$x+e^{-x}\sin y=C,$$

where $C\in \mathbb{R}$. In general, finding a multiplier that fixes a non-exact ODE like this is difficult and lies beyond the scope of this course; here it was handed to us.

Example. (Classify, then solve.) For each of the following, name the method (separable, linear or exact) before touching any algebra, then solve.
(a) $\frac{dy}{dx}=x{y}^2$
(b) $x\frac{dy}{dx}+2y=x^3$, for $x>0$
(c) $(2xy+1)+(x^2+2y)\frac{dy}{dx}=0$
(d) $\frac{dy}{dx}=e^{x+y}$

(a) The right-hand side factors as (function of $x$)$\times$(function of $y$), so this is separable (it is not linear because of the $y^2$):

$$\begin{array}{rl}\frac{1}{y^2}dy& =xdx\\ -\frac{1}{y}& =\frac{x^2}{2}+C\\ y& =\frac{-2}{x^2+C_1},\end{array}$$

where $C_1=2C$. (The constant function $y=0$ is also a solution; separation quietly divided it away.)

(b) Dividing by $x$ gives $\frac{dy}{dx}+\frac{2}{x}y=x^2$, which is the linear standard form (it is not separable since the $x^2$ term stops the right side from factorising). The integrating factor is $h(x)=e^{\int \frac{2}{x}dx}=e^{2\ln x}=x^2$, so

$$\begin{array}{rl}\frac{d}{dx}\left(x^{2}y\right)& =x^4\\ x^{2}y& =\frac{x^5}{5}+C\\ y& =\frac{x^3}{5}+\frac{C}{x^2}.\end{array}$$

(c) This is in the form $F+G\frac{dy}{dx}=0$ with $F=2xy+1$ and $G=x^2+2y$; checking the exactness test, $\frac{\partial F}{\partial y}=2x=\frac{\partial G}{\partial x}$, so it is exact. Integrating $F$ with respect to $x$ gives $H=x^{2}y+x+C_1(y)$, and comparing $\frac{\partial H}{\partial y}=x^2+C_1^{\prime }(y)$ with $G$ gives $C_1^{\prime }(y)=2y$, so $C_1(y)=y^2$. The solution is

$$x^{2}y+x+y^2=C.$$

(d) This looks like none of the types, but since $e^{x+y}=e^x{e}^y$ it is separable in disguise:

$$\begin{array}{rl}{e}^{-y}dy& =e^{x}dx\\ -e^{-y}& =e^x+C\\ y& =-\ln (C_1-e^x),\end{array}$$

where $C_1=-C$ (and we need $e^x<C_1$ for the log to make sense).
The decision order I use: try to separate first; if that fails, hunt for the linear standard form; if that also fails, rearrange into $F+G{y}^{\prime }=0$ and run the exactness test.

3.6 Solving ODEs by using a change of variable [X]

This section is MATH1241/extension [X] material only, so imma keep it brief.

Not every first order ODE is separable, linear or exact, but some can be transformed into one of these types by a suitable change of variable; you solve the transformed equation for the new variable, then convert back.

Example. Use the substitution $y(x)=x\cdot v(x)$ to solve

$$\frac{dy}{dx}=\frac{xy-y^2}{x^2}.$$

By the product rule, $\frac{dy}{dx}=v+x\frac{dv}{dx}$. Substituting (and using $y=xv$ on the right),

$$\begin{array}{rl}v+x\frac{dv}{dx}& =\frac{x(xv)-(xv{)}^2}{x^2}\\ v+x\frac{dv}{dx}& =v-v^2\\ x\frac{dv}{dx}& =-v^2\\ -\frac{1}{v^2}dv& =\frac{1}{x}dx(\text{now separable})\\ \frac{1}{v}& =\ln |x|+C\\ v& =\frac{1}{\ln |x|+C}.\end{array}$$

Since $v=\frac{y}{x}$, the general solution is

$$y=\frac{x}{\ln |x|+C},$$

where $C\in \mathbb{R}$.

Example. Solve

$$\frac{dy}{dt}+2y+y^2{t}^2{e}^{2t}=0$$

using the substitution $z=1/y$.
The $y^2$ term means this is not linear in $y$, but it will turn out linear in $z$. Differentiating $y=\frac{1}{z}$ with respect to $t$ gives $\frac{dy}{dt}=-\frac{1}{z^2}\frac{dz}{dt}$, so the ODE becomes

$$\begin{array}{rl}-\frac{1}{z^2}\frac{dz}{dt}+\frac{2}{z}+\frac{t^2{e}^{2t}}{z^2}& =0\\ \frac{dz}{dt}-2z& =t^2{e}^{2t}(\text{multiplying through by }-z^2)\\ \frac{d}{dt}\left(e^{-2t}z\right)& =t^2(\text{integrating factor }{e}^{-2t})\\ e^{-2t}z& =\frac{t^3}{3}+C_0\\ z& =e^{2t}\left(\frac{t^3}{3}+C_0\right).\end{array}$$

Since $y=1/z$, the solution is

$$y=\frac{3}{e^{2t}(t^3+C)},$$

where $C=3C_0\in \mathbb{R}$. Notice how the substitution converted a nonlinear ODE in $y$ into a first order linear ODE in $z$.

3.7 Modelling with first order ODEs

Many real-life problems can be analysed by converting them into mathematics; the theoretical framework (with all its simplifying assumptions) is called a mathematical model, and its reliability is judged by how well it predicts what actually happens. To construct one, you should

1. describe accurately the data you have,
2. decide exactly what information you want out of the model,
3. decide which variables are dependent and which are independent, and
4. describe how the dependent variables change as the independent ones vary (this is usually where the differential equation appears).

We have already built two models: Newton's law of cooling and the investment example in Section 3.4.

Mixing problems

Example. A martini is, in essence, a mixture of gin and vermouth. James Blond insists his martinis be prepared as follows: initially, 40 cc of gin are placed in a large container; then gin is poured in at 2 cc/sec while vermouth is poured in at 6 cc/sec. The mixture is constantly shaken (not stirred) and flows out at 4 cc/sec.
(a) Find an expression for the volume of vermouth in the container $t$ seconds after pouring commences.
(b) James likes his martini roughly two parts gin to three parts vermouth. How many seconds should elapse before he inserts a cocktail glass into the outflow?

(a) Let $V(t)$ denote the volume of vermouth (in cc) at time $t$; we know $V(0)=0$. The key modelling equation for any mixing problem is

$$\frac{dV}{dt}=(\text{rate of inflow})-(\text{rate of outflow}).$$

The inflow of vermouth is $6$ cc/sec. For the outflow, the total volume of liquid at time $t$ is

$$40+2t+6t-4t=40+4t,$$

so the proportion (by volume) of vermouth in the container is $\frac{V(t)}{40+4t}$, and since liquid leaves at $4$ cc/sec, vermouth leaves at $\frac{4V(t)}{40+4t}=\frac{V(t)}{10+t}$ cc/sec. Hence we obtain the IVP

$$\frac{dV}{dt}=6-\frac{V}{10+t},V(0)=0.$$

This is first order linear with integrating factor

$$h(t)=e^{\int \frac{1}{10+t}dt}=e^{\ln (10+t)}=10+t,$$

so

$$\begin{array}{rl}\frac{d}{dt}((10+t)V)& =6(10+t)\\ (10+t)V& =60t+3t^2+C.\end{array}$$

Imposing $V(0)=0$ gives $C=0$, hence

$$V(t)=\frac{3t^2+60t}{10+t}.$$

(b) Two parts gin to three parts vermouth means three-fifths of the liquid is vermouth, so we require

$$\begin{array}{rl}\frac{V(t)}{40+4t}& =\frac{3}{5}\\ \frac{3t^2+60t}{4(10+t{)}^2}& =\frac{3}{5}\\ 5(3t^2+60t)& =12(10+t{)}^2\\ 15t^2+300t& =12t^2+240t+1200\\ 3t^2+60t-1200& =0\\ t^2+20t-400& =0\\ t& =-10+\sqrt{500}\approx 12.36,\end{array}$$

where we took the positive root of the quadratic. Therefore James should insert the glass about $12.36$ seconds after mixing begins.

Example. (Mixing with a changing volume.) A tank can hold 100 litres. Initially it holds 50 litres of pure water. Brine containing 2 grams of salt per litre runs in at 3 litres per minute, and the well-stirred mixture runs out at 1 litre per minute. How much salt does the tank contain at the moment it starts to overflow?
Let $x(t)$ be the mass of salt (in grams) after $t$ minutes; $x(0)=0$. The volume is not constant here: it gains $3-1=2$ litres per minute, so the volume at time $t$ is $50+2t$ litres. Salt flows in at $2\times 3=6$ g/min and flows out at (concentration)$\times$(outflow rate) $=\frac{x}{50+2t}\times 1$ g/min. Hence

$$\frac{dx}{dt}=6-\frac{x}{50+2t},x(0)=0.$$

The most common mistake in these problems is using a constant volume in the outflow concentration; whenever the inflow and outflow rates differ, the denominator must be a function of $t$.
The ODE is linear with integrating factor

$$h(t)=e^{\int \frac{1}{50+2t}dt}=e^{\frac{1}{2}\ln (50+2t)}=\sqrt{50+2t},$$

so

$$\begin{array}{rl}\frac{d}{dt}\left(\sqrt{50+2t}x\right)& =6\sqrt{50+2t}\\ \sqrt{50+2t}x& =2(50+2t{)}^{3/2}+C\\ x& =2(50+2t)+\frac{C}{\sqrt{50+2t}}.\end{array}$$

Imposing $x(0)=0$,

$$0=100+\frac{C}{\sqrt{50}},$$

so $C=-100\sqrt{50}=-500\sqrt{2}$. The tank overflows when $50+2t=100$, i.e. at $t=25$ minutes, at which point

$$\begin{array}{rl}x(25)& =200-\frac{500\sqrt{2}}{\sqrt{100}}\\ & =200-50\sqrt{2}\\ & =50(4-\sqrt{2}).\end{array}$$

Therefore the tank contains $50(4-\sqrt{2})\approx 129.3$ grams of salt at the point of overflowing.

Population models

Suppose the city of Mathopolis initially has 3,000,000 inhabitants and (initially) grows at 2% per annum. Can we predict the population in 10, 20 or 100 years? We compare three models on this same data.

Example. (Model 1: exponential growth.) Assume the rate of change of population is proportional to the population, with constant growth rate $r$:

$$\frac{dP}{dt}=rP,P(0)=P_0.$$

This is separable:

$$\begin{array}{rl}\int \frac{1}{P}dP& =\int rdt\\ \ln P& =rt+C\\ P(t)& =P_0{e}^{rt}.\end{array}$$

(This model goes back to Thomas Malthus, 1798.) For Mathopolis, $P_0=3,000,000$ and $r=0.02$, so $P(t)=3000000e^{0.02t}$; after 100 years this predicts about 22.17 million people.
Criticisms: the model predicts indefinite growth, ignoring that resources and space are finite; and external factors (disease, disasters, wars) are ignored.

Example. (Model 2: limited growth.) Now suppose there is a critical population $P_c$ which, when exceeded, causes the population to decrease; below it the population increases. Try

$$\frac{dP}{dt}=k(P_c-P),P(0)=P_0,$$

where $k>0$ (so $\frac{dP}{dt}>0$ when $P<P_c$ and $\frac{dP}{dt}<0$ when $P>P_c$). Separating (in the case $P<P_c$),

$$\begin{array}{rl}\int \frac{1}{P_c-P}dP& =\int kdt\\ -\ln (P_c-P)& =kt+C\\ P_c-P& =A{e}^{-kt}\\ P(t)& =P_c-(P_c-P_0)e^{-kt},\end{array}$$

where imposing $P(0)=P_0$ gave $A=P_c-P_0$ (the case $P>P_c$ gives exactly the same formula). Note that $P(t)\to P_c$ as $t\to \infty$.
For Mathopolis, take $P_c=7,000,000$ (the maximum sustainable size). Since $\frac{dP}{dt}=0.02P_0=60000$ at $t=0$,

$$60000=k(7000000-3000000),$$

so $k=0.015$ and

$$P(t)=7000000-4000000e^{-0.015t}.$$

Criticism: if $P$ is close to $0$ then $\frac{dP}{dt}\approx k{P}_c$, so the model says tiny populations grow fastest, which is silly.

Example. (Model 3: the logistic model.) To fix the small-population problem, try

$$\frac{dP}{dt}=kP(P_c-P),P(0)=P_0,$$

where $k>0$; now $\frac{dP}{dt}$ is small when $P$ is small (Verhulst, 1838). Separating and using partial fractions on the left,

$$\begin{array}{rl}\int \frac{1}{P(P_c-P)}dP& =\int kdt\\ \frac{1}{P_c}\int \left(\frac{1}{P}+\frac{1}{P_c-P}\right)dP& =\int kdt\\ \frac{1}{P_c}(\ln P-\ln (P_c-P))& =kt+C\\ \ln \left(\frac{P}{P_c-P}\right)& =k{P}_{c}t+C_1\\ \frac{P}{P_c-P}& =A{e}^{k{P}_{c}t},\end{array}$$

where imposing $P(0)=P_0$ gives $A=\frac{P_0}{P_c-P_0}$. Rearranging for $P$,

$$\begin{array}{rl}P& =(P_c-P)A{e}^{k{P}_{c}t}\\ P(1+A{e}^{k{P}_{c}t})& =P_{c}A{e}^{k{P}_{c}t}\\ P& =\frac{P_{c}A{e}^{k{P}_{c}t}}{1+A{e}^{k{P}_{c}t}},\end{array}$$

and substituting in $A$ then multiplying top and bottom by $(P_c-P_0)e^{-k{P}_{c}t}$,

$$\boxed{P(t)=\frac{P_c{P}_0}{P_0+(P_c-P_0)e^{-k{P}_{c}t}}.}$$

Once again $P(t)\to P_c$ as $t\to \infty$. The resulting S-shaped curve is called a logistic curve: growth is approximately exponential at the start, then slows and approaches $P_c$ asymptotically.
For Mathopolis, $60000=k{P}_0(P_c-P_0)=k\times 3000000\times 4000000$ gives $k=5\times 10^{-9}$, so $k{P}_c=0.035$ and

$$P(t)=\frac{21000000}{3+4e^{-0.035t}}.$$

Comparing the three forecasts (in millions):

All three models still ignore external factors ($P_c$ itself can shift with technology or climate). Notice how the more realistic the model, the harder the mathematics; at the extreme end, the (partial) differential equations modelling fluid flow (the Navier–Stokes equations) are so hard that existence of smooth solutions is a Millennium problem worth one million US dollars.

3.8 Second order linear ODEs with constant coefficients

We now consider a special class of second order equations, of the form

$$\frac{d^{2}y}{d{x}^2}+a\frac{dy}{dx}+by=f(x),$$

where $a$ and $b$ are real numbers. These arise naturally in wave mechanics and predator-prey models; you have already seen $\frac{d^{2}x}{d{t}^2}+n^{2}x=0$ used to model simple harmonic motion.

The homogeneous case

We first look at the case when $f(x)\equiv 0$ (i.e. $f(x)=0$ for all $x$).

It turns out we can always solve these. The first key observation:

Proof. Let $y=A{y}_1+B{y}_2$. Then

$$\begin{array}{rl}{y}^{\prime \prime }+a{y}^{\prime }+by& =(A{y}_1+B{y}_2{)}^{\prime \prime }+a(A{y}_1+B{y}_2{)}^{\prime }+b(A{y}_1+B{y}_2)\\ & =A(y_1^{\prime \prime }+a{y}_1^{\prime }+b{y}_1)+B(y_2^{\prime \prime }+a{y}_2^{\prime }+b{y}_2)(\text{regrouping by linearity of differentiation})\\ & =A\times 0+B\times 0(\text{since }{y}_1,y_2\text{ are solutions})\\ & =0.■\end{array}$$

It can also be shown that every second order ODE has at most two linearly independent solutions (proved in second year); here two solutions are linearly independent if and only if they are not constant multiples of each other. So the entire game is: find two linearly independent solutions $y_1,y_2$, and then $y=A{y}_1+B{y}_2$ is the general solution.

To hunt for a solution, we try a function that does not change too much when differentiated, since the terms on the left-hand side need to cancel each other out to give zero. If $y=e^{\lambda x}$ (with $\lambda$ constant), then $y^{\prime }=\lambda e^{\lambda x}$ and $y^{\prime \prime }={\lambda }^2{e}^{\lambda x}$, so substituting,

$$\begin{array}{rl}{\lambda }^2{e}^{\lambda x}+a\lambda e^{\lambda x}+b{e}^{\lambda x}& =0\\ {\lambda }^2+a\lambda +b& =0(\text{dividing by }{e}^{\lambda x}\ne 0).\end{array}$$