next up previous
Next: C Progam for calculating Up: Mathieu's Equation, solution, and Previous: Hill's Method solution

Sträng's recursion formula for $ \Delta (0)$

First we note that by the symmetry of $ \Delta (0)$ , $ \gamma_{-n} = \gamma_{n}$ . Following Sträng, we define

$\displaystyle A_i = \left(\begin{array}{ccccccccc}
1&\gamma_{2i}&0&&&&&&\\
\ga...
...&&&\gamma_{2(i-1)}&1&\gamma_{2(i-1)}\\
&&&&&&0&\gamma_{2i}&1\end{array}\right)$ (2.10)

We have $ \Delta_i=det(A_i) \ni \Delta(0)=\lim_{i\rightarrow \infty}\Delta_i.$ We can decompose $ A_i$ in terms of $ A_{i-1}$ ,

$\displaystyle A_i = \left(\begin{array}{ccccc}
1&\gamma_{2i}&&&\\
\gamma_{2(i-...
...\
&.&A_{i-1}&.&\\
&.&.&.&\gamma_{2(i-1)}\\
&&&\gamma_2i&1
\end{array}\right)$ (2.11)

A Laplace decomposition yields

$\displaystyle det(A_i)=\left\vert\begin{array}{cccc}
.&.&.&\\
.&A_{i-1}&.&\\
...
...&rA_{i-1}&.&\\
&.&.&.&\gamma_{2(i-1)}\\
&&&\gamma_2i&1
\end{array}\right\vert$ (2.12)

Here $ rA_{i-1}$ represents $ A_{i-1}$ with its left most column chopped off. Again, following Sträng we define $ lA$ as the matrix A with its rightmost column removed, $ uA$ the matrix A with its lowest row removed, $ lA$ the matrix A with its upper most row removed. Ultimately, $ uldr(A_{i-1})=A_{i-2}$ , and given the symmetry involved, det$ (rd(A_{i-1}))=$det$ (ul(A_{i-1}))$ . Following this procedure we find

$\displaystyle \Delta_i=\Delta_{i-1}-2\gamma_{2i} \gamma_{2(i-1)}$det$\displaystyle (rd(A_{2(i-1)}))+(\gamma_{2i} \gamma_{2(i-1)})^2 \Delta_{i-2}$ (2.13)

We also note, similarly using Lapalcian decomposition,

$\displaystyle \Omega_i=$det$\displaystyle (ul(A_i))=$det$\displaystyle (rd(A_i))  \rightarrow  \Omega_i=det(A_{i-1})-\gamma_{2i} \gamma_{2(i-1)} \Omega_{i-2}$

so that

$\displaystyle \frac{\Delta_{i-1}-\Omega_i}{\gamma_{2i} \gamma_{2(i-1)}}=\Omega_{i-1}=$det$\displaystyle (rd(A_{i-1}))$

and

$\displaystyle \Delta_i = \Delta_{i-1} + 2(\Omega_i-\Delta_{i-1})+(\gamma_{2i} \gamma_{2(i-1)})^2\Delta_{i-2}$

$\displaystyle \frac{\Delta_i + \Delta_{i-1} -(\gamma_{2i} \gamma_{2(i-1)})^2\De...
... \Delta_{i-2} -(\gamma_{2(i-1)} \gamma_{i-2})^2\Delta_{i-3})}{2} = \Omega_{i-1}$

Plugging this into Equation 2.13,

$\displaystyle \Delta_i = (1-\gamma_{2i} \gamma_{2(i-1)})\Delta_{i-1}+\left( (\g...
...i-2} +\gamma_{2i} \gamma_{2(i-1)}(\gamma_{2(i-1)}\gamma_{2(i-2)})^2\Delta_{i-3}$

Define $ \alpha_{2i} = \gamma_{2i} \gamma_{2(i-1)}$ and $ 1-\alpha_{2i}=\beta_{2i}$ and find,

$\displaystyle \Delta_i = \beta_{2i}\Delta_{i-1}-\alpha_{2i} \beta_{2i} \Delta_{i-2} + \alpha_{2i}\alpha_{2(i-1)}^2\Delta_{i-3}$ (2.14)

We can recursively solve for $ \Delta(0)=\lim_{i\rightarrow \infty}\Delta_i$ to as much accuracy as nessessary, though the program presented below found convergence to a fair tolerance quite quickly. We first must ``seed'' the recurssion with the first three $ \Delta_i$ . This can be done by hand, though we have deferred to the kindess of our computer algebraic program Maple instead.

Maple finds,

with(linalg):
C:=matrix([[1,e6,0,0,0,0,0],[e4,1,e4,0,0,0,0],[0,e2,1,e2,0,0,0],
   [0,0,e0,1,e0,0,0],[0,0,0,e2,1,e2,0],[0,0,0,0,e4,1,e4],[0,0,0,0,0,e6,1]]):
dc:=det(C);
dc := -2*e2^2*e0*e4^2*e6+e2^2*e4^2-2*e4^2*e2*e0*e6^2+2*e2*e4^2*e6
       +e4^2*e6^2+2*e2^2*e0*e4+4*e2*e0*e6*e4-2*e2*e4-2*e6*e4-2*e2*e0+1


A:= matrix([[1,e4,0,0,0],[e2,1,e2,0,0],[0,e0,1,e0,0],
                                               [0,0,e2,1,e2],[0,0,0,e4,1]]):
da:=det(A);
da := 1-2*e2*e4-2*e2*e0+2*e2^2*e0*e4+e2^2*e4^2
 

B:=matrix([[1,e2,0],[e0,1,e0],[0,e2,1]]):
db:=det(B);
db := 1-2*e2*e0

Our program seeks to find all stable values of $ \mu $ , i.e. those that satisfy Equation 2.9 as real values (i.e. all iso-$ \mu $ for which $ \mu $ is exclusively imaginary.

Our code finds all such iso-$ \mu $ by looping through the a and q axis. If our $ \mu $ formula returns ``nan'' which is the C language's way of saying not a real number, then we set the value of $ \mu $ to zero, though of course it is only the imaginary part of $ \mu $ which is actually zero. We perform a contour plot on our data output and find the elegant avian like image of the stability region of Mathieu's equation (Figure 2.1).

Figure 2.1: Stability diagram for Mathieu's general equation
\includegraphics[width=7.5cm]{stable5.ps}

For the quadropole field, the rf linear Paul trap, we have the following stability regime (Figure 2.2). The original stability diagram is simply reflected about the x-axis as $ a \rightarrow -a$ between the two.

Figure 2.2: Stability diagram for linear rf Paul trap; stable regions are those in which the two stability diagrams intersect.
\includegraphics[width=7.5cm]{stable8.ps}

Figure 2.3: Stability diagram for linear rf Paul trap, closer view of lowest region of mutal stability (intersection).
\includegraphics[width=7.5cm]{rf2.ps}

For the chamber rf Paul trap, we recall that for the z direction we must allow for $ (a,q) \rightarrow (2a,2q)$ (Equation 1.3), and so the lowest region of stability (and the largest region at that) has a slightly different look (Figure 2.4).

Figure 2.4: Stability diagram for rf Paul trap, lowest region of stability (intersection of isobars).
\includegraphics[width=10cm]{rf4.ps}


next up previous
Next: C Progam for calculating Up: Mathieu's Equation, solution, and Previous: Hill's Method solution
tim jones 2008-07-07