Category Archives: Computer Algebra

A Proof without Calculus

Now, solve \sin(x)=x.

By inspection, x=0.

Is this the only solution?

Visually, it is difficult to tell (see Fig. 1 and Fig. 2)

Fig. 1

Fig. 2

However, we can prove that 0 is the only solution:

Fig. 3

If 0 < x \leq 1 < \frac{\pi}{2} then from Fig. 3, we have

Area of triangle OAB < Area of circular sector OAB.

That is,

\frac{1}{2}\cdot 1 \cdot \sin(x)< \frac{1}{2}\cdot 1 \cdot x.

Hence,

\forall 0< x \leq 1, \sin(x) < x\quad\quad\quad(1)

If x>1 then

\sin(x) \leq 1 \overset{x>1}{\implies} \sin(x) \leq 1< x \implies \forall x>1, \sin(x) <x\quad\quad\quad(2)

Put (1) and (2) together, we have

\forall x > 0, \sin(x) < x\quad\quad\quad(3)

If x < 0 then

-x > 0 \overset{(3)}\implies \sin(-x) < -x \overset{\sin(-x)=-\sin(x)}\implies -\sin(x) < -x \implies \sin(x) > x

i.e.,

\forall x < 0, \sin(x) > x\quad\quad\quad(4)

Therefore by (3) and (4), we conclude that

only when x = 0, \sin(x)=x.

Advertisements

From proof to simpler proof

Solve \sin(x)=2x for x.

By mere inspection, we have x=0.

Visually, it appears that 0 is the only solution (see Fig.1 or Fig. 2)

Fig. 1

Fig. 2

To show that 0 is the only solution of \sin(x)=2x analytically, let

f(x) = \sin(x)-2x

0 is a solution of \sin(x)=2x means

f(0) = 0\quad\quad\quad(1)

Suppose there is a solution

x^* \neq 0\quad\quad\quad(2)

then,

f(x^*)=\sin(x^*)-2x^*=0\quad\quad\quad(3)

Since f(x) is an function continuous and differentiable on (-\infty, +\infty),

by Lagrange’s Mean-Value Theorem (see “A Sprint to FTC“), there \exists c \in (0, x^*) such that

f(x^*)-f(0)=f'(c)(x^*-0)

\overset{(1), (3)}{\implies} 0 = f'(c) x^*\quad\quad\quad(4)

\overset{(2)}{\implies} f'(c)=0\quad\quad\quad(5).

We know

f'(x) = \cos(x)-2.

From (5), we have

\cos(c)-2=0.

i.e.,

\cos(c)=2.

This is not possible since \forall x \in R, -1 \le \cos(x) \le 1.

A simpler alternative without direct applying Lagrange’s Mean-Value Theorem is:

f(x) = \sin(x)-2x \implies f'(x) = \cos(x) - 2 \overset{-1 \le \cos(x) \le 1}{\implies} f'(x) < 0\implies

f(x) =\sin(x)-2x is a strictly decreasing function.

Since f(0) = 0, \forall x<0, f(x) > 0 and \forall x>0, f(x) < 0.

Therefore, 0 is the only solution of f(x)=0. i.e.,

0 is the only solution of \sin(x) = 2x.

Analyzing Heron’s Algorithm

Finding the value of \sqrt{a} for a>0 is equivalent to seeking a positive number x whose square yields a: x^2=a. In other words, the solution to x^2-a=0.

We can find \sqrt{4} by solving x^2-4=0. It is 2 by mere inspection. \sqrt{9} can also be obtained by solving x^2-9=0, again, by mere inspection. But ‘mere inspection’ can only go so far. For Example, what is \sqrt{15241578750190521}?

Method for finding the square root of a positive real number date back to the acient Greek and Babylonian eras. Heron’s algorithm, also known as the Babylonian method, is an algorithm named after the 1^{st} – century Greek mathematician Hero of Alexandria, It proceeds as follows:

  1. Begin with an arbitrary positive starting value x_0.
  2. Let x_{k+1} be the average of x_k and \frac{a}{x_k}
  3. Repeat step 2 until the desired accuracy is achieved.

This post offers an analysis of Heron’s algorithm. Our aim is a better understanding of the algorithm through mathematics.

Let us begin by arbitrarily choosing a number x_0>0. If x_0^2-a=0, then x_0 is \sqrt{a}, and we have guessed the exact value of the square root. Otherwise, we are in one of the following cases:

Case 1: x_0^2-a > 0 \implies x_0^2 > a \implies \sqrt{a} x_0 > a \implies \sqrt{a} > \frac{a}{x_0} \implies \frac{a}{x_0} < \sqrt{a} < x_0

Case 2: x_0^2-a <0 \implies x_0 < \sqrt{a} \implies \sqrt{a} x_0 < a \implies \sqrt{a} < \frac{a}{x_0} \implies x_0 < \sqrt{a} < \frac{a}{x_0}

Both cases indicate that \sqrt{a} lies somewhere between \frac{a}{x_0} and x_0.

Let us define e_0 as the relative error of approximating \sqrt{a} by x_0:

e_0 =\frac{x_0-\sqrt{a}}{\sqrt{a}}

The closer e_0 is to 0, the better x_0 is as an approximation of \sqrt{a}.

Since x_0>0, \sqrt{a}>0,

e_0 + 1  = \frac{x_0-\sqrt{a}}{\sqrt{a}} + 1 = \frac{x_0}{\sqrt{a}} > 0\quad\quad\quad(1)

By (1),

x_0 = \sqrt{a} (e_0+1)\quad\quad\quad(2)

Let x_1 be the mid-point of x_0 and \frac{a}{x_0}:

x_1 = \frac{1}{2} (x_0 + \frac{a}{x_0})\quad\quad\quad(3)

and, e_1 the relative error of x_1 approximating \sqrt{a},

e_1 = \frac{x_1-\sqrt{a}}{\sqrt{a}}\quad\quad\quad(4)

We have

x_1-\sqrt{a}\overset{(3)}{=}\frac{1}{2}(x_0+\frac{a}{x_0}) - \sqrt{a} \overset{(2)}{=}\frac{1}{2}(\sqrt{a}(e_0+1)+\frac{a}{\sqrt{a}(e_0+1)})-\sqrt{a}=\frac{\sqrt{a}}{2}\frac{e_0^2}{e_0+1}\overset{(1)}{>}0\quad(5)

e_1 = \frac{x_1-\sqrt{a}}{\sqrt{a}} \overset{(5)}{=} \frac{\frac{\sqrt{a}}{2}\frac{e_0^2}{e_0+1}}{\sqrt{a}} =  \frac{1}{2}\frac{e_0^2}{e_0+1}\overset{(1)}{ > }0\quad\quad\quad(6)

By (5),

x_1 > \sqrt{a} \implies \sqrt{a} x_1 > a \implies \sqrt{a} > \frac{a}{x_1}\implies \frac{\sqrt{a}}{x_1} < \sqrt{a} < x_1

i.e., \sqrt{a} lies between \frac{a}{x_1} and x_1.

We can generate more values in stages by

x_{k+1} = \frac{1}{2}(x_k + \frac{a}{x_k}), \;k=1, 2, 3, ...\quad\quad(8)

Clearly,

\forall k, x_k >0.

Let

e_k = \frac{x_k-\sqrt{a}}{\sqrt{a}}.

We have

x_k = \sqrt{a}(e_k+1)\quad\quad\quad(9)

and

e_{k+1} = \frac{x_{k+1}-\sqrt{a}}{\sqrt{a}}=\frac{\frac{1}{2}(\sqrt{a}(e_k+1) + \frac{a}{\sqrt{a}(e_k+1)})-\sqrt{a}}{\sqrt{a}}=\frac{1}{2}\frac{e_k^2}{e_k + 1}\quad\quad\quad(10)

Consequently, we can prove that

\forall k \ge 1, e_{k+1} > 0\quad\quad\quad(11)

by induction:

When k = 1, e_{1+1} \overset{(10)}{=} \frac{1}{2}\frac{e_1^2}{e_1+1} \overset{(6)} {>} 0.

Assume when k = p,

e_{p+1} =\frac{1}{2}\frac{e_p^2}{e_p^2+1} >0\quad\quad\quad(12)

When k = p + 1, e_{(p+1)+1} \overset{(10)}{=}\frac{1}{2}\frac{e_{p+1}^2}{e_{p+1} + 1}\overset{(12)}{>}0.

Since (11) implies

\forall k \ge 2, e_k > 0\quad\quad\quad(13)

(10) and (13) together implies

\forall k \ge 1, e_k > 0\quad\quad\quad(14)

It follows that

\forall k \ge 1,  0 \overset{(14)}{<} e_{k+1} \overset{(10)}{=} \frac{1}{2}\frac{e_k^2}{e_k+1} \overset{(14)}{<} \frac{1}{2}\frac{e_k^2}{e_k}=\frac{1}{2}e_k\quad\quad\quad(15)

and

\forall k \ge 1, 0 \overset{(14)}{<} e_{k+1} \overset{(10)}{=} \frac{1}{2}\frac{e_k^2}{e_k+1}\overset{(14)}{<}\frac{1}{2}\frac{e_k^2}{1}<e_k^2\quad\quad\quad(16)

(15) indicates that the error is cut at least in half at each stage after the first.

Therefore, we will reach a stage n that

0 < e_n < 1\quad\quad\quad(17)

It can be shown that

\forall k \ge 1, 0< e_{n+k} < e_n^{2^k}\quad\quad\quad(18)

by induction:

When k = 1,

0 \overset{(14)}{<} e_{n+1} \overset{(10)}{=}\frac{1}{2}\frac{e_n^2}{e_n+1} \overset{}{<} e_n^2=e_n^{2^1}\quad\quad\quad(19)

Assume when k = p,

0 < e_{n+p} < e_n^{2^p}\quad\quad\quad(20)

When k = p+1,

0 \overset{(14)}{<} e_{n+(p+1)} = e_{(n+p)+1} \overset{(16)}{<}e_{n+p}^2\overset{(20)}{<}(e_n^{2^p})^2 = e_n^{2^p\cdot 2} = e_n^{2^{p+1}}

From (17) and (18), we see that eventually, the error decreases at an exponential rate.

Illustrated below is a script that implements the Heron’s Algorithm:

/*
a   - the positive number whose square root we are seeking
x0  - the initial guess
eps - the desired accuracy
*/
square_root(a, x0, eps) := block(
    [xk],
    numer:true,
    xk:x0,
    loop,
        if abs(xk^2-a) < eps then return(xk),
        xk:(xk+a/xk)/2,
    go(loop)
)$  

Suppose we want to find \sqrt{2} and start with x_0 = 99. Fig. 2 shows the Heron’s algorithm in action.

Fig. 1

Déjà vu!

The stages of a two stage rocket have initial masses m_1 and m_2 respectively and carry a payload of mass P. Both stages have equal structure factors and equal relative exhaust speeds. If the rocket mass, m_1+m_2, is fixed, show that the condition for maximal final speed is

m_2^2 + P m_2 = P m_1.

Find the optimal ratio \frac{m_1}{m_2} when \frac{P}{m_1+m_2} = b.


According to multi-stage rocket’s flight equation (see “Viva Rocketry! Part 2“), the final speed of a two stage rocket is

v = -c \log(1-\frac{e \cdot m_1}{m_1 + m_2 +P})- c\log(1-\frac{e \cdot m_2}{m_2+P})

Let m_0 = m_1 + m_2, we have

m_1 = m_0-m_2

and,

v = -c \log(1-\frac{e \cdot (m_0-m_2)}{m_0 +P})- c\log(1-\frac{e \cdot m_2}{m_2+P})

Differentiate v with respect to m_2 gives

v' = \frac{c(e-1)e(2m_2 P-m_0 P +m_2^2)}{(P+m_2)(P-e m_2+m_2)(P+e m_2-e m_0+m_0)}= \frac{c(e-1)e(2m_2-m_0 P+m_2^2)}{(P+m_2)(P+(1-e)m_2)(P+e m_2 + (1-e)m_0)}\quad\quad\quad(1)

It follows that v'=0 implies

2m_2 P-m_0 P + m_2^2 =0.

That is, 2 m_2 P - (m_1+m_2) P + m_2^2 = (m_2-m_1) P + m_2^2=0. i.e.,

m_2^2 + P m_2 = P m_1\quad\quad\quad(2)

It is the condition for an extreme value of v. Specifically, the condition to attain a maximum (see Exercise-2)

When \frac{P}{m_1+m_2} = b, solving

\begin{cases} (m_2-m_1)P+m_2^2=0\\ \frac{P}{m_1+m_2} = b\end{cases}

yields two pairs:

\begin{cases} m_1=\frac{(\sqrt{b}\sqrt{b+1}+b+1)P}{b}\\m_2= -\frac{\sqrt{b^2+b}P+bP}{b}\end{cases}

and

\begin{cases} m_1= - \frac{(\sqrt{b}\sqrt{b+1}-b-1)P}{b}\\ m_2=\frac{\sqrt{b^2+b}P-bP}{b}\end{cases}\quad\quad\quad(3)

Only (3) is valid (see Exercise-1)

Hence

\frac{m_1}{m_2} = -\frac{(\sqrt{b}\sqrt{b+1}-b-1)P}{\sqrt{b^2+b}P-b P} = \frac{\sqrt{b+1}}{\sqrt{b}} = \sqrt{1+\frac{1}{b}}

The entire process is captured in Fig. 2.

Fig. 2


Exercise-1 Given b>0, 0<e<1, P>0, prove:

  1. - \frac{(\sqrt{b}\sqrt{b+1}-b-1)P}{b}>0
  2. \frac{\sqrt{b^2+b}P-bP}{b}>0

Exercise-2 From (1), prove the extreme value attained under (2) is a maximum.

Boosting rocket flight performance without calculus (Viva Rocketry! Part 2.1)

Fig. 1

Given

v = c\log(\frac{(a+1)(b+1)((a+1)b+1)}{(1-e+(a+1)b)((1-e)a+(a+1)b+1})

where b>0, 0<e<1, c>0, a>0, maximize v with appropriate a.

The above optimization problem is solved using calculus (see “Viva Rocketry! Part 2“). However, there is an alternative that requires only high school mathematics with the help of a Computer Algebra System (CAS). This non-calculus approach places more emphasis on problem solving through mathematical thinking, as all symbolic calculations are carried out by the CAS (e.g., see Fig. 2). It also makes a range of interesting problems readily tackled with minimum mathematical prerequisites.


The fact that

\log is a monotonic increasing function \implies v_{max} = c\log(w_{max})

where

w = \frac{(a+1)(b+1)((a+1)b+1)}{(1-e+(a+1)b)((1-e)a+(a+1)b+1)}

or

w(1-e+(a+1)b)((1-e)a+(a+1)b+1) - (a+1)(b+1)((a+1)b+1)=0\quad\quad\quad(1)

(1) can be written as

A_1 a^2 + B_1 a +C_1= 0

where

A_1 = -bew+b^2w+bw-b^2-b

B_1 = e^2w-2bew-2ew+2b^2w+3bw+w-2b^2-3b-1,

C_1 = -bew-ew+b^2w+2bw+w-b^2-2b-1.


Since A_1 = 0 means

-b e w + b^2  w + b w - b^2 - b =0.

That is

w = -\frac{b+1}{e-b-1}.

Solve

\frac{(a+1)(b+1)((a+1)b+1)}{(1-e+(a+1)b)((1-e)a+(a+1)b+1)} = -\frac{b+1}{e-b-1}

for a gives a = 0 \implies A_1 \neq 0 if a > 0.

Hence, (1) is a quadratic equation. For it to have solution, its discriminant B_1^2-4A_1C_1 must be nonnegative, i.e.,

(e^2-2be-2e+b+1)^2 w^2-2(b+1)(2be^2+e^2-2be-2e+b+1) w +(b+1)^2 \geq 0\quad(2)

Consider

(e^2-2be-2e+b+1)^2 w^2-2(b+1)(2be^2+e^2-2be-2e+b+1)w+(b+1)^2 = 0\quad(3)

If e^2-2be-2e+b+1 \neq 0, (3) is a quadratic equation.

Solving (3) yields two solutions

w_1 = -\frac{(b+1)(2\sqrt{b(b+1)}e^2-2be^2-e^2-2\sqrt{b(b+1)}e+2be+2e-b-1)}{(e^2-2be-2e+b+1)^2},

w_2=\frac{(b+1)(2\sqrt{b(b+1)}e^2+2be^2+e^2-2\sqrt{b(b+1)}e-2be-2e+b+1)}{(e^2-2be-2e+b+1)^2}.

Since 0 < e < 1,

w_1 - w_2 = -\frac{4(b+1)\sqrt{b(b+1)}(e-1)e}{(e^2-2be-2e+b+1)^2} > 0\quad(4)

(4) implies

w_1 > w_2

and, the solution to (2) is

w \leq w_2 or w \ge w_1

i.e.,

w \leq \frac{(b+1)(2\sqrt{b(b+1)}e^2+2be^2+e^2-2\sqrt{b(b+1)}e-2be-2e+b+1)}{(e^2-2be-2e+b+1)^2}\quad\quad\quad(4)

or

w \ge -\frac{(b+1)(2\sqrt{b(b+1)}e^2-2be^2-e^2-2\sqrt{b(b+1)}e+2be+2e-b-1)}{(e^2-2be-2e+b+1)^2}\quad\quad\quad(5)

We prove that (4) is true by showing (5) is false:

Consider w - w_1= 0:

\frac{(b+1)(e-1) e \cdot f(a)}{(1-e+ab+b)(a(1-e)+ab+b+1)(e^2-2be-2e+b+1)^2} = 0\quad\quad\quad(6)

where

f(a) = 2a\sqrt{b(b+1)}e^2+a^2be^2+be^2+e^2-2a^2b\sqrt{b(b+1)}

-4ab\sqrt{b(b+1)}e - 2b\sqrt{b(b+1)}e - 4a\sqrt{b(b+1)}e

-2\sqrt{b(b+1)}e-2a^2b^2e-4ab^2e-2a^2be-4abe-4be

-2e+2a^2b^2\sqrt{b(b+1)}+4ab^2\sqrt{b(b+1)} + 2b^2\sqrt{b(b+1)} +2a^2b\sqrt{b(b+1)}

+6ab\sqrt{b(b+1)} + 4b\sqrt{b(b+1)} + 2a\sqrt{b(b+1)} + 2\sqrt{b(b+1)}

+2a^2b^3 + 4ab^3 + 2b^3 +3a^2b^2 + 8ab^2 +5b^2+a^2b+4ab+4b+1.

It can be written as

A_2a^2 + B_2a + C_2\quad\quad\quad(7)

where

A_2 = be^2-2b\sqrt{b(b+1)}e-2b^2e-2be+2b^2\sqrt{b(b+1)}+2b\sqrt{b(b+1)}

+2b^3+3b^2+b,

B_2 = 2\sqrt{b(b+1)}e^2-4b\sqrt{b(b+1)}e -4\sqrt{b(b+1)}e

-4b^2e-4be+4b^2\sqrt{b(b+1)}+6b\sqrt{b(b+1)}+2\sqrt{b(b+1)}+4b^3+8b^2+4b,

C_2=be^2+e^2-2b\sqrt{b(b+1)}e-2\sqrt{b(b+1)}e-2b^2e-4be

-2e+2b^2\sqrt{b(b+1)}+4b\sqrt{b(b+1)}+2\sqrt{b(b+1)}+2b^3+5b^2+4b+1.

Since A_2 > 0 (see Exercise 1) and,

solve (7) for a yields

a = -\sqrt{1+\frac{1}{b}}.

It follows that for a > 0, f(a) > 0.

Consequently, w-w_1 is a negative quantity. i.e.,

w-w1 <  0

which tells that (5) is false.

Hence, when e^2-2be-2e+b+1 \neq 0, the global maximum w_{max} is w_2.

Solving w = w_2 for a:

\frac{(a+1)(b+1)((a+1)b+1}{(1-e+(a+1)b)((1-e)a+(a+1)b+1)} = \frac{(b+1)(2\sqrt{b(b+1)}e^2+2be^2+e^2-2\sqrt{b(b+1)}e-2be-2e+b+1)}{(e^2-2be-2e+b+1)^2},

we have

a = \sqrt{1+ \frac{1}{b}}.

Therefore,

e^2-2be-2e+b+1 \neq 0 \implies w attains maximum at a = \sqrt{1+ \frac{1}{b}}.

In fact, w attains maxima at a = \sqrt{1+\frac{1}{b}} even when e^2-2be-2e+b+1 = 0, as shown below:

Solving e^2-2be-2e+b+1 = 0 for e, we have

e_1= -\sqrt{b(b+1)}+b+1 or e_2 = \sqrt{b(b+1)} + b + 1.

Only e_1 is valid (see Exercise-2),

When e = e_1,

w(\sqrt{1+\frac{1}{b}} )- w(a) = - \frac{(b+1)g(a)}{4 \sqrt{b(b+1)} (\sqrt{b(b+1)}+ab) (a\sqrt{b(b+1)}+b+1)}\quad(8)

where

g(a) = (2a^2b+4ab+2b+2a+2)\sqrt{b(b+1)}-2a^2b^2-4ab^2-2b^2-a^2b-4ab-3b-1

Solve quadratic equation g(a) = 0 for a yields

a = \sqrt{1+\frac{1}{b}}.

The coefficient of a^2 in g(a) is 2b\sqrt{b(b+1)}-2b^2-b, a negative quantity (see Exercise-3).

The implication is that g(a) is a negative quantity when a \neq \sqrt{1 + \frac{1}{b}}.

Hence, (8) is a positive quantity, i.e.,

e^2-2be-2e+b+1 = 0, a \neq \sqrt{1+\frac{1}{b}} \implies w(\sqrt{1+\frac{1}{b}})-w(a) > 0

We therefore conclude

\forall 0 < e < 1, b > 0, w attains its maximum at a = \sqrt{1+\frac{1}{b}}.

Fig. 2


Exercise-1 Prove:0<e<1, b>0 \implies

be^2-2b\sqrt{b(b+1)}e-2b^2e-2be+2b^2\sqrt{b(b+1)}+2b\sqrt{b(b+1)}+2b^3+3b^2+b > 0

Exercise-2 Prove: b > 0 \implies 0 <-\sqrt{b(b+1)} + b +1 <1

Exercise-3 Prove: b > 0 \implies 2b\sqrt{b(b+1)}-2b^2-b < 0

Piece of Pi

A while back, we deemed that it is utterly impractical to calculate the value of \pi using the partial sum of Leibniz’s series due to its slow convergence (see “Pumpkin Pi“)

Fig. 1

As illustrated in Fig. 1, in order to determine each additional correct digit of \pi, the number of terms in the summation must increase by a factor of 10.

What we need is a fast converging series whose partial sum yields given number of correct digits with far fewer terms.

Looking back, we see that the origin of Leibniz’s series is the definite integral

\frac{\pi}{4} = \int\limits_{0}^{1}\frac{1}{1+x^2} dx\quad\quad\quad(1)

To find the needed new series, we consider a variation of (1), namely,

\frac{\pi}{6} = \int\limits_{0}^{\frac{1}{\sqrt{3}}}\frac{1}{1+x^2} dx\quad\quad\quad(2)

Given the fact (see “Pumpkin Pi“) that

\frac{1}{1+x^2} = \sum\limits_{k=1}^{n}(-1)^{k+1}x^{2k-2}+\frac{(-1)^n x^{2n}}{1+x^2}\quad\quad\quad(3)

We proceed to integrate (3) with respect to x from 0 to \frac{1}{\sqrt{3}},

\int\limits_{0}^{\frac{1}{\sqrt{3}}}\frac{1}{1+x^2} dx=\int\limits_{0}^{\frac{1}{\sqrt{3}}}\sum\limits_{k=1}^{n}(-1)^{k+1}x^{2k-2} dx + \int\limits_{0}^{\frac{1}{\sqrt{3}}}\frac{(-1)^n x^{2n}}{1+x^2}dx

= \sum\limits_{k=1}^{n}(-1)^{k+1}\int\limits_{0}^{\frac{1}{\sqrt{3}}}x^{2k-2} dx + (-1)^n\int\limits_{0}^{\frac{1}{\sqrt{3}}}\frac{x^{2n}}{1+x^2}dx

= \sum\limits_{k=1}^{n}(-1)^{k+1}\frac{x^{2k-1}}{2k-1}\bigg|_{0}^{\frac{1}{\sqrt{3}}}+ (-1)^n\int\limits_{0}^{\frac{1}{\sqrt{3}}}\frac{x^{2n}}{1+x^2}dx

As result, (2) becomes

\frac{\pi}{6}=\sum\limits_{k=1}^{n}(-1)^{k+1}\frac{(\frac{1}{\sqrt{3}})^{2k-1}}{2k-1}+(-1)^n\int\limits_{0}^{\frac{1}{\sqrt{3}}}\frac{x^{2n}}{1+x^2}dx

= \sum\limits_{k=1}^{n}(-1)^{k+1}\frac{{(\frac{1}{3}})^k}{(2k-1)\frac{1}{\sqrt{3}}} + (-1)^n\int\limits_{0}^{\frac{1}{\sqrt{3}}}\frac{x^{2n}}{1+x^2}dx

= \sqrt{3}\sum\limits_{k=1}^{n}\frac{(-1)^{k+1}}{3^k(2k-1)} + (-1)^n\int\limits_{0}^{\frac{1}{\sqrt{3}}}\frac{x^{2n}}{1+x^2}dx

i.e.,

\frac{\pi}{6} -  \sqrt{3}\sum\limits_{k=1}^{n}\frac{(-1)^{k+1}}{3^k(2k-1)} = (-1)^n\int\limits_{0}^{\frac{1}{\sqrt{3}}}\frac{x^{2n}}{1+x^2}dx\quad\quad\quad(4)

By (4),

|\frac{\pi}{6} -  \sqrt{3}\sum\limits_{k=1}^{n}\frac{(-1)^{k+1}}{3^k(2k-1)}| = |(-1)^n\int\limits_{0}^{\frac{1}{\sqrt{3}}}\frac{x^{2n}}{1+x^2}dx|=\int\limits_{0}^{\frac{1}{\sqrt{3}}}\frac{x^{2n}}{1+x^2}dx<\int\limits_{0}^{\frac{1}{\sqrt{3}}}x^{2n}dx

=\frac{x^{2n+1}}{2n+1}\bigg|_0^{\frac{1}{\sqrt{3}}}

=\frac{1}{3^n \sqrt{3} (2n+1)}

which gives

|\sqrt{3}\sum\limits_{k=1}^{n}\frac{(-1)^{k+1}}{3^k(2k-1)} - \frac{\pi}{6}| < \frac{1}{3^n \sqrt{3} (2n+1)}.

And so

-\frac{1}{3^n \sqrt{3} (2n+1)}<\sqrt{3}\sum\limits_{k=1}^{n}\frac{(-1)^{k+1}}{3^k(2k-1)}-\frac{\pi}{6}<\frac{1}{3^n \sqrt{3} (2n+1)}\quad\quad\quad(5)

Since \lim\limits_{n\rightarrow \infty}\frac{1}{3^n \sqrt{3} (2n+1)}=0, (5) implies

\lim\limits_{n\rightarrow \infty} \sqrt{3}\sum\limits_{k=1}^{n}\frac{(-1)^{k+1}}{3^k(2k-1)}-\frac{\pi}{6}= 0.

Hence,

\lim\limits_{n\rightarrow \infty} \sqrt{3}\sum\limits_{k=1}^{n}\frac{(-1)^{k+1}}{3^k(2k-1)}=\frac{\pi}{6}.

It follows that the value of \pi can be approximated by the partial sum of a new series

6\sqrt{3}\sum\limits_{k=1}^{\infty}\frac{(-1)^{k+1}}{3^k(2k-1)}

Let’s compute it with Omega CAS Explorer (see Fig. 2, 3)

Fig. 2

Fig. 2 shows the series converges quickly. The sum of the first 10 terms yields the first 6 digits!

Fig. 3

Totaling the first 100 terms of the series gives the first 49 digits of \pi (see Fig. 3)


Exercise 1. Show that \lim\limits_{n\rightarrow \infty}\frac{1}{3^n \sqrt{3} (2n+1)}=0.

Exercise 2. Can we use \frac{\pi}{3} = \int\limits_{0}^{\sqrt{3}}\frac{1}{1+x^2}dx to compute the value of \pi in a similar fashion? Explain.

Viva Rocketry! Part 2

Screen Shot 2018-12-24 at 5.13.14 AM.png

Fig. 1

A rocket with n stages is a composition of n single stage rocket (see Fig. 1) Each stage has its own casing, instruments and fuel. The nth stage houses the payload.

This image has an empty alt attribute; its file name is screen-shot-2019-01-09-at-6.24.56-pm.png

Fig. 2

The model is illustrated in Fig. 2, the i^{th} stage having initial total mass m_i and containing fuel \epsilon_i m_i (0 < \epsilon_i <1, 1 \leq i \leq n). The exhaust speed of the i^{th} stage is c_i.

The flight of multi-stage rocket starts with the 1^{st} stage fires its engine and the rocket is lifted. When all the fuel in the 1^{st} stage has been burnt, the 1^{st} stage’s casing and instruments are detached. The remaining stages of the rocket continue the flight with 2^{nd} stage’s engine ignited.

Generally, the rocket starts its i^{th} stage of flight with final velocity achieved at the end of previous stage of flight. The entire rocket is propelled by the fuel in the i^{th} casing of the rocket. When all the fuel for this stage has been burnt, the i^{th} casing is separated  from the rest of the stages. The flight of the rocket is completed if i=n. Otherwise, it enters the next stage of flight.

When all external forces are omitted, the governing equation of rocket’s i^{th} stage flight (see “Viva Rocketry! Part 1” or “An alternate derivation of ideal rocket’s flight equation (Viva Rocketry! Part 1.3)“) is

0 = m_i(t) \frac{dv_i(t)}{dt} + c_i \frac{dm_i(t)}{dt}

It can be written as

\frac{dv_i(t)}{dt} = -\frac{c_i}{m_i(t)}\frac{dm_i(t)}{dt}\quad\quad\quad(1)

Integrate (1) from t_0 to t,

\int\limits_{t_0}^{t}\frac{dv_i(t)}{dt}\;dt = -c_i \int\limits_{t_0}^{t}\frac{1}{m_i(t)}\frac{dm_(t)}{dt}\;dt

gives

v_i(t) - v_i(t_0) = -c_i(\log(m_i(t))-\log(m_i(t_0)))=-c_i\log(\frac{m_i(t)}{m_i(t_0)}).

At t=t^* when i^{th} stage’s fuel has been burnt, we have

v_i(t^*) - v_i(t_0) = -c_i\log(\frac{m_i(t^*)}{m_i(t_0)})\quad\quad\quad(2)

where

m_i(t_0) = \sum\limits_{k=i}^{n} m_k +P

and,

m_i(t^*) = \sum\limits_{k=i}^n m_k - \epsilon_i m_i + P.

Let v^*_i = v_i(t^*), \; v^*_{i-1} the velocity of rocket at the end of {i-1}^{th} stage of flight.

Since v_i(t_0) = v^*_{i-1}, (2) becomes

v^*_i - v^*_{i-1} = -c_i \log(\frac{ \sum\limits_{k=i}^n m_k - \epsilon_i m_i+P}{\sum\limits_{k = i}^{n} m_k + P}) = -c_i\log(1-\frac{\epsilon_i m_i}{\sum\limits_{k = i}^{n} m_k + P})

i.e.,

v^*_i = v^*_{i-1}-c_i\log(1-\frac{\epsilon_i m_i}{\sum\limits_{k = i}^{n} m_k + P}), \quad\quad 1\leq i \leq n, v_0=0\quad\quad\quad(3)


For a single stage rocket (n=1), (3) is

v_1^*=-c_1\log(1-\frac{\epsilon_1}{1+\beta})\quad\quad\quad(4)

In my previous post “Viva Rocketry! Part 1“, it shows that given c_1=3.0\;km\;s^{-1}, \epsilon_1=0.8 and \beta=\frac{1}{100}, (4) yields 4.7 \;km\;s^{-1}, a value far below 7.8\;km\;s^{-1}, the required speed of an earth orbiting  satellite.

But is there a value of \beta that will enable the single stage rocket to produce the speed a satellite needs? 

Let’s find out.

Differentiate (4) with respect to \beta gives

\frac{dv_1^*}{d\beta} = - \frac{c_1 \epsilon_1}{(\beta+1)^2 (1-\frac{\epsilon_1}{\beta+1})} = - \frac{c_1\epsilon_1}{(\beta+1)^2(\frac{1-\epsilon_1+\beta}{\beta+1})} < 0

since c_1, \beta are positive quantities and 0< \epsilon_1 < 1.

It means v_1^* is a monotonically decreasing function of \beta.

Moreover,

\lim\limits_{\beta  \rightarrow 0}v_1^*=\lim\limits_{\beta \rightarrow 0}-c_1\log(1-\frac{\epsilon_1}{1+\beta}) = - c_1 \log(1-\epsilon_1)\quad\quad\quad(5)

Given c_1=3.0\;km\;s^{-1}, \epsilon_1=0.8, (5) yields approximately

4.8\;km\;s^{-1}

This image has an empty alt attribute; its file name is screen-shot-2019-01-21-at-4.56.14-pm.png

Fig. 3

This upper limit implies that for the given c_1 and \epsilon_1, no value of \beta will produce a speed beyond (see Fig. 4)

screen-shot-2019-01-21-at-5.26.15-pm

Fig. 4


Let’s now turn to a two stage rocket (n=2)

From (3), we have

v_2^* = -c_1\log(1-\frac{\epsilon_1 m_1}{m_1+m_2+P}) - c_2\log(1-\frac{\epsilon_2 m_2}{m_2+P})\quad\quad\quad(6)

If c_1=c_2=c, \epsilon_1 = \epsilon_2 = \epsilon, m_1=m_2 and \frac{P}{m_1+m_2} = \beta, then

P = \beta (m_1+m_2) = 2m_1\beta = 2m_2\beta.

Consequently, 

v_2^*=-c \log(1-\frac{\epsilon}{2(1+\beta)}) - c\log(1-\frac{\epsilon}{1+2\beta})\quad\quad\quad(7)

When c=3.0\;km\;s^{-1}, \epsilon=0.8 and \beta = \frac{1}{100}

v_2^* \approx 6.1\;km\;s^{-1}

This image has an empty alt attribute; its file name is screen-shot-2019-01-21-at-10.30.21-am.png

Fig. 5

This is a considerable improvement over the single stage rocket (v^*=4.7\; km\;s^{-1}). Nevertheless, it is still short of producing the orbiting speed a satellite needs.

In fact, 

\frac{dv_2^*}{d\beta} = -\frac{2c\epsilon}{(2\beta+2)^2(1-\frac{\epsilon}{2\epsilon+2})}-\frac{2c\epsilon}{(2\beta+1)^2(1-\frac{\epsilon}{2\beta+1})}= -\frac{2c\epsilon}{(2\beta+2)^2\frac{2\beta+2-\epsilon}{2\beta+2}}-\frac{2c\epsilon}{(2\beta+1)^2\frac{2\beta+1-\epsilon}{2\beta+1}} < 0

indicates that v_2^* is a monotonically decreasing function of \beta.

In addition,

\lim\limits_{\beta \rightarrow 0} v_2^*=\lim\limits_{\beta \rightarrow 0 }-c\log(1-\frac{\epsilon}{2+2\beta})-c\log(1-\frac{\epsilon}{1+2\beta})=-c\log(1-\frac{\epsilon}{2})-c\log(1-\epsilon).

Therefore, there is an upper limit to the speed a two stage rocket can produce. When c=3.0\;km\;s^{-1}, \epsilon=0.8, the limit is approximately

6.4\;km\;s^{-1}

This image has an empty alt attribute; its file name is screen-shot-2019-01-24-at-2.07.48-pm.png

Fig. 7

In the value used above, we have taken equal stage masses, m_1 = m_2. i.e., the ratio of m_1 : m_2 = 1 : 1.

Is there a better choice for the ratio of m_1:m_2 such that a better v_2^* can be obtained?

To answer this question, let \frac{m_1}{m_2} = \alpha, we have

m_1 = \alpha m_2\quad\quad\quad(8)

Since P = \beta (m_1+m_2), by (8),

P = \beta (\alpha m_2 + m_2)\quad\quad\quad(9)

Substituting (8), (9) into (6),

v_2^* = -c\log(1-\frac{\epsilon \alpha m_2}{\alpha m_2+m_2+\beta(\alpha m_2+m_2))}) - c\log(1-\frac{\epsilon m_2}{m_2 + \beta(\alpha m_2 + m_2)})

= -c\log(1-\frac{\epsilon \alpha}{\alpha+1+\beta(\alpha+1)})-c\log(1-\frac{\epsilon}{1+\beta(\alpha+1)})

= c\log\frac{(\alpha+1)(\beta+1)((\alpha+1)\beta+1)}{(1-\epsilon+(\alpha+1)\beta)((1-\epsilon)\alpha + (\alpha+1)\beta+1)}\quad\quad\quad(10)

a function of \alpha. It can be written as

v_2^*(\alpha) = c\log(w(\alpha))

where

w(\alpha) = \frac{(\alpha+1)(\beta+1)((\alpha+1)\beta+1)}{(1-\epsilon+(\alpha+1)\beta)((1-\epsilon)\alpha + (\alpha+1)\beta+1)}

This is a composite function of \log and w.

Since \log is a monotonic increasing function (see “Introducing Lady L“), 

(v_2^*)_{max} = c\log(w_{max})

Here, (v_2^*)_{max}, w_{max} denote the maximum of v_2^* and w respectively.

To find w_{max}, we differentiate w,

\frac{dw}{d\alpha} = \frac{(\beta+1)(\alpha^2\beta-\beta-1)(\epsilon-1)\epsilon}{(\epsilon-\alpha \beta-\beta-1)^2(\alpha\epsilon-\alpha \beta-\beta-\alpha-1)^2}=\frac{(\beta+1)(\alpha^2\beta-\beta-1)(\epsilon-1)\epsilon}{(1-\epsilon+\alpha \beta+\beta)^2(\alpha(1-\epsilon)+\alpha \beta+\beta)^2}\quad\quad\quad(11)

Solving \frac{dw}{d\alpha} = 0 for \alpha gives 

\alpha = - \sqrt{1+\frac{1}{\beta}} or \sqrt{1+\frac{1}{\beta}}.

This image has an empty alt attribute; its file name is screen-shot-2019-01-29-at-11.44.15-am.png

Fig. 8

By (8), the valid solution is

\alpha = \sqrt{1+\frac{1}{\beta}}.

It shows that w attains an extreme value at \sqrt{1+\frac{1}{\beta}}.

Moreover, we observe from (11) that

\alpha < \sqrt{1+\frac{1}{\beta}} \implies \frac{dw}{d\alpha} > 0

and

 \alpha > \sqrt{1+\frac{1}{\beta}} \implies \frac{dw}{d\alpha} < 0.

i.e., w attains maximum at \alpha=\sqrt{1+\frac{1}{\beta}}.

It follows that

v_2^* attains maximum at \alpha = \sqrt{1+\frac{1}{\beta}}.

Therefore, to maximize the final speed given to the satellite, we must choose the ratio 

\frac{m_1}{m_2} = \sqrt{1+\frac{1}{\beta}}.

With \beta =\frac{1}{100}, the optimum ratio \frac{m_1}{m_2}=10.05, showing that the first stage must be about ten times large than the second.

Using this ratio and keep \epsilon=0.8, c=3.0\;km\;s^{-1} as before, (10) now gives

v_2 = 7.65\;km\;s^{-1}

a value very close to the required one.

This image has an empty alt attribute; its file name is screen-shot-2019-01-23-at-4.14.53-pm.png

Fig. 9

Setting \beta = \frac{1}{128}, we reach the goal:

v_2^* = 7.8\;km\;s^{-1}

This image has an empty alt attribute; its file name is screen-shot-2019-01-30-at-9.58.40-pm.png

Fig. 10

This image has an empty alt attribute; its file name is screen-shot-2019-01-30-at-10.00.04-pm-1.png

Fig. 11

At last, it is shown mathematically that provided the stage mass ratios (\frac{P}{m_1+m_2} and \frac{m_1}{m_2})are suitably chosen, a two stage rocket can indeed launch satellites into earth orbit.


Exercise 1. Show that \alpha < \sqrt{1+\frac{1}{\beta}} \implies \frac{dw}{d\alpha} > 0 and  \alpha > \sqrt{1+\frac{1}{\beta}} \implies \frac{dw}{d\alpha} < 0.

Exercise 2. Using the optimum \frac{m_1}{m_2} = \sqrt{1+\frac{1}{\beta}} and \epsilon=0.8, c=3.0\;km\;s^{-1}, solving (10) numerically for \beta such that v_2^* = 7.8.