A Mathematical Allegory

We have defined function y = \arcsin(x) as a set:

\{(x, y) | \sin(y) =x, \frac{-\pi}{2} \le y \le \frac{\pi}{2}\}.

By definition,

\arcsin(-1) = \frac{-\pi}{2}, \arcsin(0)=0, \arcsin(1) = \frac{\pi}{2}


\frac{d}{dx} \arcsin(x) = \frac{1}{\sqrt{1-x^2}}.

It means that \arcsin(x) is the unique solution of

\frac{dy}{dx} = \frac{1}{\sqrt{1-x^2}}\quad\quad\quad(\star)

where y(-1)=-\frac{\pi}{2}, y(0)=0 and y(1)=\frac{\pi}{2}.

To compute \arcsin(x), we solve (\star) for y(x) as follows:

Integrate \frac{dy}{dx} = \frac{1}{\sqrt{1-x^2}} from -1 to x gives

\displaystyle\int\limits_{-1}^{x}\frac{dy}{dx} \,dx=\displaystyle\int\limits_{-1}^{x}\frac{1}{\sqrt{1-\xi^2}}\; d\xi\overset{\textbf{FTC}}{\implies}y(x) - y(-1) = \displaystyle\int\limits_{-1}^{x} \frac{1}{\sqrt{1-\xi^2}}\, d\xi.


y(x) = \displaystyle\int\limits_{-1}^{x}\frac{1}{\sqrt(1-\xi^2}\,d\xi + y(-1) \overset{y(-1)=\frac{-\pi}{2}}{=} \displaystyle\int\limits_{-1}^{x}\frac{1}{\sqrt{1-\xi^2}}\,d\xi - \frac{\pi}{2}.

That is,

\arcsin(x) = \displaystyle\int\limits_{-1}^{x} \frac{1}{\sqrt{1-\xi^2}}\;d\xi - \frac{\pi}{2}.

To obtain \arcsin(x), -1 < x < 1, we numerically evaluate \int\limits_{-1}^{x}\frac{1}{\sqrt{1-\xi^2}}\,d\xi, using function ‘quad_qags’.

Fig. 1

The result is visually validated in Fig. 2.

Fig. 2

Note: ‘romberg’, another function that computes the numerical integration by Romberg’s method will not work since it evaluates \frac{1}{\sqrt{1-x^2}} at x=-1.

Fig. 3

An alternate approach is to solve \frac{dy}{dx} = \frac{1}{\sqrt{1-x^2}}, y(0)=0 as an initial-value problem of ODE using ‘rk’ , the function that implements the classic Runge-Kutta algorithm.

Fig. 4 for 0 < x < 1

Fig. 5 -1<x<0

Putting the results together, we have

Fig. 6 -1<x<1

However, we cannot solve \frac{dy}{dx} = \frac{1}{\sqrt{1-x^2}}, y(-1)=\frac{-\pi}{2} using ‘rk’:

Fig. 7

Exercise-1 Compute \arccos(x) for x \in (-1, 1).

Exercise-2 Explain why ‘rk’ cannot solve \frac{dy}{dx} = \frac{1}{\sqrt{1-x^2}}, y(-1)=\frac{-\pi}{2}.

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:

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


Problem Given

f(x) = e^x + \int\limits_{0}^{x} (t-x)f(t)\;dt\quad\quad\quad(\star)

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


From (\star), we see that

f(0) = 1;

f(x) = e^x + \int\limits_{0}^{x} t\cdot f(t) - x\cdot f(t) \;dt = e^x + \int\limits_{0}^{x} t\cdot f(t)\;dt-x\cdot \int\limits_{0}^{x}f(t)\;dt.

And so,

\frac{df(x)}{dx}=\frac{de^x}{dx} + \frac{d}{dx}\int\limits_{0}^{x}tf(t)\;dt - \frac{d}{dx}\left(x\cdot \int\limits_{0}^{x}f(t)\;dt\right)

=e^x+\frac{d}{dx}\int\limits_{0}^{x}tf(t)\;dt-\left(\int\limits_{0}^{x}f(t)\;dt + x\frac{d}{dx}\int\limits_{0}^{x}f(t)\;dt\right)

\overset{\textbf{FTC}}{=}e^x + xf(x) -\int\limits_{0}^{x} f(t)\;dt - xf(x)

= e^x - \int\limits_{0}^{x}f(t)\;dt

That is,

\frac{d}{dx}f(x)= e^x - \int\limits_{0}^{x}f(t)\;dt\implies f'(0) = 1.

It follows that



\begin{cases} f''(x)=e^x-f(x) \\f(0)=1\\f'(0)=1 \end{cases}\quad\quad\quad(\star\star)


f(x) = \frac{1}{2}(\sin(x)+\cos(x)+e^x).

Fig. 1

Notice the derivation of (\star\star) can be simplified if rgwe generalized Leibniz’s Rule (GLR, see “Deriving Generalized Leibniz’s Integral Rule”) is applied:

\frac{df(x)}{dx} = e^x + \underline{\frac{d}{dx}\int\limits_{0}^{x}(t-x)f(t)\;dt}

\overset{\textbf{GLR}}{=} e^x + \underline{(x-x)f(x)\cdot x'-(0-x)f(0)\cdot 0' + \int\limits_{0}^{x}\frac{\partial}{\partial x}(t-x)f(t)\;dt}

=  e^x + \underline{\int\limits_{0}^{x}\frac{\partial}{\partial x}(t-x)f(t)\;dt}

=e^x+\int\limits_{0}^{x}-1\cdot f(t)\;dt

= e^x-\int\limits_{0}^{x}f(t)\;dt

\implies \frac{d^2f(x)}{dx}=e^x-\frac{d}{dx}\int\limits_{0}^{x}f(t)\;dt\overset{\textbf{FTC}}{=}e^x-f(x).

Fig.2 shows that Omega CAS explorer‘s Maxima engine is both FTC and GLR aware:

Fig. 2

Exercise-1 Given:

f(x) = \int\limits_{0}^{x}t\cdot f(x-t)\;dt+\sin(x)

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

hint: Let u=x-t, t = x-u; t=0\implies u=x; t=x\implies u=0; \frac{du}{dt}=-1;

f(x) = \int\limits_{x}^{0}(x-u)\cdot f(u)\cdot (-1)\;du + \sin(x)=\int\limits_{0}^{x}(x-u)f(u)\;du+\sin(x).

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.