Category Archives: Computer Algebra

In the spirit of Archimedes’ method of exhaustion

Shown in Fig. 1 is a pile of n cylinders that approximates one half of a sphere.

Fig. 1

From Fig. 1, we see r_i^2 = r^2-h_i^2.

If \Delta V_i denotes the volume of i^{th} cylinder in the pile, then

\Delta V_i = \pi \cdot r_i^2 \cdot \frac{r}{n}

= \pi (r^2-h_i^2) \cdot \frac{r}{n}

= \pi (r^2-(i\cdot\frac{r}{n})^2)\cdot \frac{r}{n}

= \frac{\pi r^3}{n} (1-(\frac{i}{n})^2).

Let

V_{*}=\sum\limits_{i=1}^{n}\Delta V_i,

we have

V_{*}= \sum\limits_{i=1}^{n}\frac{\pi r^3}{n} (1-(\frac{i}{n})^2)

= \frac{\pi r^3}{n}\sum\limits_{i=1}^{n}(1-(\frac{i}{n})^2)

= \frac{\pi r^3}{n}(\sum\limits_{i=1}^{n}1-\sum\limits_{i=1}^{n}(\frac{i}{n})^2)

= \frac{\pi r^3}{n}(n-\frac{1}{n^2}\sum\limits_{i=1}^{n}i^2)

= \frac{\pi r^3}{n}(n-\frac{n(n+1)(2n+1)}{6n^2})\quad\quad\quad(see “Little Bird and a Recursive Generator)

= \pi r^3(1-\frac{(n+1)(2n+1)}{6n^2})

= \pi r^3(1-\frac{2n^2+3n+1}{6n^2})

= \pi r^3(1-\frac{2+\frac{3}{n}+\frac{1}{n^2}}{6}).

As n \rightarrow \infty,

V_* \rightarrow \pi r^3 (1-\frac{1}{3}) = \frac{2}{3}\pi r^3=V_{\frac{1}{2}sphere}.

Therefore, V_{sphere}= 2\cdot V_{\frac{1}{2}sphere} gives

V_{sphere} = \frac{4}{3}\pi r^3\quad\quad\quad(1-1)

Similarly, we approximate a cone

by a stack of n cylinders:

Fig. 2

Since

h_i = i\cdot\frac{h}{n}\quad\quad\quad(2-1)

and

\frac{h_i}{h} = \frac{r_i}{r},

we have

r_i = \frac{r\cdot h_i}{h}\overset{(2-1)}{=}\frac{r}{h}\cdot i\frac{h}{n}=\frac{r\cdot i}{n}.\quad\quad\quad(2-2)

If \Delta V_i denotes the i^{th} cylinder, then

\Delta V_i = \pi\cdot r_i^2\cdot\frac{h}{n}\overset{(2-2)}{=}\frac{\pi r^2 i^2 h}{n^3}.

Let

V_* = \sum\limits_{i=1}^{n}\Delta V_i,

we have

V_*= \frac{\pi r^2 h}{n^3}\sum\limits_{i=1}^{n} i^2

= \frac{\pi r h}{n^3} \cdot \frac{2n^3+3n^2+n}{6}

= \pi r^2 h (\frac{1}{3} + \frac{1}{2n} + \frac{1}{6n^2}).

Since as n \rightarrow \infty, V_*\rightarrow \frac{1}{3}\pi r^2 h, we have

V_{cone} = \frac{1}{3}\pi r^2 h.\quad\quad\quad(2-3)

We now derive the formula for a sphere’s surface area.

Let us approximate a sphere by many cones:

Fig. 3

From Fig. 3, we see that as \Delta A_i \rightarrow 0,

h \rightarrow r, \;\sum\limits_{i}\Delta A_i \rightarrow A_{sphere}.

Consequently,

\sum\limits_{i} \Delta V_i \overset{(2-3)}{=} \sum\limits_{i} \frac{1}{3}\Delta A_i h = \frac{1}{3}h \sum\limits_{i}\Delta A_i \rightarrow \frac{1}{3}rA_{sphere}

and,

\sum\limits_{i} \Delta V_i \rightarrow V_{sphere} \overset{(1-1)}{=} \frac{4}{3}\pi r^3.

It follows that

\frac{1}{3}r A_{sphere} = \frac{4}{3}\pi r^3.

Hence,

A_{sphere} = 4 \pi r^2.

Let’s Evaluate a Definite Integral Without FTC!

For the readers of “A Case of Pre-FTC Definite Integral“, this post illustrates yet another way of evaluating the following definite integral without FTC:

\int\limits_{a}^{b}x^k\;dx\quad\quad(0\le a \le b, k \in N^+)

Let us divide [0, b] into n unequal subintervals (see Fig. 1) such that

x_n = r^{n-1}b,\; ...,\; x_3=r^2b,\; x_2=rb,\; x_1=b,

i.e.,

x_i = r^{i-1}b, i = 1, 2, ..., n\quad\quad\quad(1)

in general and

r = 1-\frac{1}{n}\quad\quad\quad(2)

so that the subintervals become finer as n increases.

Fig. 1

We have, for the sum of the rectangles above the curve y=x^k,

S_r=\sum\limits_{i=1}^{n}(x^{i}-x^{i+1})\cdot(x_i)^k

= \sum\limits_{i=1}^{n}(r^{i-1}b-r^i b)\cdot(r^{i-1}b)^k

=b^{k+1}\sum\limits_{i=1}^{n}(r^{i-1}-r^i)\cdot(r^{i-1})^k

=b^{k+1}((1-r) + (r-r^2)r^k + (r^2-r^3)r^{2k} + (r^3-r^4)r^{3k} + ... + (r^{n-1}-r^n)r^{(n-1)k})

=b^{k+1}((1-r) + (1-r)r^{k+1} + (1-r)r^{2k+2} + (1-r)r^{3k+3} + ... + (1-r)r^{(n-1)k+ (n-1)})

=b^{k+1}(1-r)(1+r^{k+1} + r^{2(k+1)} + r^{3(k+1)} + ... +r^{(n-1)(k+1)})

=b^{k+1}(1-r)(1+(r^{k+1})+(r^{k+1})^2 + (r^{k+1})^3+... + (r^{k+1})^{n-1})

\overset{s=r^{k+1}}{=}b^{k+1}(1-r)(1+s+s^2+s^3+ ... + s^{n-1})

= b^{k+1}(1-r)\frac{1-s^n}{1-s}\quad\quad\quad\quad(see “Beer Theorems and Their Proofs)

=b^{k+1}(1-r)\frac{1-(r^{k+1})^n}{1-r^{k+1}}

=b^{k+1}(1-r)\frac{1-(r^{k+1})^n}{(1-r)(1+r+r^2+r^3+...+ r^k)}

=b^{k+1}\frac{1-(r^{k+1})^n}{1+r+r^2+r^3+ ... +r^k}.\quad\quad\quad(3)

Since 0<r=1-\frac{1}{n} < 1, as n\rightarrow \infty, r \rightarrow 1, (r^{k+1})^n\rightarrow 0. Consequently,

(3) \rightarrow b^{k+1}\frac{1- 0}{1 + \underbrace{1^1 + 1^2+ 1^3 + ... + 1^k}_{k\;1's}}=\frac{b^{k+1}}{1+k}.

And so,

\int\limits_{0}^{b}x^k \;dx= \frac{b^{k+1}}{k+1}.

It follows that

\int\limits_{a}^{b}x^k\;dx = \int\limits_{0}^{b}x^k\;dx - \int\limits_{0}^{a}x^k\;dx =\frac{b^{k+1}}{k+1}-\frac{a^{k+1}}{k+1},

i.e.,

For k \in N^+, \int\limits_{a}^{b}x^k\;dx = \frac{b^{k+1}-a^{k+1}}{k+1}.

Solving Kepler’s “Wine Barrel Problem” without Calculus

Shown below is a cylinder shaped wine barrel.

Fig. 1

From Fig. 1, we see that

d^2 = (\frac{h}{2})^2 + (2r)^2\implies r^2=\frac{d^2}{4}-\frac{h^2}{16}\quad\quad\quad(1)

and so,

V=\pi r^2 h \overset{(1)}{\implies} V=\pi (\frac{d^2h}{4} - \frac{h^3}{16})\quad\quad\quad(2)

Kepler’s “Wine Barrel Problem” can be stated as:

If d is fixed, what value of h gives the largest volume of V?

Kepler conducted extensive numerical studies on this problem. However, it was solved analytically only after the invention of calculus.

In the spring of 2012, while carrying out a research on solving maximization/minimization problems, I discovered the following theorem:

Theorem-1. For positive quantities a_1, a_2, ..., a_n, c_1, c_2, ..., c_n and positive rational quantities p_1, p_2, ..., p_n, if c_1a_1+c_2a_2+...+ c_na_n is a constant, then a_1^{p_1}a_2^{p_2}...a_n^{p_n} attains its maximum if \frac{c_1a_1}{p_1} = \frac{c_2a_2}{p_2} = ... = \frac{c_na_n}{p_n}.

By applying Theorem-1, the “Wine Barrel Problem” can be solved analytically without calculus at all. It is as follows:

Rewrite (2) as

V = \sqrt{16} \pi (\frac{h^2}{16})^{\frac{1}{2}}(\frac{d^2}{4} - \frac{h^2}{16})^1.\quad\quad\quad(3)

Since

\frac{h^2}{16} >0, \frac{d^2}{4} - \frac{h^2}{16} >0 and 1\cdot \frac{h^2}{16} + 1\cdot(\frac{d^2}{4} -\frac{h^2}{16}) = \frac{d^2}{4}, a constant,

by Theorem-1, when

\frac{1\cdot\frac{h^2}{16}}{\frac{1}{2}} = \frac{1\cdot(\frac{d^2}{4}-\frac{h^2}{16})}{1},\quad\quad\quad

or

\frac{\frac{h^2}{16}}{\frac{1}{2}} = \frac{d^2}{4}-\frac{h^2}{16},\quad\quad\quad(4)

V (see (3)) attains its maximum.

Solving (4) for positive h, we have

h = \frac{2}{\sqrt{3}}d.


Discovered from the same research is another theorem for solving minimization problem without calculus:

Theorem-2. For positive quantities a_1, a_2, ..., a_n, c_1, c_2, ..., c_n and positive rational quantities p_1, p_2, ..., p_n, if a_1^{p_1}a_2^{p_2}...a_k^{p_k} is a constant, then c_1a_1+c_2a_2+...+c_na_n attains its minimum if \frac{c_1a_1}{p_1} = \frac{c_2a_2}{p_2} = ... = \frac{c_na_n}{p_n}.

Let’s look at an example:

Problem: Find the minimum value of \frac{1}{\sqrt[3]{x}} + 27x for x>0.

Since for x >0, \frac{1}{\sqrt[3]{x}} >0, 27x >0 and (\frac{1}{\sqrt[3]{x}})^3\cdot 27x = 27, a constant,

by Theorem-2, when

\frac{\frac{1}{\sqrt[3]{x}}}{3} = 27x,\quad\quad\quad(5)

\frac{1}{\sqrt[3]{x}} + 27x attains its minimum.

Solving (5) for x yields

x = \frac{1}{27}.

Therefore, at x = \frac{1}{27}\approx 0.03703, \frac{1}{\sqrt[3]{x}} + 27x attains its minimum value \sqrt[3]{27}+27\cdot\frac{1}{27}= 3+1=4 (see Fig. 2).

Fig. 2

Nonetheless, neither Theorem-1 nor Theorem-2 is a silver bullet for solving max/min problems without calculus. For example,

Problem: Find the minimum value of x^3-27x for x>0.

Theorem-2 is not applicable here (see Exercise-1). To solve this problem, we proceed as follows:

From (x-3)^3 = x^3-9x^2+27x-27, we have

(x-3)^3+9x^2+27 = x^3+27x

and so,

(x-3)^3+9x^2+27-54x =x^3-27x.

That is,

(x-3)^3+9(x^2-6x+3)=x^3-27x.

Or,

x^3-27x= (x-3)^3+9(x^2-6x+9-6)

=(x-3)^3+9((x-3)^2-6)=(x-3)^3+9(x-3)^2-54

=-54+(x-3)^3+9(x-3)^2=-54+(x-3)^2(x-3+9)

=-54+(x-3)^2(x+6)

i.e.

x^3-27x = -54 +(x-3)^2(x+6).

Since x>0, x+6>0, (x-3)^2 \ge 0 \implies (x-3)^2(x+6) \ge 0,

x^3-27x = -54 + (x-3)^3(x+6) \ge -54,

with the “=” sign in “\ge” holds at x=3.

Therefore, x^3-27x attains its minimum -54 at x=3 (see FIg. 3).

Fig. 3


Exercise-1 Explain why Theorem-2 is not applicable to finding the minimum of x^3-27x for x>0.

What is the shape of a hanging rope?

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.

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


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_2and c_0. For instance, k_1=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:

An Epilogue of “Analyze This!”

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.

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.

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

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