J. Quant. Spectrosc. Radiat. Transfer. Vol. 11, pp. 589595. Pergamon Press 1971. Printed in Great Britain
A MODIFIED FEAUTRIER METHOD G. B. RYBICKI Smithsonian Astrophysical Observatory, Cambridge, Mass. 02138 AbstraetA modificationof the Feautrier method of solvingtransfer problems is presented that allows advantage to be taken of such simplificationsas isotropic scattering and complete frequencyredistribution. The method has timing and storage requirements comparable to integral equation methods, but retains the ease of formulation of differential methods. The case of a twolevelatom is treated in detail, and an appropriate generalization is sketched for the treatment of multilevelcases by means of a complete linearization technique. 1. INTRODUCTION ONE OF the most successful general methods for the numerical solution of transfer problems is that devised by FEAUTRIER31~ This method has the advantages of ease of formulation and of programming, and it has excellent stability properties that set it apart from other methods using differential equations. However, the very generality of the Feautrier method often works to its disadvantage, since certain c o m m o n simplifications, such as isotropic scattering and complete frequency redistribution, do not lead to any simplifications in the method. This accounts in part for the popularity of numerical methods using integral equations, ~25) which do take advantage of the existence of angle and frequencyindependent source functions. An interesting comparison between the integral equation method and the Feautrier method may be made for the case of line transfer in a twolevel a t o m with isotropic scattering and complete frequency redistribution. Let N be the n u m b e r of depth points, and F be the number of frequencyangle points chosen in the calculation. Then the asymptotic* timing estimate for the Feautrier method is C N F 3, while for the integral equation method it is C ' N 2 F + C"N 3. The asymptotic storage requirement is N F 2 for the Feautrier method and about 2N 2 for the integral equation method. Since F is often comparable to N in value (and can be quite a bit larger for problems involving velocity fields) the integral equation method is clearly advantageous for m a n y numerical computations. On the other hand the integral equation method is much harder to formulate and to program. The difficulties of formulation become even greater as the physical complexity of the problem increases. For example, the Feautrier method has been applied to n o n  L T E model atmospheres using the complete linearization technique of AuER and MIHALAS.(6) Although the integral equation method has not yet been completely formulated for such problems, it now appears * By this is meant the dominant terms for large N and F. The constants C, C', C",... are all roughly of the same magnitude, and may stand for different numerical values in different parts of the paper. 589
590
G.B. RYBICKI
that a comparable formulation will be quite complicated. It would seem that one must make a choice between the ease of formulation of the Feautrier method and the timing and storage advantages of the integral equation method. It is the purpose of this paper to demonstrate a simple modification of the Feautrier method that changes its timing and storage requirements to match those of the integral equation method. This modification is purely at the computational level and does not affect the ease of formulation of the method. For example, the differencing of the transfer equation and the treatment of boundary conditions is done precisely as in the usual method. "'7~ Therefore, there is a choice at the computational level of solving the equations in the usual fashion or by using the modification, whichever is most favourable for the particular problem at hand. In Section 2 the modified Feautrier method is described for the case of a twolevel atom assuming isotropic scattering and complete frequency redistribution. To facilitate the explanation, an intermediate version of the Feautrier method, suggested by the work of AUER and MIHALAS,t6) is introduced, in which the integrated mean intensity J is used as an independent variable. An "equation of constraint" then relates ? to the intensity variables. The modified Feautrier method is then derived from this intermediate version by changing the overall grouping of the variables from the usual one of depth blocks to one of frequencyangle blocks. Gaussian elimination by blocks then leads to a method with the timing and storage estimates of the integral equation method. In Section 3 it is shown how the modified Feautrier method can be generalized and applied to other transfer problems that would ordinarily admit a treatment in terms of one or more (coupled) integral equations. In particular it is indicated how the solution of a line transfer problem in a multilevel atom using a complete linearization technique may be treated by means of the modified Feautrier method. 2. TWOLEVEL ATOM The modified Feautrier method is based on the usual FEAUTRIER"~ choice of variables, which are here denoted by j(T, #, v) = ½CI(~, ~, ~) + I(~,  ~, v)] (1) h(~, ~, v) = ½EI(~,~, x )  I(~,  ~, v)].
(2)
I(z, #, v) is the specific intensity at angle (cosine) #, frequency v, and at standard optical depth z. The transfer equation may then be presented as a second order equation forj. ~
N
= js
(3)
Here X is the total opacity, normalized by the opacity used to define the zscale. For complete frequency redistribution and isotropic scattering the total source function S depends on the radiation field linearly through a single function of depth, the integrated mean intensity 1
](~) = f d#' f dv'c~(z, v')j(r, #', v'), 0
where ~ is the normalized line profile function.
(4)
A modified Feautrier method
591
These equations m a y be approximately reduced to finite terms by the introduction of various discretizations. Let z ~, 0¢ = 1. . . . N, be discrete optical depths and let (/~, v~), i   1. . . . . F, be discrete frequencyangle points. The N F values of j at these discrete points are denoted j~, 0(  1. . . . N, i = 1. . . . , F. Quadrature formulas are used to reduce the integrals in equation (4) to a sum so that the values of J at r ~ are F
Y'=
~ og,j7
(5)
i=1
The second order operator m a y be differenced and the boundary conditions introduced in the usual way. (1'7) There are three ways of proceeding from this stage, each giving rise to a different variation of the Feautrier method•
A. Usual Feautrier method The variables j~ (no subscript) will denote an Fdimensional vector with components (j~,j~ . . . . . j~)r. This gives the frequencyangle dependence o f j at a single depth point z,. Expressing the source function in terms offl via the integrated mean intensity, one obtains the overall matrix of the problem in the block tridiagonal form n I m C 1 _ A 2
J'l
0
A 3
B3C 3
0
ILll
i21
B 2 _ C 2
L2 I
i31
 AN
BN
=
IL31
(6)
L,/'J
_j J
Here A ~ and C ~ are F x F diagonal matrices, while the B ~ are full matrices• The n o n h o m o geneous terms in the equations appear in the vectors L ~. The problem is reduced to the block uppertriangular form ID
a
ID 2 ID 3
0
0
jl I
v1 I
j2 I
v 2 I!
j3 I =
v31
jNj
(7)
I'1
ivNj
by Gaussian elimination, from which the j'~ are easily determined by back substitution. The matrix I is the F x F unit matrix. The time consuming part of this solution is the computation of the N matrices D ~, each of which requires the inversion of an F x F matrix ; thus the asymptotic timing estimate of the method is C N F 3. The principal storage requirement is for storing these same N matrices for the back substitution, so that the asymptotic storage estimate is N F 2. Under the assumptions of isotropic scattering and complete frequency redistribution made here the matrices B ~ are actually the sum of a diagonal part and a matrix of rank one.
592
G.B.
RYBICKI
However, this special form for B~ does not seem to lead to any simplifications in the solution, since full matrices must still be inverted. The timing and storage estimates are the same as when general redistribution is assumed• B. Feautrier method with constraint
Instead of expressing the source function in terms of the variables j~, one may introduce 3~ as additional independent variables and regard equation (5) as an "equation of constraint," which defines a relationship between J~ andj~. This procedure is suggested by the work of AUER and MIHALAS.t6~ The N vectors ~b~, ~t = 1. . . . . N, of dimension ( F + l) are introduced, having the components (j~,j~ . . . . . j~,j~)T. These are simply the vectors j~ with the additional component J~. Then the overall matrix of the problem may again be expressed in block tridiagonal form, B 1 _
_A 2 A
C1
L1
0
B2_C
2
3
B3_C
L 2
3
t~3 l
0
(8)
L 3
LN
 AN BN .

The quantities A", B~, C ", and L ~ are of dimensionality (F+ 1) and differ from those introduced before. A" and C" are still diagonal, but now B~ has the form* x x
x
x
0
x
(9)
0 x
x
x
That is, it is diagonal except for the last row and column, which are full• The first Felements of the last column expresses the linear dependence of S on J, while the last row expresses the equation of constraint (5). This problem has the same form as before, so Gaussian elimination may be used to solve for the ~ . The special form of B ~ still does not yield any simplifications, so that the asymptotic estimates are the same as for the usual Feautrier method. C. Modified Feautrier
The particular modification that is the subject of this investigation is obtained from the preceding variation of the Feautrier method by using a different grouping of the variables. The F vectors Ji (no superscript), are defined to have the components (j~, j~ . . . . . if)r• These give the depth dependence of the particular frequencyangle group (#i, vl). Similarly the vector J with components (j1, 32 . . . . . jN)r is introduced. Then the overall matrix of the * The x's s t a n d f o r the n o n  z e r o elements.
A modified Feautrier method
593
problem may be written
7T1
0
T~.
J1
UI
T~
U2
J2
U3
I Ja
K,1 K2I Kal
(1o)
I!  •
v,
0
TF Ur
v~ v~
V~ E
ij
K,~I PI
The N × N matrices V~, U~, and E are diagonal, while each Ti has the tridiagonal form X
X
O
X
X
X
X
X
X
(11)
X
X
This different ordering might be described as a change in the overall form from depth blocks to frequencyangle blocks. The numerical values of the elements in equations (8) and (10) are unchanged, but they are to be found in different places. One notes that the "inner" and "outer" structure of the overall matrices has been interchanged, so that the block tridiagonal form of equation (8) is reflected in the tridiagonal form of the T~, while the special form of the B ~ as given in (9) now appears in the overall block form of equation (10). Equation (10) may be reduced by Gaussian elimination to the block uppertriangular form 7", K,] Ji u1 U2 K2I J2 T~ K a} U3 Ja
0
(12)
0 0
0
Tr Ue
jr
KFI
0
d
_QI
w
This is accomplished by inversion of each Ti, multiplying the ith row by this inverse and then by V~,and finally subtracting this row from the (F + 1)st row, which eliminates V~.This continues for i = 1, 2. . . . . F, until the last row reduces to the simple equation WJ = Q, (13) where F
W= E
~ V/Ti1Ui,
(14)
i=l F
Q = p
~ ViT~IKi . i=1
(15)
594
G . B . RYBICKI
Equation (13) may be solved for J by inverting IV. If needed, the vectors, Jl may be found by back substitution at a negligible cost in time, namely C ' N F ; but ordinarily the vector J will suffice. There is no need to store the overall matrix of the problem, since W is the only one needed. The quantities Wand Q are accumulated for each i in accordance with equations (14) and (15)A block of storage of size N 2 is also needed to invert each T~. The total storage required is therefore about 2N 2. The time required to invert each T~ is proportional to N 2, since it is tridiagonal. The multiplication by the diagonal matrices U i and V~are also order N 2 calculations, so that the time required to generate Wand Q is asymptotically proportional to N2F. The matrix W has no simple form and its inversion requires a time proportional to N 3. Thus the overall timing of the method is C ' N 2 F + C " N 3. These estimates are seen to be in accord with those of the integral equation method. It is interesting that each step in this modified Feautrier method corresponds to some equivalent step in the integral equation method. The inversion of a tridiagonal matrix T~ corresponds to the formal solution of the transfer equation for the ith frequencyangle group; the accumulation of terms in equation (14) corresponds to expressing the integrated mean intensity as a linear combination of intensities ; finally, the inversion of Wcorresponds to the solution of the linear system replacing the integral equation. This method was applied to a constant property, semiinfinite medium where e, the ratio of coUisional to total deexcitation rate, was taken to be 10 4. About 20 frequency angle points were taken, and about 30 depth points, which were spaced logarithmically, except for z  0. The results agreed closely with the standard ones3 2) No stability difficulties were encountered, so that it is likely that the stability of the usual Feautrier method is maintained. Similar conclusions were found by AUER (private communication), who applied the method to a problem with macroscopic velocity fields. 3. G E N E R A L I Z A T I O N S
The modification of the Feautrier method described in the preceding section can be applied to any case where the source function depends linearly on a single function of depth, sayf(z), that may in turn be expressed as a linear superposition of intensity components. The functionfthen plays the same role as the function J, and the method proceeds as before. A slightly more complicated generalization is possible when the source function depends linearly not on a single function but on L such functions fk(~), k = 1. . . . . L. It is just such a case that could be reduced to a set of L coupled integral equations by the techniques of the integral equation method. An example of such a problem would be anisotropic scattering with a phase function containing only Lterms in a Legendre polynomial expansion. L equations of constraint are required to define the fk in terms of the intensities. The matrices U i in equation (10) now are N x (NL) matrices, while the Vi becomes (NL) x N matrices. The form of these matrices is still quite simple, being composed of L diagonal matrices of dimensions N × N. The matrices E and W are (NL) × (NL) in size. Gaussian elimination now requires a time C ' N 2 F L + cIV(NL) 3 and storage (NL)2 + N 2. This is to be compared with C N F 3 time and N F 2 storage of the usual method. The usual method may be preferable when L is large. Multilevel transfer problems can also be treated using the modified Feautrier method in conjunction with the complete linearization technique of AUER and MIHALAS. t6} The
A modified Feautrier method
595
populations of the L levels may be used as independent variables, and the equations of statistical equilibrium act as the equations of constraint. Upon linearization of the equations, one obtains equations of the type just described in the preceding paragraph, except that the matrices U i are composed of L matrices of tridiagonal form rather than diagonal form; this is due to the linearization of the opacities, which couples neighbouring depth points through the differential operator. This change does not affect the asymptotic estimates, which are still C"N2FL + cIV(NL) 3 for timing and N2(L2 + 1) for storage. Since these multilevel estimates are roughly the same as for projected complete linearization treatments using integral equations, the modified Feautrier method, with its ease of formulation, appears to be an attractive alternative to such integral equation treatments. AcknowledgementsThis work was done during the tenure of a visiting fellowship at the Joint Institute for Laboratory Astrophysics in Boulder, Colorado. The author wishes to thank Dr. J. CASTORfor several useful discussions concerning the formulation of multilevel problems and for providing a fast computer subroutine for the inversion of tridiagonal matrices. REFERENCES I. 2. 3. 4. 5. 6. 7.
P. FEAUTRIER,C. r. hebd. S~anc. Acad. Sci. Paris 258, 3189 (1964). E. H. AVRETTand D. G. HUMMER,Mon. Not. R. astr. Soc. 1.30, 295 (1965). E. H. AVRETTand R. LOESER,Smithsonian Astrophysical Obs. Spec. Rpt. No. 201. Camb., Mass. (1966). E. H. AVRETTand R. LOESER,Smithsonian Astrophysical Obs. Spec. Rpt. No. 303. Camb., Mass. (1969). R. G. ATHAYand A. SKUMANICH,Annls. Astrophys. 30, 669 (1967). L. H. AUER and D. MIHALAS,Astrophys. J. 158, 641 (1969). L. H. AUER, dstrophys. J. 150, L53 (1967).