Category Archives: Computer Algebra

FTC saves the day!

Problem-1 Given

f(x) = \int\limits_{0}^{2x}f(\frac{t}{2})\;dt +\log(2)\quad\quad\quad(\star)

where f(x) is a continuous function, find f(x).



p=2x \implies \frac{dp}{dx} = 2.\quad\quad\quad(1-1)

By (\star),

\frac{df(x)}{dx} =\frac{d}{dx} \int\limits_{0}^{p} f(\frac{t}{2})\;dt + \frac{d\log(2)}{dx} = \underline{\frac{d}{dp}\left(\int\limits_{0}^{p} f(\frac{t}{2})\;dt\right)} \cdot \frac{dp}{dx}\overset{\textbf{FTC}}{=}\underline{f(\frac{p}{2})}\cdot \frac{dp}{dx}\overset{(1-1)}{=}2f(x),


\frac{df(x)}{dx} = 2f(x).\quad\quad\quad(1-2)

Moreover, we see from (\star) that

f(0) = \int\limits_{0}^{0}f(\frac{t}{2})\;dt + \log(2) = 0 + \log(2) = \log(2).\quad\quad\quad(1-3)

Solving initial-value problem

\begin{cases} \frac{df(x)}{dx} = 2f(x)\\ f(0)=\log(2)\end{cases}


f(x) = \log(2)\cdot e^{2x}.

We use Omega CAS Explorer to verify:

Fig. 1-1

Problem-2 Given

\int\limits_{0}^{1}f(u\cdot x) \;du = \frac{1}{2} f(x) +1\quad\quad\quad(\star\star)

where f(x) is a continuous function, find f(x).


Let p=u\cdot x,

u=\frac{p}{x} \implies \frac{du}{dp} = \frac{1}{x}\quad\quad\quad(2-1)

u=0\implies p=0; u=1\implies p=x.\quad\quad\quad(2-2)

\int\limits_{0}^{1}f(u\cdot x)\;du\overset{(2)}{=} \int\limits_{0}^{x}f(p)\cdot\frac{du}{dp}\cdot dp\overset{(1)}{=}\int\limits_{0}^{x}f(p)\frac{1}{x}\;dp=\frac{1}{x}\int\limits_{0}^{x}f(p)\;dp.\quad\quad\quad(2-3)

By (2-3), we express (\star\star) as

\frac{1}{x}\int\limits_{0}^{x}f(p)\;dp = \frac{1}{2}f(x)+1,


\int\limits_{0}^{x} f(p)\;dp = \frac{x}{2}f(x)+x.

It follows that

\underline{\frac{d}{dx}\left(\int\limits_{0}^{x}f(p)\;dp\right)}=\frac{d}{dx}\left(\frac{x}{2}f(x)+x\right)\overset{\textbf{FTC}}{\implies}\underline{f(x)}=\frac{1}{2}\left(f(x) + x\frac{d f(x)}{dx}\right)+1.\;(2-4)

Solving differential equation (2-4) (see Fig. 2-1) gives

f(x) = c x + 2.

Fig. 2-1

The solution is verified by Omega CAS Explorer:

Fig. 2-2

Exercise-1 Solving \begin{cases} \frac{df(x)}{dx} = 2f(x)\\ f(0)=\log(2)\end{cases} using a CAS.

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

An Epilogue to “Truth vs. Intellect”

This post illustrates an alternative of compute the approximate value of \pi.

We begin with a circle whose radius is r, and let L_{n}, L_{n+1} denotes the side’s length of regular polygon inscribed in the circle with 2^n and 2^{n+1} sides respectively, n=2, 4, ....

Fig. 1

On one hand, we see the area of \Delta ABC as

\frac{1}{2}\cdot AB\cdot BC = \frac{1}{2}\cdot AB\cdot L_{n+1}.

On the other hand, it is also

\frac{1}{2}\cdot AC\cdot BE = \frac{1}{2}\cdot 2r\cdot \frac{L_n}{2}=\frac{1}{2}\cdot r\cdot L_n.


\frac{1}{2}AB\cdot L_{n+1}= \frac{1}{2}r\cdot L_n.


AB^2\cdot L_{n+1}^2 = r^2\cdot L_n^2\quad\quad\quad(1)

where by Pythagorean theorem,

AB^2= (2r)^2 - L_{n+1}^2.\quad\quad\quad(2)

Substituting (2) into (1) gives

(4r^2-L_{n+1}^2)L_{n+1}^2 = L_n^2\implies 4r^2L_{n+1}^2 - L_{n+1}^4 = r^2 L_n^2.

That is,

L_{n+1}^4-4r^2L_{n+1}^2+r^2 L_n^2 = 0.

Let p = L_{n+1}^2, we have

p^2-4r^2 p + r^2 L_n^2=0.\quad\quad\quad(3)

Solving (3) for p yields

p = 2r^2 \pm r \sqrt{4 r^2-L_n^2}.

Since L_n^2 must be greater than L_{n+1}^2 (see Exercise 1), it must be true (see Exercise 2) that

L_{n+1}^2=2r^2 - r \sqrt{4r^2-L_n^2}.\quad\quad\quad(4)

Notice when r=\frac{1}{2}, we obtain (5) in “Truth vs. Intellect“.

With increasing n,

L_n\cdot 2^n \approx \pi\cdot 2r \implies \pi \approx \frac{L_n 2^n}{2r}.\quad\quad\quad

We can now compute the approximate value of \pi from any circle with radius r:

Fig. 2 r=2

Fig. 3 r=\frac{1}{8}

Exercise 1 Explain L_{n}^2 > L_{n+1}^2 geometrically.

Exercise 2 Show it is 2r^2-r\sqrt{4r^2-L_n^2} that represents L_{n+1}^2.

Truth vs. Intellect

It was known long ago that \pi, the ratio of the circumference to the diameter of a circle, is a constant. Nearly all people of the ancient world used number 3 for \pi. As an approximation obtained though physical measurements with limited accuracy, it is sufficient for everyday needs.

An ancient Chinese text (周髀算经,100 BC) stated that for a circle with unit diameter, the ratio is 3.

In the Bible, we find the following description of a large vessel in the courtyard of King Solomon’s temple:

He made the Sea of cast metal, circular in shape, measuring ten cubits from rim to rim and five cubits high, It took a line of thirty cubits to measure around it. (1 Kings 7:23, New International Version)

This infers a value of \pi = \frac{30}{10} = 3.

It is fairly obvious that a regular polygon with many sides is approximately a circle. Its perimeter is approximately the circumference of the circle. The more sides the polygon has, the more accurate the approximation.

To find an accurate approximation for \pi, we inscribe regular polygons in a circle of diameter 1. Let L_{n}, L_{n+1} denotes the side’s length of regular polygon with 2^n and 2^{n+1} sides respectively, n=2, 4, ....

Fig. 1

From Fig. 1, we have

\begin{cases} L_{n+1}^2 = x^2 + (\frac{1}{2} L_n)^2\quad\quad\quad(1) \\ (\frac{1}{2})^2 = (\frac{1}{2}L_n)^2 + y^2\quad\quad\;\quad(2)\\ x+y = \frac{1}{2}\;\quad\quad\quad\quad\quad\quad(3) \end{cases}

It follows that

y\overset{(2)}{=}\sqrt{(\frac{1}{2})^2-(\frac{1}{2} L_n)^2} \overset{(3)}{ \implies} x=\frac{1}{2}-\sqrt{(\frac{1}{2})^2-(\frac{1}{2}L_n)^2}.\quad\quad\quad(4)

Substituting (4) into (1) yields

L_{n+1}^2 = \left(\frac{1}{2}-\sqrt{(\frac{1}{2})^2-(\frac{1}{2}L_n)^2}\right)^2+(\frac{1}{2}L_n)^2

That is,

L_{n+1}^2 = \frac{1}{4}\left(L_n^2 + \left(1-\sqrt{1-L_n^2}\right)^2\right).

Further simplification gives

L_{n+1}^2 = \frac{1}{2}\left(1-\sqrt{1-L_n^2}\right),\quad\quad\quad(5)

Starting with an inscribed square (L_2^2 =\frac{1}{2}), we compute L_{n+1}^2 from L_{n}^2 (see Fig. 2). The perimeter of the polygon with 2^{n+1} sides is 2^{n+1} \cdot L_{n+1}.

Fig. 2


\lim\limits_{n \rightarrow \infty} 2^n \cdot L_{n} = \pi.

Exercise-1 Explain, and then make the appropriate changes:

Hint: (5) is equivalent to L_{n+1}^2 = \frac{L_n^2}{2\left(1+\sqrt{1-L_{n}^2}\right)}.

We all bleed the same color

In “Mathematical Models in Biology”, Leah Edelstein-Keshet presents a model describing the number of circulating red blood cells (RBC’s). It assumes that the spleen filters out and destroys a fraction of the cells daily while the bone marrow produces a amount proportional to the number lost on the previous day:

\begin{cases} R_{n+1} = (1-f)R_n+M_n\\ M_{n+1} = \gamma f R_n\end{cases}(1)


R_n - number of RBC’s in circulation on day n,

M_n - number of RBC’s produced by marrow on day n,

f - fraction of RBC’s removed by the spleen,

\gamma - numer of RBC’s produced per number lost.

What would be the cell count on the n^{th} day?

Observe first that (1) is equivalent to

R_{n+2} = (1-f)R_{n+1}+M_{n+1}\quad\quad\quad(2)


M_{n+1} = \gamma f  R_n.\quad\quad\quad(3)

Let n = -1,

M_0=\gamma f R_{-1} \implies R_{-1} = \frac{M_0}{\gamma f}.\quad\quad\quad(4)

Substituting (3) into (2) yields

R_{n+2} = (1-f)R_{n+1}+\gamma f R_{n}.

We proceed to solve the following initial-value problem using ‘solve_rec‘ (see “Solving Difference Equations using Omega CAS Explorer“):

\begin{cases} R_{n+2}=(1-f)R_{n+1}+\gamma f R_{n}\\ R_{0}=1, R_{-1} = \frac{1}{\gamma f}\end{cases}

Evaluate the solution with f=\frac{1}{2}, g=1, we have

R_n = \frac{4}{3} + \frac{(-1)^{n+1}2^{-n}}{3}.\quad\quad\quad(5)

Plotting (5) by ‘plot2d(4/3 + (-1)^(n+1)*2^(-n)/3, [n, 0, 10], WEB_IMAGE)’ fails (see Fig. 1) since plot2d treats (5) as a continuous function whose domain includes number such as \frac{1}{2}.

Fig. 1

Instead, a discrete plot is needed:

Fig. 2

From Fig. 2 we see that R_{n} converges to a value between 1.3 and 1.35. In fact,

\lim\limits_{n \rightarrow \infty}  \frac{4}{3} + \frac{(-1)^{n+1}2^{-n}}{3} = \frac{4}{3}\approx 1.3333....

Solving Difference Equations using Omega CAS Explorer

Maxima’s ‘solve_rec‘ solves difference equation:

Fig. 1 Geometric Sequence

It solves initial-value problem as well:

Fig. 2 Fibonacci Sequence

Mathematica has ‘RSolve‘:

and ‘RSolveValue‘:

See also “Pandora’s Box“.

Fire & Water


When a forrest fire occurs, how many fire fighter should be sent to fight the fire so that the total cost is kept minimum?


Suppose the fire starts at t=0; at t=t_1, the fire fighter arrives; the fire is extinguished later at t=t_2.

Let c_1, c_2, c_3 and x denotes the monetary damage per square footage burnt, the hourly wage per fire fighter, the one time cost per fire fighter and the number of fire fighters sent respectively. The total cost consists damage caused by the fire, the wage paid to the fire fighters and the one time cost of the fire fighters:

c(x) = c_1\cdot(total square footage burnt) +\;c_2\cdot (t_2-t_1)\cdot x + c_3\cdot x.

Notice t_2-t_1, the duration of fire fighting.

Assume the fire ignites at a single point and quickly spreads in all directions with flame velocity v_*, the growing square footage of engulfed circular area is a function of time. Namely,

b(t) = \pi r^2  = \pi (v_* t) ^2 = \pi v_*^2 t^2 \overset{k=\pi v_*^2}{=} k t^2.

Its burning rate

v_b(t) = b'(t) = 2 k t \overset{\alpha = 2 k}{=} \alpha t.\quad\quad\quad(1)

However, after the arrival of x fire fighters, \alpha is reduced by \beta x, an amount that is directly proportional to the number of fire fighters on the scene. The reduction of \alpha reflects the fire fighting efforts exerted by the crew. As a result, for t > t_1, v_b(t) declines along the line described by

v_b(t) = (\alpha-\beta x)t + d

where \alpha-\beta x <0. Or equivalently,

\beta x - \alpha >0.\quad\quad\quad(2)

Moreover, the fire is extinguished at t_2 suggests that

v_b(t_2) = 0 \implies (\alpha-\beta x) t_2 + d =0 \implies d = -(\alpha-\beta x)t_2.


\forall t_1< t \le t_2, v_b(t) = (\alpha-\beta x) t - (\alpha-\beta x) t_2=(\alpha-\beta x) (t-t_2).\quad\quad\quad(3)

Combine (1) and (3),

v_b(t) = \begin{cases} \alpha t, \quad\quad\quad\quad\quad\quad\quad\quad0\le t \le t_1\\ (\alpha-\beta x)(t-t_2),\quad\quad t_1< t \le t_2 \end {cases}.

It is further assumed that

v_b(t) is continuous at t=t_1.\quad\quad\quad(4)

We illustrate v_b(t) in Fig. 1.

Fig. 1

The fact that v_b(t) is continuous at t_1 means

\lim\limits_{t \to t_1^+}v_b(t)=(\alpha-\beta x)(t_1-t_2) = \lim\limits_{t \to t_1^-}v_b(t) = h \implies (\alpha-\beta x) (t_1-t_2)=h.

That is,

t_2-t_1=\frac{h}{\beta x - \alpha}.\quad\quad\quad(5)

The area of triangle in Fig. 1 represents b(t_2), the total square footage damaged by the fire. i.e.,

b(t_2) =\frac{1}{2}t_1 h + \frac{1}{2}(t_2-t_1)h\overset{(5)}{=}\frac{1}{2}t_1 h+ \frac{1}{2}\frac{h^2}{\beta x-\alpha}.\quad\quad\quad(6)

Consequently, the total cost

c(x)= c_1 b(t_2) + c_2 x (t_2-t_1)  + c_3 x \overset{(6), (5)}{=}\frac{c_1 t_1 h}{2} + \frac{c_1 h^2}{2(\beta x-\alpha)} + \frac{c_2 h x}{\beta x-\alpha} + c_3 x.

To minimize the total cost, we seek a value x where function c attains its minimum value c(x).

From expressing x as

x = \frac{1}{\beta}\beta x = \frac{1}{\beta}({\beta x-\alpha +\alpha} )= \frac{1}{\beta}(\beta x -\alpha) + \frac{\alpha}{\beta},

we obtain

\frac{c_2 h x}{\beta x-\alpha} = \frac{c_2 h}{\beta x-\alpha}(\frac{1}{\beta}(\beta x - \alpha) + \frac{\alpha}{\beta})=\frac{c_2 h}{\beta} + \frac{c_2 \alpha h}{\beta(\beta x-\alpha)},\quad\quad\quad(7)

c_3 x = \frac{c_3}{\beta}(\beta x-\alpha) + \frac{c_3 \alpha}{\beta}.\quad\quad\quad(8)

It follows that

c(x) = \underbrace{\frac{c_1 t_1 h}{2} + \frac{c_1 h^2}{2(\beta x-\alpha)}}_{c_1\cdot(6)} + \underbrace{\frac{c_2 h}{\beta}+\frac{c_2 \alpha h}{\beta(\beta x - \alpha)}}_{(7)} + \underbrace{\frac{c_3}{\beta}(\beta x-\alpha) + \frac{c_3 \alpha}{\beta}}_{(8)}

=\frac{c_1 t_1 h}{2} + \frac{c_2 h}{\beta} + \frac{c_3 \alpha}{\beta} + \underline{\frac{c_1 \beta  h^2+2 c_2 \alpha h}{2\beta(\beta x - \alpha)} + \frac{c_3(\beta x-\alpha)}{\beta}}.

Hence, to minimize c(x) is to find the x that minimizes \frac{c_1\beta h^2+2 c_2\alpha h }{2\beta(\beta x - \alpha)} + \frac{c_3(\beta x - \alpha)}{\beta}.

Since \frac{c_1\beta h^2+2 c_2\alpha h }{2\beta(\beta x - \alpha)} \cdot \frac{c_3(\beta x - \alpha)}{\beta}=\frac{(c_1 \beta h^2+2 c_2 \alpha h)c_3}{2\beta^2} , a constant, by the following theorem:

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

(see “Solving Kepler’s Wine Barrel Problem without Calculus“), we solve equation

\frac{c_1\beta h^2+2 c_2\alpha h }{2\beta(\beta x - \alpha)} = \frac{c_3(\beta x - \alpha)}{\beta}

for x:

Fig. 2

From Fig. 2, we see that x =\frac{\alpha}{\beta} - \frac{\sqrt{\beta c_1 c_3 h^2+2\alpha c_2 c_3 h}}{\sqrt{2}\beta c_3} \implies \beta x - \alpha < 0, contradicts (2).

However, when x =\frac{\sqrt{\beta c_1 c_3 h^2+2\alpha c_2 c_3 h}}{\sqrt{2}\beta c_3} + \frac{\alpha}{\beta}=\sqrt{\frac{\beta c_1 h^2 + 2 \alpha c_2 h}{2\beta^2 c_3}} + \frac{\alpha}{\beta}, \beta x-\alpha is a positive quantity. Therefore, it is

x \overset{h=\alpha t_1}{=} \alpha \sqrt{\frac{\beta c_1 t_1^2 + 2 c_2 t_1}{2\beta^2 c_3}} + \frac{\alpha}{\beta}

that minimizes the total cost.

That first sip of coffee in the morning

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

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

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

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

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

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

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

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

Under normal circumstances, we have

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

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

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

Fig. 2

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

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


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

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

Fig. 3 : cream first

by (3),


Otherwise (see Fig. 4),

Fig. 4: cream last


And so,

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

Fig. 5


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



from (4) , we see that

t_1-t_2 > 0;


t_1 > t_2


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

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

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

An ODE to Thanksgiving

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

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

Determine the cooking time required.

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

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

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

Fig. 1

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

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

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

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

Fig. 2

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

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

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

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

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

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

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

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

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

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

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