Introduction to Linear Algebra with Mathematica

# Preface

The Mathieu equation is an ordinary differential equation with real coefficients. Its standard form with parameters (a, q) is

$$\label{EqMathieu.1} w'' + \left( a - 2q\,\cos (2x) \right) w =0 .$$
It was introduced by the French mathematician Émile Léonard Mathieu (1835--1890) in 1868 in context of the vibrations of an elliptic membrane. Then the equation arises from the separation of variables in the Helmholtz equation
$\frac{\partial^2 u}{\partial x^2} + \frac{\partial^2 u}{\partial y^2} + k^2 u = 0$
in elliptic coordinates (ξ, η) defined by
$$\label{EqMathieu.2} x = a\,\cosh (\xi )\,\cos (\eta ) , \qquad\qquad y = a\,\sinh (\xi )\,\sin (\eta ) ,$$
where (±𝑎, 0) are the cordinates of the ellipse foci, ξ is a nonnegative real number and η ∈ [ 0, 2 π].

On the complex plane, an equivalent relationship is

$x + {\bf j} y = a\,\cosh \left( \xi + {\bf j} \eta \right) .$
These definitions correspond to ellipses and hyperbolae. The trigonometric identity
$\frac{x^2}{a^2 \cosh^2 \xi} + \frac{y^2}{a^2 \sinh^2 \xi} = \cos^2 \eta + \sin^2 \eta = 1$
shows that curves of constant ξ form ellipses, whereas the hyperbolic trigonometric identity
$\frac{x^2}{a^2 \cos^2 \eta} - \frac{y^2}{a^2 \sin^2 \eta} = \cosh^2 \xi - \sinh^2 \xi = 1 .$
\begin{align*} \nabla^2 \Psi &= \frac{1}{a^2 \left( \sinh^2 \xi + \sin^2 \eta \right)} \left( \frac{\partial^2 \Phi}{\partial \xi^2} + \frac{\partial^2 \Phi}{\partial \eta^2} \right) \\ &= \frac{2}{\cosh 2\xi - \cos 2\eta} \left( \frac{\partial^2 \Phi}{\partial \xi^2} + \frac{\partial^2 \Phi}{\partial \eta^2} \right) . \end{align*}

With[ {$Alpha] = 1}, ParametricPlot[ { \[Alpha]*Cosh[\[Mu]]*Cos[\[Theta]], \[Alpha]*Sinh[\[Mu]]*Sin[\[Theta]] }, {\[Mu], -2, 2}, {\[Theta], 0, Pi}, PlotRange -> 2*{{-1, 1}, {-1, 1}}, Mesh -> True, PlotStyle -> None ] ] Add some styling: With[ {\[Alpha] = 1}, Framed[ ParametricPlot[ { \[Alpha]*Cosh[\[Mu]]*Cos[\[Theta]], \[Alpha]*Sinh[\[Mu]]*Sin[\[Theta]] }, {\[Mu], -2, 2}, {\[Theta], 0, Pi}, PlotRange -> 2*{{-1, 1}, {-1, 1}}, Mesh -> True, PlotStyle -> None, MeshFunctions -> {#3 &, #4 &}, MeshStyle -> { Directive[Opacity[1], Thickness[0.003], Green], Directive[Opacity[1], Thickness[0.003], Blue] }, Epilog -> {Red, PointSize -> 0.02, Point[{{-\[Alpha], 0}, {\[Alpha], 0}}]}, Frame -> False, Ticks -> {None, Range[0, 10, 0.2]}, GridLines -> Automatic ], FrameStyle -> Gray, FrameMargins -> 0 ] ] To create an ellipse, fix the value of ξ: With[ {\[Alpha] = 2, \[Mu] = 2}, ParametricPlot[ {\[Alpha]*Cosh[\[Mu]]*Cos[\[Theta]], \[Alpha]*Sinh[\[Mu]]* Sin[\[Theta]]}, {\[Theta], 0, 2*Pi}, PlotRange -> 20*{{-1, 1}, {-1, 1}} ] ] Likewise, to create one of the normal curves, fix the value of η. Use Show to combine two parametric plots: With[ {\[Alpha] = 1, \[Mu]Fixed = 1}, Framed[ Show[ (* Original plot: *) ParametricPlot[ {\[Alpha]*Cosh[\[Mu]]*Cos[\[Theta]], \[Alpha]*Sinh[\[Mu]]* Sin[\[Theta]]}, {\[Mu], -2, 2}, {\[Theta], 0, Pi}, PlotRange -> 2*{{-1, 1}, {-1, 1}}, Mesh -> True, PlotStyle -> None, MeshFunctions -> {#3 &, #4 &}, MeshStyle -> { Directive[Opacity[1], Thickness[0.003], Green], Directive[Opacity[1], Thickness[0.003], Blue] }, Epilog -> {Red, PointSize -> 0.02, Point[{{-\[Alpha], 0}, {\[Alpha], 0}}]}, Frame -> False, Ticks -> {Automatic, Range[0, 10, 0.2]}, GridLines -> Automatic ] , (* Highlighted ellipse: *) ParametricPlot[ {\[Alpha]*Cosh[\[Mu]Fixed]*Cos[\[Theta]], \[Alpha]* Sinh[\[Mu]Fixed]*Sin[\[Theta]]}, {\[Theta], 0, 2*Pi}, PlotStyle -> Red ] ], FrameStyle -> Gray, FrameMargins -> 0 ] ] Calculate the radius of an ellipse The radius of each ellipse is given by the following code: radius[\[Alpha]_, \[Mu]_, \[Theta]_] := Norm[{\[Alpha]*Cosh[\[Mu]]*Cos[\[Theta]], \[Alpha]*Sinh[\[Mu]]* Sin[\[Theta]]}]; radius[\[Alpha], \[Mu], \[Theta]] Sqrt[Abs[\[Alpha] Cos[\[Theta]] Cosh[\[Mu]]]^2 + Abs[\[Alpha] Sin[\[Theta]] Sinh[\[Mu]]]^2] This expression simplifies at the points where the ellipse crosses the axes (that is, when η = 0, η = π/2, η = π or η = (3π)/2): TableForm[ Join[ {{"Eta", "Radius"}}, Map[{#, radius[\[Alpha], \[Mu], #]} &, {0, Pi/2, Pi, 3*Pi/2}] ] ]  Eta Radius 0 Abs[\[Alpha] Cosh[\[Mu]]] π/2 Abs[\[Alpha] Sinh[\[Mu]]] π Abs[\[Alpha] Cosh[\[Mu]] 3π/2 Abs[\[Alpha] Sinh[\[Mu]] radius[1, 1, 0] // N 1.54308 radius[1, 1, Pi/2] // N 1.1752 Experiment with ξ and η Change the value of ξ: With[ {\[Alpha] = 2}, Manipulate[ ParametricPlot[ { \[Alpha]*Cosh[\[Mu]]*Cos[\[Theta]], \[Alpha]*Sinh[\[Mu]]*Sin[\[Theta]] }, {\[Theta], 0, 2*Pi}, PlotRange -> 20*{{-1, 1}, {-1, 1}} ], {{\[Mu], 1, "Mu"}, 0, Pi} ] ] Change the value of θ: With[ {\[Alpha] = 2}, Manipulate[ ParametricPlot[ { \[Alpha]*Cosh[\[Mu]]*Cos[\[Theta]], \[Alpha]*Sinh[\[Mu]]*Sin[\[Theta]] }, {\[Mu], -100, 100}, PlotRange -> 20*{{-1, 1}, {-1, 1}} ], {{\[Theta], 1, "Theta"}, 0, Pi} ] ] Solving Helmholtz equation in elliptic coordinates: Rewriting the Helmholtz equation \[ \nabla^2 u + k^2 u =0$
in elliptic coordinates, we get

# Mathieu Functions

However, Mathieu's equation and its generalisations are more important than this single application would suggest. This equation is encountered in many different issues in physics, engineering and industry, including the stability of floating ships and railroad trains, the motion of charged particles in electromagnetic Paul traps, the theory of resonant inertial sensors, and many other problem (see expositions by McLachlan and Ruby).

The Mathieu equation has two types of solutions---stable and unstable, depending only on parameters 𝑎 and q and not on the initial conditions.

 Mathieu Stability Chart plots the stable and unstable zones in the plane is known as the Ince--Strutt diagram. y1[x_] = 2*Sqrt[(x*(x - 1)*(x - 4))/(3*x - 8)]; p1 = Plot[y1[x], {x, -1, 0}, PlotStyle -> Thickness[0.01], AspectRatio -> 1, Filling -> Axis, PlotRange -> {-1, 2}]; y2[x_] = .25*((((9 - 4*x)*(13 - 20*x))^{1/2}) - (9 - 4*x)); p2 = Plot[y2[x], {x, -1, 1/4}, PlotStyle -> Thickness[0.01], AspectRatio -> 1, Filling -> Top, PlotRange -> {-1, 2}]; y3[x_] = (1/4)*((9 - 4*x) - ((9 - 4*x)*(13 - 20*x))^(.5)); p3 = Plot[y3[x], {x, 1/4, 13/20}, PlotStyle -> Thickness[0.01], AspectRatio -> 1, Filling -> Top, PlotRange -> {-1, 2}]; y4[x_] = .25*((9 - 4*x) + ((9 - 4*x)*(13 - 20*x))^(.5)); p4 = Plot[y4[x], {x, 1/4, 13/20}, PlotStyle -> Thickness[0.01], AspectRatio -> 1, Filling -> Top, PlotRange -> {-1, 2}]; y5[x_] = ((2*(x - 1)*(x - 4)*(x - 9))/(x - 5))^(.5); p5 = Plot[y5[x], {x, 13/20, 1}, PlotStyle -> Thickness[0.01], AspectRatio -> 1, Filling -> Top, PlotRange -> {-1, 2}]; y6[x_] = 2*((x*(x - 1)*(x - 4))/(3*x - 8))^(.5); p6 = Plot[y6[x], {x, 1, 2}, PlotStyle -> Thickness[0.01], AspectRatio -> 1, PlotRange -> {-1, 2}, Filling -> Top]; Show[p1, p2, p3, p4, p5, p6, PlotRange -> {-1, 2}] Fig.1 Ince-Strutt diagram. Mathematica code.

It is convenient to write Mathieu equation using a slightly different notation from that in Eq. \eqref{EqMathieu.1}, a notation which will be more convenient for our further discussion, namely, as

$$\label{EqMathieu.2} \frac{{\text d}^2 u}{{\text d} t^2} + k \left( 1 - m\,\cos t \right) u = 0.$$
For the great majority of problems associated with the Mathieu equation, the crucial issue is the determination of conditions in which solutions to this equation remain bounded in the course of time, or grow indefinitely. The answer to this question is given by the well-known Ince-Strutt diagram (Figure 1), which shows the transition curves in the parameters plane. These curves divide the plane (k, m) into regions that correspond to unbounded (unstable) solutions (shaded areas in Figure 1) and stable motions.

E.Butikov showed that the transition curves of the Ince-Strutt diagram can be determined by equations:

\begin{align*} m(k) &= 2 \sqrt{\frac{k \left( k-1 \right)\left( k-4 \right)}{3k-8}} \qquad \mbox{for} \quad k < 0; \\ m(k) &= \frac{1}{4} \left( \sqrt{ \left( 9 - 4k \right) \left( 13-20k \right)} - \left( 9 -4k \right) \right) , \qquad \mbox{for} \quad k < \frac{1}{4} ; \\ m(k) &= \frac{1}{4} \left( 9 - 4k \mp \sqrt{\left( 9-4k \right)\left( 13-20k \right)} \right) , \qquad \mbox{for} \quad \frac{1}{4} < k < \frac{13}{20}; \\ m(k) &= \sqrt{\frac{2 \left( k-1 \right)\left( k-4 \right)\left( k-9 \right)}{k-5}} , \qquad \mbox{for} \quad \frac{13}{20} < k < 1 ; \\ m(k) &= 2 \sqrt{\frac{k \left( k-1 \right)\left( k-4 \right)}{3k-8}} , \qquad \mbox{for} \quad k > 1 . \end{align*}
With $$\zeta = \sin^2 z ,$$ we obtain the algebraic form of Mathieu’s equation
$\zeta \left( 1 - \zeta \right) w'' + \frac{1}{2} \left( 1 - 2\zeta \right) w' + \frac{1}{4} \left( 1 - 2q (1-2\zeta ) \right) w =0 .$
This equation has regular singularities at 0 and 1; it also has an irregular singular point at infinity.
Example 1: The idea of building traps grew out of molecular-beam physics, mass spectrometry, and particle accelerator physics. It was shown that plane electric and magnetic multipole fields are able to focus particles in two dimensions acting on the magnetic or electric dipole moment of the particles. In the electric quadrupole field the potential is quadratic in the Cartesian coordinates. The Laplace equation imposes a two dimensional model
$\Phi = \frac{\Phi_0}{2r_0^2} \left( x^2 - z^2 \right) .$
The equation of motion are
\begin{align*} \ddot{x} + \frac{e}{m r_0^2} \left( U + V \,\cos \omega t \right) x &= 0 , \\ \ddot{z} - \frac{e}{m r_0^2} \left( U + V \,\cos \omega t \right) z &= 0 , \end{align*}
where U is a dc voltage and V is an rf voltage.    ■
Example 2: Ship motion can be described by the following equation
$\left( {\bf M} + {\bf A} \right) \ddot{\phi} + D(\ddot{\phi})\,\dot{\phi} + C_{res} (\phi , \eta ) = C{ext} (\eta , \dot{\eta}, \ddot{\eta}) ,$
■
Example 3: A simple model for the motion of a conduction electron in a crystal lattice consists of the time-independent Schrödinger equation with a potential that varies periodically in space. This potential describes the effect of the forces due to the ions in the crystal lattice and the other electrons in the crystal on the motion of the electron.

Let us consider the one-dimensional version of this problem. The wavefunction of the electron then satisfies the one-dimensional Schrodinger equation in which the potential V(x) is a periodic function. Suppose that the period is a, so that V(x+a)=V(x). Then, after making the change of variables

$$u(x) = \varphi (x/a) , q(x) = (2ma^2)/h V(x/a), \lambda = 2mE/h^2 ,$$

the normalized wavefunction u(x) and energy parameter λ satisfy

$$-u'' + q(x) u = \lambda * u$$

where $$q(x+1) =q(x) .$$ We consider this Sturm--Liouville problem on the real line with a general 1-periodic potential q(x), assumed continuous.
A specific example of such an equation is the Mathieu equation:

$$-u'' +2k \,\cos (2x) u = \lambda u ,$$

where k is a real paarmeter. Its solutions, called the Mathieu functions, have been studied extensively.

Next we plot Mathieu function:

Plot3D[MathieuC[3,-3/2,x] * MathieuS[5, -3/2,y], {x, -3, 3}, {y, -3, 3}]

1. Allievi, A. and Soudack, A., Ship stability via the Mathieu equation, International Journal of Control, 1990, Volume 51, Issue 1, pp. 139--167. https://doi.org/10.1080/00207179008934054
2. Butikov, Eugene I., Analytical expressions for stability regions in the Ince–Strutt diagram of Mathieu eqution, American Journal of Physics, 86, Isuue 4, 257 (2018); https://doi.org/10.1119/1.5021895
3. Daniel, D.J., Exact solutions of Mathieu’s equation, Progress of Theoretical and Experimental Physics, 2020, Volume 2020, Issue 4, April 2020, 043A01, https://doi.org/10.1093/ptep/ptaa024
4. Jovanoski, Z. and Robinson, G., Ship Stability and Parametric Rolling, Australasian Journal of Engineering Education, Volume 15, 2009 - Issue 2, pp. 43--50. https://doi.org/10.1080/22054952.2009.11464028
5. Mathieu, É.L., Memoir on vibrations of an elliptic membrane, 1868; Mémoire sur le mouvement vibratoire d'une membrane de forme elliptique, Journal de Mathématiques Pures et Appliquées 13, 137–203 (1868);
6. McLachlan, N.W., Theory of Application of Mathieu Functions, Dover, New York, 1964.
7. Paul, W., Electromagnetic traps for charged and neutral particles, Reviews of Modern Physics, 1990, Vol. 62, No. 3, pp. 531--540. http://dx.doi.org/10.1103/RevModPhys.62.531 )
8. Rand, R.H., Lecture Notes in Nonlinear Vibrations (new edition). Published on-line by The Internet-First University Press (Cornell's digital repository), 2012.
9. Ruby, L., Applications of the Mathieu equation, American Journal of Physics, 1995, Volume 64, Issue 1 https://doi.org/10.1119/1.18290