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}

This image has an empty alt attribute; its file name is screen-shot-2022-01-07-at-5.40.02-pm.png

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

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<x<\frac{\pi}{2}, 0< \sin(x) < 1, 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<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.

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: