Snowflake and Anti-Snowflake Curves

Fig. 1

The snowflake curve made its first appearence in a 1906 paper written by the Swedish mathematician Helge von Koch. It is a closed curve of infinite perimeter that encloses a finite area.

Start with a equallateral triangle of side length $a$ and area $A_0$ , the snowflake curve is constructed iteratively (see Fig. 1). At each iteration, an outward equallateral triangle is created in the middle third of each side. The base of the newly created triangle is removed:

Fig. 2

Let $s_i$ be the number of sides the snowflak curve has at the end of $i^{th}$ iteration. From Fig. 2, we see that

$\begin{cases}s_i = 4s_{i-1}\\ s_0=3\end{cases}$ .

Fig. 3

Solving for $s_i$ (see Fig. 3) gives

$s_i = 3\cdot 4^i.$

Suppose at the end of $i^{th}$ iteration, the length of each side is $a^*_i$ , we have

$\begin{cases}a^*_i = \frac{a^*_{i-1}}{3}\\ a^*_0=a\end{cases}$ ,

therefore,

$a^*_i = (\frac{1}{3})^i a$

and, $p_i$ , the perimeter of the snowflake curve at the end of $i^{th}$ iteration is

$p_i = s_{i-1} (4\cdot a^*_i) = (3\cdot4^{i-1})\cdot (4\cdot(\frac{1}{3})^ia)=3(\frac{4}{3})^i a$ .

It follows that

$\lim\limits_{i\rightarrow \infty}p_i =\infty$

since $|\frac{4}{3}|>1$ .

Let $A^*_{i-1}, A^*_i$ be the area of equallateral triangle created at the end of $(i-1)^{th}$ and $i^{th}$ iteration respectively. Namely,

$A^*_{i-1} = \boxed{\frac{1}{2}a^*_{i-1}\sqrt{(a^*_{i-1})^2-(\frac{a^*_{i-1}}{2})^2}}$

and

$A^*_{i} = \frac{1}{2}\frac{a^*_{i-1}}{3}\sqrt{(\frac{a^*_{i-1}}{3})^2-(\frac{1}{2}\frac{a^*_{i-1}}{3})^2}= \frac{1}{9}\cdot\boxed{\frac{1}{2}a^*_{i-1}\sqrt{(a^*_{i-1})^2-(\frac{a^*_{i-1}}{2})^2}}=\frac{1}{9}A^*_{i-1}$ .

Since

$\begin{cases}A^*_i = \frac{1}{9} A^*_{i-1} \\ A^*_0=A_0\end{cases}$ ,

solving for $A^*_i$ gives

$A^*_i = (\frac{1}{9})^i A_0$ .

Hence, the area added at the end of $i^{th}$ iteration

$\Delta A_i = s_{i-1} A^*_i = (3\cdot4^{i-1})\cdot(\frac{1}{9})^{i}A_0$ .

After $n$ iterations, the total enclosed area

$A_n =A_0 + \sum\limits_{i=1}^{n}\Delta A_i = A_0 + \sum\limits_{i=1}^{n}(3\cdot4^{i-1})\cdot(\frac{1}{9})^i A_0=\frac{8}{5}A_0-\frac{3}{5}(\frac{4}{9})^{n}A_0$ .

As the number of iterations tends to infinity,

$\lim\limits_{n\rightarrow \infty}A_n = \frac{8}{5}A_0$ .

i.e., the area of the snowflake is $\frac{8}{5}$ of the area of the original triangle.

If at each iteration, the new triangles are pointed inward, the anti-snowflake is generated (see Fig. 4 ).

Fig. 4 First four iterations of anti-snowflake curve

Like the Snowflake curve, the perimeter of the anti-snowflake curve grows boundlessly, whereas its total enclosed area approaches a finite limit (see Exercise-1).

Exercise-1 Let $A_0$ be the area of the original triangle. Show that the area encloded by the anti-snowflake curve approaches $\frac{5}{2}A_0$ .

Exercise-2 An “anti-square curve” may be generated in a manner similar to that of the anti-snowflake curve (see Fig. 5). Find:

(1) The perimeter at the end of $i^{th}$ iteration

(2) The enclosed area at the end of $i^{th}$ iteration

Fig. 5 First four iterations of anti-square curve

An Edisonian Moment

Leave a reply

We are told that shortly after Edison invented the light bulb, he handed the glass section of a light bulb to one of his engineers, asking him to find the volume of the inside. This was quite a challenge to the engineer, because a light bulb is such an irregular shape. Figuring the volume of the bulb’s irregular shape was quite different from figuring the volume of a glass, or a cylinder.

Several days later, Edison passed the engineer’s desk, and asked for the volume of the bulb, but the engineer didn’t have it. He had been trying to figure the volume mathematically, and had some problems because the shape was so irregular.

Edison took the bulb from the man; filled the bulb with water; poured the water into a beaker, which measured the volume, and handed it to the amazed engineer.

Vincent A. Miller, Working with Words, Words to Work With, 2001, pp. 57-58

“Analyze This!” shows the equilibrium point $(x_*, y_*)$ of

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

can be found by solving equation

$-nk_1x^n-k_2x+c_0=0$

where $k_1>0, k_2>0, c_0>0$ .

When $n=4$ , Omega CAS Explorer‘s equation solver yields a list of $x$ ‘s (see Fig. 1).

Fig. 1

It appears that identifying $x_*$ from this list of formidable looking expressions is tedious at best and close to impossible at worst:

Fig. 2

However, by Descartes’ rule of signs ,

$\forall k_1>0, k_2>0, c_0>0, -4k_1x^4-k_2x+c_0=0$ has exactly one positive root.

It means that $x_*$ can be seen quickly from evaluating the $x$ ‘s numerically with arbitrary chosen positive values of $k_1, k_2$ and $c_0$ . For instance, $k1=1, k_2=1, c_0=1$ (see Fig. 3).

Fig. 3

Only the second value is positive:

Therefore, $x_*$ is the second on the list in Fig. 1:

It’s Magic Square!

Leave a reply

A magic square is a $n$ by $n$ square grid filled with distinct integers $1,2,..., n^2$ such that each cell contains a different integer and the sum of the integers, called magic number, is equal in each row, column and diagonal. The order of a $n$ by $n$ magic square is $n$ .

The first known example of a magic square is said to have been found on the back of a turtle by Chinese Emperor Yu in 2200 B.C. (see Fig. 1) : circular dots represents the integers 1 through 9 are arranged in a $3$ by $3$ grid. Its order and magic number are $3$ and $15$ respectively.

Fig. 1

The first magic square to appear in the Western world was the one depicted in a copperplate engraving by the German artist-mathematician Albrecht Dürer, who also managed to enter the year of engraving – $1514$ – in the two middle cells of the bottom row (see Fig. 2)

Fig. 2

Let $m$ be the magic number of a $n$ by $n$ magic square. Sum of all its rows (see “Little Bird and a Recursive Generator“)

$n\cdot m = 1+2+3+...+n^2 = \frac{n^2(n^2+1)}{2}$ .

That is,

$m = \frac{n(n^2+1)}{2}.\quad\quad\quad(1)$

Fig. 3

There is no $2$ by $2$ magic square. For if there is, it must be true (see Fig. 3) that

$a +b = a+c \implies b=c$ , contradicts “each cell contains a different integer”.

Shown in Fig. 4 is a square whose rows and columns add up to $260$ , but the diagonals do not, which makes the square only semi-magic.

Fig. 4

What’s interesting is that a chess Knight, starting its L-shaped moves from the cell marked “ $1$ ” can tour all $64$ squares in numerical order. The knight can reach the starting point from its finishing cell (marked” $64$ “) for a new round of tour.

Exercise-1 Is it true that for every $n \in \mathbb{N}\setminus\left\{2\right\}$ there exists at least one magic square?

Exercise-2 Is there a Knight’s Tour magic square?

Exercise-3 Fill in the missing numbers in the following square so that it is a magic square:

An Epilogue of “Analyze This!”

Leave a reply

In “Analyze This!“, we examined system

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

qualitatively.

Now, let us seek its equilibrium $(x_*, y_*)$ quantitatively.

In theory, one may first solve differential equation

$\frac{dx}{dt} = -nk_1x^n-k_2x+c_0$

for $x(t)$ , using a popular symbolic differential equation solver such as ‘ode2’. Then compute $x_*$ as $\lim\limits_{t \rightarrow \infty} x(t)$ , followed by $y_* = \frac{k_1}{k_2}x_*^n$ .

However,in practice, such attempt meets a deadend rather quickly (see Fig. 1).

Fig. 1

Bring in a more sophisticated solver is to no avail (see Fig. 2)

Fig. 2

An alternative is getting $x_*$ directly from polynormial equation

$-nk_1x^n-k_2x+c_0=0\quad\quad\quad(1)$

We can solve (1) for $x$ if $n \le 4$ . For example, when $n=3$ , $-3k_1 x^3-k_2x+c_0=0$ has three roots (see Fig. 3).

Fig. 3

First two roots are complex numbers. By Descartes’ rule of signs, the third root

$\frac{(\sqrt{4k_2^3+81c_0^2k_1}+9c_0\sqrt{k_1})^{\frac{2}{3}}-2^{\frac{2}{3}}k_2}{3 \cdot 2^{\frac{1}{3}}\sqrt{k_1}(\sqrt{4k_2^3+81c_0^2k_1}+9c_0\sqrt{k_1})^{\frac{1}{3}}}$

is the $x_*$ of equilibrium $(x_*, y_*)$ (see Exercise-1).

Exercise-1 Show the third root is real and positive.

Exercise-2 Obtain $x_*$ from $-4k_1x^4-k_2x+c_0=0.$

A Delightful Piece of Mathematics

Leave a reply

Prove: $\sqrt[3]{\sqrt{108} +10} - \sqrt[3]{\sqrt{108} - 10} = 2.$

We will proceed as follows:

For all $a, b$ ,

$a^3-b^3=(a-b)(a^2+ab+b^2)=(a-b)(a^2-2ab+b^2+3ab)=(a-b)((a-b)^2+3ab)$ .

That is,

$a^3-b^3 = (a-b)((a-b)^2+3ab).\quad\quad\quad(1)$

When

$a = \sqrt[3]{\sqrt{108}+10}, \quad b=\sqrt[3]{\sqrt{108}-10},$

we have

$a^3 = (\sqrt[3]{\sqrt{108}+10})^3 = \sqrt{108}+10, \quad b^3=(\sqrt[3]{\sqrt{108}-10})^3=\sqrt{108}-10.$

It means

$a^3-b^3= 20.$

We also have

$ab = \sqrt[3]{\sqrt{108}+10}\cdot\sqrt[3]{\sqrt{108}-10}=\sqrt[3]{(\sqrt{108})^2-10^2}=\sqrt[3]{8}=2.$

As a result, (1) becomes

$20 = (a-b)((a-b)^2+3\cdot 2).$

Rewrite it as

$20 = x (x^2+6)$

where $x$ denotes $a-b$ , we see that

$a-b=\sqrt[3]{\sqrt{108}+10}-\sqrt[3]{\sqrt{108}-10}$ is a positive root of $20=x(x^2+6).\quad\quad\quad(\star)$

Clearly,

$2$ is also a positive root of $20=x(x^2+6)\quad\quad\quad(\star\star)$

since $20=2\cdot(2^2+6).$

Moreover (see Exercise 1),

$20=x(x^2+6)$ has only one positive root $\quad\quad\quad(\star\star\star)$

Therefore, it must be true that

$\sqrt[3]{\sqrt{108} +10} - \sqrt[3]{\sqrt{108} - 10} = 2.$

Exercise-1 Prove: $20=x(x^2+6)$ has only one positive root.

Exercise-2 Prove: $\sqrt[3]{8+3\sqrt{21}} + \sqrt[3]{8-3\sqrt{21}} = 1$

Analyze This!

2 Replies

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

Let

$p=\frac{dx}{dt},$

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

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

$\frac{dp}{dx}=-n^2k_1x^{n-1}-k_2.$

Integrate it with respect to $x$ ,

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

Let

$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

Hence,

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.

i.e.,

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

Similary,

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

An Epilogue of “Qualitative Analysis of a Simple Chemical Reaction”

Leave a reply

In “Qualitative Analysis of a Simple Chemical Reaction“, we deduced from

$\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)\end{cases}$

a 1st order differential equation

$\frac{dx}{dt}=-2k_1x^2-k_2x+c_0.$

This differential equation can be solved to obtain $x$ first, then by (0-1),

$y = \frac{1}{2k_2}(\frac{dx}{dt}+2k_1x^2).$

As expected,

$\lim\limits_{t\rightarrow \infty} x(t) = \frac{\sqrt{k_2^2+8c_0k_1}-k_2}{4k_1}=x_*$

and

$\lim\limits_{t\rightarrow \infty} y(t) = -\frac{k_2\sqrt{k_2^2+8c0k_1}-k_2^2-4c_0k_1}{8k_1k_2}=\frac{k_1}{k_2}x_*^2=y_*.$

Exercise-1 Solve differential equation $\frac{dx}{dt}=-2k_1x^2-k_2x+c_0$ without CAS.

Hint: rewrite $\frac{dx}{dt}=-2k_1x^2-k_2x+c_0$ as $\frac{1}{2k_1x^2+k_2x-c_0}\cdot\frac{dx}{dt}=-1.$

Qualitative Analysis of a Simple Chemical Reaction

1 Reply

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

$\frac{d}{dt}(\frac{1}{2k_2}(\frac{dx}{dt}+2k_1x^2))=k_1x^2-k_2\frac{1}{2k_2}(\frac{dx}{dt}+2k_1x^2)$

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

i.e.,

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

Let

$p=\frac{dx}{dt}$ ,

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

$p\cdot\frac{dp}{dx}+p\cdot(4k_1x+k_2)=0,\quad\quad\quad(0-7)$

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

$\frac{dp}{dx}=-4k_1x-k_2.$

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

i.e.,

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

Therefore,

$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

1 Reply

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,

$\frac{x^2}{a^2}+\frac{y^2}{b^2}=1.$

Clearly,

$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

$\frac{(x'+f)^2}{a^2}+\frac{y'^2}{b^2}=1.\quad\quad\quad(2)$

Substituting $x', y'$ in (2) by

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

yields equation

$a^2r^2\sin^2(\theta)+b^2r^2\cos^2(\theta)+2b^2fr\cos(\theta)+b^2f^2-a^2b^2=0.$

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

$b^2r^2\cos^2(\theta)+a^2r^2(1-\cos^2(\theta)+2b^2fr\cos(\theta)-a^2b^2+b^2(a^2-b^2)=0.\quad\quad\quad(3)$

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

Let

$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

1 Reply

“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,