# The Homepage of Liron McLey

## Numerical method

The first thing I did was to get rid of dimensions. For each mesh point in the model I defined the following dimensionless quantities,

x \equiv r/R

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}

The dimensionless variable corresponding to the frequency is,

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

### The dimensionless equations

For the non-radial case (l \neq 0) I set the dimensionless variables as,

y_{1} = \frac{\xi_{r}}{R}

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

Turning the Cowling equations dimensionless using these variables yields,

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

\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,

y_{2}=\frac{p'}{\omega^{2}R^{2}\rho}

and so the equations become,

\frac{dy_{1}}{dx}=\frac{1}{x}\left(V_{g}-2\right)y_{1}+\frac{x}{q}\sigma^{2}Vy_{2}

\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,

\frac{y_{i}^{n+1}-y_{i}^{n}}{x^{n+1}-x^{n}}=\frac{1}{2}\Sigma_{j=1}^{2}\left(a_{ij}^{n}y_{j}^{n}+a_{ij}^{n+1}y_{j}^{n+1}\right)

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

### Mode order convention

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

n=-\Sigma_{x_{z1}}sign\left(y_{2}\frac{dy_{1}}{dx}\right)+n_{0}

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.

### The shooting method

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_{1}^{(i)}(x_{f})=C^{(o)}y_{1}^{(o)}(x_{f})

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,

\Delta\equiv y_{1}^{(i)}(x_{f})y_{2}^{(o)}(x_{f})-y_{2}^{(i)}(x_{f})y_{1}^{(o)}(x_{f})=0

### Summary

The basic function of the code is the following:

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

<< Previous     Next >>

Go back to the Computational Physics Class index page