Polar plot

The polar coordinates r and \theta can be converted to the Cartesian coordinates x and y using the trigonometry functions:

\begin{cases} x=r\cdot\cos(\theta) \\ y=r\cdot\sin(\theta)\end{cases}

It follows that a figure specified in (r, \theta) can be plotted by ‘plot2d’ as a parametric curve:

Fig. 1 r = \cos(5\theta)

It is possible to plot two or more parametric curves together:

Fig. 2 r=\theta and r=\cos(5\theta)

An alternate is the ‘draw2d’ function, it draws graphic objects created by the ‘polar’ function:

Fig. 3 r=\theta and r=\cos(5\theta)

Fig. 4 shows a graceful geometric curve that resembles a butterfly. Its equation is expressed in polar coordinates by

r = e^{\sin(\theta)} - 2\cos(4\theta)+\sin(\frac{\theta}{12})^5

Fig. 4

Let’s combine them!

It is possible to combine two or more plots into one picture.

For example, we solve the following initial-value problem

\begin{cases} \frac{dy}{dx} = (x-y)/2 \\ y(0)=1 \end{cases}\quad\quad\quad(\star)

and plot the analytic solution in Fig. 1.

Fig. 1

We can also solve (\star) numerically and plot the discrete data points:

Fig. 2

Fig. 3 is the result of combining Fig.1 and Fig. 2.


It validates the numerical solution obtained by ‘rk’: the two figures clearly overlapped.

An Alternate Solver of ODEs

Besides ‘ode2’, ‘contrib_ode’ also solves differential equations.

For example,

\frac{d^2y}{dx^2} - \frac{1+x}{x} \cdot \frac{dy}{dx}+ \frac{1}{x}\cdot y=0.

While ‘ode2’ fails:

‘contrib_ode’ succeeds:

This is an example taking from page 4 of Bender and Orszag’s “Advanced Mathematical Methods for Scientists and Engineers“. On the same page, there is another good example:

\frac{dy}{dx} = \frac{A^2}{x^4}-y^2, \quad\quad A is a constant.

Using ‘contrib_ode’, we have

It seems that ‘contrib_ode’ is a better differential equation solver than ‘ode2’:

Even though it is not perfect:

From the examples, we see the usage of ‘contrib_ode’ is the same as ‘ode2’. However, unlike ‘ode2’, ‘contrib_ode’ always return a list of solution(s). It means to solve either initial-value or boundary-value problem, the solution of the differential equation is often lifted out of this list first:

Exercise Solve the following differential equations without using a CAS:

  1. \frac{dy}{dx}= \frac{A^2}{x^4} - y^2 (hint: Riccati Equation)
  2. \frac{d^2 y}{dx^2} = \frac{y \frac{dy}{dx}} {x}


As far as I know, ‘bc2’, Maxima’s built-in function for solving two-point boundary value problem only handles the type:

\begin{cases} y'' = f(x, y, y')\\ y(a)=\alpha, y'(b) = \beta \end{cases}\quad\quad\quad(*)

For example, solving \begin{cases} y'' + y (y')^3=0 \\y(0) = 1, y(1)=3 \end{cases}:

But ‘bc2’ can not be applied to

\begin{cases} -y'' + 9 y =12 \sin(t)\\ y(0)=y(2\pi), y'(0) = y'(2\pi) \end{cases}

since it is not the type of (*). However, you can roll your own:

An error occurs on the line where the second boundary condition is specified. It attempts to differentiate the solution with respect to t under the context that t=0. i.e.,

diff(rhs(sol), 0);

which is absurd.

The correct way is to express the boundary conditions using ‘at’ instead of ‘ev’:

The following works too as the derivative is obtained before using ‘ev’:

Nonetheless, I still think using ’at’ is a better way:

“Instrumental Flying”

Consider the following set

S_1 = \{(x,y) | x=\frac{e^y-e^{-y}}{2}\}.\quad\quad\quad(1-1)

For all (x, y_1), (x, y_2) \in S_1, we have

\frac{e^{y_1}-e^{-y_1}}{2}=x, \frac{e^{y_2}-e^{-y_2}}{2}=x.

That is,

(x, y_1), (x, y_2) \in S_1 \implies \frac{e^{y_1}-e^{-y_1}}{2}- \frac{e^{y_2}-e^{-y_2}}{2}= 0.\quad\quad\quad(1-2)

Since \frac{e^{y_1}-e^{-y_1}}{2}- \frac{e^{y_2}-e^{-y_2}}{2}= 0 if and only if y_1 = y_2 (Exercise-1),

\frac{e^{y_1}-e^{-y_1}}{2}-\frac{e^{y_2}-e^{-y_2}}{2}=0 \implies y_1 = y_2.\quad\quad\quad(1-3)

From (1-2) and (1-3), we conclude:

\forall (x, y_1), (x, y_2) \in S_1 \implies y_1 = y_2.

i.e., S_1 defines a function.

Alternatively, (1-1) can be expressed as

S_1 = \{(x, y) | x = \sinh(y)\}.

It shows that S_1 is the inverse of \sinh(x). Therefore, we name the function defined by (1-1) \sinh^{-1} and write it as

y = \sinh^{-1}(x).

Let’s look at another set:

S_2 = \{(x,y)|x=\frac{e^y+e^{-y}}{2}\}.\quad\quad\quad(2-1)

For a pair (x, y_1>0) \in S_2, we have


For another pair (x, y_2=-y_1),

\frac{e^{y_2} + e^{-y_2}}{2} = \frac{e^{-y_1}+e^{-(-y_1)}}{2} = \frac{e^{-y_1} + e^{y_1}}{2} \overset{(2-2)}{=} x\implies (x, y_2=-y_1) \in S_2.

Since y_1>0, y_2 = -y_1 \neq y_1, (x, y_1), (x, y_2) \in S_1 does not implie y_1 = y_2. It means S_2 does not define a function.

However, modification of S_2 gives

S_3 = \{(x, y) | x=\frac{e^y+e^{-y}}{2}, y \ge 0\}.\quad\quad\quad(2-3)

It defines a function.

Rewrite (2-3) as

S_3 = \{(x, y) | x = \cosh(y), y \ge 0\}.\quad\quad\quad(2-4)

Then \forall (x, y_1), (x, y_2) \in S_3, we have

x=\cosh(y_1), x=\cosh(y_2)\implies \cosh(y_1) = \cosh(y_2).\quad\quad\quad(2-5)

Notice that by (2-3), y_1, y_2 \ge 0.

Cleary, (2-5) is true if and only if y_1 =y_2. For if y_1 \ne y_2, by (\star), \cosh(y_1) \ne \cosh(y_2).


(x, y_1), (x, y_2) \in S_3 \implies y_1 = y_2.

We name the function defined by (2-3) \cosh^{-1} as (2-4) shows that it is the inverse of \cosh(x).

See “Deriving Two Inverse Functions” for the explicit expressions of \sinh^{-1} and \cosh^{-1}.


t_1 \ge 0, t_2 \ge 0, t_1 \ne t_2 \implies \cosh(t_1) \ne \cosh(t_2).\quad\quad\quad(\star)

Without loss of generality, we assume that t_2 > t_1. By Lagrange’s mean-value theorem (see “A Sprint to FTC“),

\cosh(t_2) -\cos(t_1) =\cosh'(\xi) (t_2-t_1)

where \xi \in (t_1, t_2).

We have

t_2 > t_1 \implies t_2 - t_1 >0


\forall t > 0, (\cosh(t))' = (\frac{e^t+e^{-t}}{2})' = \frac{e^t-e^{-t}}{2}  \overset{t>0\implies t >-t \implies e^t  >e^{-t}}{>} 0.


\cosh(t_2) - \cosh(t_1) = \cosh'(\xi) (t_2-t_1) > 0.


\cosh(t_2) \ne \cosh(t_1).

Exercise-1 Show that \frac{e^{y_1}-e^{-y_1}}{2}- \frac{e^{y_2}-e^{-y_2}}{2}= 0 if and only if y_1 = y_2.

Deriving Two Inverse Functions

In “Instrumental Flying“, we defined function y=\sinh^{-1}(x) as

\{(x, y) | x = \frac{e^y-e^{-y}}{2}\}.\quad\quad\quad(1)

From x = \frac{e^y-e^{-y}}{2}, we obtain

(e^y)^2-2x\cdot e^y-1=0.

It means either e^y = x-\sqrt{x^2+1} or e^y = x+\sqrt{x^2+1}.

But e^y = x-\sqrt{x^2+1} suggests e^y < 0 (see Exercise-1), contradicts the fact that \forall t \in R, e^t > 0 (see “Two Peas in a Pod, Part 2“).


e^y = x+\sqrt{x^2+1} \implies y = \log(x + \sqrt{x^2+1}).


\sinh^{-1}(x) = \log(x + \sqrt{x^2+1}), \;\;x \in (-\infty, +\infty).

We also defined a non-negative valued function y = \cosh^{-1}(x):

\{(x, y) | x = \frac{e^y + e^{-y}}{2}, y \ge 0\}.\quad\quad\quad(2)

Simplifying x=\frac{e^y+e^{-y}}{2} yields

(e^y)^2-2x\cdot e^y+1=0.

It follows that either e^y = x-\sqrt{x^2-1} or e^y=x+\sqrt{x^2-1}.

For both expressions’ right-hand side to be valid requires that x \ge 1. However, when x = 2,

e^y = x-\sqrt{x^2-1} = 2 - \sqrt{3} < 1

suggests that y < 0 (see Exercise-2,3), contradicts (2).


e^y = x+\sqrt{x^2-1} \implies y = \log(x+\sqrt{x^2-1}).


\cosh^{-1}(x) = \log(x+\sqrt{x^2-1}), \;\;x \in [1, +\infty).

See also “A Relentless Pursuit“.

Exercise-1 Show that \forall x \in R, x-\sqrt{x^2+1} < 0.

Exercise-2 Without calculator or CAS, show that 2-\sqrt{3} < 1.

Exercise-3 Prove \forall t \ge 0. e^t \ge 1 (hint: see “Two Peas in a Pod, Part 2“)

Beltrami’s Identity

The Beltrami’s identity

F - y' \frac{\partial F}{\partial y'} = C

where C is a constant, is a reduced form of Euler-Lagrange equation for the special case when F does not dependent explicitly on x.

For F = F(y, y'),

\frac{dF}{dx} = \frac{\partial F}{\partial y} y' + \frac{\partial F}{\partial y'} y''.\quad\quad\quad(1)

From Euler-Lagrange equation

\frac{\partial F}{\partial y} - \frac{d}{dx}(\frac{\partial F}{\partial y'}) = 0,

we have

\frac{\partial F}{\partial y} = \frac{d}{dx}(\frac{\partial F}{\partial y'}).\quad\quad\quad(2)

Substituting (2) into (1) gives

\frac{dF}{dx} = \frac{d}{dx}(\frac{\partial F}{\partial y'}) y' + \frac{\partial F}{\partial y'} y''.\quad\quad\quad(3)

Consequently, when (3) is expressed as

\frac{dF}{dx} = \frac{d}{dx}(\frac{\partial F}{\partial y'}) y' + \frac{\partial F}{\partial y'} \frac{dy'}{dx},

it becomes clear that

\frac{dF}{dx} = \frac{d}{dx}(\frac{\partial F}{\partial y'}y').


\frac{d}{dx}(F-\frac{\partial F}{\partial y'}y') = 0.


F -  y' \frac{\partial F}{\partial y'} = C.

In “A Relentless Pursuit”, we derived differential equation

\frac{dy}{dx} = \frac{\sqrt{y^2-C_1^2}}{\pm C_1}\quad\quad\quad(4)

from Euler-Lagrange equation without taking advantage of the fact that F=y\sqrt{1+(y')^2} has no explicit dependency on x. The derivation was mostly done by a CAS.

Let’s derive (4) from Beltrami’s Identity. This time, we will not use CAS.

For F=y\sqrt{1+(y')^2}, the Beltrami’s Identity is

y\sqrt{1+(y')^2} - y'\cdot\left(y\cdot\frac{1}{2}\frac{2y'}{\sqrt{1+(y')^2}}\right) = C.

It simplifies to

\frac{y}{\sqrt{1+(y')^2}} = C.

Further simplification yields

C^2(y')^2 = y^2-C^2.

For C \ne 0,

(y')^2 = \frac{y^2-C^2}{C^2}.


\frac{dy}{dx} = \pm \sqrt{\frac{y^2-C^2}{C^2}}=\frac{\sqrt{y^2-C^2}}{\pm C}.

Prequel to “A Relentless Pursuit”

Fig. 1

Illustrated in Fig. 2 are two circular hoops of unit radius, centered on a common x-axis and a distance 2a apart. There is also a soap films extends between the two hoops, taking the form of a surface of revolution about the x-axis. If gravity is negligible, the film takes up a state of stable, equilibrium in which its surface area is a minimum.

Fig. 2

Our problem is to find the function y(x), satisfying the boundary conditions

y(-a) = y(a) = 1,\quad\quad\quad(1)

which makes the surface area


a minimum.


F(x,y, y') = 2\pi y \sqrt{1+(y')^2}.

We have

\frac{\partial F}{\partial y} = 2\pi \sqrt{1+(y')^2}


\frac{\partial F}{\partial y'} = 2 \pi y \cdot\frac{1}{2}\left(1+(y')^2\right)^{-\frac{1}{2}}\cdot 2y'=\frac{2 \pi y y'}{\sqrt{1+(y')^2}}.

The Euler-Lagrange equation

\frac{\partial F}{\partial y} - \frac{d}{dx}\left(\frac{\partial F}{\partial y'}\right) = 0


\sqrt{1+(y')^2} - \frac{d}{dx}\left(\frac{y y'}{\sqrt{1+(y')^2}}\right) = 0.

Fig. 3

Using Omega CAS Explorer (see Fig. 3), it can be simplified to:

y \frac{d^2 y}{dx^2}- \left(\frac{dy}{dx}\right)^2=1.

This is the differential equation solved in “A Relentless Pursuit” whose solution is

y = C_1\cdot \cosh(\frac{x+C_2}{C_1}).

We must then find C_1 and C_2 subject to the boundary condition (1), i.e.,

C_1\cdot \cosh(\frac{a+C_2}{C_1}) = C_1\cdot\cosh(\frac{-a+C_2}{C_1})\implies \cosh(\frac{a+C_2}{C_1}) = \cosh(\frac{-a+C_2}{C_1}).

The fact that \cosh is an even function implies either

a+C_2 = -a+C_2\quad\quad\quad(3)


a+C_2 = -(-a+C_2).\quad\quad\quad(4)

While (3) is clearly false since it claims for all a >0, a = -a, (4) gives


And so, the solution to boundary-value problem

\begin{cases} y \frac{d^2 y}{dx^2}- \left(\frac{dy}{dx}\right)^2=1,\\ y(-a)=y(a)=1\end{cases}\quad\quad\quad(5)


y = C_1\cdot \cosh(\frac{x}{C_1}).\quad\quad\quad(6)

To determine C_1, we deduce the following equation from the boundary conditions that y=1 at x=\pm a:

C_1\cdot \cosh(\frac{a}{C_1}) = 1.\quad\quad\quad(7)

This is a transcendental equation for C_1 that can not be solved explicitly. Nonetheless, we can examine it qualitatively.


\mu = \frac{a}{C_1}

and express (7) as

\cosh(\mu) = \frac{\mu}{a}.\quad\quad\quad(8)

Fig. 4

A plot of (8)’s two sides in Fig. 4 shows that for sufficient small a, the curves z = \cosh(\mu) and z = \frac{\mu}{a} will intersect. However, as a increases, z=\frac{\mu}{a}, a line whose slope is \frac{1}{a} rotates clockwise towards \mu-axis. The curves will not intersect if a is too large. The critical case is when a=a^*, the curves touch at a single point, so that

\cosh(\mu) = \frac{\mu}{a^*}\quad\quad\quad(9)

and y=\frac{\mu}{a} is the tangent line of z=\cosh(\mu), i.e.,

\sinh(\mu) = \frac{1}{a^*}.\quad\quad\quad(10)

Dividing (9) by (10) yields

\coth(\mu) = \mu. \quad\quad\quad(11)

What the mathematical model (5) predicts then is, as we gradually move the rings apart, the soap film breaks when the distance between the two rings reaches 2a^*, and for a > a^*, there is no more soap film surface connects the two rings. This is confirmed by an experiment (see Fig. 1).

We compute the value of a^*, the maximum value of a that supports a minimum area soap film surface as follows.

Fig. 5

Solving (11) for \mu numerically (see Fig. 5), we obtain

\mu = 1.1997.

By (10), the corresponding value of

a^* = \frac{1}{\sinh(\mu)} = \frac{1}{\sinh(1.1997)} = 0.6627.

We also compute the surface area of the soap film from (2) and (6) (see Fig. 6). Namely,

A = 2\pi \displaystyle\int\limits_{-a}^{a} C_1 \cosh\left(\frac{x}{C_1}\right) \sqrt{1+\left(\frac{d}{dx}C_1\cosh\left(\frac{x}{C_1}\right)\right)^2}\;dx =  \pi C_1^2\left(\sinh\left(\frac{2a}{C_1}\right) + \frac{2a}{C_1}\right).

Fig. 6

Exercise-1 Given a=\frac{1}{2}, solve (7) numerically for C_1.

Exercise-2 Without using a CAS, find the surface area of the soap film from (2) and (6).

An Epilogue to “A Relentless Pursuit”



differential equation (1) in “A Relentless Pursuit” can be expressed as

\frac{dp}{dy}-\frac{1}{y}p = \frac{1}{y}p^{-1}.

This is Bernoulli’s equation

\frac{dp}{dy} + f(y) p = g(y) p^{\alpha}

with f(y) = -\frac{1}{y}, g(y) = \frac{1}{y} and \alpha = -1 (see “Meeting Mr. Bernoulli“).


p^{1-\alpha} = e^{(\alpha-1)\displaystyle\int f(y)\;dy}\left((1-\alpha)\displaystyle\int e^{-(\alpha-1)\displaystyle\int f(y)\;dy} g(y)\;dy + C\right).\quad\quad\quad(1)

Substitute f(y), g(y) and \alpha into (1) gives

p^{1-(-1)} = e^{(-1-1)\displaystyle\int -\frac{1}{y}\; dy}\left((1-(-1))\displaystyle\int e^{-(-1-1)\displaystyle\int -\frac{1}{y}\;dy}\frac{1}{y}\;dy+C\right).


p^2 = Cy^2-1

where it must be true that C>0. Therefore,

p^2=\frac{1}{C_1^2}y^2-1, C_1>0\implies p^2 = \frac{y^2-C_1^2}{C_1^2} \implies p \overset{C_1>0}{=} \frac{\sqrt{y^2-C_1^2}}{\pm C_1}.

That is,

\frac{dy}{dx} = \frac{\sqrt{y^2-C_1^2}}{\pm C_1}.\quad\quad\quad(2)

There is yet another way to obtain (2):

Since equation (1) in “A Relentless Pursuit“, namely,

yy''-(y')^2 = 1 \implies 1+(y')^2 = y y''\quad\quad\quad(3)


(1+(y')^2)'= 2y' y''.\quad\quad\quad(4)

(4)/(3) yields

\frac{(1+(y')^2)'}{1+(y')^2} = \frac{2y'}{y}.\quad\quad\quad(5)

Integrate (5) with respect to y, we have

\frac{1}{2}\log(1+(y')^2) = \log(y) + C \implies \log(\frac{\sqrt{1+(y')^2}}{y}) = C.


\frac{\sqrt{1+(y')^2}}{y}= e^C=\frac{1}{C_1}

where C_1>0.

It follows that

\frac{dy}{dx} = \frac{\sqrt{y^2-C_1^2}}{\pm C_1}.