The Homepage of Liron McLey

Snow on the way to Lake Louise, AB, Canada

Homework files     Useful links     Project


Background     The Equations     Data feed     Numerical method     Compiling and running     Validation of the code


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:


<< Previous     Next >>


Go back to the Computational Physics Class index page