A mathematical view of our approach

0.65![clip]diag_th_model/rad_shells.eps

In this subsection, we explain briefly how the gaseous model is solved numerically. We concentrate on aspects of the method not exposed in previous work. This description is therefore complementary to Sec. 2.2 of GS94. The algorithm used is a partially implicit Newton-Raphson-Henyey iterative scheme (HenyeyEtAl59, see also KW94, Sec. 11.2).

Putting aside the bounding conditions, the set of equations to be solved are Eqs. 2.22 to 2.25. In practice, however, the equations are rewritten using the logarithm of all positive quantities as dependant functions. As explained in GS94, this greatly improves energy conservation. Formally, one may write this system as follows

where the $x\ensuremath{^{(i)}}$ are the local quantities defining the state of the cluster, i.e.

with $\ensuremath{N_{\rm eq}}=7$ in the present application.

To be solved numerically, this set of coupled partial differential equations have to be discretised according to time and radius. Let us first consider time stepping. Let $\Delta t$ be the time step. Assume we know the solution $\underline{x}(t-\Delta t)$ at time $t-\Delta t$ and want to compute $\underline{x}(t)$. For the sake of numerical stability, a partially implicit scheme is used. We adopt the shorthand notations $x\ensuremath{^{(i)}}\equiv x\ensuremath{^{(i)}}(t)$ and $y\ensuremath{^{(i)}}\equiv x\ensuremath{^{(i)}}(t-\Delta t)$. Time derivation is replaced by finite differences, x t t^-1(x -y ). In the terms $f\ensuremath{^{(i)}}$, we replace the $x\ensuremath{^{(j)}}$ by $\tilde{x}\ensuremath{^{(j)}}$ which are values intermediate between $y\ensuremath{^{(j)}}$ and $x\ensuremath{^{(j)}}$, $\tilde{x}\ensuremath{^{(j)}}= \zeta x\ensuremath{^{(j)}}+ (1-\zeta)y\ensuremath{^{(j)}}$, with $\zeta= 0.55$ for stability purpose GS94.

Spatial discretisation is done by defining all quantities (at a given time) on a radial mesh, $\{r\ensuremath{_{1}}, r\ensuremath{_{2}},\dotso r\ensuremath{_{N_{\rm r}}}\}$ with $r\ensuremath{_{1}}=0$ and $r\ensuremath{_{\ensuremath{N_{\rm r}}}}=r_{\rm max}$. A staggered mesh is implemented. While values of $r$, $u$, $v_{\rm t}$, $v_{\rm r}$ and $M_{\rm r}$ are defined at the boundaries of the mesh cells, $\rho$, $p_{\rm t}$ and $p_{\rm r}$ are defined at the centre of each cell, see Fig.2.2. When the value of a ``boundary'' quantity is needed at the centre of a cell, or vice-versa, one resort to simple averaging, i.e. $\hat{b}\ensuremath{_{k}}= 0.5(b\ensuremath{_{k-1}}+b\ensuremath{_{k}})$, $\hat{c}\ensuremath{_{k}}= 0.5(c\ensuremath{_{k}}+c\ensuremath{_{k+1}})$, if $b$ and $c$ are is border- and centre-defined quantities, and $\hat{b}$, $\hat{c}$ their centre- and border-interpolations, respectively. For all runs presented here, $\ensuremath{N_{\rm r}}=300$ and the points $r\ensuremath{_{2}}$ to $r_{\rm max}$ are logarithmically equidistant with $r_{\rm max}=10^4\,{\rm pc}$ and $r\ensuremath{_{2}}\simeq 1.7\times 10^{-6}\,{\rm pc}$. Let us adopt the notation $x\ensuremath{^{(j)}}\ensuremath{_{k}}$ for the value of $x\ensuremath{^{(j)}}$ at position $r\ensuremath{_{k}}$ (or $\hat{r}\ensuremath{_{k}}$) and $\Delta r\ensuremath{_{k}}\equiv r\ensuremath{_{k}}-r\ensuremath{_{k-1}}$. Then, radial derivatives in the terms $f\ensuremath{^{(i)}}$ are approximated by finite differences, x r #<2692#>x -#<2698#>x r if the derivative has to be evaluated at a point where $x\ensuremath{_{k}}$ is defined (centre or border of a cell), or x r #<2710#>x -#<2716#>x r = x -x 2r otherwise. As an exception we use upstream differencing in $\partial{u}/\partial{r}$ for the second equation in set 2.22, i.e. the difference quotient is displaced by half a mesh point upstream to improve stability.

By making the substitutions for ${\partial x\ensuremath{^{(j)}}}/{\partial t}$ and ${\partial x\ensuremath{^{(j)}}}/{\partial r}$ in the set of differential equations 2.72, one obtains, at each mesh point $r\ensuremath{_{k}}$, a set of $\ensuremath{N_{\rm eq}}$ non-linear algebraic equations linking the new values to be determined, $\underline{x}\ensuremath{_{k-1}}$ and $\underline{x}\ensuremath{_{k}}$, to the ``old'' ones, $\underline{y}\ensuremath{_{k-1}}$ and $\underline{y}\ensuremath{_{k}}$, which are known,

Note that the structure of the equations is the same at all mesh points, except $k=1$ and $k=\ensuremath{N_{\rm r}}$. In particular, terms $k-1$ do not appear in ${\cal F}\ensuremath{^{(i)}}\ensuremath{_{1}}$. Also, one has to keep in mind that only the $\underline{x}\ensuremath{_{k-1}}$ and $\underline{x}\ensuremath{_{k}}$ are unknown; the $\underline{y}\ensuremath{_{k-1}}$ and $\underline{y}\ensuremath{_{k}}$ play the role of fixed parameters in these equation (as do the $\Delta r\ensuremath{_{k}}$). If one defines a $(\ensuremath{N_{\rm eq}}\times \ensuremath{N_{\rm r}})$-dimension vector ${{\cal X}^*}$ whose component $\ensuremath{N_{\rm eq}}(k-1)+i$ is $x\ensuremath{^{(i)}}\ensuremath{_{k}}$, one can write the system of $\ensuremath{N_{\rm eq}}\times \ensuremath{N_{\rm r}}$ equations as ${{\cal F}^*}({{\cal X}^*})=0$, i.e. F^*(X^*) (

) = (

) with X^* (

).

The system is solved iteratively using Newton-Raphson scheme. If ${{\cal X}^*}\ensuremath{_{[ m]}}$ is the approximation to the solution of Eq. 2.78 after iteration $m$, with ${{\cal F}^*}\ensuremath{_{[ m]}}\equiv {{\cal F}^*}({{\cal X}^*}\ensuremath{_{[ m]}})\ne 0$, the solution is refined through the relation

X^* = X^* - ( )^-1 F^*

where $(\ensuremath{{\partial{{\cal F}^*}}/{\partial {{\cal X}^*}}})^{-1}$ is the inverse of the matrix of derivatives. The latter, of dimension $(\ensuremath{N_{\rm eq}}\, \ensuremath{N_{\rm r}})\times(\ensuremath{N_{\rm eq}}\, \ensuremath{N_{\rm r}})$, has the following structure

F^*X^* =

. In this diagram, each square is a $\ensuremath{N_{\rm eq}}\times \ensuremath{N_{\rm eq}}$ sub-matrix. For $2\le k \le \ensuremath{N_{\rm r}}-1$, lines ${\ensuremath{N_{\rm eq}}}k-6$ to ${\ensuremath{N_{\rm eq}}}k$ of $\ensuremath{{\partial{{\cal F}^*}}/{\partial {{\cal X}^*}}}$ are composed of a group of 3 such ${\ensuremath{N_{\rm eq}}}\times {\ensuremath{N_{\rm eq}}}$ matrices, $\ensuremath{{\square_{-}}}\ensuremath{_{k}}, \ensuremath{{\blacksquare}}\ensuremath{_{k}}, \ensuremath{{\square_{+}}}\ensuremath{_{k}}$ that span columns ${\ensuremath{N_{\rm eq}}}k-13$ to ${\ensuremath{N_{\rm eq}}}k+{\ensuremath{N_{\rm eq}}}$, while the rest is composed of zeros,

We can see this more explicitely,

More schematically, to have an overview of the structure,

The Heyney method is a way to take advantage of the special structure of matrix $\ensuremath{{\partial{{\cal F}^*}}/{\partial {{\cal X}^*}}}$ to solve system (2.79) efficiently, with a number of operation scaling like rather than as would be the case if one uses a general-purpose matrix inversion scheme. Setting and , Eq.(2.79) is equivalent to ( ) W^* = B^* with the unknown vector. We further decompose vectors and into -dimensional sub-vectors, each one representing the values at a given mesh points, W^* =

. Then, the system (2.83) can be written as a set of coupled -dimensional vector equations,

The algorithm operates in two passes. First, going from $k=1$ to , one defines recursively a sequence of -vectors and -matrices through

is not defined. In the second pass, the values of the unknown are computed, climbing back from $k=\ensuremath{N_{\rm r}}$ to , with

Note that, with this algorithm, only matrices have to be inverted. We use Gauss elimination for this purpose because this venerable technique proves to be robust enough to properly deal with the kind of badly conditioned matrices that often appear in this application.

The initial model for the Newton-Raphson algorithm is given by the structure of the cluster at the previous time, One iterates until the following convergence criteria are met. Let us set . Then, the condition for logarithmic quantities is _i=1... 1 _k=1... (x )^2 < _1, with . For velocities ($u$, , ), one checks _i=1... 1 _k=1... ( x x +_1 w )^2 < _2, with and . Generally, two iterations are sufficient to reach convergence.

0.95