Shooting method for finding global eigenmodes

After using Fourier spectrum expansion and taking the inner product over $\theta$, Eq. (190) can be written as the following system of ordinary differential equations:

$$\frac{d}{d \psi} \left(\begin{array}{c} \overline{P}_1^{(1)} (\ps... ... \xi_{\psi}^{(1)} (\psi)\\ \vdots\\ \xi_{\psi}^{(L)} (\psi) \end{array}\right),$$ (251)

where $L$ is the total number of the poloidal harmonics included in the Fourier expansion, the matrix elements $A_{\alpha \beta}^{(i j)}$ are functions of $\psi$ and $\overline{\omega}^2$. Next, we specify the boundary condition for the system. Note that equations system (251) has $2 L$ first-order differential equations, for which we need to specify $2 L$ boundary conditions to make the solution unique. The geometry determines that the radial displacement at the magnetic axis must be zero, i.e.,

$$\xi_{\psi}^{(l)} (\psi = \psi_0) = 0, \ensuremath{\operatorname{For}} l = 1, 2, \ldots, L.$$ (252)

We consider only the modes that vanish at the plasma boundary, for which we have the following boundary conditions:

$$\xi_{\psi}^{(l)} (\psi = \psi_{\ensuremath{\operatorname{LCFS}}}) = 0, \ensuremath{\operatorname{For}} l = 1, 2, \ldots, L$$ (253)

Now Eqs. (252) and (253) provide $2 L$ boundary conditions, half of which are at the boundary $\psi = \psi_0$ and half are at the boundary $\psi = \psi_{\ensuremath{\operatorname{LCFS}}}$. Therefore equations system (251) along with the boundary conditions Eqs. (252) and (253) constitutes a standard two-points boundary problem[5]. Note, however, that we are solving a eigenvalue problem, for which there is an additional equation for $\overline{\omega}^2$:

$$\frac{d \overline{\omega}^2}{d \psi} = 0.$$ (254)

This increases the number of equations by one and so we need one additional boundary condition. Note that, by eliminating all $\overline{P}_1^{(l)}$, equations system (251) can be written as a system of second-order differential equations for $\xi_{\psi}^{(l)}$. Further note that the unknown functions $\xi_{\psi}^{(l)}$ satisfy homogeneous equations and homogeneous boundary conditions, which indicates that if $\xi_{\psi}^{(l)}$ with $l = 1, 2, \ldots L$ are solutions, then $c \xi_{\psi}^{(l)}$ are also solutions to the original equations, where $c$ is a constant. Therefore the value of the derivative of $d \xi_{\psi}^{(l)} / d \psi$ at the boundary have one degree of freedom. Due to this fact, one of the derivatives $d \xi_{\psi}^{(1)} / d \psi$, $d \xi_{\psi}^{(2)} / d \psi$,..., $d \xi_{\psi}^{(L)} / d \psi$ at $\psi_0$ can be set to be a nonzero value. For example, setting the value of $d \xi_{\psi}^{(1)} / d \psi$ at $\psi_0$ to be $0.5$ and making use of $\xi_{\psi}^{(l)} = 0$ at $\psi_0$, we obtain

$$A_{21}^{(11)} (\psi_0) \overline{P}_1^{(1)} (\psi_0) + A_{21}^{(1... ...\psi_0) + \ldots + A_{21}^{(1 L)} (\psi_0) \overline{P}_1^{(L)} (\psi_0) = 0.5,$$ (255)

which can be solved to give

$$\overline{P}_1^{(L)} (\psi_0) = \frac{1}{A_{21}^{(1 L)} (\psi_0)}... ...(1)} (\psi_0) - A_{21}^{(12)} (\psi_0) \overline{P}_1^{(2)} (\psi_0) - \ldots],$$ (256)

which provides the additional boundary condition we need. In the present version of my code, for convenience, I directly set the value of $\overline{P}_1^{(L)} (\psi_0)$ to a small value, instead of using Eq. (256). The following sketch map describes the function $\mathbf{F} (\mathbf{X})$ for which we need to find roots in the shooting process.

$$\mathbf{X}= \left(\begin{array}{c} \overline{P}_1^{(1)} (\psi_0)\... ...\ \xi_{\psi}^{(L)} (\psi_{\ensuremath{\operatorname{LCFS}}}) \end{array}\right)$$ (257)