B \equiv q/x^{3}\;\;\; q=m/M

V \equiv -\frac{1}{\Gamma_{1}}\frac{dln\, p}{dln\, r}=\frac{Gm\rho}{\Gamma_{1}pr}

\Gamma_{1} \equiv \Gamma_{1}

A \equiv \frac{1}{\Gamma_{1}}\frac{dln\, p}{dln\, r}-\frac{dln\,\rho}{dln\, r}

\sigma^{2}=\frac{R^{3}}{GM}\omega^{2}

y_{2} = \frac{l\left(l+1\right)}{R}\xi_{h}

Turning the Cowling equations dimensionless using these variables yields,

\frac{dy_{2}}{dx}=\frac{1}{x}\left(l(l+1)-\eta A\right)y_{1}+\frac{1}{x}\left(A-1\right)y_{2}

where \eta=l(l+1)B/\sigma^{2}. For the radial case, y_{2} is replaced by,

and so the equations become,

\frac{dy_{2}}{dx}=\left(1-\frac{BA}{\sigma^{2}}\right)y_{1}+\frac{1}{x}Ay_{2}

The code uses the second order centered differences method to solve the differential equations. Hence the equations above are replaced by the difference equations,

where a_{ij}^{n}=a_{ij}\left(x^{n}\right), y_{j}^{n}=y_{j}\left(x^{n}\right) and x^{n} are the points that construct the mesh of the star. Hence, the solution at x^{n+1} is determined by that at x^{n}. Note that the equation is symmetric under the exchange y_{i}^{n+1}\leftrightarrow y_{i}^{n} and x^{n+1}\leftrightarrow x^{n}.

The code finds the mode order according to the definition made by Scuflaire (1974). The order is defined as follows,

where x_{z1} are the zeros in y_{1} (excluding the center) and sign is the sign function. The value of n_{0} depends on the way the solution behaves near the inner boundary : if y_{1} and y_{2} have a similar sign then n_{0}=0, otherwise n_{0}=1. For the cases intended to be run by this code (a complete stellar model that includes the center together with the Cowling approximation), it follows that n_{0}=1 for radial oscillations and n_{0}=0 for non-radial oscillation. Hence, this defines n>0 for p modes, n=0 for f modes and n<0 for g modes.

For the second order system, the differential equations have a unique (up to normalization) solution satisfying the inner boundary conditions, y_{i}^{(i)}
, and a unique solution satisfying the outer boundary condition, y_{i}^{(o)}.
These can be obtained by numerical integration. The final solutions are then y_{i}=C^{(i)}y_{i}^{(i)}=C^{(o)}y_{i}^{(o)}.
The eigenvalue is obtained by demanding that the solutions agree at a suitable matching point, x=x_{f}.
Hence,

C^{(i)}y_{2}^{(i)}(x_{f})=C^{(o)}y_{2}^{(o)}(x_{f})

A non trivial solution \left(C^{(i)},C^{(o)}\right) for these equations exists only when the determinant vanishes,

The basic function of the code is the following:

- Read the data file

- Select parameters to solve for - select a value of l and a value of n
for which to find the eigenfrequency \omega
and the eigenfunctions \xi_{i}.

- Convert everything into its dimensionless form

- Use the shooting method - numerically integrate the equations for some \sigma^{2}, once for the interior boundary condition and once for
the surface boundary condition. Then evaluate whether \sigma^{2} is an eigenvalue or not - does \Delta=0?

- If \sigma^{2} is an eigenvalue, use the eigenfunction to determine its order.

- If the order is larger or smaller than n, keep looking in the appropriate direction.

- If it is equal to n, print an output file containing the data.

Go back to the Computational Physics Class index page