Chebyshev spectral-SN method for the neutron transport equation

SN CHEBYSHEV SPECTRAL- METHOD FOR THE NEUTRON TRANSPORT EQUATION M. ASADZADEH1 AND A. KADEM2 Abstra t. We study onvergen e of a ombined spe tral and (SN ) dis rete ordinates approximation for a multidimensional, steady state, linear transport problem with isotropi s attering. The pro edure is based on expansion of the angular ux in a trun ated series of Chebyshev polynomials in spatial variables that results in the transformation of the multidimensional problems into a set of one-dimensional problems. The onvergen e of this approa h is studied in the ontext of the dis rete-ordinates equations based on a spe ial quadrature rule for the s attering integral. The dis rete-ordinates and quadrature errors are expanded in trun ated series of Chebyshev polynomials of degree  L, and the onvergen e is derived assuming L  t 4s where t and s are totaland s attering ross-se tions, respe tively. 1. Introdu tion In this note we develop spe tral approximations for two and three dimensional, steady state, linear transport equation with isotropi s attering, in bounded domains. The pro edure is based on the expansion of the angular ux in a trun ated series of Chebyshev polynomials in the spatial variables. We study the onvergen e of this method in two dimensional ase, where we use a spe ial quadrature rule to dis retize in the angular variables, approximating the s alar ux. The similarity of the spe tral method to the nite element method is evident: the bases fun tions have a onstant norm and the pro edure is to represent the approximate solution as a linear ombination of nite number of basis fun tions (trun ated series of Chebyshev polynomials) and then use a variational formulation. The main di eren e is that: the nite element bases fun tions are lo ally supported, whereas the hebyshev polynomials are global fun tions. See also [6℄ for further details. In [16℄ this approa h, with no onvergen e rate analysis, is onsidered for a trun ated series of general orthogonal polynomials. The detailed study in [16℄ is arried out for the Legendre polynomials, where an index mix aused that a signi ant drift term is argued to be of lower order and therefore its ontribution is not in luded in the estimates. We apply this pro edure using Chebyshev polynomials with,e.g., the advantage of having onstant weighted-L2 norms, and give a full onvergen e study in luding estimates of the ontribution from the whole drift term. The nal estimation is via an inverse iterative/indu tion argument, based on an estimate derived from some elementary properties of Chebyshev polynomials in Appendix I. In our knowledge onvergen e rate analysis, in this setting, is not onsidered in the literature. 1991 Mathemati s Subje t Classi ation. 65N35, 65D32, 82D75, 40A10, 41A50. Key words and phrases. Convergen e analysis, linear transport equation, isotropi s attering, Chebyshev spe tral method, dis rete-ordinates method. 1 2 M. ASADZADEH 1 AND A. KADEM 2 Related problems, in di erent settings, are studied in the nu lear engineering literature, see, e.g., referen es in Vilhena et al in [16℄. Barros and Larsen [4℄ arried out a spe tral nodal method for ertain dis rete-ordinates problems. Chebyshev spe tral methods for radiative transfer problems are studied, e.g., by Kim and Ishimaru in [11℄ and by Kim and Mos oso in [12℄. In, e.g., astrophysi al aspe ts, spe tral methods are onsidered for relativisti gravitation and gravitational radiation by Bonazzola el al in [6℄. A multi-domain spe tral method is studied by Grand lement et al [10℄, for s alar and ve torial Poisson equations. C++ software library, developed for multi-domain, is available in publi domain (GPL), . For more detailed study on Chebyshev spe tral method and also approximations by the spe tral methods we refer the reader to monographs by Boyd [7℄ and Bernardi and Maday [5℄. An outline of this paper is as follows: In Se tion 2 we derive the trun ated spe tral equations in 2 dimensions. In Se tion 3 we prove that a ertain weighted-L norm for the error in the dis rete-ordinates approximation of the spe tral solution is dominated by that of a quadrature approximation. In Se tion 4 we onstru t a spe ial quadrature rule and derive onvergen e rates for the quadrature error. Combining the results of Se tions 3 and 4, we on lude the onvergen e of the dis rete-ordinates for the spe tral method. Appendix I is devoted to ertain properties of the Chebyshev polynomials, that are frequently used in the paper, and also the proof of a ru ial estimate used in the approximation of the ontribution from the drift term. Finally in Appendix II we derive the spe tral equations in a three dimensional setting. http://www.lorene.obspm.fr 2 2. The two-dimensional spe tral solution Consider the two-dimensional linear, steady state, transport equation given by (2.1)   x = p Z (x; ; ) + 1 1 1 Z 2 s 0 ; 0 ( 0 2  (x; ; ) + t (x; ; ) os  y ! ; ) (x; 0 ; 0 )d0 d0 + S (x; ; ); in the re tangular domain = fx := (x; y) : 1  x  1; 1  y  1g and the dire tions in D = f(; ) : 1    1; 0    2g. Here (x; ; ) is the angular ux, t and s denote the total- and the di erential ross se tions, respe tively, (s (0 ; 0 ! ; ) des ribes the s attering from an assumed pre- ollision angular oordinates (0 ; 0 ) to a post- ollision oordinates (; )), and S is the sour e term. See [14℄ for the details. Note that, in the ase of one-speed neutron transport equation; taking the angular variable in a dis , this problem would orresponds to a three dimensional ase with all fun tions being onstant in the azimuthal dire tion of the z variable. In this way the a tual spatial domain may be assumed to be a ylinder with the rossse tion and the axial symmetry in z . Then D will orrespond to the proje tion of the points on the unit sphere (the \speed") onto the unit dis (whi h oin ides with D.) See, [1℄ for the details. Given the fun tions f (y; ; ) and f (x; ; ), des ribing the in ident ux, we seek for a solution of (2.1) subje t to the following boundary onditions: 1 2 CHEBYSHEV SPECTRAL-S N METHOD FOR THE NEUTRON TRANSPORT 3 For 0    2, let  0 <   1; (2.2) (x = 1; y; ; ) = f (y; ; 0);; x x== 11;; 1   < 0: For 1 <  < 1, let  0 < os   1; (2.3) (x; y = 1; ; ) = f (y; ; 0);; y y== 11;; 1  os  < 0: Expanding the angular ux (x; ; ) in terms of the Chebyshev polynomials, in the y variable, leads to 1 2 (2.4) (x; ; ) = I X i=0 ( ) () i x; ;  Ti y : Below we determine the rst omponent, i.e., (x; ; ) expli itly, whereas the other omponents, i (x; ; ); i = 1; : : : I , will appear as the unknowns in I one dimensional transport equations: We start to determine (x; ; ), by inserting (2.4) into the boundary onditions (2.3) at y = 1, to nd that: I X ( 1)i i (x; ; ); 0 < os   1; (2.5) (x; ; ) = f (x; ; ) 0 0 0 (2.6) 2 0 (x; ; ) = i=1 I X i=1 ( ) i x; ;  ; 1  os  < 0: where 1  x  1, 1 <  < 1, and we have used the fa t that for the Chebyshev polynomials T (x)  1, Ti (1)  1 and Ti ( 1) = ( 1)i . See Appendix I. If we now insert from (2.4) into (2.1), multiply the resulting equation by pTk yy2 , k = 1; : : : ; I , and integrate over y we nd that the omponents k (x; ; ), k = 1; ::::; I; satisfy the following I one-dimensional transport equations:   k (x; ; ) + t k (x; ; ) x Z Z  (2.7) = s (0 ; 0 ! ; ) k (x; 0 ; 0 )d0 d0 + Gk (x; ; ): 0 ( ) 1 1 1 2 0 The same pro edure with the boundary ondition (2.2) at x = 1, and (2.4) yields I X (2.8) ( 1; y; ; ) = f (y; ; ) = i ( 1; ; )Ti (y ): 1 i=0 Now multiply (2.8) by p ; = 1; : : : ; I , and integrating over y we nd that 2 Z f (y; ; ) pTk (y) dy: (2.9) ( 1 ; ;  ) = k  1 y Similarly, (note the sign of  below), the boundary ondition at x = 1 is written as I X 0 <   1: (2.10) i (1; ; )Ti (y ) = 0; Tk (y) k 1 y2 1 1 i=0 1 2 M. ASADZADEH1 AND A. KADEM2 4 Multiplying (2.10) by pTk yy2 ; k = 1; : : : ; I and integrating over y, we get (2.11) 0 <   1; 0    2: k (1; ; ) = 0; We an easily he k that Gk in (2.7) is written as ( ) 1 (2.12) Gk (x; ; ) = Sk (x; ; ) where (2.13) Aki = 2 Z 1 1 p1 2 os  XI Ak i=k+1 i ( k x; ;  ) Tk (y ) d ( Ti (y )) dy dy 1 y2 p and 2 Z S (x; y; ; ) pTk (y) dy: (2.14) Sk (x; ; ) =  1 y Note that the solutions to the one dimensional problems given through the equations (2.7){(2.14) de ne the omponents k (x; ; ) for k = I; : : : ; 1; in this de reasing order to avoid the oupling of the equations. On e this is done, the angular ux is ompletely determined by (2.4). Here, we have used the onvention PIi I : : : = 0. Hen e, the staring GI (x; ; )  SI (x; ; ). Note also that although the solution, developed in here, rely on spe i boundary onditions the pro edure is quite general in the sense that the expression for the rst omponent, (x; ; ); keeps the information from the boundary onditions in the y variable, while the other omponents are derived based on the boundary onditions in x. 3. Convergen e of the spe tral solution In the sequel we fo us on the two dimensional, steady state linear transport pro ess with isotropi s attering, i.e., s (0 ; 0 ! ; )  s = onstant. For this problem we show, using a weighted-L norm, onvergen e of the spe tral solution de ned for the spatial variables. More spe i ally we show that: in a ertain weighted- L norm, the (trun ated) dis rete ordinates approximation error for the spe tral solution is dominated by that of a spe ial quadrature error. The study of onvergen e of this quadrature approximation is the matter of the next se tion. Assuming isotropi s attering, the equation (2.1) is written as p   ( x; ; ) + 1  os  (x; ; ) + t (x; ; )  x y (3.1) Z Z 0 0 0 0 = s (x;  ;  )d d + S (x; ; ) for x 2 := f(x; y) : 1  x  1; 1  y  1g,  2 [ 1; 1℄ and  2 [0; 2℄. The study of the problem with the anisotropi s attering is a rather involved task. See, e.g., [3℄ for an approa h involving anisotropi s attering. Consider now the dis rete ordinates (SN ) approximation of the equation (3.1): for m = 1; :::; M , let 1 2 1 = +1 0 2 2 2 1 1 (3.2) where (3.3) m 2 0 X M   ( x ) +  ( x ) +  ( x ) =  !n n (x) + Sm (x); m m m t m s x y n=1 m p = 1 2m os m ; SN CHEBYSHEV SPECTRAL- METHOD FOR THE NEUTRON TRANSPORT 5 and m (x) := m (x; y) is the angular ux in the dire tions de ned by m and m and asso iated with the quadrature weights !m . Finally Sm (x) is the orresponding inhomogeneous sour e term de ned in the dis rete dire tion (m ; m ) 2 [ 1; 1℄ . We assume a quadrature mesh (m ; m ) 6= (0; 0),   <  < :::: < M ; (3.4)  <  < ::: < M; satisfying the following onditions: 2 1 2 1 (3.5) !m 2 M X  4=M; m=1 !m  4; = 1; :::; M: m Further, we assume that the dis rete-ordinates equation (3.2) satisfy the same boundary onditions, in the dis rete dire tions, as the ontinuous one, i.e., (3.1) (as stated is Se tion 2). We shall prove that, under ertain assumptions, the solution of the equation (3.2) would onverge to that of the equation (3.1) as M ! 1. To this approa h we de ne the error in the approximate ux by (3.6) "m (x) = (x; m ; m ) m = 1; :::; M; m (x); and the trun ation error in the quadrature formula as (3.7)  (x) = Z Z 2 1 1 (x; 0 ; 0 )d0 d0 0 M X n=1 !n (x; n ; n ): Subtra ting the dis rete ordinates equation (3.2) from the ontinuous equation (3.1) in the dis rete dire tions, we obtain, for ea h m = 1; : : : ; M , an equation relating the dis rete ordinates approximation error to the quadrature error, viz, M X (3.8)  "m(x) +  "m(x) +  " (x) =  ! " (x) +   (x): m m x t m y s n n n=1 s We expand both the approximation and the quadrature errors in a trun ated series of Chebyshev polynomials in y, L X "m (x) = "lm (x)Tl (y ); (3.9) l=0 (3.10) and de ne the l (3.11)  (x) = L X  l (x)Tl (y ) l=0 moments of the errors by M 2 Æ Z X l; "l = !m ("lm (x))  th 0 1 2 1=2 1 Æl;0 m=1 Z 1 2 dx 1=2 ; ( l (x)) dx : =  Remark. Note that (3.9) and (3.10) involve further, trun ated, approximations of  (x), in (3.7) and the solution "m (x) of (3.6). We keep using the same notation as before the trun ations. Also, despite the re ent trun ation in y, we use equalities in (3.9), (3.10), as well as in the subsequent relations below. (3.12) l 2 1 M. ASADZADEH1 AND A. KADEM2 6 The main result of this paper is as follows: Theorem 3.1. L = O( ) = 4 l = 0; 1; : : : ; L " ! 0; M ! 1: In the remaining part of this se tion we show that, for !  4=M; m = 1; : : : ; M , the L norm of the trun ated spe tral error " , ounted in a reverse order on l = L; L 1; : : : ; 0, is dominated by that of the quadrature error  . The next se tion is devoted to proof of the following result: Theorem 3.2. !  4=M; m = 1; ; : : : ; M 2 L (; ) ! 0; M ! 1:  To prepare for the proof of the Theorem 3.1, we substitute (3.9) and (3.10) into the equation (3.8) to get X X X d" (x) dT T (y ) +  (y) +  " (x)T (y) " (x)  Let , where s, tr l then for , as m l 2 l For , if m l L (3.13) m l dx l=0 as L l m L l l m m = M X s !n L X n=1 Tj (y ) 1 y2 (3.14) l=0 =2 where j "ln (x)Tl (y ) + s L X  l (x)Tl (y ): l=0 , j = 0; : : : ; L and integrating over y yields X d"j (x)  m m +m Æj;0 dx (3.15) l l=0 l=0 L 2 l m t dy l=0 Multiplying (3.13) by p , then 1 Z (l) = j (l)" (x) + 2 l m  t "jm (x) Æj;0 X  s !n "jn (x) + Æj;0 n=1 M 2  s  j (x); Æj;0 dTl (y)
pTj (y) 2 dy: dy 1 y 1 1 Finally, we multiply the equation (3.14) by " (x) and integrate over x to obtain Z Z X  d" (x)  dx +  " (x) ( l) " (x)" (x)dx 2 Æ dx j m 1 j;0 (3.16) m 1  t Æj;0 +2 =2 X  s !n Æj;0 n=1 M Z Z m j  "jm (x) 2 1 m 1 j m j m "jm (x)"jn (x)dx + m j m l m dx 2 Now we rewrite the rst term in equation (3.16) as Z d" (x)   (3.17)  " (x) dx = dx 2 (" (1)) 1 1 l=0 1 1 1 1 L j m j m j m  s Æj;0 2 Z 1 1 "jm (x) j (x)dx: (" ( 1)) j m 2  : CHEBYSHEV SPECTRAL-S N METHOD FOR THE NEUTRON TRANSPORT 7 Note that m [("jm (1))2 ("jm ( 1))2 ℄ > 0. Indeed, for m > 0, using the boundary ondition "m ( 1; y ) = 0 and the identity (3.18) " (x) = j m 2 Æj;0  Z 1 1 "m (x; y )Tj (y ) p 1 1 dy; y2 we nd that "jm ( 1) = 0. The same is valid for x = 1, when m < 0. Consequently, 2 (3.19) X Æj;0 m  L j (l) Z 1 1 l=0  Z "jm (x)"lm (x)dx + t Z M X !n s 1  "jm (x) 1 2 " (x)" (x)dx + s j m 1 n=1 1 dx Z 1 j n 1 "jm (x) j (x)dx: To pro eed we multiply the inequality (3.19) by !m and sum over m to obtain Z M X 1 t !m " (x) j m 1 m=1
2 dx  s + s (3.20) Z " 1 1 Z " 1 1 2 Æj;0  M X #2 M X !m " (x) m=1 M X # !m "jm (x)  j (x)dx m=1 " !m m m=1 dx j m L X j (l) Z 1 # " (x)" (x)dx j m l m 1 l=0 := I + II + III: The ru ial part is now to estimate the -term III using the elementary properties of the Chebyshev polynomials. We start with the simpler terms I and II : Lemma 3.3. With !m  4=M; m = 1; : : : ; M , jI j 4 2 Æ p jII j  4 2 Æ j;0 (3.21) "j s Proof. j;0 j = 0; : : : ; L, 2 "j s we have, for j We use the elementary relation (a1 + a2 + : : : + aM )2  M (a 2 1 + a22 + : : : + a2M ); to write " (3.22) M X #2 !m " (x) j m M m=1 max   1 m M j! j M X m  !m "jm (x) 2 : m=1 Integrating (3.22) over x and using !m  4=M we get (3.23) Z 1 1 " M X m=1 #2 !m " (x) j m dx  4 Z 1 M X 1 m=1  !m "jm (x) 2 dx; that 8 M. ASADZADEH1 AND A. KADEM2 and hen e the rst estimate follows re alling (3.11). As for the se ond estimate, applying the Cau hy-S hwarz inequality, (3.23), (3.11) and (3.12) we get 1 Z " 1 M X m =1 # !m " (x)  j (x)dx j m 0  (3.24) p  Z " 1 1 M X 11=2 #2 !m "jm (x) =1 1 X m Z M 4 1 =1 p  dxA  !m "jm (x) 2 Z !1=2 dx 1 1 r  2 m  4 k"j kk j k; 2 Æj;0  1=2 2  j (x) dx  k j k Æj;0  whi h gives the desired estimate for II and the proof is omplete. Next using the Proposition 5.1 from the Appendix I we estimate the ontribution from the term III and derive the following key estimate: Proposition 3.4. For k" (3.25) L k k k = 0; 1; 2; : : : ; L,  j ( 1)j +k 1 k X  =0 we have the re ursive estimates  j )k"L (L Hen e, in parti ular the starting estimate, for p k" k  (3.26) L 4s  k = 0, j k+ p 4s  k L k k: is: k k: L With these estimates we an now easily prove our main result: Proof of Theorem 3.1. Proposition 3.4 and Theorem 3.2 give the desired result. Proof of Proposition 3.4. By the Proposition 5.1 (see Appendix I) we have that (3.27) j whereas for j < l,  (3.28)  j (l ) = (l) = 0; 0; l; for j for for  l; j+l j+l even odd. Therefor if we start with j = L , then j (L) = 0 and hen e (3.20) ombined with the de nition (3.11) and Lemma 3.3 yields (3.29) t  "L 2  4  "L 2 p + 4s  "L L : 2 2 2 Now rearranging the terms and re alling that  := t 4s we obtain (3.26). The proof of (3.25) is a reversed indu tive argument as follows: For j = L 1 we have that j (L) = L 1 (L) = L , whereas L 1 (l) = 0, for l < L. Hen e, using (3.27) we get (3.30) L X =0 l l j (l )"m (x) = s L X =0 l L 1(l)" (x) = l m L 1(L)" (x) = L" (x): L m L m CHEBYSHEV SPECTRAL-S N METHOD FOR THE NEUTRON TRANSPORT 9 Thus using the Cau hy-S hwarz inequality jIII j = Z  2 L 1 X Z M X 1 1 1 (x) dx 2 !m ["L (x)℄ dx m r "L M X 1  L !m ["m 1 (x)℄ 2 1=2 dx 1 L "L = "L (3.32) "L 2 2 1  4  " (3.33) L "L 2 is even), L X (3.34) L 2( L L  2 1 j t =  2L 1 The same pro edure applied to j+L  +  4s j "L 2 + L or equivalently using the notation here s p " L = l)"lm (x) = 2( L = L L 1 "L 1, we get "L 1 ; 1 4s , L + p L) j( 1) and 1)"m "L 1 1 L : 1 = L l) = 0, L) = 0, l<L 2( 2( L L 2( L 4s 2 yields 1) = (L : 1 Inserting in (3.20) and using also (3.11) and Lemma 3.3, with   m=1 "L 2 1=2 m=1 1 t  m=1 Z 2 (x) dx l=0 L m !m "L (x)"m m m  1 M m  2L L l)"lm (x)"m 1( L m=1 1 1 r L X 1 !m m  2L(max j j) (3.31) Z  M X Æj;0  2 (x) = (L for 1)"m 1 L (note that 1. Thus (x); l=0 so that, as in the previous step  "L (3.35) 2  2(L 1) Similarly sin e for j = (L 2) and L L 3; we have l) = 0 for l < L 3( L X (3.36) l)"lm (x) = 3( L L + p 4s L) = L , 3( L L 3( L 2 : 1) = 0, L L 3( 2) = 2, we get L 3( 2)"m L 2 (x) + L L)"L m (x) 3( l=0 = 2(L 2)"m L 2 (x) + 2L"m (x); L whi h using the same pro edure as before yields (3.37) L  "L 3  2L "L + 2(L 2) "L 2 + p 4s L 3 : Now the formula (3.25) is proved by an indu tion argument.  4. The quadrature rule and Proof of Theorem 3.2 In this se tion we onstru t a spe ial quadrature mesh satisfying the in (3.5) and prove the Theorem 3.2 in this setting. remaining step in the proof of the Theorem 3.1 and onditions This would provide us the omplete the onvergen e 10 M. ASADZADEH1 AND A. KADEM2 analysis. We also derive onvergen e rates for the quadrature error (3.7) where we identify the angular domain (4.1) D = f(; ) : 1    1; 0    2 g ; by n o ~ := (; ) : 1  ;   1;  = p1 2 os  : (4.2) D Then the quadrature ( ubature) rule, for the multiple integral in (3.1) an be onstru ted using (4.2) as in (3.7), see [9℄. To derive onvergen e rates, below we onstru t an equivalent rule, dire tly dis retizing D given by (4.1), and with a somewhat general features: Z 2 Z 1 X (4.3) ( ; ; )dd  !kj ( ; k ; j ); 0 x 1 x
where
:= f(k ; j ); k = 1; :::; K and j = 1; :::; J; J  K g  D is a M = JK , dis rete set of points in D onsisting of the Gauss quadrature points k 2 [ 1; 1℄ 2 j asso iated with the equally spa ed j = J ; j = 1; :::; J; and weights !kj = Ak Wj where Wj = 2J , j = 1; :::; J and Ak are given below. Thus the error in (4.3) an be split into two de oupled quadrature errors: Z jeM ( )j :=  (4.4) + 2 Z 1 0 2 Z Z 0 := x 1 Ak =1 2 Z 0 K X ( ; ; )d x 1 K X k 1 X ( ; ; )dd 2 h Z k ( x; k ;  0 jeK [ (x; )℄j d + =1 )d =1 ( !kj Ak J X j K X k
=1 Ak jeJ x; k ; j ( x; k ;  ) d ( Wj ) x; k ; j [ ( x; k ) i )℄j; with the obvious notations for the two quadrature errors: (4.5) eJ [ ( ; )℄ := (4.6) eK x [ ( ; )℄ := Z 2 0 Z x 1 1 ( x; ;  J X )d j =1 K X ( ; ; )d x k =1 Wj ( x; ; j ); Ak ( x; k ;  ): Below we derive error estimates for the quadrature rules (4.5) and (4.6), with onvergen e rates with respe t to the assumed regularity of in  and . eJ [ ℄ (4.5) J  r (x;; ) j 2 [0; 2 ℄ [0; 2℄  r Z 2  r ( ; ; ) d; (4.7) jeJ [ ℄j  Cr optimal Lemma 4.1. Let points denote the error in . Suppose that , with is integrable on x 0 Jr where Cr is independent of J equally spa ed quadrature and . r , then CHEBYSHEV SPECTRAL-SN METHOD FOR THE NEUTRON TRANSPORT Lemma 4.2. Let eK [ ℄ denote the error in imation of the integral of integrable on Cs  , then jeK [ ℄j  KCss (4.8) where [ 1; 1℄ on is independent of K Z 2 [ 1; 1℄ K -point Gaussian quadrature approx- . Suppose that 1 s 1 and (x; ; )
(1 s 11  s (x;; ) 2 s=2 s (1 ) is 2 s=2 d; ) . We postpone the proofs of these lemmas and rst derive the proof of Theorem 3.2 from them. For the transport equation (3.1), in polygonal domains, the regularity requirements in the lemmas 4.1 and 4.2 are proved for r = s = 1 in [1℄:  2 Lw  2 L [0; 2 ℄ ~ Proposition 4.3. w ~ := (1 2)1=2 1 1 [ 1; 1℄    (x) (4.3)   (4.9) k kL2 ( )  C J1 + K1 kgkH 1 ( ) ; ~ = R 11 R02 g (3.1) g = s ~ + S H 1( ) L2 Now we are ready to derive our nal error estimate: We multiply (3.10) by pTk1(yy)2 ; k = 0; : : : ; L, integrate over y 2 [ 1; 1℄ and use the Cau hy-S hwarz inequality to get for l = 0; : : : ; L, 2 Æl;0 Z 1  (x) pTl (y) dy  l (x) =  1 y2 1 Z 1 1=2  Z 1 1=2 dy dy  (x)2 p  2 Æl;0 Tl (y )2 p (4.10) 1 y2 1 y2 1 1 Z 1=2 1  = 2 Æl;0  (x)2 p dy 2 : 1 y 1 Now re alling (3.12) it follows that Z Z 1 1=2 2 Æl;0  1 dy l dx  C k kL2 ( ) : (4.11)    (x)2 p  1 y2 1 1 Combining with (4.9), re alling also M  J 1=2  K 1=2 we get the desired result.  and Let Then for the quadrature error where is the right hand side of is the usual , where of the approximation with , i.e., . we have, -based Sobolev spa e of order one on , and . Proof of Theorem 3.2. The onvergen e rate in Lemmas 4.1 and 4.2, as well as the rates in Proposition 4.3, an be improved up to the optimal order O(J 2 " )  O(K 2 " ), " arbitrarily small, for the neutron transport equation, in polygonal domains using, e.g., a post pro essing pro edure f. Asadzadeh [2℄. Now it remains to verify the estimates in Lemmas 4.1-4.2. We may assume that is 2-periodi in  and in the quadrature formula Remark. Proof of Lemma 4.1. (4.12) Z 0 2 (x; ; )d  J X j =1 Wj (x; ; j ); 12 M. ASADZADEH1 AND A. KADEM2 approximate by trigonometri polynomials in . Then we an easily he k that: no matter how we hoose the quadrature points j and weights Wj , the formula (4.12) an not be exa t for trigonometri polynomials of degree J , (see, e.g., [13℄ for the details). It turns out that the highest degree of pre ision J 1 is a hieved just for our simplest quadrature formula: equally spa ed nodes j = 2Jj and onstant weights Wj = 2J ; j = 1; 2; : : : ; J . Thus we have Z 2 J   X (j 1) 2 : (4.13) ()d  2 0 J J j =1 We an easily verify that (4.13) is exa t for the fun tions eimx ; m = 0; 1; : : : ; J 1. Further a trigonometri polynomial of degree J , with the Fourier series expansion ( )  a20 + (4.14) TJ x J X j =1 (aj os jx + bj sin jx); having 2J + 1 degrees of freedom (a0 ; aj ; bj ; j = 1; : : : ; J ) orresponds to an algebrai polynomials of degree 2J . Thus (4.13) is exa t for algebrai polynomials of degree 2J 1, so that for 2 C (r) [0; 2℄; r = 2J , ( is 2J times ontinuously di erentiable in ), using Taylor expansion up to degree 2J 1, in both sides of (4.12), we obtain the desired result.  Lemma 4.2 is a spe ial ase of the a lassi al result due to DeVore and S ott (Theorem 3 in [8℄, Proposition 4.4 below): Consider, for positive integer s, the fun tion spa e (4.15) 2 Yws := fu 2 L1lo (℄ 1; 1[) : kukw;s < 1g with w being a weight fun tion and (4.16) kukw;s = Z 1 1 [ju()j + ju(s) ()j(1 2 s w  d; )℄ ( ) where u(s) is interpreted as a weak derivative. Proposition 4.4 (DeVore and S ott). Let eK [ ℄ denote the error in K -point Gaussian quadrature approximation of the integral of on ; . Suppose that [ 1 1℄ (1 (weak derivative) is integrable on [ 1; 1℄, i.e., 2 Y1s , where s s is any positive integer su h that 1  s  2K . Then p Z 1 n 1 2 s ; (1 2 )s o d;  s (x; ; ) min (4.17) je [ ℄j  C  2 )s s (x;;) K s 1 s where Cs is independent of K and K . This follows, evidently, from the Proposition 4.4.  Below we review a pro edure, based on analyzing the Peano kernel for the quadrature error (4.6), and establish the bound (4.8) for s = 1, see [1℄ or [8℄. This would suÆ es to justify the use of Proposition 4.3. The full proof of (4.8), or (4.17), for s  1 is treated as in [8℄. Consider the Gauss quadrature rule Proof of Lemma 4.2. (4.18) Z 1 1 (x; ; )d  K X k=1 Ak (x; k ; ); CHEBYSHEV SPECTRAL-SN METHOD FOR THE NEUTRON TRANSPORT where (4.19) k := os k ; k 13 i h 2 (22kK +1)1 ; 2K2k+ 1 ; ; k = 1; : : : ; K; are zeros of Legendre polynomials and (4.20) Ak := Z1Y x xl dx; k = 1; : : : ; K; 1 l6=k xk xl are the integrals of the asso iated Lagrange interpolation polynomials. Now using the Peano kernel theorem we an write (4.21) eK [ ℄ = Z1 1 ( ) 0 ( ) d; where ( ) = eK [H ℄, j j  1 and H is the Heaviside fun tion   < ; (4.22) H () := 01;;   : It follows that X X (4.23) ( ) = 1  Ak = Ak  1: k > k < By the Chebyshev-Markov-Stieltjes ( f. [17℄ p. 50) inequality we have 1 + k  (4.24) k X i=1 Ai  1 + k+1 ; k = 1; : : : ; K: Thus with 1 = 0 < 1 < : : : < K < K +1 = 1 we get for k = 1; : : : ; K that (4.25) k 1 k  (k )  0  (k +)  k+1 k : Sin e  vanishes on ea h interval [k 1 ; k ℄ and has the slope one almost everywhere, we have (4.26) maxfj()j :  2 [k 1 ; k ℄g  k k 1 ; k = 1; : : : ; K: To bound k k 1 , we de ne Ik := [ k 1 ; k ℄, then (4.27) k k 1 = os k 1 ( k os k = Z k fsin k 1 ) max 2Ik k 1 sin d g  23K max fsin g: 2Ik Now sin e (sin )= is de reasing in [0; ℄, using (4.19) we get   k   sin k 1  sin k 1  4 sin (4.28) sin  k 1 k 1 k 1; 2 Ik ; for k = 2; : : : ; K . By the symmetry properties of j ( f. [17℄) we also get (4.29) sin  4 sin k ; 2 Ik ; k = 2; : : : ; K: Thus for k = 2; : : : ; K , p (4.30) max f sin g  4 min f sin g  4 min f 1 os2 g: 2I 2I 2I k k k 14 M. ASADZADEH 1 AND A. KADEM 2 Hen e, ombining (4.27) and (4.30), and using (4.19) we have for k = 2; : : : ; K , p p 6 6 (4.31)   = min f 1 os g = min f 1  g: k k 2 K 2Ik Thus, by (4.26), for  2 [ ; N ℄, 1 1 p K 2Ik 2 (4.32) j()j  6 K1  : The orresponding estimate for  2 [ 1;  ℄ and  2 [N ; 1℄ is (see [8℄): p  1  (4.33) j()j  p : 2K Summing up we have shown Z p 
1  d: (4.34) je [ ℄j  6 2 1 2 1 K K 1  2 This proves (4.8) for s = 1. For further details we refer to [1℄ and [8℄. Referen es [1℄ Asadzadeh, M., Analysis of Fully Dis rete S heme for Neutron Transport in Two-Dimensional Geometry. SIAM J. Numer. Anal. 23 (1986), pp 543-561. , Lp and eigenvalue error estimates for dis rete ordinates method for two-dimensional [2℄ neutron transport, SIAM J. Numer. Anal. 26 (1989), pp 66{87. , The Fokker-Plan k Operator as an Asymptoti Limit in Anisotropi Media. Math. [3℄ Comput. Modelling, 35 (2002) pp 1119-1133. [4℄ Barros, R. C. and Larsen E. W., A Spe tral Nodal Method for One-Group X; Y {Geometry Dis rete Ordinates Problems. Nu l. S i. and Eng., 111 (1992). [5℄ Bernardi, C. and Maday, Y., Approximations spe trales de problemes aux limites elliptiques, Springer-Verlag, Paris, (1992). [6℄ Bonazzola, S., Gourgoulhon, E, and Mar k, J.-A., Spe tral methods. Relativisti gravitation and gravitational radiation, Comb. Contemp. Astrophys., Cambridge University Press, (1977). [7℄ Boyd, John. P., Chebyshev and Fourier Spe tral Methods, Se ond Edition, Dover Publi ation, New York, (2001). [8℄ DeVore, R. A. and S ott, L. R., Error bounds for Gaussian quadrature and weighted-L1 polynomial approximation, SIAM J. Numer. Anal., 21 (1984), pp 400{412. [9℄ Engels, H., Numeri al Quadrature and Cubature, A ademi Press, London, (1980). [10℄ Grand lement, P., Bonazzola, S., Gourgoulhon, E, and Mar k, J.-A., A multidomain spe tral method for s alar and ve torial Poisson equations with non ompa t sour es, J. Comp. Phys. 170 (2001), pp. 231{260. [11℄ Kim, Arnold. D. and Ishimaru, Akira, A Chebyshev Spe tral Method for Radiative Transfer Equations Applied to Ele tromagneti Medium, Wave Propagation and S attering in Dis rete Random J. Comp. Phy. 152 (1999), pp. 264{280. [12℄ Kim, Arnold. D. and Mos oso, Miguel, Chebyshev Spe tral Methods for Radiative Transfer, SIAM J. S i. Comput., 23 (2002), pp. 2075{2095. [13℄ Krylov, V. E., Approximate al ulation of integrals, Translated by Stroud, A. H., Ma Millan, New York, London, (1962). [14℄ Lewis, E. E. and Miller, W. F. Jr. Computational Methods of Neutron Transport, John Wiley & Sons, New York, (1984). [15℄ Rivlin, T. J., The Chebyshev Polynomials.John Wiley & Sons, New York, (1974). [16℄ Vilhena, M. T., Bari hello, L. B., Zabadal, J. R., Segatto, C. F., Cardona, A. V., and Pazos, R. P., Solution to the multidimensional linear transport equation by the spe tral method, Progress in Nu lear Energy, 35 (1999), pp 275-291. [17℄ Szego, G., Orthogonal Polynomials, AMS Colloquium Publi ations 23, Ameri an Mathemati al So iety, New York, (1957). CHEBYSHEV SPECTRAL- N METHOD FOR THE NEUTRON TRANSPORT S 5. 15 Appendix I: Elementary properties of Chebyshev Polynomials Chebyshev polynomials are weighted orthogonal polynomials de ned by (5.1) Tn (x) = os(n ar os(x)); with the weight fun tion w(x) = p11 x2 . Thus Chebeshev polynomials are a sublass of Ja obi polynomials, where the Ja obi weights wJ = (1 + x)a (1 x)b , a; b > 1 are restri ted to a = b = 1=2. It follows that 8 < T (x)T (x)w(x) dx = : 1 Z1 (5.2) i j Hen e kT k (5.3) i w = i 6= j i = j 6= 0 i = j = 0: 0 =2   ; 2 Æi;0 i = 0; 1; : : : : Tn (x) is a polynomial of degree n, orthogonal to all polynomials of degree  n 1: On di erentiating Tn (x) = os n with respe t to x(= os ) we obtain a polynomial of degree n 1 alled the Chebyshev polynomials of se ond kind: sin n 1 ; Un 1 = Tn0 (x) = n sin (5.4) x = os : Further we an easily verify the following properties (see [15℄ for the details): For even (odd) n only even (odd) powers of x o ur in Tn (x). (5.5) (5.6) Tn ( x) = ( 1)n Tn (x): 1 U (x) + T2 (x) + T4 (x) + : : : + T2k (x) = 2k ; 2 2 k = 0; 1; : : : ; U (x) T1 (x) + T3 (x) + : : : + T2k+1 (x) = 2k+1 ; k = 0; 1; : : : ; : 2 Below we formulate and prove the property that has been essential in deriving the basi estimate in se tion 3 (Proposition 3.4.): (5.7) Proposition 5.1. Let (5.8) j (l) := Z1 d T (y ) Tl (y )
p j 2 dy; dy 1 y 1 we have that (5.9) and for (5.10) j (l) = 0; j < l, j (l) =  0; l; for j+l j+l j  l; even odd. M. ASADZADEH1 AND A. KADEM2 16 The rst assertion is a trivial onsequen e of the fa t that Tj is orthogonal to all polynomials of degree  j 1. As for the se ond assertion we note that Tl0 (x) = lUl (x): Thus if l is odd then l 1 is even, say l 1 = 2k, hen e using (5.6) Z h1 i Tj (y) dy + T (x) + T (x) + : : : + Tl (x)
p j (l ) = 2l 2 1 y (5.11)  odd, i.e., j + l even = 02;l  ; jj even, i.e., j + l odd: Æ 0 The ase l is even is treated similarly and using (5.7) and the proof is omplete.  6. Appendix II: The three-dimensional spe tral solution We extend now the approa h presented in Se tion 2 to the transport pro ess in three dimensions,   p    ( x; ; ) + 1  os  ( x; ; ) + sin  ( x; ; )  x y z (6.1) Z Z  + t (x; ; ) = s (0 ; 0 ! ; ) (x; 0 ; 0 )d0 d0 + S (x; ; ); where we assume that the spatial variable x := (x; y; z ) varies in the ubi domain := f(x; y; z ) : 1  x; y; z  1g; and (x; ; ) := (x; y; z; ; ) is the angular ux in the dire tions de ned by  2 [ 1; 1℄ and  2 [0; 2℄. We seek for a solution of (6.1) satisfying the following boundary onditions: For the boundary terms in x; for0    2, 0 <   1; (6.2) (x = 1; y; z; ; ) = f (y; z; ; 0);; x x== 11;; 1   < 0; For the boundary terms in y and for 1 <  < 1,  0 < os   1; (6.3) (x; y = 1; z; ; ) = f (x; z; ; 0);; y y== 11;; 1  os  < 0: Finally, for the boundary terms in z ; for 1 <  < 1, (6.4) (x; y; z = 1; ; ) = f (x; y; ; 0);; z z== 11;;  0<<2; : Here we assume that f (y; z; ; ), f (x; z; ; ) and f (x; y; ; ) are given fun tions. Expanding the angular ux (x; y; z; ; ) in a trun ated series of Chebyshev polynomials Ti (y)and Rj (z ) leads to Proof. 1 1 2 1 2 4 1 2 j; 2 1 1 2 0 1 2 3 1 (6.5) 2 (x; y; z; ; ) = 3 I X J X i=0 j =0 i;j (x; ; )Ti (y)Rj (z ): We repeat the pro edure in Se tion 2 and insert (x; y; z; ; ) given by (6.5) into the boundary onditions in (6.3), for y = 1. Multiplying the resulting expressions by pRj zz2 and integrating over z , we get the omponents ;j (x; ; ), j = 0; : : : ; J : ( ) 0 1 (6.6) j ;j (x; ;  ) = f (x; ;  ) 0 2 I X i=1 ( 1)i i;j (x; ; ); 0 < os   1; CHEBYSHEV SPECTRAL-SN METHOD FOR THE NEUTRON TRANSPORT and (6.7) 0;j X I (x; ; ) = i;j (x; ; ); 17 1  os  < 0: i=1 Similarly, we substitute (x; y; z; ; ) from (6.5) into the boundary onditions for z = 1, multiply the resulting expressions by p , i = 1; : : : ; I and integrate over 2 y , to de ne the omponents (x; ; ), i = 0; : : : ; I : For 1  x  1; 1 <  < 1; Ti (y ) 1 y i;0 (6.8) i;0 X J (x; ; ) = f (x; ; ) i ( 1) 3 j i;j (x; ; ); 0   < ; j =1 (6.9) i;0 X J (x; ; ) = i;j (x; ; ); <  2; j =1 where (6.10) f2 x; ;  and (6.11) f3i x; ;  j ( )=2 Z Æ0;j 1  1 Z ( f2 x; z; ;  ) pR (z ) 1 z j 2 dz pT (y) dy: )  1 y To determine the omponents (x; ; ), i = 1; : : : ; I , and j = 1; : : : ; J , we substitute (x; ; ), from (6.5) into (6.1) and the boundary onditions for x = 1. Multiplying the resulting expressions by p 2  p 2 ; and integrating over y and z we obtain I  J one-dimensional transport problems, viz   (x; ; ) +  (x; ; ) = G (x; ; ) x Z Z (6.12) +  ( ;  ! ; ) (x;  ;  )d d ; ( )= 2 Æi;0 1 1 ( i f3 x; y; ;  2 i;j Rj (z ) Ti (y ) 1 1 y z i;j t i;j i;j 1 1 0 0 0 s 1 0 0 0 i;j 1 with the boundary onditions (6.13) ( 1; ; ) = f (; ); where 4 Z Z p T (y)R (z ) f (y; z; ; )dzdy; (6.14) f (; ) =  (1 y )(1 z ) and (6.15) (1; ; ) = 0; for 0 <   1; and 0    2. Finally G (x; ; ) = S (x; ; ) 2 3 X X p (6.16) 1   4 os  A (x; ; ) + sin  B (x; ; )5 ; i;j 1 i;j i;j 1 2 1 1 1 1 i j 2 2 1 i;j i;j i;j I 2 J k i k=i+1 l j k;j l=j +1 i;l 18 M. ASADZADEH1 AND A. KADEM2 with (6.17) Si;j Z1Z1 4 (x; ; ) = 2 Z 2 = 1 1 T (y )R (z ) p S ( ; ; )dzdy; 1 (1 y2 )(1 z 2 ) i j x ( ( )) pT (y) 2 dy 1 y 1 Z 2 1 d (R (z )) pR (z ) dz: (6.19) B =  1 z2 1 dz Now, starting from the solution of the problem given by equations (6.12){(6.19) for (x; ; ), we then solve the problems for thePother omponents, in the deP reasing order in i and j . Re all that = +1 : : : = = +1 : : :  0: Hen e, solving I  J one-dimensional problems, the angular ux (x; ; ), is now ompletely determined through (6.5). (6.18) A k i  d Tk y dy i j l j l I ;J J I i I j J If we have to deal with a di erent type of boundary onditions, we have to keep in mind that the rst omponents 0 (x; ; ) and 0 (x; ; ) are determined from the boundary onditions for z and y and the other ones, (x; ; ) for i = 1; : : : ; I and j = 1; : : : ; J will satisfy one-dimensional transport problems subje t to the same type of boundary onditions of the original problem in the variable x. Remark: i; ;j i;j The resear h of the rst author is supported by the INTAS grant Ref. Nr 04-77-6902. A knowledgments: 1 Department of Mathemati s, Chalmers University of Te hnology, SE-412 96 Gothen- burg, Sweden. E-mail address : mohammadmath. halmers.se 2 Department of Mathemati E-mail address : s, Fa ulty of S ien e, University of Setif 19000, Algeria. abdelouahabkyahoo.fr