Category Archives: Computer Algebra

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.

Advertisements

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.

Does gravity matter ? (Viva Rocketry! Part 1.2)

p152.jpg

A single rocket expels its propellant at a constant rate k.

Assuming constant gravity is the only external force, show that the equation of motion is

(p+m_0-kt)\frac{dv}{dt}=ck-(p+m_0-kt)g

where v is the rocket’s speed, c the speed of the propellant relative to the rocket, p the payload mass, and m_0 the initial rocket mass.

If the rocket burn is continuous, show that the burn time is \frac{\epsilon m_0}{k} and deduce that the final speed given to the payload is

v=-c \log(1-\frac{\epsilon m_0}{m_0+p}) - \frac{g\epsilon m_0}{k}

where 1-\epsilon is the structural factor of the rocket.

Estimate the percentage reduction in the predicted final speed due to the inclusion of the gravity term if

\epsilon=0.8, \frac{p}{m_0}=\frac{1}{100}, c=3.0\;km\;s^{-1}, m_0=10^5\;kg, and k=5 \times 10^3\;kg\;s^{-1}.

Find an expression for the height reached by the rocket during the burn and estimate its value using the data above.


Let’s recall the governing equation of rocket’s flight derived in “Viva Rocketry! Part 1“, namely,

F = m\frac{dv}{dt} + u\frac{dm}{dt}.

In the present context, m = m_0 + p- k t. It implies that

\frac{dm}{dt}=-k

and,

F=-mg=-(m_0+p-kt)g.

With u = c, we have

-(p+m_0-kt)g = (p+m_0-kt)\frac{dv}{dt}+c\cdot(-k),

i.e.,

(p+m_0-kt)\frac{dv}{dt}=ck-(p+m_0-kt)g

or

\frac{dv}{dt}=\frac{ck}{m_0+p-kt}-g.

The structural factor 1-\epsilon indicates the amount of fuel is \epsilon m_0. Since the fuel is burnt at a constant rate k, it must be true that at burnt out time T,

\epsilon m_0 = kT.

Therefore,

T=\frac{\epsilon m_0}{k}.

The solution to initial-value problem

\begin{cases} \frac{dv}{dt}=\frac{ck}{m_0+p-kt}-g\\ v(0) = 0 \end{cases}

tells the speed of the rocket during its flight while fuel is burnt (see Fig. 1):

v = -c \log(1-\frac{kt}{m_0+p})-gt\quad\quad\quad(1)

Screen Shot 2018-12-03 at 4.48.48 PM.png

Fig. 1

Evaluate (1) at burnt out time gives the final speed of the payload:

v_1=-c \log(1-\frac{\epsilon m_0}{m_0+p}) - \frac{g\epsilon m_0}{k}\quad\quad\quad(2)

Notice the first term of (2) is the burnt out velocity without gravity (see “Viva Rocketry! Part 1“)

It follows that the percentage reduction in the predicted final speed due to the inclusion of gravity is

\frac{\frac{g \epsilon m_0} {k}}{-c \log(1-{\epsilon \over {1+\frac{p}{m_0}}})}\quad\quad\quad(3)

Using the given values which are typical, the estimated value of (3) (see Fig. 2) is

0.003\%.

Screen Shot 2018-12-03 at 3.43.03 PM.png

Fig. 2

This shows the results obtained without taking gravity into consideration can be regarded as a reasonable approximation and the characteristics of rocket flight indicated in “Viva Rocketry! Part 1” are valid.

Since v = \frac{dy}{dt}, (1) can be written as

\frac{dy}{dt} = -c \log(1-\frac{kt}{m_0+p})-gt

To find the distance travelled while the fuel is burnt, we solve yet another initial-value problem:

\begin{cases}\frac{dy}{dt} = -c \log(1-\frac{kt}{m_0+p})-gt \\ y(0) = 0 \end{cases}

Screen Shot 2018-12-03 at 4.54.24 PM.png

Fig. 3

The solution (see Fig. 3) is

y= -\frac{1}{2}g t^2 + ct - c\cdot (t-\frac{m_0+p}{k}) \cdot \log(1-\frac{kt}{m_0+p}).

Hence, the height reached at the burnt out time t=\frac{\epsilon m_0}{k} is

h = -\frac{g\epsilon^2 m_0^2}{2k^2}+\frac{c\epsilon m_0}{k}+\frac{c}{k}\cdot (p+(1-\epsilon)m_0) \cdot \log(1-\frac{\epsilon m_0}{m_0+p}).

Using the given values, we estimate that h \approx 27 \; km (see Fig. 4)

Screen Shot 2018-12-03 at 3.34.13 PM.png

Fig. 4


Exercise 1: Find the distance the rocket travelled while the fuel is burnt by solving the following initial-value problem:

\begin{cases}\frac{d^2y}{dt^2} =\frac{ck}{p+m_0-kt}-g \\ y(0) = 0, y'(0)=0 \end{cases}

Thunderbolt (Viva Rocketry! Part 1.1)

2612512.jpg

Fig. 1

Shown in Fig. 1 is an experimental car propelled by a rocket motor. The drag force (air resistance) is given by R = \beta v^2. The initial mass of the car, which includes fuel of mass m_f, is m_0. The rocket motor is burning fuel at the rate of q with an exhaust velocity of u relative to the car.  The car is at rest at t=0. Show that the velocity of the car is given by, for 0 \leq t \leq T,

v(t) = \mu \cdot \frac{1-({m \over m_0})^{\frac{2\beta \mu}{q}}}{1+({m \over m_0})^{\frac{2\beta \mu}{q}}},

where m=m_0-qt, \mu^2=\frac{qu}{\beta}, and T=\frac{m_f}{q} is the time when the fuel is burnt out.


We have derived the governing equation of rocket flight in “Viva Rocketry! Part 1“, namely,

F = m \frac{dv}{dt} + u \frac{dm}{dt}\quad\quad\quad(1)

From m=m_0-qt,  we have

\frac{dm}{dt} = -q.

Apply air resistance R=\beta v^2 as the external force, (1) becomes

-\beta v^2 = (m_0 - q t) \frac{dv}{dt} - u q.

And the car is at rest initially implies

v(0)=0.

It follows that the motion of the car can be modeled by an initial-value problem

\begin{cases} -\beta v^2 = (m_0 - q t) \frac{dv}{dt} - u q \\ v(0) = 0 \end{cases}\quad\quad\quad(2)

It suffices to show that the given v(t) is the solution to this initial-value problem:

Screen Shot 2018-11-28 at 9.28.46 AM.png

Fig. 2

An alternative is obtaining the stated v(t) through solving (2).

Screen Shot 2018-11-27 at 4.40.18 PM.png

Fig. 3

The fact that m = m_0 -qt, (-1)^{2 \frac{\sqrt{b}\sqrt{u}}{\sqrt{q}}} = 1 simplifies the result considerably,

\frac{\sqrt q \sqrt u}{\sqrt \beta}\cdot \frac{m_0^{\frac{2 \sqrt \beta \sqrt u}{\sqrt q}}-m^{\frac{2 \sqrt \beta \sqrt u}{\sqrt q}}}{m_0^{\frac{2 \sqrt \beta \sqrt u}{\sqrt q}}-m^{\frac{2 \sqrt \beta \sqrt u}{\sqrt q}}}\quad\quad\quad(3)

Divide both the numerator and denominator of (3) by m_0^{\frac {2 \sqrt{\beta} \sqrt{u}}{\sqrt{q}}} then yields

\frac{\sqrt q \sqrt u}{\sqrt \beta}\cdot \frac{1-(\frac{m}{m_0})^{\frac{2 \sqrt \beta \sqrt u}{\sqrt q}}}{1-(\frac{m}{m_0})^{\frac{2 \sqrt \beta \sqrt u}{\sqrt q}}}

which is equivalent to the given v(t) since \mu^2=\frac{qu}{\beta}.

At time t=T, the fuel is burnt out. It means

m_0-m_f = m_0 - qT.

Therefore,

T = \frac{m_f}{q}


Exercise 1: Solve (2) manually.

Hint: The differential equation of (2) can be written as \frac{1}{uq - \beta v^2} \frac{dv}{dt} = \frac{1}{m_0 - q t}.

Exercise 2: For m_0=900\;kg, m_f=450\;kg, q=15\;kg/sec, u=500\;m/sec, \beta=0.3, what is the burnout velocity of the car?

Viva Rocketry! Part 1

maxresdefault.jpg

In this post, we will first look at the main characteristics of rocket flight, and then examine the feasibility of launching a satellite as the payload of a rocket into an orbit above the earth.

A rocket accelerates itself by ejecting part of its mass with high velocity.

Screen Shot 2018-11-09 at 9.37.30 AM.png

Fig. 1

Fig. 1 shows a moving rocket. At time t+\Delta t, the mass \Delta m leaves the rocket in opposite direction. As a result, the rocket is being propelled away with an increased speed.

Let

m(\square), m_{\square} – the mass of rocket at time {\square}

\vec{v}_{\square} – the velocity of rocket at time \square

v(\square), v_{\square} – the magnitude of \vec{v}_{\square}

\vec{v}^*_{t+\Delta t} – the velocity of ejected mass \Delta m at t + \Delta t

v^*_{t+\Delta t} – the magnitude of \vec{v}^*_{t+\Delta t}

u – the magnitude of \Delta m‘s velocity relative to the rocket when it is ejected. It is time invariant.


From Fig. 1, we have

\Delta m = m_t - m_{t + \Delta t},

\vec{v}_t = v_{t},

\vec{v}_{t + \Delta t} = v_{t + \Delta t}

and most notably, the relationship between v^*_{t+\Delta t}, v_{t+\Delta t} and u (see “A Thought Experiment on Velocities”):

 v^*_{t+\Delta t} = u - v_{t + \Delta t}.

It follows that

\vec{v}^*_{t+\Delta t} = -v^*_{t+\Delta t} = v_{t + \Delta t} - u,

momentum at time t: \vec{p}(t) = m_t \vec{v}_t = m_t v_t

and,

momentum at time t+\Delta t\vec{p}(t+\Delta t) = m_{t+\Delta t}\vec{v}_{t+\Delta t} + {\Delta m} \vec{v}^*_{t+\Delta t}=m_{t+\Delta t}\vec{v}_{t+\Delta t} + (m_t - m_{t+\Delta t}) \vec{v}^*_{t+\Delta t}= m_{t+\Delta t}v_{t+\Delta t} + (m_t - m_{t+\Delta t})(v_{t+\Delta t}-u).

Consequently, change of momentum in \Delta t is \vec{p}(t+\Delta t)- \vec{p}(t) = m_t (v_{t + \Delta t} - v_t) + u (m_{t + \Delta t} - m_t).

Apply Newton’s second law of motion to the whole system,

\vec{F}= {d \over dt} \vec{p}(t)

= \lim\limits_{\Delta t \rightarrow 0} {{\vec{p}(t+\Delta t) - \vec{p}(t)} \over \Delta t}

= \lim\limits_{\Delta t \rightarrow 0} { {m_t (v_{t + \Delta t} - v_t) + u (m_{t + \Delta t} - m_t)} \over {\Delta t} }

= \lim\limits_{\Delta t \rightarrow 0} {m_t {{v_{t + \Delta t} - v_t} \over {\Delta t}} + {u {{m_{t + \Delta t} - m_t} \over {\Delta t}}}}

= m_t \lim\limits_{\Delta t \rightarrow 0}{(v_{t+\Delta t} - v_t) \over {\Delta t}} + u \lim\limits_{\Delta t \rightarrow 0}{(m_{t +\Delta t} - m_t) \over \Delta t}

That is,

\vec{F} = m(t) {d \over dt} v(t) + u {d \over dt} m(t)

where \vec{F} is the sum of external forces acting on the system.

To get an overall picture of the rocket flight, we will neglect all external forces.

Without any external force, \vec{F} = 0. Therefore

0 = m(t) {d \over dt} v(t) + u {d \over dt} m(t)

i.e.,

{d \over dt} v(t) = -{u \over m(t)} {d \over dt} m(t)\quad\quad\quad(1)

That fact that u, m(t) in (1) are positive quantities shows as the rocket loses mass ({d \over dt} m(t) < 0), its velocity increases ({d \over dt} v(t) > 0)

Integrate (1) with respect to t,

\int {d \over dt} v(t)\;dt = -u \int {1 \over m(t)} {d \over dt} m(t)\;dt

gives

v(t) = -u \log(m(t)) + c

where c is the constant of integration.

At t = 0, v(0)=0, m(0) = m_1 + P where m_1 is the initial rocket mass (liquid or solid fuel + casing and instruments, exclude payload) and P the payload.

It means c = u \log(m_1+P).

As a result,

v(t) = -u \log(m(t)) + u \log(m_1+P)

= -u (\log(m(t) - \log(m_1+P))

= -u \log({m(t) \over m_1+P})

i.e.,

v(t) = -u \log({m(t) \over m_1+P})\quad\quad\quad(2)

Since m_1 is divided into two parts, the initial fuel mass \epsilon m_1 (0 < \epsilon < 1), and the casing and instruments of mass (1-\epsilon)m_1, m(0) can be written as

m(0) = \epsilon m_1 + ( 1 - \epsilon) m_1 + P

When all the fuel has burnt out at t_1,

m(t_1) = (1 - \epsilon)m_1 + P

By (2), the rocket’s final speed at t_1

v(t_1) = -u \log({m(t_1) \over {m_1+P}})

= -u \log({(1-\epsilon)m_1+P  \over {m_1 + P}})

= -u \log({{m_1 + P -\epsilon m_1} \over {m_1+P}})

= -u \log(1-{{\epsilon m_1} \over {m_1+P}})

= -u \log(1-{\epsilon \over {1 + {P \over m_1}}})

= -u \log(1-{\epsilon \over {1 + \beta}})

where \beta = {P \over m_1}.

In other words,

v(t_1) =-u \log(1-{\epsilon \over {1 + \beta}})\quad\quad\quad(3)

Hence, the final speed depends on three parameters

u, \epsilon and \beta

Typically,

u = 3.0\;km\;s^{-1}, \epsilon = 0.8 and \beta = 1/100.

Using these values, (3) gives

Screen Shot 2018-11-05 at 9.21.14 PM.png

v_1 = 4.7\;km\;s^{-1}\quad\quad\quad(4)

This is an upper estimate to the typical final speed a single stage rocket can give to its payload. Neglected external forces such as gravity and air resistance would have reduced this speed.

With (4) in mind, let’s find out whether a satellite can be put into earth’s orbit as the payload of a single stage rocket.

We need to determine the speed that a satellite needs to have in order to stay in a circular orbit of height h above the earth, as illustrated in Fig. 2.

Screen Shot 2018-11-08 at 3.47.09 PM.png

Fig. 2

By Newton’s inverse square law of attraction, The gravitational pull on satellite with mass m_{s} is

{\gamma \; m_{s} M_{\oplus} \over (R_{\oplus} + h)^2}

where universal gravitational constant \gamma = 6.67 \times 10^{-11}, the earth’s mass M_{\oplus} = 5.9722 \times 10^{24}\; kg, and the earth’s radius R_{\oplus} = 6371\;km.

For a satellite to circle around the earth with a velocity of magnitude v, it must be true that

{\gamma \; m_{s} M_{\oplus} \over (R_{\oplus} + h)^2} = {m_{s} v^2 \over (R_{\oplus}+h) }

 i.e,

v = \sqrt{\gamma \; M_{\oplus} \over (R_{\oplus}+h)}

On a typical orbit, h = 100\;km above earth’s surface,

Screen Shot 2018-11-02 at 10.04.41 AM.png

v = 7.8\;km\cdot s^{-1}

This is far in excess of (4), the value obtained from a single stage rocket.

The implication is that a typical single stage rocket cannot serve as the launching vehicle of satellite orbiting around earth.

We will turn to multi-stage rocket in “Viva Rocketry! Part 2”.

 

What moves fast, will slow down, Part One

rockwel1.jpg

This post aims to explain mathematically how populations change.

Our first attempt is based on ideas put forward by Thomas Malthus’ article “An Essay on the Principle of Population” published in 1798.

Let p(t) denotes total population at time t.

Assume in a small interval \Delta t, births and deaths are proportional to p(t) and \Delta t. i.e.

births = a \cdot p(t) \Delta t

deaths = b \cdot p(t) \Delta t

where a, b are constants.

It follows that the change of total population during time interval \Delta t is

p(t+\Delta t) - p(t) = a\cdot p(t)\Delta t - b \cdot p(t)\Delta t = r\cdot p(t)\Delta t

where r = a - b.

Dividing by \Delta t and taking the limit as \Delta t \rightarrow 0, we have

\lim\limits_{\Delta \rightarrow 0} {p(t+\Delta t) - p(t) \over \Delta t} = r \cdot p(t)

which is

{d \over dt} p(t) = r \cdot p(t)\quad\quad\quad(1)

a first order differential equation.

Since (1) can be written as

{1 \over p(t)} {d \over dt} p(t) = r,

integrate with respect to t; i.e.

\int {1 \over p(t)}{d \over dt} p(t)dt = \int {r} dt

leads to

\log p(t) = r\cdot t + c

where c is the constant of integration.

If at t=0, p(0) = p_0, we have

c = \log p_0

and so

p(t) = p_0 e^{r\cdot t}\quad\quad\quad(2)

The result of our first attempt shows that the behavior of the population depends on the sign of constant r. We have exponential growth if r > 0, exponential decay if r < 0 and no change if r = 0.

The world population has been on a upward trend ever since such data is collected (see “World Population by Year“)

Qualitatively, our model (2) with r>0 indicates this trend. However, it also predicts the world population would grow exponentially without limit. And that, is most unlikely to occur, since there are so many limitation factors to growth: lack of food, insufficient energy, overcrowding, disease and war.

Therefore, it is doubtful that model (1) is The One.


Our second attempt makes a modification to (1). It takes the limitation factors into consideration by replacing constant r in (1) with a function r(t). Namely,

r(t) = \gamma - \alpha \cdot p(t)\quad\quad\quad(3)

where \gamma and \alpha are both positive constants.

Replace r in (1) with (3),

 {d \over dt} p(t) = (\gamma - \alpha \cdot p(t)) p(t) = \gamma (1 - {p(t) \over {\gamma \over \alpha}}) p(t)\quad\quad\quad(4)

Since r(t) is a monotonic decreasing function, it shows as population grows, the growth slows down due to the limitation factors.

Let p_{\infty} = {\gamma \over \alpha},

{d \over dt} p(t) = \gamma (1- {p(t) \over p_{\infty}}) p(t)\quad\quad\quad(5)

This is the Logistic Differential Equation.

Written differently as

{d \over dt} p(t) - \gamma \cdot p(t) = -{\gamma \over p_{\infty}} p(t)^2,

the Logistic Differential Equation is also a Bernoulli’s equation (see “Meeting Mr. Bernoulli“)

Let’s understand (5) geometrically without solving it.

Two constant functions, p(t) = 0 or p_{\infty} are solutions of (5), since

{d \over dt} 0 = \gamma (1-{0\over p_{\infty}})  0 = 0

and

{d \over dt} p_{\infty} = \gamma (1-{p_{\infty} \over {p_{\infty}}}) p_{\infty} = 0.

Plot p(t) vs. {d \over dt} p(t) in Fig. 1, the points, 0 and p_{\infty}, are where the curve of {d \over dt} p(t) intersects the axis of p(t).

Screen Shot 2018-10-23 at 9.42.20 AM.png

Fig. 1

At point A where p(t) > p_{\infty}, since {d \over dt} p(t) < 0, p(t) will decrease; i.e., A moves left toward p_{\infty}.

Similarly, at point B where p(t)  < p_{\infty}, {d \over dt} p(t) > 0 implies that p(t) will increase and B moves right toward p_{\infty}.

The model equation can also tell the manner in which p(t) approaches p_{\infty}.

Let p = p(t),

{d^2 \over dt^2} p(t) = {d \over dt}({d \over dt} p)

= {d \over dp} ({d \over dt}p) \cdot {d \over dt} p

= {d \over d p}(\gamma (1-{p \over p_{\infty}})p)\cdot {d \over dt }p

= \gamma(1 - {2 p\over p_{\infty}})\cdot \gamma (1-{p \over p_{\infty}})p

= \gamma^2 p ({{2 p} \over p_{\infty}} -1)({p \over p_{\infty}}-1)

As an equation with unknown p,  \gamma^2 p ({{2 p} \over p_{\infty}} -1)({p \over p_{\infty}}-1)=0 has three zeros:

0, {p_{\infty} \over 2} and p_{\infty}.

Therefore,

{d^2 \over dt^2}p > 0 if p > p_{\infty},

{d^2 \over dt^2} p < 0 if {p_{\infty} \over 2} < p < p_{\infty}

and

{d^2 \over dt^2} p > 0 if p < {p_{\infty} \over 2}.

Consequently p(t), the solution of initial-value problem

\begin{cases} {d \over dt} p(t) = \gamma (1-{p(t) \over p_{\infty}}) p(t) \\  p(0)=p_0 \end{cases}\quad\quad(6)

where p_0 \neq 0, p_{\infty} behaves in the manner illustrated in Fig. 2.

Screen Shot 2018-10-23 at 9.50.25 AM.png

Fig. 2

If  p_0 > p_{\infty}, p(t) approaches p_{\infty} on a concave curve. Otherwise, when {p_{\infty} \over 2} \leq p_0 < p_{\infty}, p(t) moves along a convex curve.  For p_0 < {p_{\infty} \over 2}, the curve is concave first.  It turns convex after p(t) reaches {p_{\infty} \over 2}.


Next, let’s solve the initial-value problem analytically for p_0 \neq 0, p_{\infty}.

Instead of using the result from “Meeting Mr. Bernoulli“, we will start from scratch.

At t where p(t) \neq 0, p_{\infty},  we re-write (5) as

{1 \over p(t)(1-{p(t) \over p_{\infty}}) }{d \over dt} p(t) = \gamma.

Expressed in partial fraction,

({1 \over p(t)} + {{1 \over p_{\infty}} \over {1-{p(t) \over p_{\infty}}}}) {d \over dt} p(t) = \gamma.

Integrate it with respect to t,

\int ({1 \over p(t)} + {{1 \over p_{\infty}} \over {1-{p(t) \over p_{\infty}}}}) {d \over dt} p(t) dt = \int \gamma dt

gives

\log p(t)  - \log (1-{p(t) \over p_{\infty}}) = \gamma t + c

where c is the constant of integration.

i.e.,

\log {p(t) \over {1-{p(t) \over p_{\infty}}}} = \gamma t + c.

Since p(0) = p_0, we have

c = {\log {p_{0} \over {1-{p_0 \over p_{\infty}}}}}

and so

\log ({{p(t) \over {1-{p(t) \over p_{\infty}}}} \cdot {{1-{p_0 \over p_{\infty}}}\over p_0}} )=\gamma t.

Hence,

{{p(t) \over {1 - {p(t) \over p_\infty}}}= {{p_0 \cdot e^{\gamma t}} \over {1-{p_0 \over p_\infty}}}}.

Solving for p(t) gives

p(t) = { p_{\infty} \over {1+({p_{\infty} \over p_0}-1)e^{-\gamma \cdot t}}}\quad\quad\quad(7)

We proceed to show that (7) expresses the value of p(t), the solution to (6) where p_0 \neq 0, p_{\infty}, for all t ‘s (see Fig.3)

Screen Shot 2018-10-25 at 4.29.47 PM.png

Fig. 3

From (7), we have

 \lim\limits_{t \rightarrow \infty} p(t) = p_{\infty}.

It validates Fig. 1.

(7) also indicates that none of the curves in Fig. 2 touch horizontal line p(t) = p_{\infty}.

If this is not the case, then there exists at least one instance of t where p(t) = p_{\infty}; i.e.,

 {p_{\infty} \over {1+({p_{\infty} \over p_0}-1)e^{-\gamma \cdot t}}} = p_{\infty}.

It follows that

{({p_{\infty} \over {p_0}} - 1) e^{-\gamma t}} = 0

Since {e^{-\gamma t}} > 0 (see “Two Peas in a Pod, Part 2“), it must be true that

p_0 = p_{\infty}.

But this contradicts the fact that (7) is the solution of the initial-value problem (6) where p_0 \neq 0,p_\infty.

Reflected in Fig.1 is that A and B will not become p_{\infty}. They only move ever closer to it.

Last, but not the least,

{\lim \limits_{t \rightarrow \infty}} {d \over dt} p(t) =  \gamma (1-{{ \lim\limits_{t \rightarrow \infty} p(t)} \over p_{\infty}}) {\lim\limits_{t \rightarrow \infty} p(t)} = \gamma (1 - {p_{\infty} \over p_{\infty}}) p_{\infty} = 0.

Hence the title of this post.

You say, “y” I say, “y[x]”

220px-Wolfram_Language_Logo_2016.svg.png

 

You see things; and you say “Why?”

But I dream things that never were; and I say “Why not?”

George Bernard Shaw in Back to Methuselah

 

The Wolfram Language function DSolve and NDSolve can solve differential equations.

Let’s look at a few examples.

Example 1 Solving an ODE symbolically. The solution, a function, is evaluated at a given point.

Screen Shot 2018-10-06 at 8.55.21 PM.png

Example 2 Solving an ODE symbolically. Redefine a function and evaluate it at a given point.

Screen Shot 2018-10-06 at 8.52.38 PM.png

Example 3 Solving an ODE initial-value problem symbolically. Get the value at a given point from the symbolic solution.

Screen Shot 2018-10-06 at 8.44.42 PM.png

Example 4 Solving an ODE initial-value problem numerically. Get the value at a given point from the numerical solution.

Screen Shot 2018-10-06 at 8.57.59 PM.png

Regarding whether to specify ‘y‘ or ‘y[x]‘ in DSolve, the only decent explanation I can find is in Stephen Wolfram’s book “The Mathematica Book”. This is straight from horse’s mouth:

When you ask DSolve to get you a solution for y[x], the rule it returns specify how to replace y[x] in any expression. However, these rules do not specify how to replace objects such as y'[x]. If you want to manipulate solutions that you get from DSolve, you will often find it better to ask for solutions for y, rather than y[x].

He then proceeds to give an illustration:

Screen Shot 2018-10-07 at 12.30.46 AM.png

Had you started with DSolve[y'[x]==x+y[x], y[x], x], the result would be

Screen Shot 2018-10-07 at 12.34.12 AM.png

As expected, only y[x] is replaced.