Analyze This!

Consider the following chemical reaction

\underbrace{A + A + ... + A}_{n} \underset{k_2}{\stackrel{k_1}{\rightleftharpoons}} A_n

where n molecule of A combine reversibly to form A_n and, k_1, k_2 are the reaction rates.

If x, y are the concentrations of A, A_n respectively, then according to the Law of Mass Action, the reaction is governed by

\begin{cases}\frac{dx}{dt}=n k_2 y - n k_1 x^n \quad\quad(0-1)\\ \frac{dy}{dt}=k_1 x^n-k_2 y\quad\quad\quad(0-2)\\x(0)=x_0, y(0)=y_0\;\quad(0-3)\end{cases}

Without solving this initial-value problem quantitatively, the future state of system can be predicted through qualitatively analyzing how the value of (x, y) changes over the course of time.

To this end, we solve (0-1) for y first:

y=\frac{1}{n k_2}(\frac{dx}{dt} +n k_1 x^n).

Substitute it in (0-2),

\frac{d}{dt} (\frac{1}{n k_2}(\frac{dx}{dt} +n k_1 x^n)) =k_1 x^n -k_2 \cdot \frac{1}{n k_2}(\frac{dx}{dt} +n k_1 x^n).

It simplifies to

\frac{d^2x}{dt^2} + (n^2 k_1 x^{n-1} + k_2)\frac{dx}{dt} =0.



we have


Substituting p, p\frac{dp}{dx} for \frac{dx}{dt}, \frac{d^2x}{dt^2} respectively in (1-1) gives

p\frac{dp}{dx}+(n^2 k_1 x^{n-1}+k_2)p=0.

It means p=0 or


Integrate it with respect to x,

p = \frac{dx}{dt}=-n k_1 x^n - k_2 x +c_0.


f(x) = -n k_1 x^n - k_2 x +c_0,

we have

\frac{df(x)}{dx} = -n^2 k_1 x^{n-1} - k_2 < 0 \implies f(x) = -n k_1 x^{n-1} - k_2 x + c_0 is a monotonically decreasing function.

In addition, Descartes’ rule of signs reveals that

f(x)=0 has exactly one real positive root.

By definition, this root is the x_* in an equilibrium point (x_*, y_*) .

Fig. 1


As time advances, x\uparrow if x_0 < x_*. Otherwise (x_0>x_*), x\downarrow \quad\quad\quad(1-1)

Dividing (0-2) by (0-1) yields

\frac{dy}{dx} = -\frac{1}{n}.

That is,

y=-\frac{1}{n} x + c_1.

By (0-3),

c_1 = y_0 + \frac{1}{n}x_0.

And so,

y=-\frac{1}{n} x + y_0 + \frac{1}{n}x_0.

Since y is a line with a negative slope,

y is a monotonically decreasing function of x.\quad\quad\quad(1-2)

Moreover, from (0-1) and (0-2), we see that

\forall x > 0, (x, \frac{k_1}{k_2}x^n) is an equilibrium point.


All points on the curve y = \frac{k_1}{k_2}x^n in the first quadrant of x-y plane are equilibriums \quad\quad\quad(1-3).

Based on (1-1), (1-2) and (1-3), for a initial state (x_0, y_0),

x_0 < x_* \implies x\uparrow, y\downarrow.


x_0 > x_* \implies x\downarrow, y\uparrow.

Fig. 2

A phase portrait of the system is shown in Fig. 3.

Fig. 3

It shows that (x, y) on the trajectory approaches the equilibrium point (x_*, y_*) over the course of time. Namely, the system is asymptotically stable.

Qualitative Analysis of a Simple Chemical Reaction

We will study a simple chemical reaction described by

A + A \underset{k_2}{\stackrel{k_1}{\rightleftharpoons}} A_2

where two molecules of A are combined reversibly to form A_2 and, k_1, k_2 are the reaction rates.

If x is the concentration of A, y of A_2, then according to the Law of Mass Action,

\begin{cases}\frac{1}{2}\cdot\frac{dx}{dt}=k_2 y -k_1 x^2 \\ \frac{dy}{dt}=k_1 x^2-k_2 y\\x(0)=x_0, y(0)=y_0\end{cases}

or equivalently,

\begin{cases}\frac{dx}{dt}=2 k_2 y - 2 k_1 x^2 \quad\quad(0-1)\\ \frac{dy}{dt}=k_1 x^2-k_2 y\quad\quad\quad(0-2)\\x(0)=x_0, y(0)=y_0\;\quad(0-3)\end{cases}

We seek first the equilibrium points that represent the steady state of the system. They are the constant solutions where \frac{dx}{dt}=0 and \frac{dy}{dt}=0, simultaneously.

From \frac{dx}{dt}=2 k_2 y - 2 k_1 x^2=0 and \frac{dy}{dt}=k_1 x^2-k_2 y=0, it is apparent that

\forall  x \ge 0, (x_*, y_*) = (x, \frac{k_1}{k_2}x^2)\quad\quad\quad(0-4)

is an equilibrium point.

To find the value of x_*, we solve for y from (0-1),

y = \frac{1}{2k_2}(\frac{dx}{dt}+2 k_1x^2).\quad\quad\quad(0-5)

Substitute it in (0-2),


\implies  \frac{d}{dt}(\frac{dx}{dt}+2k_1x^2)=2k_2k_1x^2-k_2(\frac{dx}{dt}+2k_1x^2),


\frac{d^2x}{dt^2} + \frac{dx}{dt}(4k_1x+k_2)=0.\quad\quad\quad(0-6)

This is a 2nd order nonlinear differential equaion. Since it has no direct dependence on t, we can reduce its order by appropriate substitution of its first order derivative.



we have

\frac{d^2x}{dt^2} = \frac{d}{dt}(\frac{dx}{dt})=\frac{dp}{dt}=\frac{dp}{dx}\frac{dx}{dt}=\frac{dx}{dt}\frac{dp}{dx}=p\cdot\frac{dp}{dx}

so that (0-6) is reduced to


a 1st order differential eqution. It follows that either p=\frac{dx}{dt}=0 or \frac{dp}{dx}+(4k_1x+k_2)=0.

The second case gives


Integrate it with respect to x,

p=\frac{dx}{dt}=-2k_1x^2-k_2x + c_0.\quad\quad\quad(0-8)

Hence, the equilibrium points of (0-1) and (0-2) can be obtained by solving a quadratic equation

-2k_1x^2-k_2x + c_0=0.

Notice in order to have x_* \ge 0 as a solution, c_0 must be non-negative .

Fig. 1

The valid solution is

x_* = \frac{\sqrt{k_2^2+8c_0k_1}-k_2}{4k_1}.

Fig. 2

By (0-4),

y_* =\frac{k_1}{k_2}x_*^2=\frac{(\sqrt{k_2^2+8c_0k_1}-k2)^2}{16k_1 k_2}.

and so, the equilibrium point is

(x_*, y_*) = (\frac{\sqrt{k_2^2+8c_0k_1}-k_2}{4k_1}, \frac{(\sqrt{k_2^2+8c_0k_1}-k2)^2}{16k_1 k_2}).

Next, we turn our attentions to the phase-plan trajectories that describe the paths traced out by the (x, y) pairs over the course of time, depending on the initial values.

For (x, y) \ne (x_*, y_*), \frac{dx}{dt} \ne 0. Dividing (0-2) by (0-1) yields

\frac{\frac{dy}{dt}}{\frac{dx}{dt}}= \frac{k_1 x^2-k_2 y}{2 k_2 y- 2 k_1 x^2}=-\frac{1}{2}


\frac{dy}{dx} = -\frac{1}{2}.

Integrating it with respect to x,

y = -\frac{1}{2} x + c_1.

By (0-3),

c_1= y_0 +\frac{1}{2}x_0.


y = -\frac{1}{2}x +y_0 + \frac{1}{2}x_0.\quad\quad\quad(1-1)

Moreover, by (0-5)

y=\frac{1}{2k_2}(\frac{dx}{dt} + 2 k_1 x^2) \overset{(0-8)}{=} \frac{1}{2k_1}(-2k_1 x^2-k_2 x + x_0 +2 k_1 x^2)=\frac{1}{2k_2}(-k_2x +x_0).

As a result,

y_0 = \frac{1}{2k_2}(-k_2 x_0 + x_0).

Substitute y_0 in (1-1), we have

y = -\frac{1}{2} x + \frac{1}{2k_2}(-k_2x_0+x_0) + \frac{1}{2}x_0.

This is the trajectory of the system. Clearly, all trajectories are monotonically decreaseing lines.

At last, let us examine how the system behaves in the long run.

If x < x_* then \frac{dx}{dt}>0 (see Fig. 2) and x will increase. As a result, y will decrease. Similarly, if x > x_*, \frac{dx}{dt}<0 ensures that x will decrease. Consequently, y will increase.

Fig. 3 Trajectories and Equilibriums

It is evident that as time t advances, (x, y) on the trajectory approaches the equilibrium point (x_*, y_*).

A phase portrait of the system is illustrated in Fig. 4.

Fig. 4

It shows that the system is asymptotically stable.

An Ellipse in Its Polar Form

In this appendix to my previous post “From Dancing Planet to Kepler’s Laws“, we derive the polar form for an ellipse that has a rectangular coordinate system’s origin as one of its foci.

Fig. 1

We start with the ellipse shown in Fig. 1. Namely,



f^2=a^2-b^2\implies f< a.\quad\quad\quad(1)

After shifting the origin O to the right by f, the ellipse has the new origin O' as one of its foci (Fig. 2).

Fig. 2

Since x = x' + f, y=y', the ellipse in x'O'y' is


Substituting x', y' in (2) by

x'=r\cos(\theta), y'= r\sin(\theta)

yields equation


Replacing \sin^2(\theta), f^2 by 1-\cos^2(\theta), a^2-b^2 respectively, the equation becomes


Fig. 3

Solving (3) for r (see Fig. 3) gives

r=\frac{b^2}{f cos(\theta)+a} or r=\frac{b^2}{f \cos(\theta)-a}.

The first solution

r=\frac{b^2}{f\cos(\theta)+a} \implies r= \frac{\frac{b^2}{a}}{\frac{f}{a}\cos(\theta)+1}.


p=\frac{b^2}{a}, e=\frac{f}{a},

we have

r = \frac{p}{1+e\cdot\cos(\theta)}.

The second solution is not valid since it suggests that r < 0:

\cos(\theta)<1 \implies f\cdot\cos(\theta) < f \overset{(1)}{\implies} f\cos(\theta)<a \implies f\cos(\theta)-a<0.

From Dancing Planet to Kepler’s Laws

This most beautiful system of the sun, planets, and comets, could only proceed from the counsel and dominion of an intelligent powerful Being” Sir. Issac Newton

When I was seven years old, I had the notion that all planets dance around the sun along a wavy orbit (see Fig. 1).

Fig. 1

Many years later, I took on a challenge to show mathematically the orbit of my ‘dancing planet’ . This post is a long overdue report of my journey.

Shown in Fig. 2 is the sun and a planet in a x-y-z coordinate system. The sun is at the origin. The moving planet’s position is being described by x=x(t), y=y(t), z=z(t).

Fig. 2 r=\sqrt{x^2+y^2+z^2}, F=G\frac{Mm}{r^2}, F_z=-F\cos(c)=-F\cdot\frac{z}{r}

According to Newton’s theory, the gravitational force sun exerts on the planet is

F=-G\cdot M \cdot m \cdot \frac{1}{r^2}(\frac{x}{r},\frac{y}{r}, \frac{z}{r})=-\mu\cdot m \cdot \frac{1}{r^3}\cdot(x, y, z)

where G is the gravitational constant, M, m the mass of the sun and planet respectively. \mu = G\cdot M.

By Newton’s second law of motion,

\frac{d^2x}{dt^2} = -\mu\frac{x}{r^3},\quad\quad\quad(0-1)

\frac{d^2y}{dt^2} = -\mu\frac{y}{r^3},\quad\quad\quad(0-2)

\frac{d^2z}{dt^2} = -\mu\frac{z}{r^3}.\quad\quad\quad(0-3)

y \cdot (0-3) - z \cdot (0-2) yields

y\frac{d^2z}{dt^2}-z\frac{d^2y}{dt^2} = -\mu\frac{yz}{r^3}+ \mu\frac{yz}{r^3}=0.


y\frac{d^2z}{dt^2}-z\frac{d^2y}{dt^2} = \frac{dy}{dt}\frac{dz}{dt}+y\frac{d^2z}{dt^2}-\frac{dz}{dt}\frac{dy}{dt}-z\frac{d^2y}{dt^2}=\frac{d}{dt}(y\frac{dz}{dt}-z\frac{dy}{dt}),

it must be true that

\frac{d}{dt}(y\frac{dz}{dt}-z\frac{dy}{dt}) = 0.



where A is a constant.


z\frac{dx}{dt}-x\frac{dz}{dt}= B,\quad\quad\quad(0-5)

x\frac{dy}{dt}-y\frac{dx}{dt}= C\quad\quad\quad(0-6)

where B,C are constants.


Ax=xy\frac{dz}{dt} - xz\frac{dy}{dt},

By=yz\frac{dx}{dt} - xy\frac{dz}{dt},



Ax + By +Cz=0.\quad\quad\quad(0-7)

If C \ne 0 then by the following well known theorem in Analytic Geometry:

If A, B, C and D are constants and A, B, and C are not all zero, then the graph of the equation Ax+By+Cz+D=0 is a plane“,

(0-7) represents a plane in the x-y-z coordinate system.

For C=0, we have


It means


where k is a constant. Simply put,

y=k x.

Hence, (0-7) still represents a plane in the x-y-z coordinate system (see Fig. 3(a)).

Fig. 3

The implication is that the planet moves around the sun on a plane (see Fig. 4).

Fig. 4

By rotating the axes so that the orbit of the planet is on the x-y plane where z \equiv 0 (see Fig. 3), we simplify the equations (0-1)-(0-3) to

\begin{cases} \frac{d^2x}{dt^2}=-\mu\frac{x}{r^3} \\ \frac{d^2y}{dt^2}=-\mu\frac{y}{r^3}\end{cases}.\quad\quad\quad(1-1)

It follows that

\frac{d}{dt}((\frac{dx}{dt})^2 + (\frac{dy}{dt})^2)

= 2\frac{dx}{dt}\cdot\frac{d^2x}{dt^2} + 2 \frac{dy}{dt}\cdot\frac{d^2y}{dt^2}

\overset{(1-1)}{=}2\frac{dx}{dt}\cdot(-\mu\frac{x}{r^3})+ 2\frac{dy}{dt}\cdot(-\mu\frac{y}{r^3})

= -\frac{\mu}{r^3}\cdot(2x\frac{dx}{dt}+2y\frac{dy}{dt})

=  -\frac{\mu}{r^3}\cdot\frac{d(x^2+y^2)}{dt}

= -\frac{\mu}{r^3}\cdot\frac{dr^2}{dt}

= -\frac{\mu}{r^3} \cdot 2r \cdot \frac{dr}{dt}

= -\frac{2\mu}{r^2} \cdot \frac{dr}{dt}.


\frac{d}{dt}((\frac{dx}{dt})^2 + (\frac{dy}{dt})^2) = -\frac{2\mu}{r^2} \cdot \frac{dr}{dt}.

Integrate with respect to t,

(\frac{dx}{dt})^2+(\frac{dy}{dt})^2 = \frac{2\mu}{r} + c_1\quad\quad\quad(1-2)

where c_1 is a constant.

We can also re-write (0-6) as


where c_2 is a constant.

Using polar coordinates

\begin{cases} x= r\cos(\theta) \\ y=r\sin(\theta) \end{cases},

Fig. 5

we obtain from (1-2) and (1-3) (see Fig. 5):

(\frac{dr}{dt})^2 + (r\frac{d\theta}{dt})^2-\frac{2\mu}{r} = c_1,\quad\quad\quad(1-4)

r^2\frac{d\theta}{dt} = c_2.\quad\quad\quad(1-5)

If the speed of planet at time t is v then from Fig. 6,

Fig. 6

v = \lim\limits_{\Delta t \rightarrow 0}\frac{\Delta l}{\Delta t} = \lim\limits_{\Delta t\rightarrow 0}\frac{l_r}{\Delta t}\overset{l_r=r\Delta \theta}{=}\lim\limits_{\Delta t \rightarrow 0}\frac{r\cdot \Delta \theta}{\Delta t}=r\cdot\lim\limits_{\Delta t\rightarrow 0}\frac{\Delta \theta}{\Delta t}=r\cdot\frac{d\theta}{dt}


v = r\frac{d\theta}{dt}.\quad\quad\quad(1-6)

Suppose at t=0, the planet is at the greatest distance from the sun with r=r_0, \theta=0 and speed v_0. Then the fact that r attains maximum at t=0 implies (\frac{dr}{dt})_{t=0}=0. Therefore, by (1-4) and (1-5),

(\frac{dr}{dt})^2_{t=0} + (r\frac{d\theta}{dt})^2_{t=0}-\frac{s\mu}{r} = 0+ v_0^2-\frac{2\mu}{r}=c_1,

r (r\frac{d\theta}{dt})_{t=0}=r_0v_0=c_2.



c_2=v_0 r_0.\quad\quad\quad(1-8)

We can now express (1-4) and (1-5) as:

\frac{dr}{dt} = \pm \sqrt{c_1+\frac{2\mu}{r}-\frac{c_2^2}{r^2}},\quad\quad\quad(1-9)

\frac{d\theta}{dt} = \frac{c_2}{r^2}.\quad\quad\quad(1-10)


\rho = \frac{c_2}{r}\quad\quad\quad(1-11)


\frac{d\rho}{dr} = -\frac{c_2}{r^2},\quad\quad\quad(1-12)


By chain rule,

\frac{d\theta}{dt} = \frac{d\theta}{d\rho}\cdot\frac{d\rho}{dr}\cdot\frac{dr}{dt}.


\frac{d\theta}{d\rho} = \frac{\frac{d\theta}{dt}}{ \frac{d\rho}{dr} \cdot \frac{dr}{dt}}

\overset{(1-10), (1-12), (1-9)}{=} \frac{\frac{c_2}{r^2}}{ (-\frac{c_2}{r^2})\cdot(\pm\sqrt{c_1+\frac{2\mu}{r}-\frac{c_2^2}{r^2}}) }

\overset{(1-11)}{=} \mp\frac{1}{\sqrt{c_1-\rho^2+2\mu(\frac{\rho}{c_2})}}

= \mp\frac{1}{\sqrt{c_1+(\frac{\mu}{c_1})^2-\rho^2+2\mu(\frac{\rho}{c_2}) -(\frac{\mu}{c_2})^2}}

= \mp\frac{1}{\sqrt{c_1+(\frac{\mu}{c_1})^2-(\rho^2-2\mu(\frac{\rho}{c_2}) +(\frac{\mu}{c_2})^2)}}

= \mp\frac{1}{\sqrt{c_1+(\frac{\mu}{c_1})^2-(\rho-\frac{\mu}{c_2})^2}}.

That is,

\frac{d\theta}{d\rho} = \mp\frac{1}{\sqrt{c_1+(\frac{\mu}{c_1})^2-(\rho-\frac{\mu}{c_2})^2}}.\quad\quad\quad(1-14)



we let

\lambda = \sqrt{c_1 + (\frac{\mu}{c_2})^2}=\sqrt{(v_0-\frac{\mu}{v_0r_0})^2}=|v_0-\frac{\mu}{v_0r_0}|.

Notice that \lambda \ge 0.

By doing so, (1-14) can be expressed as

\frac{d\theta}{d\rho} =\mp \frac{1}{\sqrt{\lambda^2-(\rho-\frac{\mu}{c_2})^2}} .

Take the first case,

\frac{d\theta}{d\rho} = -\frac{1}{\sqrt{\lambda^2-(\rho-\frac{\mu}{c_2})^2}} .

Integrate it with respect to \rho gives

\theta + c = \arccos(\frac{\rho-\frac{\mu}{c_2}}{\lambda})

where c is a constant.

When r=r_0, \theta=0,

c = \arccos(1)=0 or c = \arccos(-1) = \pi.

And so,

\theta = \arccos(\frac{\rho-\frac{\mu}{c_2}}{\lambda}) or \theta+\pi = \arccos(\frac{\rho-\frac{\mu}{c_2}}{\lambda}).

For c = 0,

\lambda\cos(\theta) = \rho-\frac{\mu}{c_2}.

By (1-11), it is

\frac{c_2}{r}-\frac{\mu}{c_2} = \lambda \cos(\theta).\quad\quad\quad(1-15)

Fig. 7

Solving (1-15) for r yields

r=\frac{c_2^2}{c_2 \lambda \cos(\theta)+\mu}=\frac{\frac{c_2^2}{\mu}}{\frac{c_2}{\mu}\lambda \cos(\theta)+1}\overset{p=\frac{c_2^2}{\mu}, e=\frac{c_2 \lambda}{\mu}}{=}\frac{p}{e \cos(\theta) + 1}.


r = \frac{p}{e \cos(\theta) + 1}.\quad\quad\quad(1-16)

Studies in Analytic Geometry show that for an orbit expressed by (1-16), there are four cases to consider depend on the value of e:

We can rule out parabolic and hyperbolic orbit immediately for they are not periodic. Given the fact that a circle is a special case of an ellipse, it is fair to say:

The orbit of a planet is an ellipse with the Sun at one of the two foci.

In fact, this is what Kepler stated as his first law of planetary motion.

Fig. 8

For c=\pi,

\theta + \pi = \arccos(\frac{\rho-\frac{\mu}{c_2}}{\lambda})

from which we obtain

r=\frac{c_2^2}{c_2 \lambda \cos(\theta+\pi)+\mu}=\frac{\frac{c_2^2}{\mu}}{\frac{c_2}{\mu}\lambda\cos(\theta+\pi)+1}\overset{p=\frac{c_2^2}{\mu}, e=\frac{c_2 \lambda}{\mu}}{=}\frac{p}{e \cos(\theta+\pi) + 1}.\quad\quad(1-17)

This is an ellipse. Namely, the result of rotating (1-16) by hundred eighty degrees or assuming r attains its minimum at t=0.

The second case

\frac{d\theta}{d\rho} = +\frac{1}{\sqrt{\lambda^2-(\rho-\frac{\mu}{c_2})^2}}

can be written as

-\frac{d\theta}{d\rho} = -\frac{1}{\sqrt{\lambda^2-(\rho-\frac{\mu}{c_2})^2}} .

Integrate it with respect to \rho yields

-\theta + c = \arccos(\frac{\rho-\frac{\mu}{c_2}}{\lambda})

from which we can obtain (1-16) and (1-17) again.

Fig. 9

Over the time duration \Delta t, the area a line joining the sun and a planet sweeps an area A (see Fig. 9).

A = \int\limits_{t}^{t+\Delta t}\frac{1}{2}r\cdot v\;dt \overset{(1-6)}{=} \int\limits_{t}^{t+\Delta t}\frac{1}{2}r\cdot r\frac{d\theta}{dt}\;dt=\int\limits_{t}^{t+\Delta t}\frac{1}{2}r^2\frac{d\theta}{dt}\;dt\overset{(1-5)}{=}\int\limits_{t}^{t+\Delta t}\frac{1}{2}c_2\;dt = \frac{1}{2}c_2\Delta t.

It means

A = \frac{1}{2}c_2\Delta t\quad\quad\quad(2-1)

or that

\frac{A}{\Delta t} = \frac{1}{2}c_2

is a constant. Therefore,

A line joining the Sun and a planet sweeps out equal areas during equal intervals of time.

This is Kepler’s second law. It suggests that the speed of the planet increases as it nears the sun and decreases as it recedes from the sun (see Fig. 10).

Fig. 10

Furthermore, over the interval T, the period of the planet’s revolution around the sun, the line joining the sun and the planet sweeps the enire interior of the planet’s elliptical orbit with semi-major axis a and semi-minor axis b. Since the area enlosed by such orbit is \pi ab (see “Evaluate a Definite Integral without FTC“), setting \Delta t in (2-1) to T gives

{\frac{1}{2}c_2 T} = {\pi a b} \implies  {\frac{1}{4}c_2^2 T^2}={\pi^2 a^2 b^2} \implies T^2 = \frac{4\pi^2 a^2 b^2}{c_2^2} \implies \frac{T^2}{a^3} = \frac{4\pi^2b^2}{c_2^2a}. (3-1)

While we have p = \frac{c_2^2}{\mu} in (1-16), it is also true that for such ellipse, p=\frac{b^2}{a} (see “An Ellipse in Its Polar Form“). Hence,

\frac{b^2}{a}=\frac{c_2^2}{\mu}\implies c_2^2=\frac{\mu b^2}{a}.\quad\quad\quad(3-2)

Substituting (3-2) for c_2^2 in (3-1),

\frac{T^2}{a^3} = \frac{4\pi^2 b^2}{(\frac{\mu b^2}{a})a}=\frac{4\pi^2}{\mu} \overset{\mu=GM}{=}\frac{4\pi^2}{GM}.\quad\quad\quad(3-3)

Thus emerges Kepler’s third law of planetary motion:

The square of the orbital period of a planet is directly proportional to the cube of the semi-major axis of its orbit.

Established by (3-3) is the fact that the proportionality constant is the same for all planets orbiting around the sun.

This image has an empty alt attribute; its file name is main-qimg-e7446d7decf4bdc83f8b712109470213.jpg

A Joint Work with David Deng on CAS-Aided Analysis of Rocket Flight Performance Optimization

Research on rocket flight performance has shown that typical single-stage rockets cannot serve as the carrier vehicle for launching satellite into orbit. Instead, multi-stage rockets are used in practice with two-stage rockets being the most common. The jettisoning of stages allows decreasing the mass of the remaining rocket in order for it to accelerate rapidly till reaching its desired velocity and height.

Optimizing flight performance is a non-trivial problem in the field of rocketry. This post examines a two-stage rocket flight performance through rigorous mathematical analysis. A Computer Algebra System (CAS) is employed to carry out the symbolic computations in the process. CAS has been proven to be an efficient tool in carry out laborious mathematical calculations for decades. This post reports on the process and the results of using Omega CAS explorer, a Maxima based CAS to solve this complex problem.

A two-stage rocket consists of a payload P propelled by two stages of masses m_1 (first stage) and m_2(second stage), both with structure factor 1-e. The exhaust speed of the first stage is c_1, and of second stage c_2. The initial total mass, m_1 + m_2 is fixed. The ratio b = \frac{P}{m_1+m_2} is small.

Based on Tsiolkovsky’s equation, we derived the multi-stage rocket flight equation [1]. For a two-stage rocket, the final velocity can be calculated from the following:

v = -c_1\log(1-\frac{em_1}{m_1+m_2+P}) - c_2\log(1-\frac{em_2}{m_2+P})\quad\quad\quad(1)

Let a = \frac{m_2}{m_1+m_2}, so that (1) becomes

v=-c_1\log(1-\frac{e(1-a)}{1+b}) - c_2\log(1-\frac{ea}{a+b})\quad\quad\quad(2)

where 0<a<1, b>0, 0<e<1, c_1>0, c_2>0.

We seek an appropriate value of a that maximizes v.

Consider v as a function of a, its derivative \frac{d}{da}v is computed (see Fig. 1)

Fig. 1

We have \frac{d}{da}v = -\frac{e(abc_2e+abc_1e+a^2c_1e+b^2c_2+bc_2-b^2c_1-2abc_1-a^2c_1)}{(b+a)(ae-b-a)(ae-e+b+1)}.

That is, \frac{d}{da}v = \frac{e(abc_2e+abc_1e+a^2c_1e+b^2c_2+bc_2-b^2c_1-2abc_1-a^2c_1)}{(b+a)(b+a(1-e))(ae+b+1-e)}.

Fig. 2

As shown in Fig. 2, \frac{d}{da}v can be expressed as

\frac{d}{da}v = \frac{e(A_2a^2+B_2a+C_2)}{(b+a)(b+a(1-e))(ae+b+1-e)}\quad\quad\quad(3)

Notice that A_2 = c_1(e-1) = -c_1(1-e) < 0.

Solving \frac{d}{da}v = 0 for a gives two solutions a_1, a_2 (see Fig. 3)

Fig. 3

We rewrite the expression under the square root in a_1 and a_2 as a quadratic function of e: Ae^2+Be+ Cand compute B^2-4AC (see Fig. 4)

Fig. 4

If c_1 \ne c_2, B^2-4AC< 0. It implies that Ae^2+Be+C is positive since A>0. When c_1=c_2, Ae^2+Be+C where A>0, 0<e<1 is still positive since as a result of B^2-4AC=0, the zero point of function Ax^2+Bx+C is \frac{-B}{2A} = \frac{(8b^2+8b)c_1^2}{2(b^2c_1^2+(2b^2+4b)c_1^2+b^2c_1^2)} = 1.

The expression under the square root is positive means both a_1 and a_2 are real-valued and a_1-a_2 >0 (see Fig. 5), i.e., a_1 > a_2.

Fig. 5

From (3) where A_2  < 0, we deduce the following:

For all a \ge 0, if a>a_1 then \frac{d}{da}v(a) < 0 \quad\quad\quad(\star)

For all a \ge 0, if a_2<a<a_1 then \frac{d}{da}v(a) >0 \quad\quad\quad(\star\star)

For all a \ge 0, if a<a_2 then \frac{d}{da}v(a) < 0 \quad\quad\quad(\star\star\star)

Fig. 6

Moreover, from Fig. 6,

\frac{d}{da}v(0)\cdot \frac{d}{da}v(1) = -\frac{e^2(c_1e+bc_2-bc_1-c_1)(c_2e-bc_2-c_2+bc_1)}{b(b+1)(e-b-1)^2}=-\frac{e^2(c_1e+bc_2-bc_1-c_1)(c_2e-bc_2-c_2+bc_1)}{b(b+1)(b+1-e)^2}.

Since the expression in the numerator of \frac{d}{da}v(0)\cdot \frac{d}{da}v(1), namely

(c_1e + bc_2-bc_1-c_1)(c_2e-bc_2-c_2 +bc_1)

= (c_1(e-1-b) + bc_2)(c_2(e-1-b)+bc_1)


\ge c_1c_2(e-1-b)^2 + b^2c_1c_2+2bc_1c_2(e-1-b)



It follows that

\frac{d}{da}v(0) \cdot \frac{d}{da}v(1) < 0\quad\quad\quad(4)

The implication is that \frac{d}{da}v has at least one zero point between 0 and 1.

However, if both a_1 and a_2, the two known zero points of \frac{d}{da}v are between 0 and 1, by (\star) and (\star\star\star), \frac{d}{da}v(0)\cdot  \frac{d}{da}v(1) must be positive, which contradicts (4). Therefore, \frac{d}{da}v must have only one zero point between 0 and 1.

We will proceed to show that the only zero lies between 0 and 1 is a_1.

There are two cases to consider.

Case 1 (c_1 \le c_2) \frac{d}{da} v(0)=\frac{e(bc_2(1-e)+b^2(c_2-c1))}{b^2(b+1-e)} >0 since b>0, c_2>0, 0<e<1 and c_2-c_1>0. But this contradicts (\star\star\star). Therefore, a_2 must not be positive.

Case 2 (c_1 >  c_2) The denominator of a_2 is negative since c_1 > 0, e < 1. However, (-bc_2-bc_1)e + 2bc_1, the terms not under the square root in the numerator of a_2 can be expressed as bc_1(1-e) + b(c_1-c_2e). This is a positive expression since c_1 >  c_2 > 0, 0<e<1 implies that c_1-c_2e >  c_1 - c_1e = c_1(1-e) >0. Therefore, a_2< 0.

The fact that only a_1 lies between 0 and 1, together with (\star) and (\star\star) proves that a_1 is where the global maximum of v occurs.

a_1 can be simplified to a Taylor series expansion (see Fig. 7)

Fig. 7

The result -\frac{b((c_2+c_1)e-2c_1)}{2c_1e-2c_1} - \frac{\sqrt{b}\sqrt{c_2}}{\sqrt{c_1}} produced by CAS can be written as -\sqrt{\frac{bc_2}{c_1}} + O(b). However, it is incorrect as \lim\limits_{b\rightarrow 0} \frac{-\sqrt{\frac{bc_2}{c_1}}+O(b)}{\sqrt{b}} = -\sqrt{\frac{c_2}{c_1}} < 0 would suggest that a_1 is a negative quantity when b is small.

To obtain a correct Taylor series expansion for a_1, we rewrite a_1 as \sqrt{b^2D+bE}+bF first where




Its first order Taylor series is then computed (see Fig. 8)

Fig. 8

The first term of the result can be written as O(b). Bring the value of E into the result, we have:

a_1= \sqrt{bE} + O(b) = \sqrt{\frac{bc_2}{c_1}} + O(b)\quad\quad\quad(6)

To compute v(a_1) from (6) , we substitute \sqrt{Db^2+Eb} + Fb for a in (2) and compute its Taylor series expansion about b=0 (see Fig. 9 )

Fig. 9

Writing its first term as O(b) and substituting the value E yields:

v = -(c_1+c_2)\log(1-e)+\frac{\sqrt{b}(c_1E^{1/2}+c_2E^{-1/2})}{e-1} + O(b)

= -(c_1+c_2)\log(1-e) + \frac{\sqrt{b}(c_1e\sqrt{\frac{c_2}{c_1}}+c_2e\sqrt{\frac{c_1}{c_2}})}{e-1}+ O(b)

= -(c_1+c_2)\log(1-e)-2e\frac{\sqrt{c_1c_2b}}{1-e} + O(b)

It is positive when b is small.

We have shown the time-saving, error-reduction advantages of using CAS to aid manipulation of complex mathematical expressions. On the other hand, we also caution that just as is the cases with any software system, CAS may contain software bugs that need to be detected and weeded out with a well- trained mathematical mind.


[1] M. Xue, Viva Rockettry! Part 2

[2] Omega: A Computer Algebra System Explorer

Round Two on Finite Difference Approximations of Derivatives

There is another way to obtain the results stated in “Finite Difference Approximations of Derivatives“.

Let y'_i, y''_i denotes y'(x_i) and y''(x_i) respectively and,

h=x_{i+1}-x_i=x_i-x_{i-1} >0.

We define

y'_i = \alpha_1 y_{i-1} + \alpha_2 y_{i+1}\quad\quad\quad(1).

By Taylor’s expansion around x_i,

y_{i-1} \approx y_i + y'_i(-h)+\frac{y''_i}{2!}(-h)^2\quad\quad\quad(1-1)

y_{i+1} \approx y_i + y'_ih+\frac{y''_i}{2!}h^2\quad\quad\quad(1-2)

Substituting (1-1), (1-2) into (1),

y'_i \approx \alpha_1 y_i - \alpha y'_ih +\alpha_1\frac{y''_i}{2}h^2 + \alpha_2 y_i +\alpha_2y'_i h +\alpha_2\frac{y''_i}{2}h^2.

That is,

y'_i \approx (\alpha_1+\alpha_2) y_i +(\alpha_2-\alpha_1)h y'_i + (\alpha_1+\alpha_2)\frac{y''_i}{2}h^2.

It follows that

\begin{cases} \alpha_1+\alpha_2=0 \\ (\alpha_2-\alpha_1)h=1\\ \alpha_1+\alpha_2=0 \end{cases}\quad\quad\quad(1-3)

Fig. 1

Solving (1-3) for \alpha_1, \alpha_2 (see Fig. 1) yields

\alpha_1=-\frac{1}{2h}, \alpha_2=\frac{1}{2h}.


y'_i \approx -\frac{1}{2h} y_{i-1} + \frac{1}{2h} y_{i+1} = \frac{y_{i+1}-y_{i-1}}{2h}


\boxed{y'_i \approx \frac{y_{i+1}-y_i}{2h}}

Now, let

y''_i = \alpha_1y_{i-1}+\alpha_2 y_i + \alpha_3 y_{i+1}.


y_{i-1} \approx  y_i - y'_i h + \frac{y''_i}{2!} h^2


y_{i+1} \approx y_i +y'_i h + \frac{y''_i}{2!}h^2,

we have,

y''_i \approx \alpha_1 y_i - \alpha_1 y'_i h + \frac{\alpha_1}{2}y''_i h^2 + \alpha_2 y_i + \alpha_3 y_i + \alpha_3 y'_i h + \frac{\alpha_3}{2}y''_i h^2.

It leads to

\begin{cases} \alpha_1+\alpha_2+\alpha_3=0\\-\alpha_1+\alpha_3=0\\(\frac{\alpha_1}{2}+\frac{\alpha_3}{2})h^2=1 \end{cases}

Fig. 2

whose solution (see Fig. 2) is

\alpha_1 = \frac{1}{h^2}, \alpha_2 = -\frac{2}{h^2}, \alpha_3=\frac{1}{h^2}.


y''_i \approx \frac{1}{h^2}y_{i-1}-\frac{2}{h^2}y_i +\frac{1}{h^2}y_{i+1}=\frac{y_{i-1}-2y_i+y_{i+1}}{h^2}


\boxed{ y''_i \approx \frac{y_{i+1}-2y_i + y_{i-1}}{h^2}}

A Primer on Mathematical Epidemiology

In Memory of Dr. Li WenLiang (1986-2020)

This post is an introduction to deterministic models of infectious diseases and their Computer Algebra-Aided mathematical analysis.

We assume the followings for the simplistic SI model:

(A1-1) The population N under consideration remains a constant.

(A1-2) The population is divided into two categories: the infectious and the susceptible. Their percentages are denoted by i(t) and s(t) respectively. At t=0, i(0) = i_0.

(A1-3) The infectious’ unit time encounters with other individuals is \lambda. Upon an encounter with the infectious, the susceptible becomes infected.

When a infectious host have \lambda encounters with the population, \lambda\cdot s(t) susceptible become infected. There are N\cdot i(t) infectious in total at time t. It means that within any time interval [t, t+\Delta t], the infectious will increase by N\cdot i(t) \cdot \lambda s(t)\cdot \Delta t. i.e.,

N\cdot i(t+\Delta t) -N\cdot i(t) = N\cdot i(t)\cdot\lambda s(t)\cdot\Delta t.

Cancelling out the N‘s,

\frac{i(t+\Delta t)-i(t)}{\Delta t} = \lambda\cdot i(t) s(t)

and so

\lim\limits_{\Delta t\rightarrow 0}\frac{i(t+\Delta t)-i(t)}{\Delta t} = \lambda i(t) s(t).

That is,

\frac{d i(t)}{dt} = \lambda\cdot  i(t) s(t).

Deduce further from (A1-2) (i+s=1) is that

\begin{cases} \frac{d i}{dt} = \lambda i (1-i) \\ i(0) = i_0 \end{cases}\quad\quad\quad(1-1)

Let’s examine (1-1) qualitatively first.

We see that the SI model has two critical points:

i_1^* = 0, \quad i_2^* = 1.


\forall i \in (0, 1), \frac{d i}{dt} = i(1-i) > 0 \implies \lim\limits_{t \rightarrow \infty} i(t)=1.

This indicates that in the presence of any initial infectious hosts, the entire population will be infected in time. The rate of infection is at its peak when i = 0.5.

The qualitative results can be verified quantitatively by Omega CAS Explorer.

Fig. 1-1

From Fig. 1-1, we see that

i(t) = \frac{1}{1-(1-\frac{1}{i_0})e^{-\lambda t}}.


\lim\limits_{t \rightarrow \infty} i(t) = \lim\limits_{t \rightarrow \infty} \frac{1}{1-(1-\frac{1}{i_0})e^{-\lambda t}} = 1.

Fig. 1-2 confirms that the higher the number of initial infectious hosts(i_0), the sooner the entire population becomes infected (i = 1)

Fig. 1-2

The SI model does not take into consideration any medical practice in combating the spread of infectious disease. It is pessimistic and unrealistic.

An improved model is the SIR model. The assumptions are

(A2-1) See (A1-1)

(A2-2) See (A1-2)

(A2-3) See (A1-3)


(A2-4) Number of individuals recovered from the disease in unit time is \mu. The recovered are without immunity. They can be infected again.

By (A2-1) – (A2-4), the modified model is

N \frac{d i}{dt} = \lambda N s i - \mu N i


\begin{cases} \frac{d i}{dt} = \lambda i \cdot(1-i)-\mu\cdot i \\ i(0) = i_0 \end{cases}


\sigma = \frac{\lambda}{\mu},

we have

\begin{cases} \frac{d i}{dt} = -\lambda i \cdot(i - (1-\frac{1}{\sigma}))\\ i(0) = i_0 \end{cases}\quad\quad\quad(2-1)

The new model has two critical points:

i_1^*=0, \quad i_2^* = 1-\frac{1}{\sigma}.


\frac{d^2i}{dt^2} =\frac{d}{dt}(\frac{di}{dt})=\frac{d}{di}(-\lambda i (i-(1-\frac{1}{\sigma})))\cdot \frac{di}{dt}=-2\lambda (i-\frac{1}{2}(1-\frac{1}{\sigma}))\cdot \frac{di}{dt}


\frac{d^2i}{dt^2} = 2\lambda^2 i\cdot (i-\frac{1}{2}(1-\frac{1}{\sigma}))\cdot (i-(1-\frac{1}{\sigma}))

Without solving (2-1), we extract from it the following qualitative behavior:

Case-1 \sigma > 1

[2-1-1] \lim\limits_{t \rightarrow \infty} i(t) = 1-\frac{1}{\sigma}

[2-1-2] \forall i>1-\frac{1}{\sigma}, \frac{di}{dt} < 0.

[2-1-3] \forall i < 1-\frac{1}{\sigma}, \frac{di}{dt} > 0.

[2-1-4] i > 1-\frac{1}{\sigma}, \frac{d^2i}{dt^2}  > 0.

[2-1-5] \forall (\frac{1}{2}(1-\frac{1}{\sigma}) < i < (1-\frac{1}{\sigma})), \frac{d^2i}{dt^2} < 0.

[2-1-6] \forall(0 < i < \frac{1}{2}(1-\frac{1}{\sigma})), \frac{d^2i}{dt^2} > 0.

Case-2 \sigma < 1

[2-2-1] \lim\limits_{t\rightarrow \infty} i(t) = 0.

[2-2-2] \forall i, \frac{di}{dt} < 0.

[2-2-3] \forall i, \frac{d^2i}{dt^2} >0.

Case-3 \sigma = 1

[2-3-1] \lim\limits_{t\rightarrow \infty} i(t) = 0.

[2-3-2] \forall i, \frac{di}{dt} < 0.

[2-3-3] \forall i, \frac{d^2i}{dt^2} >0.

The cases are illustrated by solving (2-1) analytically using Omega CAS Explorer (see Fig. 2-1,2-2,2-3)

Fig. 2-1 \sigma> 1, i approaches 1-\frac{1}{\sigma} asymptotically.

Fig. 2-2 \sigma < 1, i approaches 0 asymptotically.

Fig. 2-3 \sigma = 1, i approaches 0 asymptotically.

Fig. 2-4 \sigma > 1, i‘s monotonicity depends on i_0.

Fig. 2-4 shows that for \sigma > 1, if \frac{1}{2}(1-\frac{1}{\sigma})< i_0 < 1-\frac{1}{\sigma}, then i increases on a convex curve. Otherwise, \forall i_0 < \frac{1}{2}(1-\frac{1}{\sigma}), i increases on a concave curve first. The curve turns convex after i reaches \frac{1}{2}(1-\frac{1}{\sigma}). However, \forall i_0 > (1-\frac{1}{\sigma}), i monotonically decreases along a concave curve.

Fig. 2-5 \sigma < 1, i monotonically decrease.

Fig 2-5 illustrates the case \sigma >1.

We also have:

Fig. 2-6 \sigma =1, i monotonically decrease.

From these results we may draw the following conclusion:

If \sigma >1, the monotonicity of i(t) depends on the level of i_0. Otherwise (\sigma = \frac{\lambda}{\mu} \le 1), i(t) will decrease and approach to 0 since the rate of recovery from medical treatment \mu is at least on par with the rate of infection \lambda.

This model is only valid for modeling infectious disease with no immunity such as common cold, dysentery. Those who recovered from such disease become the susceptible and can be infected again.

However, for many disease such as smallpox, measles, the recovered is immunized and therefore, falls in a category that is neither infectious nor susceptible. To model this type of disease, a new mathematical model is needed.

Enter the Kermack-McKendrick model of infectious disease with immunity.

There are three assumptions:

(A3-1) The total population N does not change.

(A3-2) Let i(t),  s(t) and r(t) denote the percentage of the infectious, susceptible and recovered respectively. At t=0, i(0) = i_0, s(0) = s_0, r(0) = r_0 = 0. The recovered are the individuals who have been infected and then recovered from the disease. They will not be infected again or transmit the infection to others.

(A3-3) \lambda is the unit time number of encounters with the infectious, \mu the unit time recoveries from the disease.

For the recovered, we have

N\frac{dr}{dt} = \mu N  i.


\begin{cases} \frac{d i(t)}{dt} = \lambda s(t) i(t) - \mu i(t),\;i(0) = i_0\quad\quad\; (3-1-1)\\ \frac{d s(t)}{dt} = -\lambda s(t)i(t),\;s(0) = s_0\quad\quad\quad\quad\;\;(3-1-2) \\ \frac{d r(t)}{dt} = \mu i(t),\;r(0)=r_0=0\quad\quad\quad\quad\;\;\;\;(3-1-3)\\i_0, s_0 \in (1, 0), i_0+s_0=1\quad\quad\quad\quad\quad\;\;\;\;(3-1-4)\end{cases}

This system of differential equations appears to defy any attempts to obtain an analytic solution (i.e., no solution can be expressed in terms of known function).

Numerical treatments for two sets of given \lambda, \mu and (s_0, i_0) are depicted in Fig. 3-1 and Fig. 3-2.

Fig. 3-1 \lambda = 0.25, \mu=0.1, i_0= 0.15, s_0=0.85

Fig. 3-2 \lambda = 0.25, \mu=0.1, i_0=0.75, s_0=0.25

However, it is only the rigorous analysis in general terms gives the correct interpretations and insights into the model.

To this end, we let

\rho = \frac{\mu}{\lambda}.

For \frac{ds}{dt} \ne 0 (i.e., i \ne 0, s \ne 0),

\frac{\frac{di}{dt}}{\frac{ds}{dt}} = \frac{di}{ds} = -1+\rho\frac{1}{s}


\frac{di}{ds} = -1+\rho\frac{1}{s}.

It has the following qualitatives:

[3-1] \forall s < \rho, \frac{di}{ds} > 0 \implies i increases.

[3-2] s = \rho, \frac{di}{ds} = 0\implies i reaches its maximum.

[3-3] \forall s>\rho, \frac{di}{ds} < 0\implies i decreases.

The analytical solution to

\begin{cases} \frac{di}{ds} = -1+\frac{\rho}{s}\\ i(s_0) = i_0\end{cases}\quad\quad\quad(3-2)

(see Fig. 3-3) is

i(s) = \rho\log(\frac{s}{s_0}) -s +i_0+s_0 \overset{[3-1-4]}{=} \rho\log(\frac{s}{s_0}) -s +1\quad\quad\quad(3-3)

Fig. 3-3

Notice that

\forall s, (s, 0) is a critical point of (3-1-1) and (3-1-2).


all points on the s-axis of the s-i phase plane are critical points of (3-1-1) and (3-1-2).

By a theorem of qualitative theory of ordinary differential equations (see Fred Brauer and John Nohel: The Qualitative Theory of Ordinary Differential Equations, p. 192, Lemma 5.2),

\forall t>0, i(t)>0\quad\quad\quad(3-4)

Moreover, let s' = s_0e^{-\frac{1}{\rho}},

i(s')=1-s' + \rho\log(\frac{s'}{s_0}) <1 + \rho\log(\frac{s'}{s_0}) = 1+ \rho\log(\frac{s_0e^{-\frac{1}{\rho}}}{s_0})=1+\rho\cdot(-\frac{1}{\rho})=0


i(s') < 0\quad\quad\quad(3-5)

Since s' = s_0e^{-\frac{1}{\rho}} = \frac{s_0}{e^{\frac{1}{\rho}}} \overset{\frac{1}{\rho}>0\implies e^{\frac{1}{\rho}}>e^0=1}{<} s_0, together, (3-5) and i(s_0) = i_0 >0 implies

\exists s_{\infty} : s' < s_{\infty}<s_0 \ni 1-s_{\infty} + \rho \log(\frac{s_{\infty}}{s_0})=0


\exists s_{\infty} : 0 \overset{s' = s_0e^{-\frac{1}{\rho}}> 0}{<} s_{\infty} < s_0 \ni 1-s_{\infty} + \rho\log(\frac{s_{\infty}}{s_0})=0.

Clearly, (s_{\infty}, 0) is a critical point of (3-1-1) and (3-1-2). Lemma 5.2 thus ensures

\forall t>0, s(t) > s_{\infty} > 0\quad\quad\quad(3-6)

It follows that

\forall t>0, \frac{ds(t)}{dt} = -\lambda i(t)\cdot s(t) \overset{(3-4), (3-6)}{<} 0\quad\quad\quad(3-7)

To the list ([3-1]-[3-3]) , we now add:

[3-4] \forall t_2> t_1>0, s(t_2) \overset{(3-7)}{<} s(t_1)

[3-5] \lim\limits_{t \rightarrow \infty} i(t) = 0,  \lim\limits_{t \rightarrow \infty} s(t) = s_{\infty}

[3-6] \forall t>0, \frac{dr(t)}{dt} = \mu i(t) \overset{(3-4)}{>}0

And so, for all (s_0, i_0) : 0 < s_0 <1, 0<i_0 < 1,

if \rho < s_0 <1 then

\forall  s_+: \rho \le s_+ < s_0, i_+ = i(s_+) \overset{[3-1], [3-4], (3-1-4)}{>} i(s_0)=i_0.


i_{max} \overset{[3-2]}{=} i(\rho) \overset{(3-3)}{=} \rho\log(\frac{\rho}{s_0})-\rho+1\overset{s_0 > \rho \implies \frac{\rho}{s_0} < 1 \implies \log(\frac{\rho}{s_0}) < 0}{<}1\quad\quad\quad(3-8),

we have

\forall s_+: 0 < s_+ \le \rho,  i_+ = i(s_+) \le i_{max} < 1.

Fig. 3-4

If s_0=\rho then,

i_+ = i(s_+ < \rho) \overset{[3-2], [3-4]}{<} i_{max} \overset{(3-8)}{<} 1.

Fig. 3-5

If 0 <s_0 < \rho then

i_+=i(s_+ < s_0) \overset{[3-1], [3-4]}{<} i(s_0)=i_0 < 1.

Fig. 3-6

In fact, for all finite t>0,

\forall s_+ : 0<s_+ <s_0, i_+ = i(s_+) <1.

Thus, the orbits of (3-1-1) and (3-1-2) have the form illustrated in Fig. 3-7.

Fig. 3-7

For example,

Fig. 3-8 \rho=0.3

What we see is that as time t advances, (i(t), s(t)) moves along the curve (3-3) in the direction of decreasing s. Consequently, if s_0 is less than \rho, then i(t) decreases monotonically to 0, and s(t) decreases to s_{\infty}. Therefore, if a small group of infectious i_0 is introduced into the population with the susceptibles s_0, with s_0 < \rho, the disease will die out rapidly. On the other hand, if s_0 is greater than \rho, then i(t) increases as s(t) decreases to \rho. Only after attaining its maximum at s=\rho, i(t) starts to decrease when the number of susceptibles falls below the threshold value \rho.

We therefore conclude:

An epidemic will occur only if the number of susceptibles in a population exceeds the threshold value \rho.

It means a larger \rho is preferred.

To increase \rho = \frac{\mu}{\lambda}, the recovery rate \mu is boosted through adequate medical care and treatment. Meanwhile, \lambda is reduced by quarantine and social distancing.

In addition to increase \rho, we can also decrease s_0 through immunizing the population.

If the number of susceptibles s_0 is initially greater than, but close to the threshold value \rho:

\rho < s_0 \approx \rho\quad\quad\quad(3-9),

and i_0 is very small compared to s_0:

i_0 \ll s_0\quad\quad\quad(3-10),

we can estimate the number of individuals who ultimately contracted the disease.

From (3-1-4), we have


and \lim\limits_{t\rightarrow \infty}(s(t) + i(t) + r(t)) yields

s_{\infty} + 0 + r_{\infty} = 1\quad\quad\quad(3-12)


s_0-s_{\infty} + i_0 - r_{\infty} =0 \overset {(3-10)}{\implies} s_0 - s_{\infty} - r_{\infty} = 0


r_{\infty} = s_0 - s_{\infty}\quad\quad\quad(3-13)

Given (3-9), we deduce from it that

1< \frac{s_0}{\rho} \approx 1 \implies 0 < \frac{s_0}{\rho}-1 \approx 0 \implies 0 <\frac{s_0-\rho}{\rho} \approx 0\implies 0 <\frac{s_0-\rho}{\rho} \ll 1


\rho < s_0 \approx \rho \implies 0<\frac{s_0-\rho}{\rho} \ll 1\quad\quad\quad(3-14)

(3-1-2) / (3-1-3) gives

\frac{\frac{ds}{dt}}{\frac{dr}{dt}}=\frac{ds}{dr} = -\frac{\lambda}{\mu}s=-\frac{1}{\rho}s


\begin{cases} \frac{ds}{dr} =-\frac{1}{\rho}s\\s(r_0=0)=s_0\end{cases}


s = s_0 e^{-\frac{r}{\rho}}\quad\quad\quad(3-15)

After substituting (3-15) in (3-1-3),

\frac{dr}{dt} = \mu(1-r-s) \overset{(3-14)}{=}\mu(1 - r -s_0 e^{-\frac{r}{\rho}})\quad\quad\quad(3-16)

In view of the fact that

\frac{r}{\rho} \overset{[3-6]}{<} \frac{s_0-s_{\infty}}{\rho} \overset{(3-9)}{\approx} \frac{s_0-s_{\infty}}{s_0}<1\quad\quad\quad(3-17),

we approximate the term e^{-\frac{r}{\rho}} in (3-16) with a Taylor expression up to the second order (see Fig. 3-9)

Fig. 3-9

The result is an approximation of equation (3-16):




It can be solved analytically (see Fig. 3-10).

Fig. 3-10

As a result,

\lim\limits_{t \rightarrow \infty} r(t) = 2\rho(1-\frac{\rho}{s_0})\quad\quad\quad(3-18)

It follows from

2\rho(1-\frac{\rho}{s_0}) = 2\rho (\frac{s_0-\rho}{s_0}) =2(s_0-\rho)(\frac{\rho}{s_0-\rho+\rho})=2(s_0-\rho)(\frac{1}{1+\frac{s_0-\rho}{\rho}}) \overset{(3-13)}{\approx} 2(s_0-\rho)

that (3-18) yields

r_{\infty} \approx 2(s_0-\rho)

Namely, the size of the epidemic is roughly 2(s_0-\rho). Consequently, by (3-13),

s_{\infty} = s_0 - 2(s_0-\rho).

The above analysis leads to the following threshold theorem of epidemiology:

(a) An epidemic occurs if and only if s_0 exceeds the threshold \rho.

(b) If \rho<s_0 \approx \rho and i_0 \ll s_0, then after the epidemic, the number of susceptible individuals is reduced by an amount approximately 2(s_0-\rho), namely, s_{\infty} \approx s_0-2(s_0-\rho).

We can also obtain (b) without solving for r(t):

From (3-3), as t\rightarrow \infty,

0 =\rho \log(\frac{s_{\infty}}{s_0})-s_{\infty} + i_0 + s_0

\overset{i_0 \ll s_0}{=} s_0 - s_{\infty} + \rho\log(\frac{s_0-s_0+s_{\infty}}{s_0})

= s_0 -s_{\infty} + \rho\log(\frac{s_0-(s_0-s_{\infty})}{s_0})

= s_0 -s_{\infty} + \rho\log(1-\frac{s_0-s_{\infty}}{s_0})

When s_0-s_{\infty} is small compared to s_0 (see (3-17)), we approximate \log(1-\frac{s_0-s_{\infty}}{s_0}) with a truncated Taylor series after two terms. Namely,

\log(1-\frac{s_0-s_{\infty}}{s_0}) = -\frac{s_0-s_{\infty}}{s_0} - \frac{1}{2}(\frac{s_0-s_{\infty}}{s_0})^2+ ...


0 = s_0-s_{\infty}+ \rho(-\frac{s_0-s_{\infty}}{s_0}-\frac{1}{2}(\frac{s_0-s_{\infty}}{s_0})^2)

=(s_0-s_{\infty})(1-\frac{\rho}{s_0} - \frac{\rho}{2s_0^2}(s_0-s_{\infty}))

Solving for s_0-s_{\infty} yields

s_0-s_{\infty} = 2s_0(\frac{s_0}{\rho}-1)

= 2(s_0-\rho+\rho)(\frac{s_0-\rho}{\rho})

= 2(s_0-\rho)(\frac{s_0-\rho+\rho}{\rho})

= 2(s_0-\rho)(1+\frac{s_0-\rho}{\rho})

\overset{(3-14)}{\approx} 2(s_0-\rho)

Exercise-1 For the Kermack-McKendrick model, show that \forall t>0,

  1. s(t)+i(t)+r(t)=1.
  2. 0<r(t)<1.