next up previous
Next: Formal solution Up: Method Previous: Radiative transfer equation

The operator splitting method

The mean intensity $J$ is obtained from the source function $S$ by a formal solution of the RTE which is symbolically written using the $\Lambda $-operator $\Lambda $ as

\begin{displaymath}
J = \Lambda S.
\end{displaymath} (2)

The source function is given by $S=(1-\epsilon )J + \epsilon B$, where $\epsilon $ denotes the thermal coupling parameter and $B$ is Planck's function.

The $\Lambda $-iteration method, i.e. to solve Eq. 2 by a fixed-point iteration scheme of the form

$\displaystyle { J_{\rm new}}= \Lambda {S_{\rm old}}, \quad
{S_{\rm new}}= (1-\epsilon ){ J_{\rm new}}+ \epsilon B ,$     (3)

fails in the case of large optical depths and small $\epsilon $. Here, ${S_{\rm old}}$ is the current estimate for the source function $S$ and ${S_{\rm new}}$ is new, improved, extimate of $S$ for the next iteration. The failure to converge of the $\Lambda $-iteration is caused by the fact that the largest eigenvalue of the amplification matrix is approximately (Mihalas et al., 1975) $\lambda _{\rm max} \approx (1-\epsilon )(1-T^{-1})$, where $T$ is the optical thickness of the medium. For small $\epsilon $ and large $T$, this is very close to unity and, therefore, the convergence rate of the $\Lambda $-iteration is very poor. A physical description of this effect can be found in Mihalas (1980).

The idea of the ALI or operator splitting (OS) method is to reduce the eigenvalues of the amplification matrix in the iteration scheme (Cannon, 1973) by introducing an approximate $\Lambda $-operator (ALO) $\ifmmode{\Lambda^*}\else\hbox{$\Lambda^*$}\fi $ and to split $\Lambda $ according to

\begin{displaymath}
\Lambda = \ifmmode{\Lambda^*}\else\hbox{$\Lambda^*$}\fi +(\Lambda -\ifmmode{\Lambda^*}\else\hbox{$\Lambda^*$}\fi )
\end{displaymath} (4)

and rewrite Eq. 3 as
\begin{displaymath}
{ J_{\rm new}}= \ifmmode{\Lambda^*}\else\hbox{$\Lambda^*$}\...
...\ifmmode{\Lambda^*}\else\hbox{$\Lambda^*$}\fi ){S_{\rm old}}.
\end{displaymath} (5)

This relation can be written as Hamann (1987)
\begin{displaymath}
\left[1-\ifmmode{\Lambda^*}\else\hbox{$\Lambda^*$}\fi (1-\e...
...bda^*}\else\hbox{$\Lambda^*$}\fi (1-\epsilon ){ J_{\rm old}},
\end{displaymath} (6)

where ${ J_{\rm fs}}=\Lambda {S_{\rm old}}$ and ${ J_{\rm old}}$ is the current estimate of the mean intensity $J$. Equation 6 is solved to get the new values of ${ J}$ which is then used to compute the new source function for the next iteration cycle.

Mathematically, the OS method belongs to the same family of iterative methods as the Jacobi or the Gauss-Seidel methods (Golub & Van Loan, 1989). These methods have the general form

\begin{displaymath}
M x^{k+1} = Nx^{k} + b
\end{displaymath} (7)

for the iterative solution of a linear system $Ax=b$ where the system matrix $A$ is split according to $A=M-N$. In the case of the OS method we have $M=1-\ifmmode{\Lambda^*}\else\hbox{$\Lambda^*$}\fi (1-\epsilon )$ and, accordingly, $N=(\Lambda -\ifmmode{\Lambda^*}\else\hbox{$\Lambda^*$}\fi )(1-\epsilon )$ for the system matrix $A=1-\Lambda (1-\epsilon )$. The convergence of the iterations depends on the spectral radius, $\rho(G)$, of the iteration matrix $G=M^{-1}N$. For convergence the condition $\rho(G)<1$ must be fulfilled, this puts a restriction on the choice of $\ifmmode{\Lambda^*}\else\hbox{$\Lambda^*$}\fi $. In general, the iterations will converge faster for a smaller spectral radius. To achieve a significant improvement compared to the $\Lambda $-iteration, the operator $\ifmmode{\Lambda^*}\else\hbox{$\Lambda^*$}\fi $ is constructed so that the eigenvalues of the iteration matrix $G$ are much smaller than unity, resulting in swift convergence. Using parts of the exact $\Lambda $ matrix (e.g., its diagonal or a tri-diagonal form) will optimally reduce the eigenvalues of the $G$. The calculation and the structure of $\ifmmode{\Lambda^*}\else\hbox{$\Lambda^*$}\fi $ should be simple in order to make the construction of the linear system in Eq. 6 fast. For example, the choice $\ifmmode{\Lambda^*}\else\hbox{$\Lambda^*$}\fi =\Lambda $ is best in view of the convergence rate (it is equivalent to a direct solution by matrix inversion) but the explicit construction of $\Lambda $ is more time consuming than the construction of a simpler $\ifmmode{\Lambda^*}\else\hbox{$\Lambda^*$}\fi $. The solution of the system Eq. 6 in terms of linear algebra, using modern linear algebra packages such as, e.g., LAPACK (Anderson et al., 1992), is so fast that its CPU time can be neglected for the small number of variables encountered in 1D problems (typically the number of discrete shells is about 50). However, for 2D or 3D problems the size of $\Lambda $ gets very large due to the much larger number of grid points as compared to the 1D case. Matrix inversions, which are necessary to solve Eq. 6 directly, therefore become extremely time consuming. This makes the direct solution of Eq. 6 more CPU intensive even for $\ifmmode{\Lambda^*}\else\hbox{$\Lambda^*$}\fi $'s of moderate bandwidth, except for the trivial case of a diagonal $\ifmmode{\Lambda^*}\else\hbox{$\Lambda^*$}\fi $. Different methods like modified conjugate gradient methods (Turek, 1993) may be effective for these 2D or 3D problems.

The CPU time required for the solution of the RTE using the OS method depends on several factors: (a) the time required for a formal solution and the computation of ${ J_{\rm fs}}$, (b) the time needed to construct $\ifmmode{\Lambda^*}\else\hbox{$\Lambda^*$}\fi $, (c) the time required for the solution of Eq. 6, and (d) the number of iterations required for convergence to the prescribed accuracy. Points (a), (b) and (c) depend mostly on the number of spatial points, and can be assumed to be fixed for any given configuration. However, the number of iterations required to convergence depends strongly on the bandwidth of $\ifmmode{\Lambda^*}\else\hbox{$\Lambda^*$}\fi $. This indicates, that there is an optimum bandwidth of the $\ifmmode{\Lambda^*}\else\hbox{$\Lambda^*$}\fi $-operator which will result in the shortest possible CPU time needed for the solution of the RTE, see Hauschildt et al. (1994).


next up previous
Next: Formal solution Up: Method Previous: Radiative transfer equation
Peter Hauschildt 2006-01-09