# How to solve a quartic equation

Begin with a general quartic equation

$x^4+bx^3+cx^2+dx+e=0,\quad\quad\quad(*)$

depress it using the substitution

$x = y - \frac{b}{4},$

creates a new quartic equation (without the $y^3$ term) in $y:$

$y^4+qy^2+ry+s=0\quad\quad\quad(**)$

where

$\begin{cases}q=c-\frac{3b^2}{8}\\r=d-\frac{bc}{2}+\frac{b^3}{8}\\s=e-\frac{bd}{4}+\frac{b^2c}{16}-\frac{3b^4}{256}\end{cases}$

When $r=0,$ i.e., $d-\frac{bc}{2}+\frac{b^3}{8}=0,$ the quartic equation in $y$ becomes a quadratic equation in $y^2:$

$(y^2)^2+qy^2+s=0.$

It follows that

$y^2 = \frac{-q\pm \sqrt{q^2-4s}}{2} \implies y=\pm\sqrt{\frac{-q\pm\sqrt{q^2-4s}}{2}}\implies x=y-\frac{b}{2}.$

For example, to solve quartic equation $x^4+2x^3-3x^2-4x+4=0$, depress it using the substitution $x = y-\frac{1}{2},$ we obtain

$y^4-\frac{9}{2}y^2+\frac{81}{16}=0.$

This is a quadratic equation in $y^2:$

$(y^2)^2-\frac{9}{2}y^2+\frac{81}{16}=0.$

Therefore,

$y^2 = \frac{\frac{9}{2} \pm \sqrt{(\frac{9}{4})^2-4\cdot 1\cdot \frac{81}{16}}}{2}=\frac{9}{4}\implies y = \pm\sqrt{\frac{9}{4}} = \pm\frac{3}{2}.$

Consequently,

$x_1 = \frac{3}{2}-\frac{1}{2} = 1, x_2=-\frac{3}{2}-\frac{1}{2} = -2.$

In general, we solve the depressed quartic equation as follows:

Substituting $u+v+w$ for $y$ in the left side of (**) yields

$y^4+qy^2+ry+s$

By $\begin {cases} y^2 = (u+v+w)^2=u^2+v^2+w^2+2(uv+uw+vw),\\ y^4 = (y^2)^2 = \left((u^2+v^2+w^2)+2(uv+uw+vw)\right)^2 \\= (u^2+v^2+w^2)^2+4(u^2+v^2+w^2)(uv+uw+vw) + 4(uv+uw+vw)^2 \\= (u^2+v^2+w^2)^2+2(u^2+v^2+w^2)\cdot 2(uv+uw+vw) + 4(uv+uw+vw)^2\end {cases}$

$= \underbrace{(u^2+v^2+w^2)^2 + 4(u^2+v^2+w^2)(uv+uw+vw)+4(uv+uw+vw)^2}_{y^4}+ q\underbrace{(u^2+v^2+w^2+2(uv+uw+vw))}_{y^2}+ r\underbrace{(u+v+w)}_{y}+s$

$= (u^2+v^2+w^2)^2+4(u^2+v^2+w^2)(uv+uw+vw)+4(uv+uw+vw)^2 + +q(u^2+v^2+w^2)+2q(uv+vw+uw) + r(u+v+w)+s$

$= (u^2+v^2+w^2) +q(u^2+v^2+w^2) + 4(u^2+v^2+w^2)(uv+uw+vw) + 2q(uv+uw+vw) +4(uv+uw+vw)^2 + r(u+v+w)+s$

$= (u^2+v^2+w^2) + q(u^2+v^2+w^2) + 2(uv+uw+vw)(2(u^2+v^2+w^2) + q) +4(u^2v^2+u^2w^2+v^2w^2 + 2u^2vw + 2v^2uw +2w^2uv) + q(u^2+v^2+w^2)+r(u+v+w)+s$

$= (u^2+v^2+w^2) + q(u^2+v^2+w^2) + 2(uv+uw+vw)(2(u^2+v^2+w^2) + q) +4(u^2v^2+u^2w^2+v^2w^2) + 8u^2vw + 8v^2uw +8w^2uv) + q(u^2+v^2+w^2)+r(u+v+w)+s$

$= u^2+v^2+w^2) + q(u^2+v^2+w^2) + 2(uv+uw+vw)(2(u^2+v^2+w^2) + q) +4(u^2v^2+u^2w^2+v^2w^2) + 8u^2vw + 8v^2uw +8w^2uv + r(u+v+w)+s$

$=(u^2+v^2+w^2) + q(u^2+v^2+w^2) + 2(uv+uw+vw)(2(u^2+v^2+w^2) + q) +4(u^2v^2+u^2w^2+v^2w^2) + 8uvw(u + v + w) + r(u+v+w)+s$

$= (u^2+v^2+w^2) + q(u^2+v^2+w^2) + 2(uv+uw+vw)(\underline{2(u^2+v^2+w^2) + q}) +4(u^2v^2+u^2w^2+v^2w^2) + (\underline{8uvw+r})(u + v + w) +s$

If $\begin{cases} 2(u^2+v^2+w^2) + q=0 \\ 8uvw + r=0\end{cases}\quad\quad\quad(1)$

$=(u^2+v^2+w^2)^2+q(u^2+v^2+w^2) + 4(u^2v^2+u^2w^2+v^2w^2)+s$

Since $(1) \implies \begin{cases} u^2+v^2+w^2 = -\frac{q}{2}\\ uvw = -\frac{r}{8}\end{cases}\quad\quad\quad(2)$

$= (-\frac{q}{2})^2 + q \cdot (-\frac{q}{2}) + 4(u^2v^2+u^2w^2+v^2w^2)+s$

$= \frac{q^2}{4} - \frac{q^2}{2} + 4(u^2v^2+u^2w^2+v^2w^2)+s$

$= -\frac{q^2}{4} + s + 4(u^2v^2+u^2w^2+v^2w^2)$

$= -\frac{q^2-4s}{4} + 4(\underline{u^2v^2+u^2w^2+v^2w^2})$

If $u^2v^2 + u^2w^2+v^2w^2 = \frac{q^2-4s}{16}\quad\quad\quad(3)$

$= 0.$

It means that $u+v+w$ is a solution of $(**)$ if

$\begin{cases} u^2+v^2+w^2=-\frac{q}{2} \\ u^2v^2+u^2v^2+v^2w^2=\frac{q^2-4s}{16}\\ uvw=-\frac{r}{8}\end{cases}\quad\quad\quad(4-1, 4-2, 4-3)$

Moreover, squaring (4-3) gives

$\begin{cases} u^2+v^2+w^2=-\frac{q}{2} \\ u^2v^2+u^2v^2+v^2w^2=\frac{q^2-4s}{16}\\ u^2v^2w^2=\frac{r^2}{64}\end{cases}\quad\quad\quad(5-1, 5-2, 5-3)$

By Vieta’s theorem, $u^2, v^2, w^2$ satisfying (5-1, 5-2, 5-3) are the three solutions of cubic equation

$z^3+\frac{q}{2}z^2+\frac{q^2-4s}{16}z-\frac{r^2}{64}= 0.$

(See “How to solve a cubic equation“)

Suppose the three solutions are $z_1, z_2, z_3,$ we have

$u^2=z_1 \implies u = \pm\sqrt{z_1}, v^2 = z_2 \implies v = \pm\sqrt{z_2}, w^2=z_3 \implies w=\pm \sqrt{z_3}.$

Clearly, there are eight combinations of $u,v,w:$

$\begin{cases} u=\sqrt{z_1} \\ v= \sqrt{z_2} \\w=\sqrt{z_3} \end{cases} \begin{cases} u=\sqrt{z_1} \\ v= \sqrt{z_2} \\w=-\sqrt{z_3} \end{cases} \begin{cases} u=\sqrt{z_1} \\ v= -\sqrt{z_2} \\w=\sqrt{z_3} \end{cases} \begin{cases} u=\sqrt{z_1} \\ v= -\sqrt{z_2} \\w=-\sqrt{z_3} \end{cases}$

$\begin{cases} u=-\sqrt{z_1} \\ v= \sqrt{z_2} \\w=\sqrt{z_3} \end{cases} \begin{cases} u=-\sqrt{z_1} \\ v= \sqrt{z_2} \\w=-\sqrt{z_3} \end{cases} \begin{cases} u=-\sqrt{z_1} \\ v= -\sqrt{z_2} \\w=\sqrt{z_3} \end{cases} \begin{cases} u=-\sqrt{z_1} \\ v= -\sqrt{z_2} \\w=-\sqrt{z_3} \end{cases}$

Among them, only four are valid due to constraint (4-3) placed on the product $u\cdot v\cdot w.$

From (5-3), we see $uvw = \sqrt{z_1}\sqrt{z_2}\sqrt{z_3}= \pm\frac{r}{8}.$

If $\sqrt{z1}\sqrt{z_2}\sqrt{z_3} = +\frac{r}{8},$ the valid ones are:

$\begin{cases} u=\sqrt{z_1} \\ v= \sqrt{z_2} \\w=-\sqrt{z_3} \end{cases}$ since $uvw=\sqrt{z_1}\cdot \sqrt{z_2}\cdot -\sqrt{z_3}=-\sqrt{z_1}\sqrt{z_2}\sqrt{z_3} = -\frac{r}{8}$

$\begin{cases} u=\sqrt{z_1} \\ v= -\sqrt{z_2} \\w=\sqrt{z_3} \end{cases}$ since $uvw=\sqrt{z_1}\cdot -\sqrt{z_2}\cdot \sqrt{z_3}=-\sqrt{z_1}\sqrt{z_2}\sqrt{z_3} = -\frac{r}{8}$

$\begin{cases} u=-\sqrt{z_1} \\ v= \sqrt{z_2} \\w=\sqrt{z_3} \end{cases}$ since $uvw=-\sqrt{z_1}\cdot \sqrt{z_2}\cdot \sqrt{z_3}=-\sqrt{z_1}\sqrt{z_2}\sqrt{z_3} = -\frac{r}{8}$

$\begin{cases} u=-\sqrt{z_1} \\ v=-\sqrt{z_2} \\w=-\sqrt{z_3} \end{cases}$ since $uvw=-\sqrt{z_1}\cdot -\sqrt{z_2}\cdot -\sqrt{z_3}=-\sqrt{z_1}\sqrt{z_2}\sqrt{z_3} = -\frac{r}{8}$

Otherwise $\sqrt{z_1}\sqrt{z_2}\sqrt{z_3}=-\frac{r}{8}$, we have

$\begin{cases} u=\sqrt{z_1} \\ v= \sqrt{z_2} \\w=\sqrt{z_3} \end{cases} \begin{cases} u=\sqrt{z_1} \\ v= -\sqrt{z_2} \\w=-\sqrt{z_3} \end{cases} \begin{cases} u=-\sqrt{z_1} \\ v= \sqrt{z_2} \\w=-\sqrt{z_3} \end{cases} \begin{cases} u=-\sqrt{z_1} \\ v= -\sqrt{z_2} \\w=\sqrt{z_3} \end{cases}$

Consequently, a solution to the general quartic equation (*) is

$x = u+v+w - \frac{b}{4}.$

There are four such solutions.

Exercise-1 Show that (5-3) $\implies uvw = \sqrt{z_1}\sqrt{z_2}\sqrt{z_3}= \pm\frac{r}{8}.$

# How to solve a cubic equation

Consider the general cubic equation

$x^3+bx^2+cx+d=0.\quad\quad\quad(*)$

Let

$y = x + \omega \implies x = y - \omega.\quad\quad\quad(1)$

Substituting (1) into (*), we have

$(y-\omega)^3 + b(y-\omega)^2+c(y-\omega) + d = 0.$

This is a new cubic equation in $y$:

$y^3 + (b-3\omega)y^2 + (3\omega^2-2b\omega+c)y - \omega^3+b\omega^2-c\omega+d=0\quad\quad\quad(2)$

Let

$\omega = \frac{b}{3}\quad\quad\quad(3)$

so that (2) becomes

$y^3 + (c-\frac{b^2}{3})y+d-\frac{bc}{3}+\frac{2b^3}{27}=0.$

That is

$y^3+py=q\quad\quad\quad(4)$

where $p = c-\frac{b^2}{3}, q = -d+\frac{bc}{3}-\frac{2b^3}{27}$ (see Fig. 1).

Thus, (4) is the so called depressed cubic equation. A solution is readily available through Cardano’s formula (see “Through the Mind’s Eye“).

Once we know $y$, (1) gives

$x = y - \frac{b}{3}.$

Fig. 1

# Wallis’ Pi

There is a remarkable expression for the number $\pi$ as an infinite product. Starting with definite integral $\int\limits_{0}^{\frac{\pi}{2}} \sin^{m}(x)\;dx, m=0,1,2,3,4...$, we derive it as follows:

$\int\limits_{0}^{\frac{\pi}{2}} \sin^m(x)\;dx$

$=\int\limits_{0}^{\frac{\pi}{2}}\sin^{m-1}(x)\cdot\sin(x)\;dx$

$= \int\limits_{0}^{\frac{\pi}{2}} \sin^{m-1}(x)\cdot(-\cos(x))'\;dx$

By $\int\limits_{a}^{b} u\cdot v'\;dx = u\cdot v \bigg|_{b}^{a}- \int\limits_{b}^{a} u' \cdot v\;dx,$

$= \sin^{m-1}(x)\cdot(-\cos(x))\bigg|_{0}^{\frac{\pi}{2}} - \int\limits_{0}^{\frac{\pi}{2}}(\sin^{m-1}(x))'\cdot(-\cos(x))\;dx$

$= 0 - \int\limits_{0}^{\frac{\pi}{2}}(m-1)\sin^{m-2}(x)\cos(x)\cdot(-\cos(x))\;dx$

$= +\int\limits_{0}^{\frac{\pi}{2}}(m-1)\sin^{m-2}(x)\cdot\cos^2(x)\;dx$

$\overset{\cos^2(x) = 1-\sin^2(x)}{=}\int\limits_{0}^{\frac{\pi}{2}}(m-1)\sin^{m-2}(x)\cdot (1-\sin^{2}(x))\;dx$

$= \int\limits_{0}^{\frac{\pi}{2}}(m-1)\sin^{m-2}(x) - (m-1)\sin^{m}(x)\;dx$

$= (1-m)\int\limits_{0}^{\frac{\pi}{2}}\sin^{m}(x)\;dx + (m-1)\int\limits_{0}^{\frac{\pi}{2}}\sin^{m}(x)\;dx.$

That is,

$m\int\limits_{0}^{\frac{\pi}{2}}\sin^{m}(x)\;dx = (m-1)\int\limits_{0}^{\frac{\pi}{2}}\sin^{m-2}(x)\;dx.$

And so,

$\int\limits_{0}^{\frac{\pi}{2}}\sin^{m}(x) \; dx= \frac{m-1}{m}\int\limits_{0}^{\frac{\pi}{2}}\sin^{m-2}(x)\;dx.\quad\quad\quad(1)$

By repeated application of (1) we have the following values for $I_{m} = \int\limits_{0}^{\frac{\pi}{2}} \sin^{m}(x)\;dx$:

For even $m$,

$I_{2n} = \left(\prod\limits_{k=0}^{n-1}\frac{2n-2k-1}{2n-2k}\right)I_0\overset{I_0 = \int\limits_{0}^{\frac{\pi}{2}}\;dx = \frac{\pi}{2}}{\implies}I_{2n}= \left(\prod\limits_{k=0}^{n-1}\frac{2n-2k-1}{2n-2k}\right)\cdot \frac{\pi}{2}.$

Similarly, for odd $m$,

$I_{2n+1} = \left(\prod\limits_{k=0}^{n-1}\frac{2n+1-2k-1}{2n+1-2k}\right)\cdot I_{1}\overset{I_1 = \int\limits_{0}^{\frac{\pi}{2}}\sin^{1}(x)\;dx = 1}{\implies} I_{2n+1}= \prod\limits_{k=0}^{n-1}\frac{2n-2k}{2n-2k+1}.$

i.e.,

$I_{2n} = \frac{2n-1}{2n}\cdot\frac{2n-3}{2n-2}\cdot\frac{2n-5}{2n-4}\cdot\frac{2n-7}{2n-6}\cdot...\cdot\frac{5}{6}\cdot\frac{3}{4}\cdot\frac{1}{2}\cdot\frac{\pi}{2}\quad\quad\quad(2)$

$I_{2n+1} =\frac{2n}{2n+1}\cdot\frac{2n-2}{2n-1}\cdot\frac{2n-4}{2n-3}\cdot\frac{2n-6}{2n-5}\cdot ...\cdot\frac{6}{7}\cdot\frac{4}{5}\cdot\frac{2}{3}\quad\quad\quad\quad\quad(3)$

Since for $0 we have $\sin^{2n-1}(x) > \sin^{2n}(x) > \sin^{2n+1}(x).$ It means

$\int\limits_{0}^{\frac{\pi}{2}}\sin^{2n-1}(x) > \int\limits_{0}^{\frac{\pi}{2}}\sin^{2n}(x)\;dx > \int\limits_{0}^{\frac{\pi}{2}}\sin^{2n+1}(x)\;dx>0$

or,

$I_{2n-1} > I_{2n} > I_{2n+1}>0.$

Hence,

$\frac{I_{2n-1}}{I_{2n+1}}>\frac{I_{2n}}{I_{2n+1}} > 1.\quad\quad\quad(4)$

Moreover, we have

$I_{2n+1} \overset{(1)}{=} \frac{(2n+1)-1}{2n+1}I_{(2n+1)-2} = \frac{2n}{2n+1} I_{2n-1},$

so that

$\frac{I_{2n-1}}{I_{2n+1}} = \frac{2n+1}{2n}.\quad\quad\quad(5)$

And,

$\frac{I_{2n}}{I_{2n+1}} \overset{(2), (3)}{=} \frac{(2n+1)\cdot(2n-1)^2\cdot (2n-3)^2...7^2\cdot 5^2\cdot 3^2}{(2n)^2\cdot (2n-2)^2\cdot (2n-4)^2\cdot ...\cdot 6^2\cdot 4^2\cdot 2^2}\cdot \frac{\pi}{2}$

$= \frac{(2n)^2\cdot (2n-1)^2\cdot (2n-2)^2\cdot (2n-3)^2\cdot (2n-4)^2\cdot...\cdot 7^2\cdot 6^2\cdot 5^2\cdot 4^2\cdot 3^2\cdot 2^2}{(2n)^4(2n-2)^4\cdot (2n-4)^4...6^4\cdot 4^4\cdot 2^4}\cdot (2n+1)\cdot\frac{\pi}{2}$

$= \frac{\left((2n)\cdot (2n-1)\cdot (2n-2)\cdot (2n-3)\cdot (2n-4)\cdot...\cdot 7\cdot 6\cdot 5\cdot 4\cdot 3\cdot 2\right)^2}{(2^4\cdot n^4)\cdot (2^4\cdot(n-1)^4) \cdot (2^4\cdot (n-2)^4)\cdot...\cdot (2^4\cdot 3^4)\cdot (2^4\cdot 2^4)\cdot 2^4}\cdot (2n+1)\cdot\frac{\pi}{2}$

$= \frac{((2n)!)^2}{\underbrace{2^4\cdot 2^4\cdot ...\cdot 2^4}_{n 2^4s}\cdot (n\cdot(n-1)\cdot(n-2)\cdot ...\cdot 3\cdot 2\cdot 1)^4}(2n+1)\cdot\frac{\pi}{2}$

$= \frac{((2n)!)^2\cdot (2n+1)}{2^{4n}\cdot(n!)^4}\cdot \frac{\pi}{2}$

gives

$\frac{I_{2n}}{I_{2n+1}} = \frac{((2n)!)^2(2n+1)}{2^{4n}(n!)^4}\cdot\frac{\pi}{2}.\quad\quad\quad(6)$

Substituting (5) and (6) into (4) yields

$\frac{2n+1}{2n} >\frac{((2n)!)^2(2n+1)}{2^{4n}(n!)^4}\cdot\frac{\pi}{2}>1.$

Since $\lim\limits_{n \rightarrow \infty} \frac{2n+1}{2n} = 1, \lim\limits_{n \rightarrow \infty}1 = 1,$ by Squeeze Theorem for Sequences,

$\lim\limits_{n \rightarrow \infty}\frac{((2n)!)^4(2n+1)}{2^{4n}(n!)^4} \cdot \frac{\pi}{2}= 1\implies \lim\limits_{n \rightarrow \infty}\frac{((2n)!)^2(2n+1)}{2^{4n}(n!)^4}=\frac{2}{\pi}.$

Consequently,

$\lim\limits_{n \rightarrow \infty} \frac{2^{4n}(n!)^4}{((2n)!)^2(2n+1)}= \frac{\pi}{2},$

i.e.,

$\pi = 2\cdot\lim\limits_{n \rightarrow \infty}\frac{2^{4n}(n!)^4}{((2n)!)^2(2n+1)}.$

This is Wallis’ product representation for $\pi,$ named after John Wallis who discovered it in 1665.

Maxima knows Wallis’ Pi:

Fig. 1

So does Mathematica:

Fig. 2

Its convergence to $\pi$ is illustrated in Fig. 3:

Fig. 3

# A Cautionary Tale of Compute Inverse Trigonometric Functions

From $\arcsin(x)$‘s definition:

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

we have

$\arcsin(0) = 0, \arcsin(\frac{1}{2})= \frac{\pi}{6},\arcsin(\frac{\sqrt{3}}{2}) = \frac{\pi}{3}, \arcsin(1) = \frac{\pi}{2}$

and,

$\arcsin(-\frac{1}{2})= -\frac{\pi}{6},\arcsin(-\frac{\sqrt{3}}{2}) = -\frac{\pi}{3}, \arcsin(-1) = -\frac{\pi}{2}.$

To obtain other values of $\arcsin(x)$, we may simply solve

$\sin(y) - x = 0$

for $y$ where $-1 \le x \le 1.$

For example (see Fig. 1), solving $\sin(y) = \frac{1}{\sqrt{2}}$ for $y$ gives $y \approx \frac{\pi}{4}.$ It is in agreement with the fact that $\arcsin(\frac{1}{\sqrt{2}})= \frac{\pi}{4}.$

Fig. 1

In Fig. 2, we compute $\arcsin(x)$ from repeatedly solving $\sin(y) = x$ for $y$ where

$x=-1+i \cdot\frac{1-(-1)}{n}, i=1,2,...n.$

Fig. 2

Since $\sin(y)$ is a periodic function, $\sin(y) = x$ has infinitely many solutions. It is possible that the solution obtained by Newton’s method lies outside of $[-\frac{\pi}{2}, \frac{\pi}{2}]$, the range of $\arcsin(x)$ by its definition. Such solution cannot be considered the value of $\arcsin(x).$

Fig. 3

Exercise-1 Compute $\arccos(x)$ by solving $\cos(y) = x$ for $y.$

Exercise-2 Explain Fig. 3.

# A Mathematical Allegory

We have defined function $\arcsin$ 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}$

and

$\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.$

Therefore,

$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

Putting the results together, we have

Fig. 6 $-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}$.

# Oasis

An oasis awaits

The above image is created by Omega CAS Explorer:

The governing equations are:

$x_{n+1} = \rho + c_2(x_n\cos(t_n)-y_n\sin(t_n),$

$y_{n+1} = c_2(x_n\sin(t_n) + y_n\cos(t_n))$

where $t_n = c_1-\frac{c_3}{1+x_n^2 + y_n^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.

Fig.3

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

# DIY

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: