A Gift That Keeps On Giving

February 15, 2023February 17, 2023 / Michael Xue / Leave a comment

We see from “Seek-Lock-Strike!” Again that given the missile’s position $(x, y),$

$\frac{dx}{dy}=\frac{a+v_a t-x}{b-y}$

where $x$ and $y$ are themselves functions of time $t.$

It means

$\frac{dx}{dy} = \frac{\frac{dx}{dt}}{\frac{dy}{dt}}=\frac{a+v_a t-x}{b-y}\implies \frac{dx}{dt} = \frac{a+v_a t-x}{b-y}\cdot\frac{dy}{dt}.$

That is, let $\kappa(x,y,t) = \frac{a+v_a t-x}{b-y},$

$\frac{dx}{dt}=\kappa\cdot \frac{dy}{dt}.\quad\quad\quad(1)$

We also have (see “Seek-Lock-Strike!”)

$v_m = \sqrt{(\frac{dx}{dt})^2+(\frac{dy}{dt})^2}\implies v_m = \sqrt{1+\kappa^2}|\frac{dy}{dt}| .$

Since $\frac{dy}{dt} >0,$

$\frac{dy}{dt} = \frac{v_m}{\sqrt{1+\kappa^2}}.\quad\quad\quad(2)$

Substitute (2) into (1) yields

$\frac{dx}{dt} = \frac{v_m\cdot \kappa}{\sqrt{1+\kappa^2}}.\quad\quad\quad(3)$

It follows that $(x(t), y(t))$ , the position of the missile satisfies the initial-value problem

$\begin{cases} \frac{dx}{dt} = \frac{\kappa\cdot v_m}{\sqrt{1+\kappa^2}} \\ \frac{dy}{dt} = \frac{v_m}{\sqrt{1+\kappa^2}} \\x(0)=0, y(0)=0\end{cases}\quad\quad\quad(4)$

To obtain the missile’s trajectory, we solve (4) numerically using the Runge-Kutta algorithm. It integrates (4) from $t=0$ to $t= t_*$ (see “Seek-Lock-Strike!”).

Fig. 1 $a=100\;m, b=3000\;m, v_a=1500\;ms^{-1},v_m=2000\;ms^{-1}$

The missile strike is illustrated in Fig. 1 and 2.

Fig. 2 $a=100\;m, b=3000\;m, v_a=1500\;ms^{-1},v_m=2000\;ms^{-1}$

Fig. 3 $a=-1200\;m, b=3000\;m, v_a=1500\;ms^{-1},v_m=2000\;ms^{-1}$

The trajectories shown are much smoother than those in “Seek-Lock-Strike!” Animated.

“Seek-Lock-Strike!” Animated

February 1, 2023February 15, 2023 / Michael Xue / Leave a comment

In “Seek-Lock-Strike!” Again, we obtained the missile’s trajectory. Namely,

$x = \frac{1}{2} \left(\frac{(b-y)^{r+1}}{b^r \cdot f \cdot (r+1)}-\frac{b^r \cdot f \cdot (b-y)^{1-r}}{1-r}\right) -k\quad\quad\quad(*)$

where

$f = \frac{a}{b}+\sqrt{1+(\frac{a}{b})^2},$

$k = \frac{b}{2}\left(\frac{1}{f \cdot (r+1)}-\frac{f}{1-r}\right).$

Since the fighter jet maintains its altitude ( $y= b$ ), the missile must strike it at $(x_*, b)$ . Setting $y=b$ in $(*)$ gives $x_* = -k.$

Hence, we can plot $(*):$

Fig. 1 $a = 100\;m, b=3000\;m, v_a=1500\;ms^{-1}, v_m=2000\;ms^{-1}$

We can also illustrate “Seek-Lock-Strike” in an animation:

Fig. 2 $a = 100\;m, b=3000\;m, v_a=1500\;ms^{-1}, v_m=2000\;ms^{-1}$

Fig. 3 $a = -2500\;m, b=3000\;m, v_a=1500\;ms^{-1}, v_m=2000\;ms^{-1}$

“Seek-Lock-Strike!” Again

January 23, 2023February 9, 2023 / Michael Xue / 1 Comment

We can derive a different governing equation for the missile in “Seek-Lock-Strike!“.

Fig. 1

Looking from a different viewpoint (Fig. 1), we see

$\frac{dx}{dy} = \frac{a+v_at-x}{b-y}.\quad\quad\quad(1)$

Solving (1) for $t$ ,

$t = -\frac{\frac{dx}{dy}y-b\frac{dx}{dy}-x+a}{v_a}.\quad\quad\quad(2)$

We also have

$v_m t = \int\limits_{0}^{y}\sqrt{1+(\frac{dx}{dy})^2} \implies t = \frac{\int\limits_{0}^{y}\sqrt{1+(\frac{dx}{dy})^2}\;dy}{v_m}.\quad\quad\quad(3)$

Equate (1) and (2) gives

$-\frac{\frac{dx}{dy}y-b\frac{dx}{dy}-x+a}{v_a}-\frac{\int\limits_{0}^{y}\sqrt{1+(\frac{dx}{dy})^2}\;dy}{v_m} = 0.\quad\quad\quad(4)$

The governing eqaution emerges after differentiate (4) with respect to $x:$

$-\frac{d^2x}{dy^2}y+b\frac{d^2x}{dy^2}-\frac{v_a\sqrt{1+(\frac{dx}{dy})^2}}{v_m} = 0.\quad\quad\quad(5)$

We let $p = \frac{dx}{dy}$ so $\frac{d^2x}{dy^2} = \frac{d}{dy}\left(\frac{dx}{dy}\right)=\frac{dp}{dy}$ and express (5) as

$\frac{dp}{dy}(b-y) -r\sqrt{1+(\frac{dx}{dy})^2} = 0\quad\quad\quad(*)$

where $r = \frac{v_a}{v_m}.$

Fig. 2

Using Omega CAS Explorer, we compute the missile’s striking time $t_*$ (see Fig. 3). It agrees with the result obtained previously.

Fig. 3

Exercise-1 Obtain the missile’s trajectory from (*).

“Seek-Lock-Strike!” Simplified

January 12, 2023January 19, 2023 / Michael Xue / Leave a comment

There is an easier way to derive the governing equation ((5), “Seek-Lock-Strike!“) for the missile.

Solving

$\frac{dy}{dx}(a+v_at-x) = b-y$

for $t,$ we have

$t = \frac{(x-a)\frac{dy}{dx}-y+b}{v_a\frac{dy}{dx}}.\quad\quad\quad(1)$

From

$v_m t = \int\limits_{0}^{x}\sqrt{1+(\frac{dy}{dx})^2}\;dx,$

we also have

$t = \frac{\int\limits_{0}^{x}\sqrt{1+(\frac{dy}{dx})^2}\;dx}{v_m}.\quad\quad\quad(2)$

Equate (1) and (2) gives

$\frac{\int\limits_{0}^{x}\sqrt{1+(\frac{dy}{dx})^2}\;dx}{v_m}-\frac{(x-a)\frac{dy}{dx}-y+b}{v_a\frac{dy}{dx}}=0.\quad\quad\quad(3)$

Differentiate (3) with repect to $x,$ we obtain

$(bv_m-v_my)\frac{d^2y}{dx^2} + v_a(\frac{dy}{dx})^2\sqrt{1+(\frac{dy}{dx})^2} = 0.\quad\quad\quad(4)$

(see Fig. 1)

Fig. 1

Let $r=\frac{v_a}{v_m}, p=\frac{dy}{dx} \implies \frac{d^2y}{dx^2} = p\frac{dp}{dy},$ (4) bcomes

$(b-y)p\frac{dp}{dy} + r p^2 \sqrt{1+p^2} = 0.$

Since $p = \frac{dy}{dx} \ne 0$ , dividing $p$ through yields

$(b-y)\frac{dp}{dy} + rp\sqrt{1+p^2}=0,$

the governing equation for the missile.

Seek-Lock-Strike!

November 22, 2022February 7, 2023 / Michael Xue / Leave a comment

A guided missile is launched to destroy a fighter jet (Fig. 1).

Fig. 1

We introduce a coordinate axes such that at $t=0,$ the missile is at origin $(0, 0)$ and the jet at $(a,b)$ . The jet flies parallel to the x-axis with constant speed $v_a$ . The missile has locked onto the jet so it is always pointing at the jet as it moves. Its speed is $v_m.$

Find the time and position the missile strikes its target.

“The missile has locked onto the jet so it is always pointing at the jet as it moves” means that the tangent to the missile’s path at any point $(x, y)$ will pass through the position of the jet. The equation of the tangent is

$\frac{dy}{dx}\cdot(a + v_a t -x) = b-y.\quad\quad\quad(1)$

In addition,

$(\frac{dx}{dt})^2 + (\frac{dy}{dt})^2=v_m^2.\quad\quad\quad(2)$

Notice $x, y$ are functions of time $t : x=x(t), y=y(t).$

Differentiate (1) with respecte to $t:$

$\frac{d}{dt}(\frac{dy}{dx})\cdot(a+v_at-x) + \frac{dy}{dx}\cdot\frac{d}{dt}(a+v_at-x) = -\frac{dy}{dt},$

$\frac{d}{dt}(\frac{dy}{dx}) = \frac{d}{dx}(\frac{dy}{dx})\cdot\frac{dx}{dt}=\frac{d^2y}{dx^2}\cdot\frac{dx}{dt}$

$\left(\frac{d^2y}{dx^2}\cdot\frac{dx}{dt}\right)\cdot(a+v_at-x) + \frac{dy}{dx}\cdot(v_a -\frac{dx}{dt}) = -\frac{dy}{dt},$

$\left(\frac{d^2y}{dx^2}\cdot\frac{dx}{dt}\right)\cdot(a+v_at-x) + \frac{dy}{dx}\cdot v_a -\underbrace{\frac{dy}{dx}\cdot\frac{dx}{dt}}_{\frac{dy}{dt}}= -\frac{dy}{dt},$

$\frac{d^2y}{dx^2}\cdot\frac{dx}{dt}\cdot(a+v_at-x) + \frac{dy}{dx}\cdot v_a = 0,$

$\frac{d^2y}{dx^2}= \frac{d}{dx}(\frac{dy}{dx})=\frac{d}{dy}(\frac{dy}{dx})\cdot\frac{dy}{dx}$

$\frac{d}{dy}(\frac{dy}{dx})\cdot\frac{dy}{dx}\cdot\frac{dx}{dt}(a+v_at-x) + \frac{dy}{dx}v_a=0,$

$\frac{d}{dy}(\frac{dy}{dx})\cdot\frac{dx}{dt}\cdot\underline{\frac{dy}{dx}(a+v_at-x)} + \frac{dy}{dx}v_a=0.$

By (1), substituting $b-y$ for $\frac{dy}{dx}\cdot(a + v_a t -x)$ ,

$\frac{d}{dy}(\frac{dy}{dx})\cdot\frac{dx}{dt}\cdot(b-y) + \frac{dy}{dx}v_a=0.\quad\quad\quad(3)$

Let $\frac{dy}{dt} = \frac{dy}{dx}\cdot\frac{dx}{dt},$ we express (2) as

$(\frac{dx}{dt})^2 + (\frac{dy}{dx}\cdot\frac{dx}{dt})^2=v_m^2,$

$(\frac{dx}{dt})^2\cdot(1+(\frac{dy}{dx})^2)= v_m^2.$

Suppose $a>0$ , which implies that $\frac{dx}{dt} > 0$ (see Exercise-4). Solving for $\frac{dx}{dt}$ gives

$\frac{dx}{dt} = \frac{v_m}{\sqrt{1+(\frac{dy}{dx})^2}}.\quad\quad\quad(4)$

Submitting (4) into (3),

$\frac{d}{dy}(\frac{dy}{dx})\cdot\frac{v_m}{\sqrt{1+(\frac{dy}{dx})^2}}\cdot(b-y) + \frac{dy}{dx}\cdot v_a=0,$

$\frac{d}{dy}(\frac{dy}{dx})\cdot(b-y) + \frac{v_a}{v_m}\cdot\frac{dy}{dx}\cdot\sqrt{1+(\frac{dy}{dx})^2}=0,\quad\quad\quad(5)$

Let $p = \frac{dy}{dx}, r=\frac{v_a}{v_m},$

(5) becomes

$\frac{dp}{dy}(b-y) + r\cdot p\cdot\sqrt{1+p^2}=0.$

We solve this non-linear differential equation as follows:

For $y < b$ ,

$\frac{1}{p\sqrt{1+p^2}}\frac{dp}{dy} = \frac{-r}{b-y},$

$\int (\frac{\sqrt{1+p^2}}{p}-\frac{p}{\sqrt{1+p^2}})\frac{dp}{dy} dy = \int \frac{-r}{b-y} dy,$

$\int\frac{\sqrt{1+p^2}}{p}dp-\int\frac{p}{\sqrt{1+p^2}}dp = r\log(b-y)+k_2.\quad\quad\quad(6)$

Since $\int\frac{\sqrt{1+p^2}}{p}dp=\sqrt{1+p^2} - \mathrm{arcsinh}\left(\frac{1}{|p|}\right)$ (See “I vs. CAS“), $\int\frac{p}{\sqrt{1+p^2}}dp=\sqrt{1+p^2},$

(6) gives

$-\mathrm{arcsinh}\left(\frac{1}{|p|}\right) = r\log(b-y)+k_1.\quad\quad\quad(7)$

At $t=0, y=0, p=\frac{dy}{dx}=\frac{b}{a},$

(7) yields

$k_1 = -\mathrm{arcsinh}\left(\frac{a}{b}\right) - r\log(b).$

Moreover, since $p = \frac{dy}{dx} > 0$ , we have

$\frac{dy}{dx} = \frac{1}{-\sinh(r\log(b-y)+k_1)},$

$-\sinh(r\log(b-y)+k_1)\cdot \frac{dy}{dx} = 1,$

$\displaystyle\int -\sinh(r\log(b-y) + k_1) \cdot \frac{dy}{dx}\;dx = \displaystyle \int 1\; dx = x +k_2.$

Using Omega CAS Explorer:

we obtain

$\frac{1}{2}\left(\frac{e^{k_1}(b-y)^{r+1}}{r+1} - \frac{e^{-k_1}(b-y)^{1-r}}{1-r}\right) = x + k_2.\quad\quad\quad(8)$

Since $y=0, x=0$ , (8) gives

$k_2 = \frac{1}{2}\left(\frac{e^{k_1}b^{r+1}}{r+1} - \frac{e^{-k_1}b^{1-r}}{1-r}\right).\quad\quad\quad(9)$

Suppose when $t = t_*,$ the missile hits the target. Then the striking coordinates $(x_*, y_*)$ are the same as that of the fighter jet, i.e.,

$x_* = a + v_a t_*, \;\;y_*=b.\quad\quad\quad(10)$

Hence,

$\lim\limits_{t \rightarrow t_*} x = x_* \overset{(10)}{=} a+v_a t_*, \;\;\lim\limits_{t \rightarrow t_*} y = y_* \overset{(10)}{=}b \implies \lim\limits_{t \rightarrow t_*} b-y = 0.\quad\quad\quad(11)$

For $r \le 1$ (i.e., $v_a \le v_m$ ), as $t \rightarrow t_*,$ (8) gives (see Exercise-1)

$0 = a +v_a t_* + k_2 \implies -k_2-a = v_a t_*.$

That is,

$-\underbrace{\frac{1}{2}\left(\frac{e^{k_1}b^{1+r}}{1+r}-\frac{e^{-k_1}b^{1-r}}{1-r}\right)}_{(9):\;\;k_2}-a=v_a t_*.$

By $(\star)$ (see below),

$-\frac{1}{2}\left(\frac{1}{b^r\left(\frac{a}{b}+\sqrt{1+(\frac{a}{b})^2}\right)}\cdot\frac{b^{1+r}}{1+r}-b^r\left(\frac{a}{b}+\sqrt{1+(\frac{a}{b})^2}\right)\cdot\frac{b^{1-r}}{1-r}\right) - a = v_a t_*,$

$-\frac{1}{2}\left(\frac{1}{\left(\frac{a}{b}+\sqrt{1+(\frac{a}{b})^2}\right)}\cdot\frac{b}{1+r}-\left(\frac{a}{b}+\sqrt{1+(\frac{a}{b})^2}\right)\cdot\frac{b}{1-r}\right) - a = v_a t_*,$

$\frac{\frac{-b}{\left(\frac{a}{b}+\sqrt{1+(\frac{a}{b})^2}\right)(1+r)} + \frac{\left(\frac{a}{b}+\sqrt{1+(\frac{a}{b})^2}\right)b}{1-r}}{2} -a = v_a t_*,$

$\frac{\frac{-b\left(\frac{a}{b}-\sqrt{1+(\frac{a}{b})^2}\right)}{-1(1+r)} + \frac{\left(\frac{a}{b}+\sqrt{1+(\frac{a}{b})^2}\right)b}{1-r}}{2} -a = v_a t_*,$

$\frac{\frac{b\left(\frac{a}{b}-\sqrt{1+(\frac{a}{b})^2}\right)}{1+r} + \frac{\left(\frac{a}{b}+\sqrt{1+(\frac{a}{b})^2}\right)b}{1-r}}{2} -a = v_a t_*,$

$\frac{b\left(\frac{a}{b}-\sqrt{1+(\frac{a}{b})^2}\right)}{2(1+r)} + \frac{\left(\frac{a}{b}+\sqrt{1+(\frac{a}{b})^2}\right)b}{2(1-r)} -a = v_a t_*,$

$\frac{a-\sqrt{b^2+a^2}}{2(1+r)} + \frac{a+\sqrt{b^2+a^2}}{2(1-r)} -a = v_a t_*,$

$\frac{(1-r)(a-\sqrt{b^2+a^2}) +(1+r)(a+\sqrt{b^2+a^2})-2a(1-r^2)}{2(1+r)(1-r)} = v_a t_*,$

$\frac{a-\sqrt{b^2+a^2}-ar+r\sqrt{b^2+a^2}+ar+r\sqrt{b^2+a^2}+a+\sqrt{b^2+a^2}-2a + 2ar^2}{2(1-r^2)}=v_a t_*,$

$\frac{2r\sqrt{b^2+a^2} + 2ar^2}{2(1-r^2)} = v_a t_*,$

$\frac{ar^2+\sqrt{b^2+a^2}r}{1-r^2} = v_a t_*,$

$\frac{r(ar + \sqrt{b^2+a^2})}{v_a(1-r^2)} = t_*.$

Since $r = \frac{v_a}{v_m}$ , we have

$t_* = \frac{v_a}{v_m}\cdot\frac{ar + \sqrt{b^2+a^2}}{v_a(1-r^2)}=\frac{1}{v_m}\cdot\frac{ar + \sqrt{b^2+a^2}}{1-r^2}.$

Namely,

$t_* = \frac{\sqrt{b^2+a^2} + a\cdot r}{(1-r^2)\cdot v_m}, \quad\quad r=\frac{v_a}{v_m}.$

It follows that the guided missile strikes the fighter jet at $(a+r\cdot\frac{\sqrt{b^2+a^2}+a\cdot r}{1-r^2}, b).$

$e^{k_1} = e^{-\mathrm{arcsinh}(\frac{a}{b})-r\log(b)}$

$= \frac{1}{e^{\mathrm{arcsinh}(\frac{a}{b})+r\log(b)}}$

$= \frac{1}{e^{\mathrm{arcsinh}(\frac{a}{b})}\cdot e^{r\log(b)}}$

$\mathrm{arcsinh}(\blacksquare) = \log\left(\blacksquare + \sqrt{1+\blacksquare^2}\right)$ (see “Deriving Two Inverse Functions“)

$= \frac{1}{e^{\log\left(\frac{a}{b}+\sqrt{1+(\frac{a}{b})^2}\right)}\cdot e^{r\log(b)}}$

$e^{r\log(b)}=e^{\log(b^r)}=b^r$ (see “Introducing Lady L” and “Two Peas in a Pod, Part 3“)

$= \frac{1}{\left(\frac{a}{b}+\sqrt{1+(\frac{a}{b})^2}\right)b^r}.$

i.e.,

$e^{k_1} = \frac{1}{b^r\left(\frac{a}{b}+\sqrt{1+(\frac{a}{b})^2}\right)}.\quad\quad\quad(\star)$

Exercise-1 Show that for $r < 1$ , (8) gives $0 = a +v_a t_* + k_2 \implies -k_2-a = v_a t_*.$ (Hint: (10))

Exercise-2 For $r < 1$ , what is the total distance traveled by the missile when it strikes the fighter jet? (hint: Don’t make things harder than they are)

Exercise-3 Show that if $r \ge 1$ (i.e., $v_a \ge v_m$ ), the missile will not strike the fighter jet. Explain.

Excercise-4 Show that $a>0 \implies \frac{dx}{dt} >0$ .

Excercise-5 For $a <0$ , find the time and position the missile strikes its target .

Prequel to “A Relentless Pursuit”

July 15, 2021August 18, 2021 / Michael Xue / Leave a comment

Fig. 1

Illustrated in Fig. 2 are two circular hoops of unit radius, centered on a common x-axis and a distance $2a$ apart. There is also a soap films extends between the two hoops, taking the form of a surface of revolution about the x-axis. If gravity is negligible, the film takes up a state of stable, equilibrium in which its surface area is a minimum.

Fig. 2

Our problem is to find the function $y(x)$ , satisfying the boundary conditions

$y(-a) = y(a) = 1,\quad\quad\quad(1)$

which makes the surface area

$A=2\pi\displaystyle\int\limits_{-a}^{a}y\sqrt{1+(y')^2}\;dx\quad\quad\quad(2)$

a minimum.

Let

$F(x,y, y') = 2\pi y \sqrt{1+(y')^2}.$

We have

$\frac{\partial F}{\partial y} = 2\pi \sqrt{1+(y')^2}$

and

$\frac{\partial F}{\partial y'} = 2 \pi y \cdot\frac{1}{2}\left(1+(y')^2\right)^{-\frac{1}{2}}\cdot 2y'=\frac{2 \pi y y'}{\sqrt{1+(y')^2}}.$

The Euler-Lagrange equation

$\frac{\partial F}{\partial y} - \frac{d}{dx}\left(\frac{\partial F}{\partial y'}\right) = 0$

becomes

$\sqrt{1+(y')^2} - \frac{d}{dx}\left(\frac{y y'}{\sqrt{1+(y')^2}}\right) = 0.$

Fig. 3

Using Omega CAS Explorer (see Fig. 3), it can be simplified to:

$y \frac{d^2 y}{dx^2}- \left(\frac{dy}{dx}\right)^2=1.$

This is the differential equation solved in “A Relentless Pursuit” whose solution is

$y = C_1\cdot \cosh(\frac{x+C_2}{C_1}).$

We must then find $C_1$ and $C_2$ subject to the boundary condition (1), i.e.,

$C_1\cdot \cosh(\frac{a+C_2}{C_1}) = C_1\cdot\cosh(\frac{-a+C_2}{C_1})\implies \cosh(\frac{a+C_2}{C_1}) = \cosh(\frac{-a+C_2}{C_1}).$

The fact that $\cosh$ is an even function implies either

$a+C_2 = -a+C_2\quad\quad\quad(3)$

$a+C_2 = -(-a+C_2).\quad\quad\quad(4)$

While (3) is clearly false since it claims for all $a >0, a = -a$ , (4) gives

$C_2=0.$

And so, the solution to boundary-value problem

$\begin{cases} y \frac{d^2 y}{dx^2}- \left(\frac{dy}{dx}\right)^2=1,\\ y(-a)=y(a)=1\end{cases}\quad\quad\quad(5)$

$y = C_1\cdot \cosh(\frac{x}{C_1}).\quad\quad\quad(6)$

To determine $C_1$ , we deduce the following equation from the boundary conditions that $y=1$ at $x=\pm a:$

$C_1\cdot \cosh(\frac{a}{C_1}) = 1.\quad\quad\quad(7)$

This is a transcendental equation for $C_1$ that can not be solved explicitly. Nonetheless, we can examine it qualitatively.

Let

$\mu = \frac{a}{C_1}$

and express (7) as

$\cosh(\mu) = \frac{\mu}{a}.\quad\quad\quad(8)$

Fig. 4

A plot of (8)’s two sides in Fig. 4 shows that for sufficient small $a$ , the curves $z = \cosh(\mu)$ and $z = \frac{\mu}{a}$ will intersect. However, as $a$ increases, $z=\frac{\mu}{a}$ , a line whose slope is $\frac{1}{a}$ rotates clockwise towards $\mu$ -axis. The curves will not intersect if $a$ is too large. The critical case is when $a=a^*$ , the curves touch at a single point, so that

$\cosh(\mu) = \frac{\mu}{a^*}\quad\quad\quad(9)$

and $y=\frac{\mu}{a}$ is the tangent line of $z=\cosh(\mu),$ i.e.,

$\sinh(\mu) = \frac{1}{a^*}.\quad\quad\quad(10)$

Dividing (9) by (10) yields

$\coth(\mu) = \mu. \quad\quad\quad(11)$

What the mathematical model (5) predicts then is, as we gradually move the rings apart, the soap film breaks when the distance between the two rings reaches $2a^*$ , and for $a > a^*$ , there is no more soap film surface connects the two rings. This is confirmed by an experiment (see Fig. 1).

We compute the value of $a^*$ , the maximum value of $a$ that supports a minimum area soap film surface as follows.

Fig. 5

Solving (11) for $\mu$ numerically (see Fig. 5), we obtain

$\mu = 1.1997.$

By (10), the corresponding value of

$a^* = \frac{1}{\sinh(\mu)} = \frac{1}{\sinh(1.1997)} = 0.6627$ .

We also compute the surface area of the soap film from (2) and (6) (see Fig. 6). Namely,

$A = 2\pi \displaystyle\int\limits_{-a}^{a} C_1 \cosh\left(\frac{x}{C_1}\right) \sqrt{1+\left(\frac{d}{dx}C_1\cosh\left(\frac{x}{C_1}\right)\right)^2}\;dx = \pi C_1^2\left(\sinh\left(\frac{2a}{C_1}\right) + \frac{2a}{C_1}\right).$

Fig. 6

Exercise-1 Given $a=\frac{1}{2}$ , solve (7) numerically for $C_1.$

Exercise-2 Without using a CAS, find the surface area of the soap film from (2) and (6).

That first sip of coffee in the morning

December 31, 2020January 2, 2021 / Michael Xue / Leave a comment

It is a good idea to enjoy a cup of coffee before starting a busy day.

Suppose the coffee fresh out of the pot with temperature $\alpha^{\circ} C$ is too hot, we can immediately add cream to reduce the temperature by $\delta^{\circ} C$ instantly, then wait for the coffee to cool down naturally to $\omega^{\circ} C$ before sipping it comfortably. We can also wait until the temperature of the coffee drops to $(\omega+\delta)^{\circ} C$ first, then add the cream to further reduce it instantly to $\omega^{\circ} C$ .

Typically, $\alpha = 90, \omega = 75$ , and $\delta = 5$ .

If we are in a hurry and want to wait the shortest possible time, should the cream be added right after the coffee is made, or should we wait for a while before adding the cream?

The heat flow from the hot water to the surrounding air obeys Newton’s cooling and heating law, described by the following ordinary differential equation:

$\frac{d}{dt}\theta(t) = k (E-\theta(t))$

where $\theta(t)$ , a function of time $t$ , is the temperature of the water, $E$ is the temperature of its surroundings, and $k>0$ is a constant depends on the heat transfer mechanism, the contact are with the surroundings, and the thermal properties of the water.

Fig. 1 a place where Newton’s law breaks down

Under normal circumstances, we have

$E \ll \omega < \omega+\delta < \alpha-\delta < \alpha\quad\quad\quad(1)$

Based on Newton’s law, the mathematical model of coffee cooling is:

$\begin{cases} \frac{d}{dt}\theta(t) = k (E-\theta(t)) \\ \theta(0)= \theta_0\end{cases}\quad\quad\quad(2)$

Fig. 2

Solving (2), an initial-value problem (see Fig. 2) gives

$\theta(t) = E + e^{-kt}(\theta_0 - E).$

Therefore,

$t = \frac{\log\left(\frac{E-\theta_0}{E-\theta(t)}\right)}{k}=\frac{\log\left(\frac{\theta_0-E}{\theta(t)-E}\right)}{k}.\quad\quad\quad(3)$

If cream is added immediately (see Fig. 3),

Fig. 3 : cream first

by (3),

$t_1=\frac{\log\left(\frac{(\alpha-\delta)-E}{\omega-E}\right)}{k}.$

Otherwise (see Fig. 4),

Fig. 4: cream last

$t_2=\frac{\log\left(\frac{\alpha-E}{(\omega+\delta)-E}\right)}{k}.$

And so,

$t_1-t_2 = \frac{\log\left(\frac{(\alpha-\delta)-E}{\omega-E}\right)}{k}- \frac{\log\left(\frac{\alpha-E}{(\omega+\delta)-E}\right)}{k}=\frac{1}{k}\left(\log\left(\frac{(\alpha-\delta)-E}{\omega-E}\right)-\log\left(\frac{\alpha-E}{(\omega+\delta)-E}\right)\right)\quad(4)$

Fig. 5

Since

$\frac{(\alpha-\delta)-E}{\omega-E}- \frac{\alpha-E}{(\omega+\delta)-E}=\frac{(\alpha-\delta-E)(\omega+\delta-E)-(\omega-E)(\alpha-E)}{(\omega-E)(\omega+\delta-E)}=\frac{-\delta(\omega+\delta-\alpha)}{(\omega-E)(\omega+\delta-E)}\overset{(1)}{>0}$

implies

$\log\left(\frac{(\alpha-\delta)-E}{\omega-E}\right)-\log\left(\frac{\alpha-E}{(\omega+\delta)-E}\right)>0,\quad\quad\quad(5)$

from (4) , we see that

$t_1-t_2 > 0;$

i.e.,

$t_1 > t_2$

Hence,

If we are in a hurry and want to wait the shortest possible time, we should wait for a while before adding the cream!

Exercise-1 Solve (2) without using a CAS.

Exercise-2 Show that $\frac{\alpha-\delta-E}{\omega-E}\cdot\frac{\omega+\delta-E}{\alpha-E} >1.$

An ODE to Thanksgiving

December 8, 2020March 7, 2023 / Michael Xue / Leave a comment

A turkey is taken from the refrigerator at ${\theta_0}^{\circ} C$ and placed in an oven preheated to $E^{\circ} C$ and kept at that temperature; after $t_1$ minutes the internal temperature of the turkey has risen to ${\theta_1}^{\circ} C$ . The fowl is ready to be taken out when its internal temperature reaches ${\theta_2}^{\circ} C$ .

Typically, ${\theta_0} = 2, E=200, t_1=30, \theta_1=16, \theta_2=88$ .

Determine the cooking time required.

According to Newton’s law of heating and cooling (see “Convective heat transfer“), the rate of heat gain or loss of an object is directly proportional to the difference in the temperatures between the object and its surroundings. This law is best described by the following ODE (Ordinary Differential Equation):

$\frac{d}{dt}{\theta(t)} = k\cdot(E-\theta(t)),\quad\quad\quad(1)$

where $\theta(t), E$ are the temperatures of the object and its surroundings respectively. $k > 0$ is the constant of proportionality.

Fig. 1

We see that (1) has a critical point $\theta^* = E$ . Fig. 1 illustrates the fact that depending on its initial temperature, an object either heats up or cools down, trending towards $E$ in both cases.

We formulate the problem as a system of differential-algebraic equations:

$\begin{cases} \frac{d}{dt}{\theta(t)} = k\cdot(E-\theta(t)) \\ \theta(0)=\theta_0\\ \theta(t_1)=\theta_1 \\ \theta(\boxed{t_2}) = \theta_2\end{cases}(2)$

To find the required cooking time, we solve (2) for $t_2$ (see Fig. 2).

Fig. 2

Using Omega CAS Explorer, the typical cooking time is found to be approximately $4$ hours $(3.88 \cdots\approx 4)$

Luise Lange of Woodrow Wilson Junior College once wrote (see ” A Century of Calculus, Part I”, p. 50):

“In many calculus texts problems are formulated too one-sidedly in terms of particular, numerical data rather than in general terms. While pedagogically it may be wise to begin a new type of problem with some numerical examples, it is only the general formulation, and the interpretation of the answer in general terms, which can give insight into the functional relation between the given and the derived data.”

I agree with her wholeheartedly! On encountering a mathematical modeling problem stated with numerical values, I prefer to re-state it using symbols first. Then solve the problem symbolically. The numerical values are substituted for the symbols at the very end.

This post is a case in point, as the problem is re-formulated from page 1005 of Jan Gullberg’s “Mathematics From the Birth of Numbers”:

Exercise-1 Solving (2) without using a CAS.

Exercise-2 Given $\theta_0< \theta_1 < E$ , show that

$k = \frac{t\log(\frac{E-\theta_0}{E-\theta_1})}{t_1} > 0.$

Exercise-3 Given $\theta_0 < \theta_1< E$ , verify that $\lim\limits_{t\rightarrow \infty} \theta(t) = E$ from

$\theta(t) = -Ee^{ -\frac{t\log(\frac{E-\theta_0}{E-\theta_1})}{t_1} } + \theta_0 e^{-\frac{t\log(\frac{E-\theta_0}{E-\theta_1})}{t_1}} + E.$

Exercise-4 A slice is cut from a loaf of bread fresh from the oven at $180^{\circ} C$ and placed in a room with a constant temperature of $20^{\circ} C$ . After 1 minute, the temperature of the slice is $140^{\circ} C$ . When has the slice of bread cooled to $32^{\circ} C$ ?

What is the shape of a hanging rope?

August 14, 2020August 31, 2020 / Michael Xue / Leave a comment

Question: What is the shape of a flexible rope hanging from nails at each end and sagging under the gravity?

Answer:

First, observe that no matter how the rope hangs, it will have a lowest point $A$ (see Fig. 1)

Fig. 1

It follows that the hanging rope can be placed in a coordinate system whose origin coincides with the lowest point $A$ and the tangent to the rope at $A$ is horizontal:

Fig. 2

At $A$ , the rope to its left exerts a horizontal force. This force (or tension), denoted by $T_0$ , is a constant:

Fig. 3

Shown in Fig. 3 also is an arbitrary point $B$ with coordinates $(x, y)$ on the rope. The tension at $B$ , denoted by $T_1$ , is along the tangent to the rope curve. $\theta$ is the angle $T_1$ makes with the horizontal.

Since the section of the rope from $A$ to $B$ is stationary, the net force acting on it must be zero. Namely, the sum of the horizontal force, and the sum of the vertical force, must each be zero:

$\begin{cases}T_1cos(\theta)=T_0\quad\quad\quad(1)\\ T_1\sin(\theta) = \rho gs\;\;\quad\quad(2)\end{cases}$

where $\rho$ is the hanging rope’s mass density and $s$ its length from $A$ to $B$ .

Dividing (2) by (1), we have

$\frac{T\sin(\theta)}{T\cos(\theta)} = \tan(\theta) = \frac{\rho g}{T_0}s\overset{k=\frac{\rho g}{T_0}}{\implies} \tan(\theta) = ks.\quad\quad\quad(3)$

Since

$\tan(\theta) = \frac{dy}{dx}$ , the slope of the curve at $B$ ,

and

$s = \int\limits_{0}^{x}\sqrt{1+(\frac{dy}{dx})^2}$ ,

we rewrite (3) as

$\frac{dy}{dx} = k \int\limits_{0}^{x}\sqrt{1+(\frac{dy}{dx})^2}\;dx$

and so,

$\frac{d^2y}{dx^2}=k\cdot \frac{d}{dx}(\int\limits_{0}^{x}\sqrt{1+(\frac{dy}{dx})^2}\;dx)=k\sqrt{1+(\frac{dy}{dx})^2}$

i.e.,

$\frac{d^2y}{dx^2}=k\sqrt{1+(\frac{dy}{dx})^2}.\quad\quad\quad(4)$

To solve (4), let

$p = \frac{dy}{dx}$ .

We have

$\frac{dp}{dx} = k\sqrt{1+p^2} \implies \frac{1}{\sqrt{1+p^2}}\frac{dp}{dx}=k.\quad\quad\quad(5)$

Integrate (5) with respect to $x$ gives

$\log(p+\sqrt{1+p^2}) = kx + C_1\overset{p(0)=y'(0)=0}{\implies} C_1 = 0.$

i.e.,

$\log(p+\sqrt{1+p^2}) = kx.\quad\quad\quad(6)$

Solving (6) for $p$ yields

$p = \frac{dy}{dx} =\sinh(kx).\quad\quad\quad(7)$

Integrate (7) with respect to $x$ ,

$y = \frac{1}{k} \cosh(kx) + C_2\overset{y(0)=0,\cosh(0)=1}{\implies}C_2=-\frac{1}{k}$ .

Hence,

$y = \frac{1}{k}\cosh(kx)-\frac{1}{k}$ .

Essentially, it is the hyperbolic cosine function that describes the shape of a hanging rope.

Exercise-1 Show that $\int \frac{1}{\sqrt{1+p^2}} dp = \log(p + \sqrt{1+p^2})$ .

Exercise-2 Solve $\log(p+\sqrt{1+p^2}) = kx$ for $p$ .

From Dancing Planet to Kepler’s Laws

May 12, 2020January 5, 2023 / Michael Xue / 2 Comments

“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$ .

Since

$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$ .

i.e.

$y\frac{dz}{dt}-z\frac{dy}{dt}=A\quad\quad\quad(0-4)$

where $A$ is a constant.

Similarly,

$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.

Consequently,

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

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

$Cz=xz\frac{dy}{dt}-yz\frac{dx}{dt}$ .

Hence

$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

$\frac{d}{dt}(\frac{y}{x})=\frac{x\frac{dy}{dt}-y\frac{dx}{dt}}{x^2}\overset{(0-6)}{=}\frac{C}{x^2}=\frac{0}{x^2}=0$ .

It means

$\frac{y}{x}=k$

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}$ .

i.e.,

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

$x\frac{dy}{dt}-y\frac{dx}{dt}=c_2\quad\quad\quad(1-3)$

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}$

gives

$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.$

i.e.,

$c_1=v_0^2-\frac{2\mu}{r_0},\quad\quad\quad(1-7)$

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

Let

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

then

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

$r=\frac{c_2}{\rho}.\quad\quad\quad(1-13)$

By chain rule,

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

Thus,

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

Since

$c_1+(\frac{\mu}{c_2})^2\overset{(1-7)}{=}v_0^2-\frac{2\mu}{r_0}+(\frac{\mu}{v_0r_0})^2=(v_0-\frac{\mu}{v_0r_0})^2$ ,

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}$ .

i.e.,

$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