THE DISCRETE ORDINATES METHOD FOR

THE NEUTRON TRANSPORT EQUATION

IN AN INFINITE CYLINDRICAL DOMAIN^{∗}

Mohammad Asadzadeh, Peter Kumlin and Stig Larsson

Abstract. We prove a regularity result for a Fredholm integral equation with weakly singular kernel, arising in connection with the neutron transport equation in an infinite cylindrical domain. The theorem states that the solution has almost two derivatives in L_{1}, and is proved using Besov space techniques. This result is applied in the error analysis of the discrete ordinates method for the numerical solution of the neutron transport equation. We derive an error estimate in the L_{1}-norm for the scalar flux, and as a consequence, we obtain an error bound for the critical eigenvalue.

1. Introduction. We consider the numerical solution of the steady state isotropic monoenergetic neutron transport equation in a cylindrical domain in R^{3} with a polygonal cross section Ω. The restriction to the mono-energetic case means that we assume that the velocity space is the unit sphere S^{2} ⊂ R^{3}. The cylindrical symmetry reduces the problem to R^{2} by projection along the axis of the cylinder. Thus we study the neutron transport equation in a bounded polygonal domain Ω ⊂ R^{2} with the velocity space equal to the unit disc D ⊂ R^{2}.

We analyze a semidiscrete numerical method, the discrete ordinates method, involving the discretization of a weighted integral over D, described in polar coordinates, by means of an N-point Gaussian quadrature rule in the radial variable, and a uniform M-point quadrature rule in the angular variable.

For this method we give an L_{1} error estimate for the scalar flux of order N^{−4} + M^{−2+ϵ}, as well as an error bound of the same order for the critical eigenvalue. In order to prove these estimates we apply an error bound for weighted-L_{1} polynomial interpolation due to De Vore and Scott [7], together with a new L_{1} regularity result for the scalar flux. Loosely speaking, the latter result states that the scalar flux, which is the solution of a Fredholm integral equation with weakly singular kernel, has almost two derivatives in L_{1}. This is proved using Besov space techniques. The limited regularity of the solution strongly affects the error analysis. Although our main concern is L_{1}, we also derive some results in L_{∞}

Typeset by AMS-TEX

1991 Mathematics Subject Classification. Primary 65N15, 65N30.

Key words and phrases. Neutron transport equation, Gauss quadrature, error estimate, weighted-L_{1} interpolation, Fredholm integral equation, weakly singular kernel, regularity, Besov space, interpolation space, critical eigenvalue.

∗Published in Math. Models Methods Appl. Sci. 2 (1992), 317–338.

(regularity and error bounds), which are needed in a duality argument in the proof of the error estimate for the critical eigenvalue.

Problems of this type have been studied in various settings by several authors. The slab geometry, Ω ⊂ R^{1} and velocity space [−1, 1], was considered by Pitkäranta and Scott [11], who proved L_{p} and eigenvalue estimates for both semidiscrete and fully discrete schemes. Two-dimensional geometry, Ω ⊂ R^{2} and velocities in the unit circle S^{1}, was considered by Johnson and Pitkäranta [8] and Asadzadeh [2]. In [8] semidiscrete and fully discrete schemes were analyzed in L_{2}, whereas [2] contains L_{p} and eigenvalue estimates for the discrete ordinates method.

In Asadzadeh [3] the discrete ordinates method was studied in L_{2} in a fully threedimensional setting, Ω ⊂ R^{3} and velocity space S^{2}. In Asadzadeh [1] the geometry is the same as in the present work, but the analysis takes place in L_{2}, this norm being more suitable for spatial discretization based on the finite element method. Due to the limited regularity of the exact solution, the error bound in [1] for the angular discretization was O(N^{−1} + M^{−1}). By using the L_{1} norm together with our new regularity result (and a better analysis of the Gauss quadrature), we are able to improve this to O(N^{−4} +M^{−2+ϵ}). Thus, due to regularity limitations, an error estimate in the L_{1} norm for functions yields a better error bound for the critical eigenvalue. Moreover, the L_{1} norm is the most relevant norm from a physical point of view, since the scalar flux represents a particle density.

The outline of the paper is as follows: In Section 2 we describe the continuous problem that we want to solve, and in Section 3 we formulate the semidiscrete approximation and state our main result. Section 4 is devoted to an analysis of the regularity of the scalar flux. Section 5 contains the error analysis for the scalar flux, and in Section 6 the results obtained in the previous section are used to prove an error bound for the critical eigenvalue.

2. The continuous problem. We consider the following model problem for one-velocity neutron transport in an infinite cylinder Ω̃ ⊂ R^{3} with boundary Γ̃:

(2.1)  \begin{matrix}μ · ∇u(x, μ) + u(x, μ) = λ ∫_{S^{2}} u(x, η) dη + f(x),  (x, μ) ∈ Ω̃ × S^{2}, \\ u(x, μ) = 0,  (x, μ) ∈ Γ̃_{μ}^{−} × S^{2}. \end{matrix}

Here λ is a positive parameter and u(x, μ) is the density of neutrons at the point x ∈ Ω̃ flowing in the direction μ ∈ S^{2} = {μ ∈ R^{3} : |μ| = 1}. The boundary condition is specified on the inflow boundary

(2.2)  Γ̃_{μ}^{−} = {x ∈ Γ̃ : μ · n̂(x) < 0},

where n̂(x) is the unit outward normal.

We assume that the cross-section of the cylinder Ω̃ is a bounded convex polygonal domain Ω ⊂ R^{2} with boundary Γ. Assuming also that the source term is constant along the axial direction of the cylinder, we may project the integro-differential equation (2.1) onto the cross-section Ω:

(2.3) \begin{matrix}μ · ∇u(x, μ) + u(x, μ) = λ ∫_{D} u(x, η) (1 − |η|^{2})^{−1/2} dη + f(x),  (x, μ) ∈ Ω × D, \\ u(x, μ) = 0,  (x, μ) ∈ Γ_{μ}^{−} × D, \end{matrix}

where now the velocities μ vary over the unit disc D = {μ ∈ R^{2} : |μ| ≤ 1} and Γ_{μ}^{−} is defined analogously to (2.2).

Let us introduce the so-called scalar flux U defined by

(2.4)  U(x) = ∫_{D} u(x, μ) (1 − |μ|^{2})^{−1/2} dμ.

In order to derive an equation for U, we consider the following hyperbolic partial differential equation: given g = g(x) find v(x, μ) such that

(2.5)  \begin{aligned}μ · ∇v + v & = g,  (x, μ) ∈ Ω × D, \\ v & = 0,  (x, μ) ∈ Γ_{μ}^{−} × D. \end{aligned}

The solution is

(2.6)  v(x, μ) = (T_{μ}g)(x) = ∫_{0}^{d(x,μ)/|μ|} e^{−s}g(x − sμ) ds,

where d(x, μ) is the distance from x to the inflow boundary in the direction −μ:

d(x, μ) = inf\{s > 0 : x − s \frac{μ}{|μ|} ∉ Ω\}.

Using (2.3) and (2.6) we obtain the following equation for U:

(2.7)  (I − λT)U = T f,

where

(2.8)  T g(x) = ∫_{D} T_{μ}g(x) (1 − |μ|^{2})^{−1/2} dμ.

We show in Section 4 below that T is a compact operator on L_{1}(Ω) by rewriting it as an integral operator with weakly singular kernel. The Fredholm integral equation (2.7) thus has a unique solution U ∈ L_{1}(Ω) for any f ∈ L_{1}(Ω) as long as λ^{−1} is not an eigenvalue of T. The solution u of (2.3) can then be computed from (2.6) with g = λU + f.

3. The semidiscrete problem. Our aim in this section is to formulate a quadrature rule for the integral in (2.4) or (2.8) using a certain set of quadrature points Δ = {μ_{1}, …, μ_{n}} ⊂ D and weights ω_{μ}, μ ∈ Δ, thereby obtaining a semidiscrete approximation U_{n} of the scalar flux U. Thus, if

(3.1)  ∫_{D} v(μ) (1 − |μ|^{2})^{−1/2} dμ ≈ ∑_{μ∈Δ} v(μ)ω_{μ} ≡ Q[v(μ)]

is this quadrature rule, then we define U_{n} by

(3.2)  (I − λT_{n})U_{n} = T_{n}f,

where

T_{n}g(x) = ∑_{μ∈Δ} ω_{μ}T_{μ}g(x).

Using polar coordinates μ = rμ̂(φ), μ̂(φ) = (cos φ, sin φ), we want to define two quadrature rules

(3.3)  ∫_{0}^{1} v(r)\frac{r dr}{\sqrt{1 − r^{2}}} ≈ ∑_{k=1}^{N} A_{k}v(r_{k}) ≡ Q_{r}[v(r)],

(3.4)  ∫_{0}^{2π} v(φ) dφ ≈ ∑_{j=1}^{M} W_{j}v(φ_{j}) ≡ Q_{φ}[v(φ)],

and the quadrature rule in (3.1) will then be

Q[v(μ)] = Q_{φ}[Q_{r}[v(rμ̂(φ))]]

with

Δ = {r_{k}μ̂(φ_{j})}_{k=1}^{N}_{j=1}^{M}, ω_{kj} = A_{k}W_{j}.

For convenience we introduce the notation

n = NM, n̄ = min(N, M).

Thus n is the total number of quadrature points, and if n̄ is large then both N and M are large.

For the angular integral in (3.4) we choose the trapezoidal rule:

(3.5)  φ_{j} = \frac{2πj}{M} , W_{j} = \frac{2π}{M} ,  j = 1, …, M.

For the radial integral in (3.3) we choose the Gauss rule associated with the measure dα(r) = (1 − r^{2})^{−1/2}r dr:

(3.6)  r_{k} = sin θ_{k}, A_{k} = α(s_{k}) − α(s_{k−1}),  k = 1, …, N,

where α(r) = −\sqrt{1 − r^{2}}, and θ_{k}, s_{k} are certain points satisfying

(3.7)  \begin{matrix}θ_{k} ∈ [\frac{(2k − 1)π}{4N + 2} , \frac{2kπ}{4N + 2}] , \\ s_{k} ∈ (r_{k}, r_{k+1}), s_{0} = 0, s_{N} = 1, \end{matrix}

see Szegő [13], pp. 47–50.

In Lemma 5.7 below we demonstrate that equation (3.2) has a unique solution U_{n} ∈ L_{1}(Ω) for any f ∈ L_{1}(Ω), provided that λ^{−1} is not an eigenvalue of T and that n̄ is large. The main result of this paper is the following error bound for U_{n}. Its proof will be given in Section 5.

Theorem 3.1. Assume that λ^{−1} is not an eigenvalue of T. For each ϵ > 0 there is C such that, for large n̄,

∥U − U_{n}∥_{L_{1}(Ω)} ≤ C (\frac{1}{N^{4}} + \frac{1}{M^{2−ϵ}} ) (∥U∥_{L_{1}(Ω)} + ∥f∥_{W_{1}^{1}(Ω)}).

4. Regularity. In this section we study the mapping properties of the operator T defined in (2.8), which can be rewritten as

(4.1)  T g(x) = ∫_{Ω} ∫0^{\frac{π}{2}} \frac{exp(−|x − y|/cos t)}{|x − y|} dt g(y) dy,

where Ω is our bounded polygonal domain in R^{2} and g ∈ L_{1}(Ω). Hence T can be viewed as an integral operator with a weakly singular kernel. We recall that T is associated with the scalar flux U(x) via the Fredholm integral equation

(4.2)  U = λT U + T f.

Our main result in this section is the following estimate:

Theorem 4.1. Let f ∈ W_{1}^{1}(Ω) and let U ∈ L_{1}(Ω) be a solution of (4.2). Then for all θ ∈ [0, 1) and all q ∈ [1, ∞] there is a constant C = C(λ, θ, q) such that

∥U∥_{(L_{1}(Ω),W_{1}^{2}(Ω))_{θ,q}} ≤ C (∥U∥_{L_{1}(Ω)} + ∥f∥_{W_{1}^{1}(Ω)}).

Here (L_{1}(Ω), W_{1}^{2}(Ω))_{θ,q} denotes the real interpolation space between L_{1}(Ω) and W_{1}^{2}(Ω), see Bergh and Löfström [5]. The value of the exponent q is irrelevant for our discussion and we let q ∈ [1, ∞] be fixed throughout this section. Loosely speaking this result means that U has almost two derivatives in L_{1}. Since there is no restriction on λ, the result is also applicable to eigenfunctions of T. Furthermore we prove:

Theorem 4.2. The operator T is a bounded linear operator from L_{1}(Ω) into B_{1}^{s,q}(Ω) for all s ∈ [0, 1), and from L_{∞}(Ω) into (L_{∞}(Ω), W_{∞}^{1}(Ω))_{θ,q} for all θ ∈ [0, 1).

The Banach spaces B_{1}^{s,q}(Ω) appearing in the theorem are Besov spaces. We recall their definitions below. Employing the fact that B_{1}^{s,q}(Ω) is compactly imbedded in L_{1}(Ω) for s > 0, we conclude:

Corollary. T is a compact operator on L_{1}(Ω).

We remark that the regularity of solutions to (4.1) was also investigated by Pitkäranta [10] in terms of weighted Hölder spaces instead of Sobolev spaces as in the present work.

The main steps in the proof of Theorems 4.1 and 4.2 are the following:

(1) identifying T as a convolution operator on R^{2};

(2) imbeddings between interpolation spaces (W_{1}^{k}, W_{1}^{l})_{θ,q} and Besov spaces;

(3) multilinear interpolation;

(4) a bootstrap argument.

Step 1. Set h_{τ} (x) = exp(−|x|/τ)/|x| for x ∈ R^{2}, τ > 0 and note that

(4.3)  T g(x) = ∫_{Ω} ∫0^{\frac{π}{2}} \frac{exp(−|x − y|/cos t)}{|x − y|} dt g(y) dy = ∫0^{\frac{π}{2}} [h_{cos t} ∗ (χ_{Ω} g)](x) dt,

where χ_{Ω} is the characteristic function of the set Ω ⊂ R^{2} and ∗ denotes convolution in R^{2}. The functions T g(x) and h_{τ} (x) are defined for all x ∈ R^{2} and in the sequel we shall think of all functions in L_{1}(Ω) as being defined on R^{2} (extended by 0 on R^{2} \ Ω unless otherwise stated). We write L_{1} = L_{1}(R^{2}), B_{1}^{s,q} = B_{1}^{s,q}(R^{2}), etc. In particular, by (W_{1}^{k}, W_{1}^{l})_{θ,q} we mean (W_{1}^{k}(R^{2}), W_{1}^{l}(R^{2}))_{θ,q} although the final result refers to (W_{1}^{k}(Ω), W_{1}^{l}(Ω))_{θ,q}. The justification for this is the existence of continuous linear extension operators

(4.4)  W_{1}^{k}(Ω) → W_{1}^{k}(R^{2}),  k ∈ N,

(see Stein [12]), which are defined for bounded Lipschitz domains, and in particular, for the domain Ω above.

Step 2. All explicit norm calculations will be done in Besov spaces. The reason for this is that they have nice intrinsic norm definitions

∥f∥_{B_{1}^{s,q}} = ∑_{|α|≤[s]} ∥D^{α}f∥_{L_{1}} + ∑_{|α|=[s]}(∫_{0}^{1}(t^{[s]−s}ω(D^{α}f)(t))^{q} \frac{dt}{t})^{1/q},

for q ∈ [1, ∞] and noninteger s > 0, where [s] denotes the integer part of s and

ω(h)(t) = sup_{|η|≤t} ∥h(· + η) − h∥_{L_{1}},

and with the usual modification for q = ∞, see Bergh and Löfström [5]. Moreover, they are interpolation spaces, i.e.,

(B_{1}^{s_{1},q_{1}}, B_{1}^{s_{2},q_{2}})_{θ,q} = B_{1}^{θs_{2}+(1−θ)s_{1},q},

for 0 ≤ θ ≤ 1, s_{1}, s_{2} > 0, and they approximate the Sobolev spaces W_{1}^{k} in the sense that

B_{1}^{k+ϵ,q} ⊂ W_{1}^{k} ⊂ B_{1}^{k−ϵ,q},

for small ϵ > 0, k ∈ Z_{+}, see Nikolskii [9]. This trivially implies the imbeddings

(4.5)  B_{1}^{θl+(1−θ)k+ϵ,q_{1}} ⊂ (W_{1}^{k}, W_{1}^{l})_{θ,q} ⊂ B_{1}^{θl+(1−θ)k−ϵ,q_{2}},

for small ϵ > 0, 0 ≤ θ ≤ 1, k, l ∈ N. We also note the Sobolev imbedding

(4.6)  B_{1}^{s,q} ⊂ L_{p},  for 1 − \frac{s}{2} < \frac{1}{p}.

In the argument below we need the following lemmas.

THE NEUTRON TRANSPORT EQUATION 7

Lemma 4.3. We have

sup_{τ∈(0,1]} ∥h_{τ} ∥_{B_{1}^{s,q}} < ∞, for s ∈ [0, 1).

Proof. Recall that h_{τ} (x) = exp(−|x|/τ)/|x|. Hence we have

∥h_{τ} ∥_{L_{1}} = τ∥h_{1}∥_{L_{1}} ≤ C, for τ ∈ (0, 1],

and similarly

ω(h_{τ})(t) = sup_{|η|≤t} ∥h_{τ} (· + η) − h_{τ} ∥_{L_{1}} = τ sup_{|η|≤t/τ} ∥h_{1}(· + η) − h_{1}∥_{L_{1}} = τω(h_{1})(t/τ),

so that

(4.7)  (∫_{0}^{1}(t^{−s}ω(h_{τ})(t))^{q} \frac{dt}{t})^{1/q} = τ^{1−s}(∫_{0}^{1/τ} (t^{−s}ω(h_{1})(t))^{q} \frac{dt}{t})^{1/q}.

We shall show that

(4.8)  ω(h_{1})(t) ≤ C min(t(1 + | log t|), 1).

Then the expression in (4.7) is bounded by

τ^{1−s}(∫_{0}^{1/τ} (t^{−s}ω(h_{1})(t))^{q} \frac{dt}{t})^{1/q} ≤ C_{s,q},

for 0 < s < 1, and the desired result follows. For the proof of (4.8) we write

\begin{aligned}h_{1}(x + η) − h_{1}(x) & = e^{−|x+η|}(\frac{1}{|x + η|} − \frac{1}{|x|}) + (e^{−|x+η|} − e^{−|x|})\frac{1}{|x|} \\ & ≡ g_{1}(x, η) + g_{2}(x, η). \end{aligned}

Setting, for fixed η,

R = R(η) = max(1, 3|η|, \frac{1}{|η|}),

we have

∥g_{1}(·, η)∥_{L_{1}} = ∫_{A_{1}} |g_{1}(x, η)| dx + ∫_{A_{2}} |g_{1}(x, η)| dx,

with

A_{1} = {x ∈ R^{2} : 2|η| ≤ |x| ≤ R}, A_{2} = R^{2} \ A_{1}.

Straightforward calculations give

∫_{A_{1}} |g_{1}(x, η)| dx ≤ C|η|(1 + |log |η||)

and

∫_{A_{2}} |g_{1}(x, η)| dx ≤ C|η|.

Moreover, we obtain

\begin{aligned}∥g_{2}(·, η)∥_{L_{1}} & ≤ ∫_{|x|≤R} \frac{1}{|x|}|e^{−|x+η|} − e^{−|x|}| dx + C|η| \\ & ≤ C ∫_{|x|≤R} \frac{e^{−|x|}}{|x|} min(1, e^{|η|} − 1) dx + C|η| ≤ C|η|. \end{aligned}

Together these estimates prove (4.8). □ Lemma 4.4. Let s ∈ (0, 1) and let Ω be the polygonal domain above. Then the following estimates are true:

∥χ_{Ω} g∥_{B_{1}^{s,q}} ≤ C∥g∥_{B_{1}^{2s+ϵ,q}}, ϵ > 0, (4.9) ∥χ_{Ω} g∥_{B_{1}^{s,q}} ≤ C∥g∥_{W_{1}^{1}}, (4.10)

where C depends only on s, q, ϵ and s, q, respectively. Proof. To prove the first estimate we note that

ω(χ_{Ω} g)(t) ≤ ω(g)(t) + sup_{|η|≤t} ∫_{R^{2}} |g(x)| |χ_{Ω} (x + η) − χ_{Ω} (x)| dx.

Hölder’s inequality applied to the second term on the right-hand side together with the Sobolev imbedding (4.6) yields (\frac{1}{p′} ≡ s + δ = 1 − \frac{1}{p})

\begin{aligned}ω(χ_{Ω} g)(t) & ≤ ω(g)(t) + ∥g∥_{L_{p}} sup_{|η|≤t}(∫_{R^{2}} |χ_{Ω} (x + η) − χ_{Ω} (x)|^{p′} dx)^{1/p′} \\ & ≤ ω(g)(t) + Ct^{1/p′}∥g∥_{L_{p}} ≤ ω(g)(t) + Ct^{s+δ}∥g∥_{B_{1}^{2s+ϵ,q}}, \end{aligned}

for any ϵ > 0, 0 < δ < ϵ/2. This proves (4.9).

The estimate (4.10) follows from

(4.11)  ∥f∥_{B_{1}^{s,q}} ≤ C(∥f∥_{L_{1}} + ∑_{i=1}^{2} ∥∂_{x_{i}}f∥_{B_{1}^{s−1,q}})

applied to f = χ_{Ω} g and the estimate

(4.12)  ∥∂_{x_{i}}(χ_{Ω} g)∥_{B_{1}^{s−1,q}} ≤ C∥g∥_{W_{1}^{1}},  i = 1, 2.

Proof of (4.12): For g ∈ W_{1}^{1} we have by the divergence theorem that

(4.13)  ∂_{x_{i}}(χ_{Ω} g) = χ_{Ω} ∂_{x_{i}}g − gn_{i} dS,  in D′(R^{2}),

where n = (n_{1}, n_{2}) is the unit outer normal on Γ and dS denotes the arclength on Γ. Here we note that gn_{i} dS is well-defined in D′(R^{2}), since g has an L_{1}-trace on Γ. The first term on the right-hand side is an L_{1}-function. For the second term we use an equivalent form of the Besov norm,

∥f∥_{B_{1}^{s,q}} = (∑_{j=0}^{∞}(2^{sj}∥F^{−1}(ϕ_{j}Ff)∥_{L_{1}})^{q})^{1/q},

where s ∈ R, F denotes the Fourier transform and

ϕ_{j}(ξ) = ϕ(2^{−j}ξ), j = 1, 2, …,

ϕ_{0}(ξ) = 1 − ∑_{j=1}^{∞} ϕ_{j}(ξ),

where ϕ ∈ C_{0}^{∞}(R^{2}) satisfies

i) supp ϕ = \{ξ : \frac{1}{2} ≤ |ξ| ≤ 2\},

ii) ϕ > 0, in \{ξ : \frac{1}{2} < |ξ| < 2\},

iii) ∑_{j=−∞}^{∞} ϕ(2^{−j}ξ) = 1, for ξ ≠ 0.

To prove

∥gn_{i} dS∥_{B_{1}^{s−1,q}} ≤ C∥g∥_{W_{1}^{1}}

it suffices to show

∥F^{−1}(ϕ_{j}F(gn_{i} dS))∥_{L_{1}} ≤ C∥g∥_{W_{1}^{1}}, j = 0, 1, 2, …,

with C independent of j. The latter bound follows from

∥F^{−1}(ϕ_{j}F(gn_{i} dS))∥_{L_{1}} = ∥ ∫_{Γ} 2^{2j}(F^{−1}ϕ)(2^{j}(· − y)) (gn_{i} dS)(y)∥_{L_{1}}
≤ ∫_{R^{2}} ∫_{Γ} 2^{2j}|F^{−1}ϕ(2^{j}(x − y))| |g(y)| dS(y) dx
≤ ∫_{R^{2}} |F^{−1}ϕ(x)| dx ∫_{Γ} |g(y)| dS(y)
≤ C∥g∥_{L_{1}(Γ)} ≤ C∥g∥_{W_{1}^{1}}, j = 1, 2, …,

where in the last step we used the trace estimate. The case j = 0 follows in a similar way and completes the proof of (4.12).

Proof of (4.11): For i = 1, 2 set

ψ_{i}(ξ) = ξ_{i}ϕ(ξ),

ψ_{i,j}(ξ) = ψ_{i}(2^{−j}ξ), j = 1, 2, …,

ψ_{i,0}(ξ) = ξ_{i} − ∑_{j=1}^{∞} ψ_{i,j}(ξ).

With this notation we have

∥∂_{x_{i}}f∥_{B_{1}^{s−1,q}} = (∑_{j=0}^{∞}(2^{sj}∥F^{−1}(ψ_{i,j}Ff)∥_{L_{1}})^{q})^{1/q}, i = 1, 2,

and so it suffices to show

(4.14)  ∥F^{−1}(ϕ_{j}Ff)∥_{L_{1}} ≤ C ∑_{i=1}^{2} ∥F^{−1}(ψ_{i,j}Ff)∥_{L_{1}},  j = 1, 2, …,

and

(4.15)  ∥F^{−1}(ϕ_{0}Ff)∥_{L_{1}} ≤ C∥f∥_{L_{1}}.

Inequality (4.14) is equivalent to (not the same f)

∥F^{−1}(ϕFf)∥_{L_{1}} ≤ C ∑_{i=1}^{2} ∥F^{−1}(ψ_{i}Ff)∥_{L_{1}}.

Let χ_{i} ∈ C_{0}^{∞}(R^{2}), χ_{i} ≥ 0, be supported away from the ξ_{i}-axis for i = 1, 2, and such that χ_{1} + χ_{2} = 1 on supp ϕ. We obtain

∥F^{−1}(ϕFf)∥_{L_{1}} ≤ ∑_{i=1}^{2} ∥F^{−1}(\frac{1}{ξ_{i}} χ_{i}ψ_{i}Ff)∥_{L_{1}}

= ∑_{i=1}^{2} ∥F^{−1}(\frac{1}{ξ_{i}} χ_{i}) ∗ F^{−1}(ψ_{i}Ff)∥_{L_{1}}

≤ C ∑_{i=1}^{2} ∥F^{−1}(ψ_{i}Ff)∥_{L_{1}}.

Finally for (4.15) we find

∥F^{−1}(ϕ_{0}Ff)∥_{L_{1}} = ∥(F^{−1}ϕ_{0}) ∗ f∥_{L_{1}} ≤ C∥f∥_{L_{1}},

and (4.11) is proved. □

Remark. Actually the right-hand side of of (4.11) is an equivalent norm for B_{1}^{s,q}.

Step 3. Consider the bilinear mapping (f_{1}, f_{2}) ↦ f_{1} ∗ f_{2}. We observe that D^{α}(f_{1} ∗ f_{2}) = (D^{α_{1}}f_{1}) ∗ (D^{α_{2}}f_{2}) for any multi-index α = α_{1} + α_{2}. This yields Lemma 4.5. The following estimates hold true:

(4.16)  ∥f_{1} ∗ f_{2}∥_{B_{1}^{s,q}} ≤ C∥f_{1}∥_{B_{1}^{s_{1},q}} ∥f_{2}∥_{B_{1}^{s_{2},q}}, 0 < s < s_{1} + s_{2},
(4.17)  ∥f_{1} ∗ f_{2}∥_{(L_{∞},W_{∞}^{k})_{θ,q}} ≤ C∥f_{1}∥_{(L_{∞},W_{∞}^{k})_{θ,q}} ∥f_{2}∥_{L_{∞}}, k ≥ 0, 0 ≤ θ ≤ 1.

Proof. For fixed s, s_{1}, s_{2} satisfying the hypotheses in (4.16) we choose k > s_{1}, l > s_{2} and θ̄ ∈ (0, 1) such that kθ̄ < s_{1}, lθ̄ < s_{2}, (k + l)θ̄ > s. Multilinear interpolation applied to

∥f_{1} ∗ f_{2}∥_{W_{1}^{k+l}} ≤ ∥f_{1}∥_{W_{1}^{k}} ∥f_{2}∥_{W_{1}^{l}},
∥f_{1} ∗ f_{2}∥_{L_{1}} ≤ ∥f_{1}∥_{L_{1}} ∥f_{2}∥_{L_{1}},

yields

∥f_{1} ∗ f_{2}∥_{(L_{1},W_{1}^{k+l})_{θ,q}} ≤ C∥f_{1}∥_{(L_{1},W_{1}^{k})_{θ,q}}∥f_{2}∥_{(L_{1},W_{1}^{l})_{θ,q}},

for θ ∈ [0, 1]. An application of (4.5) now proves (4.16). The remaining estimate (4.17) is proved analogously. □

We can now prove Theorem 4.2. Let s ∈ [0, 1) and choose s_{1} ∈ (s, 1). Using the representation (4.3) and Lemma 4.5 together with Lemma 4.3 we obtain

∥T g∥_{B_{1}^{s,q}(Ω)} ≤ ∥T g∥_{B_{1}^{s,q}} ≤ ∫0^{\frac{π}{2}} ∥h_{cos t} ∗ (χ_{Ω} g)∥_{B_{1}^{s,q}} dt

≤ C ∫0^{\frac{π}{2}} ∥h_{cos t}∥_{B_{1}^{s_{1},q}} dt ∥χ_{Ω} g∥_{L_{1}} ≤ C∥g∥_{L_{1}(Ω)},

(4.18)

which is the first claim of Theorem 4.2. The second assertion is proved in the same way using (4.17) instead of (4.16).

Step 4. We now prove Theorem 4.1 by means of Theorem 4.2 and a bootstrap argument. Consider the integral equation (4.2), i.e.,

U = λT U + T f.

Assume that f ∈ W_{1}^{1}(Ω) and that U ∈ L_{1}(Ω) is a solution of (4.2). We shall show that

(4.19)  ∥U∥_{B_{1}^{s,q}(Ω)} ≤ C(∥U∥_{L_{1}(Ω)} + ∥f∥_{W_{1}^{1}(Ω)}),  for 0 ≤ s < 2,

from which the desired bound follows in view of the imbedding (4.5). If 0 ≤ s < 1, then (4.19) follows directly from (4.18). Thus let 1 ≤ s < 2. Let Ũ denote the extension of U to R^{2} by λT U + T f. Extend f ∈ W_{1}^{1}(Ω) to f̃ ∈ W_{1}^{1}(R^{2}) as in (4.4). Using again the representation (4.3) together with the bounds in Lemmas 4.4 and 4.3 we obtain

\begin{aligned}∥Ũ∥_{B_{1}^{s,q}} & ≤ C ∫0^{\frac{π}{2}} ∥h_{cos t}∥_{B_{1}^{1−ϵ,q}} dt ∥χ_{Ω} (λŨ + f̃)∥_{B_{1}^{s−1+ϵ′,q}} \\ & ≤ C(∥χ_{Ω} Ũ∥_{B_{1}^{s−1+ϵ′,q}} + ∥χ_{Ω} f̃∥_{B_{1}^{s−1+ϵ′,q}}) \\ & ≤ C(∥Ũ∥_{B_{1}^{2(s−1+ϵ′)+ϵ′′,q}} + ∥f̃∥_{W_{1}^{1}}), \end{aligned} (4.20)

for small ϵ, ϵ′, ϵ′′ > 0. Choosing ϵ′, ϵ′′ such that

s − (2(s − 1 + ϵ′) + ϵ′′) > \frac{2 − s}{2},

and repeating the argument in (4.20) a finite number of times we arrive at

∥Ũ∥_{B_{1}^{s,q}} ≤ C(∥Ũ∥_{B_{1}^{s_{1},q}} + ∥f̃∥_{W_{1}^{1}}),

with 0 < s_{1} < 1. Since ∥U∥_{B_{1}^{s,q}(Ω)} ≤ ∥Ũ∥_{B_{1}^{s,q}}, and since (4.18) shows

∥Ũ∥_{B_{1}^{s_{1},q}} ≤ C(∥U∥_{L_{1}(Ω)} + ∥f∥_{L_{1}(Ω)}),

and by (4.4)

∥f̃∥_{W_{1}^{1}} ≤ C∥f∥_{W_{1}^{1}(Ω)},

this completes the proof of (4.19) and Theorem 4.1 is proved.

5. Error analysis. The aim of this section is to prove Theorem 3.1. The proof is based on the observation that, since

(I − λT)U = T f, (I − λT_{n})U_{n} = T_{n}f,

we have

(I − λT_{n})(U − U_{n}) = (T − T_{n})(λU + f).

We show in Lemma 5.7 below that (I − λT_{n})^{−1} is uniformly bounded on L_{1}(Ω). This implies that

(5.1)  ∥U − U_{n}∥_{L_{1}(Ω)} ≤ C∥(T − T_{n})(λU + f)∥_{L_{1}(Ω)}.

We therefore need to estimate (T − T_{n})g, which can be viewed as the quadrature error

(T − T_{n})g(x) = ∫_{D} v(x, μ) \frac{dμ}{\sqrt{1 − |μ|^{2}}} − Q[v(x, μ)], with v(x, μ) = T_{μ}g(x).

This will also be the key to the proof of the boundedness of (I − λT_{n})^{−1}.

Using polar coordinates μ = rμ̂(φ), μ̂(φ) = (cos φ, sin φ), we split the quadrature error as follows

(5.2) \begin{aligned}(T − T_{n})g(x) & = ∫_{0}^{2π} (∫_{0}^{1} v(x, rμ̂(φ))\frac{r}{\sqrt{1 − r^{2}}} dr − Q_{r}[v(x, rμ̂(φ))]) dφ \\ & + Q_{r} [∫_{0}^{2π} v(x, rμ̂(φ)) dφ − Q_{φ}[v(x, rμ̂(φ))]] . \end{aligned}

We estimate the latter two terms in the following sequence of lemmas.

Lemma 5.1. Let Q_{φ} be defined by the trapezoidal rule (3.5). Then

|∫_{0}^{2π} u(φ) dφ − Q_{φ}[u(φ)]| ≤ \frac{C}{M^{s}} ∫_{0}^{2π} |u^{(s)}(φ)| dφ, s = 1, 2.

Proof. This is a standard result. □

Lemma 5.2. Let Q_{r} be defined by the Gauss rule (3.6). Then

|∫_{0}^{1} u(r)\frac{r}{\sqrt{1 − r^{2}}} dr − Q_{r}[u(r)]| ≤ \frac{C_{s}}{N^{s}} ∫_{0}^{1} |u^{(s)}(r)|(r(1 − r))^{\frac{s−1}{2}} r dr, s = 1, 2, ….

Proof. We modify the argument of De Vore and Scott [7] so as to handle an integral with weight w(r) = r/\sqrt{1 − r^{2}}. By Taylor’s formula and using the fact that N-point Gauss quadrature is exact for polynomials of degree ≤ 2N − 1, we can represent the error as

e[u] = ∫_{0}^{1} u(r)\frac{r}{\sqrt{1 − r^{2}}} dr − Q_{r}[u(r)] = ∫_{0}^{1} u′(r)E(r) dr,

where

E(r) = ∫_{r}^{1} \frac{t}{\sqrt{1 − t^{2}}} dt − ∑_{r_{k}>r} A_{k}.

We show below that

(5.3)  |E(r)| ≤ \frac{C}{N} r.

Hence

|e[u]| = inf_{P ∈P_{2N−1}}|e[u − P]| ≤ \frac{C}{N} inf_{P ∈P_{2N−2}} ∫_{0}^{1} |u′(r) − P (r)| r dr,

and our task is now reduced to a problem of weighted-L_{1} polynomial approximation. Applying a general result from [7] we obtain

inf_{P ∈P_{2N−2}} ∫_{0}^{1} |u′(r) − P (r)| r dr ≤ \frac{C_{s}}{N^{s−1}} ∫_{0}^{1} |u^{(s)}(r)|(r(1 − r))^{\frac{s−1}{2}} r dr,

which proves the desired result.

It remains to prove (5.3). Consider first the situation that r ∈ [0, r_{N}]. We may assume, without loss of generality, that r_{l−1} < r < r_{l}. Then, by (3.6), we have

E(r) = \sqrt{1 − r^{2}} − ∑_{r_{k}>r} A_{k} = \sqrt{1 − r^{2}} − α(s_{N}) + α(s_{l−1}) = \sqrt{1 − r^{2}} − \sqrt{1 − s_{l−1}^{2}},

and, since r_{l−1} < s_{l−1} < r_{l} , by Taylor expansion we get

|E(r)| ≤ \sqrt{1 − r_{l−1}^{2}} − \sqrt{1 − r_{l}^{2}} = cos θ_{l−1} − cos θ_{l} ≤ \frac{C}{N} r.

Finally, for r ∈ [r_{N}, 1] we have

E(r) = \sqrt{1 − r^{2}} ≤ \sqrt{1 − r_{N}^{2}} ≤ cos \frac{(2N − 1)π}{4N + 2} ≤ \frac{C}{N} . □

We now need to consider the regularity of v(x, μ) with respect to r and φ.

Lemma 5.3. Let S_{i}, i = 1, ..., P, denote the sides of the polygonal domain Ω and let τ̂_{i} = (cos ψ_{i}, sin ψ_{i}) be a tangent to S_{i}. Let v(·, μ) = T_{μ}g with μ = rμ̂(φ), μ̂(φ) = (cos φ, sin φ). Then for 0 ≤ r ≤ 1, 0 ≤ φ ≤ 2π, we have

(5.4)  ∥\frac{∂v(·, μ)}{∂φ}∥_{L_{1}(Ω)} ≤ C∥g∥_{W_{1}^{1}(Ω)},

(5.5)  ∥\frac{∂v(·, μ)}{∂φ}∥_{L_{∞}(Ω)} ≤ \frac{C}{min_{1≤i≤P} | sin(φ − ψ_{i})|} ∥g∥_{W_{∞}^{1}(Ω)},

(5.6)  ∥\frac{∂^{2}v(·, μ)}{∂φ^{2}} ∥_{L_{1}(Ω)} ≤ \frac{C}{min_{1≤i≤P} | sin(φ − ψ_{i})|} ∥g∥_{W_{1}^{2}(Ω)},

(5.7)  ∥\frac{∂^{2}v(·, μ)}{∂r^{2}} ∥_{L_{∞}(Ω)} ≤ \frac{C}{r^{2}} ∥g∥_{L_{∞}(Ω)}.

(5.8)  ∥\frac{∂^{4}v(·, μ)}{∂r^{4}} ∥_{L_{1}(Ω)} ≤ \frac{C}{r^{3}} ∥g∥_{W_{1}^{1}(Ω)},

Proof. In order to prove (5.4) we recall from (2.6) that v(x, μ) = T_{μ}g(x) is given by

(5.9)  v(x, μ) = ∫_{0}^{d/r} e^{−s}g(x − srμ̂) ds,

where d = d(x, μ), μ = rμ̂, μ̂ = (cos φ, sin φ). Making the transformation of variables σ = sr we find that

(5.10)  v(x, μ) = \frac{1}{r} ∫_{0}^{d} e^{−σ/r}g(x − σμ̂) dσ,

so that

\frac{∂v}{∂φ} = \frac{1}{r} \frac{∂d}{∂φ}e^{−d/r}g(x − dμ̂) − ∫_{0}^{d} \frac{σ}{r} e^{−σ/r} \frac{∂μ̂}{∂φ} · ∇g(x − σμ̂) dσ.

We now divide Ω into strips

Ω_{μ,i} = {x ∈ Ω : x − dμ̂ ∈ S_{i}},

where the S_{i} are the sides of Ω. In each Ω_{μ,i} we use a transformation of coordinates

x = P_{i} + ξ_{1}μ̂ + ξ_{2}τ̂_{i}, 0 ≤ ξ_{1} ≤ B_{i}(ξ_{2}), 0 ≤ ξ_{2} ≤ L_{i},

where P_{i} is an endpoint of S_{i} and τ̂_{i} = (cos ψ_{i}, sin ψ_{i}) is a tangent of S_{i}. Thus d = ξ_{1} and we compute that

\frac{∂d}{∂φ} = \frac{∂ξ_{1}}{∂φ} = −\frac{cos(φ − ψ_{i})}{sin(φ − ψ_{i})}ξ_{1}.

Since the area element is dx = | sin(φ − ψ_{i})| dξ we thus have

∥\frac{∂v(·, μ)}{∂φ}∥_{L_{1}(Ω_{μ,i})}

≤ ∫_{0}^{L_{i}} ∫_{0}^{B_{i}(ξ_{2})} |\frac{cos(φ − ψ_{i})}{sin(φ − ψ_{i})}| \frac{ξ_{1}}{r} e^{−ξ_{1}/r} |g(P_{i} + ξ_{2}τ̂_{i})| | sin(φ − ψ_{i})| dξ_{1}dξ_{2}

+ ∫_{0}^{L_{i}} ∫_{0}^{B_{i}(ξ_{2})} ∫_{0}^{ξ_{1}} \frac{σ}{r} e^{−σ/r} |\frac{∂μ̂}{∂φ}| |∇g(P_{i} + (ξ_{1} − σ)μ̂ + ξ_{2}τ̂_{i})| dσ | sin(φ − ψ_{i})| dξ_{1}dξ_{2}

≤ C (∥g∥_{L_{1}(S_{i})} + ∥g∥_{W_{1}^{1}(Ω_{μ,i})}) .

Using also the trace theorem and summing over i we obtain (5.4). The maximum norm bound (5.5) is proved in a similar way.

In order to prove (5.6) we differentiate v(x, μ) once more:

\frac{∂^{2}v}{∂φ^{2}} = [\frac{1}{r} \frac{∂^{2}ξ_{1}}{∂φ^{2}} − \frac{1}{r^{2}} (\frac{∂ξ_{1}}{∂φ} )^{2}] e^{−ξ_{1}/r}g(P_{i} + ξ_{2}τ̂_{i})

+ \frac{1}{r} \frac{∂ξ_{1}}{∂φ} e^{−ξ_{1}/r} \frac{∂ξ_{2}}{∂φ} τ̂_{i} · ∇g(P_{i} + ξ_{2}τ̂_{i})

− \frac{ξ_{1}}{r} \frac{∂ξ_{1}}{∂φ} e^{−ξ_{1}/r} \frac{∂μ̂}{∂φ} · ∇g(P_{i} + ξ_{2}τ̂_{i})

− ∫_{0}^{ξ_{1}} \frac{σ}{r} e^{−σ/r} \frac{∂^{2}μ̂}{∂φ^{2}} · ∇g(P_{i} + (ξ_{1} − σ)μ̂ + ξ_{2}τ̂_{i}) dσ

+ ∫_{0}^{ξ_{1}} \frac{σ^{2}}{r} e^{−σ/r} (\frac{∂μ̂}{∂φ} · ∇)^{2} g(P_{i} + (ξ_{1} − σ)μ̂ + ξ_{2}τ̂_{i}) dσ.

One easily computes that

\frac{∂ξ_{2}}{∂φ} = \frac{ξ_{1}}{sin(φ − ψ_{i})}, \frac{∂^{2}ξ_{1}}{∂φ^{2}} = \frac{1 + cos^{2}(φ − ψ_{i})}{sin^{2}(φ − ψ_{i})}ξ_{1}.

Hence

∥\frac{∂^{2}v(·, μ)}{∂φ^{2}} ∥_{L_{1}(Ω_{μ,i})} ≤ C ∫_{0}^{L_{i}} ∫_{0}^{B_{i}(ξ_{2})} \frac{(\frac{ξ_{1}}{r} + (\frac{ξ_{1}}{r})^{2}) e^{−ξ_{1}/r}}{| sin(φ − ψ_{i})|} |g(P_{i} + ξ_{2}τ̂_{i})| dξ_{1}dξ_{2}

+ C ∫_{0}^{L_{i}} ∫_{0}^{B_{i}(ξ_{2})} (\frac{\frac{ξ_{1}}{r} e^{−ξ_{1}/r}}{| sin(φ − ψ_{i})|} + \frac{ξ_{1}^{2}}{r} e^{−ξ_{1}/r}) |∇g(P_{i} + ξ_{2}τ̂_{i})| dξ_{1}dξ_{2}

+ C ∫_{0}^{L_{i}} ∫_{0}^{B_{i}(ξ_{2})} ∫_{0}^{ξ_{1}}(\frac{σ}{r} e^{−σ/r}|∇g(P_{i} + (ξ_{1} − σ)μ̂ + ξ_{2}τ̂_{i})|

+ \frac{σ^{2}}{r} e^{−σ/r}|D^{2}g(P_{i} + (ξ_{1} − σ)μ̂ + ξ_{2}τ̂_{i})|) dσ dξ_{1}dξ_{2}

≤ \frac{C}{| sin(φ − ψ_{i})|} (∥g∥_{L_{1}(S_{i})} + ∥∇g∥_{L_{1}(S_{i})}) + C∥g∥_{W_{1}^{2}(Ω_{μ,i})}

≤ \frac{C}{| sin(φ − ψ_{i})|}∥g∥_{W_{1}^{2}(Ω)} ≤ \frac{C}{min_{j} | sin(φ − ψ_{j})|}∥g∥_{W_{1}^{2}(Ω)},

where |D^{2}g| denotes the sum of the absolute values of all second order partial derivatives, and where we have used the trace theorem. Summation with respect to i now leads to (5.6).

We now turn to the proof of (5.8). We want to take four derivatives of v with respect to r, placing as few derivatives as possible on g. To achieve this we first differentiate (5.10) three times, using the fact that d is independent of r. With P (s) denoting various polynomials, not necessarily the same at each occurrence, we obtain

\frac{∂^{3}v}{∂r^{3}} = r^{−4} ∫_{0}^{d} P (σ/r)e^{−σ/r}g(x − σμ̂) dσ = r^{−3} ∫_{0}^{d/r} P (s)e^{−s} g(x − srμ̂) ds.

Differentiating once more we get

\begin{aligned}\frac{∂^{4}v}{∂r^{4}} & = r^{−4}P (\frac{d}{r}) e^{−d/r} g(x − dμ̂) + r^{−4} ∫_{0}^{d/r} P (s) e^{−s} g(x − srμ̂) ds \\ & + r^{−3} ∫_{0}^{d/r} P (s) e^{−s} μ̂ · ∇g(x − srμ̂) ds, \end{aligned}

or, after integration by parts in the middle term,

\frac{∂^{4}v}{∂r^{4}} = r^{−4}P (\frac{d}{r}) e^{−d/r} g(x − dμ̂) + r^{−3} ∫_{0}^{d/r} P (s) e^{−s} μ̂ · ∇g(x − srμ̂) ds ≡ I_{1} + I_{2}.

Here

\begin{aligned}∥I_{1}∥_{L_{1}(Ω_{μ,i})} & ≤ r^{−4} ∫_{0}^{L_{i}} ∫_{0}^{B_{i}(ξ_{2})} |P (\frac{ξ_{1}}{r} )| e^{−ξ_{1}/r} |g(P_{i} + ξ_{2}τ̂_{i})| dξ_{1} dξ_{2} \\ & ≤ r^{−4} ∫_{0}^{∞} |P (\frac{ξ_{1}}{r} )| e^{−ξ_{1}/r} dξ_{1} ∥g∥_{L_{1}(S_{i})} ≤ C r^{−3}∥g∥_{L_{1}(Γ)}, \end{aligned}

and

\begin{aligned}∥I_{2}∥_{L_{1}(Ω_{μ,i})} & ≤ r^{−3} ∫_{0}^{L_{i}} ∫_{0}^{B_{i}(ξ_{2})} ∫_{0}^{ξ_{1}/r} |P (ξ_{1})| e^{−ξ_{1}}|∇g(P_{i} + (ξ_{1} − s)μ̂ + ξ_{2}τ̂_{i})| ds dξ_{1} dξ_{2} \\ & ≤ C r^{−3}∥g∥_{W_{1}^{1}(Ω)}. \end{aligned}

Using also the trace theorem and summing over i we obtain (5.8). The maximum norm bound (5.7) is proved in a similar way. □

Lemma 5.4. (a) Let Q_{r} be defined by the Gauss rule (3.6). Then

∥∫_{0}^{1} v(·, μ)\frac{r}{\sqrt{1 − r^{2}}} dr − Q_{r}[v(·, μ)]∥_{L_{1}(Ω)} ≤ \frac{C}{N^{4}} ∥g∥_{W_{1}^{1}(Ω)},

∥∫_{0}^{1} v(·, μ)\frac{r}{\sqrt{1 − r^{2}}} dr − Q_{r}[v(·, μ)]∥_{L_{∞}(Ω)} ≤ \frac{C}{N^{2}} ∥g∥_{L_{∞}(Ω)},

independently of φ.

(b) Let Q_{φ} be defined by the trapezoidal rule (3.5). Then

∥∫_{0}^{2π} v(·, μ) dφ − Q_{φ}[v(·, μ)]∥_{L_{1}(Ω)} ≤ \frac{C}{M} ∥g∥_{W_{1}^{1}(Ω)},

∥∫_{0}^{2π} v(·, μ) dφ − Q_{φ}[v(·, μ)]∥_{L_{1}(Ω)} ≤ \frac{C}{M^{2}} log M ∥g∥_{W_{1}^{2}(Ω)},

∥∫_{0}^{2π} v(·, μ) dφ − Q_{φ}[v(·, μ)]∥_{L_{∞}(Ω)} ≤ \frac{C}{M} log M ∥g∥_{W_{∞}^{1}(Ω)},

independently of r.

Proof. Part (a) and the first estimate of part (b) are immediate consequences of Lemmas 5.2, 5.3, and Lemmas 5.1, 5.3, respectively. For the second estimate of part (b) we have to pay special attention to angles φ that are close to an angle of direction ψ_{i} of a side S_{i}. For this purpose we define intervals

I_{k} = [\frac{2π(k − 1)}{M}, \frac{2πk}{M} ] = [φ_{k−1}, φ_{k}], k = 1, ..., M,

and for each ψ_{i} we let J_{i} be the union of the interval I_{k} that contains ψ_{i} and the adjacent interval I_{k±1} closest to ψ_{i}. Further, we define

S_{0} = ∪_{i=1}^{P}J_{i},

where the union is taken over all the sides of Ω, and

S = [0, 2π], S_{1} = \overline{S \ S_{0}}.

Corresponding to the division ∫_{S} = ∫_{S_{0}} + ∫_{S_{1}} of the integral we divide the trapezoidal rule over S into the trapezoidal rule over S_{0} and S_{1}. For this we have the weights

W̃_{j} = \{ \begin{matrix}W_{j} = 2π/M, & if φ_{j} ∉ S_{0} ∩ S_{1}, \\ W_{j}/2 = π/M, & if φ_{j} ∈ S_{0} ∩ S_{1}. \end{matrix}

Using the standard error estimates for the trapezoidal rule from Lemma 5.1 together with Lemma 5.3 we now have

∥∫_{0}^{2π} v(·, μ) dφ − Q_{φ}[v(·, μ)]∥_{L_{1}(Ω)} ≤ ∥∫_{S_{0}} v(·, rμ̂(φ)) dφ − ∑_{φ_{j}∈S_{0}} W̃_{j}v(·, rμ̂(φ_{j}))∥_{L_{1}(Ω)}

+ ∥∫_{S_{1}} v(·, rμ̂(φ)) dφ − ∑_{φ_{j}∈S_{1}} W̃_{j}v(·, rμ̂(φ_{j}))∥_{L_{1}(Ω)}

≤ \frac{C}{M} ∫_{S_{0}} ∥\frac{∂v(·, rμ̂(φ))}{∂φ}∥_{L_{1}(Ω)} dφ + \frac{C}{M^{2}} ∫_{S_{1}} ∥\frac{∂^{2}v(·, rμ̂(φ))}{∂φ^{2}} ∥_{L_{1}(Ω)} dφ

≤ \frac{C}{M} ∫_{S_{0}} dφ ∥g∥_{W_{1}^{1}(Ω)} + \frac{C}{M^{2}} ∫_{S_{1}} \frac{1}{min_{i} | sin(φ − ψ_{i})|} dφ ∥g∥_{W_{1}^{2}(Ω)}

≤ \frac{C}{M^{2}} ∥g∥_{W_{1}^{1}(Ω)} + \frac{C}{M^{2}} log M ∥g∥_{W_{1}^{2}(Ω)} ≤ \frac{C}{M^{2}} log M ∥g∥_{W_{1}^{2}(Ω)}.

The maximum norm bound is proved in a similar way. □

Lemma 5.5. We have

∥(T − T_{n})g∥_{L_{1}(Ω)} ≤ C (\frac{1}{N^{4}} + \frac{1}{M} ) ∥g∥_{W_{1}^{1}(Ω)},

∥(T − T_{n})g∥_{L_{1}(Ω)} ≤ C \frac{1}{N^{4}} ∥g∥_{W_{1}^{1}(Ω)} + C \frac{1}{M^{2}} log M ∥g∥_{W_{1}^{2}(Ω)}, (5.11)

∥(T − T_{n})g∥_{L_{∞}(Ω)} ≤ C \frac{1}{N^{2}} ∥g∥_{L_{∞}(Ω)} + C \frac{1}{M} log M ∥g∥_{W_{∞}^{1}(Ω)}. (5.12)

Proof. This is an immediate consequence of the previous lemma and (5.2). □

We can now prove the boundedness of (I−λT_{n})^{−1}. We need the following simple lemma, the proof of which can be found in Johnson and Pitkäranta [8].

Lemma 5.6. Let T : L_{1}(Ω) → L_{1}(Ω) be a bounded linear operator such that for some C

∥(I − λT)g∥_{L_{1}(Ω)} ≥ C∥g∥_{L_{1}(Ω)} ∀g ∈ L_{1}(Ω),

and let {T_{n}}_{n=1}^{∞} be a bounded sequence of continuous linear operators on L_{1}(Ω) such that for some positive integer m

(5.13)  ∥(T − T_{n})T_{n}^{m}∥ → 0  as n → ∞.

Then there is C_{1} such that for n large

(5.14)  ∥(I − λT_{n})g∥_{L_{1}(Ω)} ≥ C_{1}∥g∥_{L_{1}(Ω)}  ∀g ∈ L_{1}(Ω).

Lemma 5.7. Assume that λ^{−1} is not an eigenvalue of T. Then the operator I − λT_{n} is invertible on L_{1}(Ω) for n̄ large and, moreover, its inverse is uniformly bounded.

Proof. We only sketch the proof, which is a simple modification of the corresponding argument in [1]. The idea is to show that (5.13) holds with m = 2. In order to do this one splits T_{n}^{2} into two terms:

(5.15)  T_{n}^{2} = ∑_{(μ,ν)∈Δ^{2}} ω_{μ}ω_{ν}T_{μ}T_{ν} = ∑_{(μ,ν)∈Δ_{ϵ}^{2}} ω_{μ}ω_{ν}T_{μ}T_{ν} + ∑_{(μ,ν)∈Δ^{2}\Δ_{ϵ}^{2}} ω_{μ}ω_{ν}T_{μ}T_{ν},

where the first sum on the right is over all velocities μ, ν ∈ Δ, (1) whose directions are bounded away from each other;

(2) whose directions are bounded away from the directions of the sides of Ω; and

(3) whose magnitudes are bounded away from 0; closeness being quantified by the positive parameter ϵ. One can then show that

(5.16)  ∥T_{μ}T_{ν}g∥_{W_{1}^{1}(Ω)} ≤ Cϵ^{−a}∥g∥_{L_{1}(Ω)},  ∀(μ, ν) ∈ Δ_{ϵ}^{2},

for some a > 0. The proof of this is based on direct estimates on the formula (5.9). For example one has

(5.17)  ∥μ · ∇T_{μ}g∥_{L_{1}(Ω)} + ∥T_{μ}g∥_{L_{1}(Ω)} + ∫_{Γ} |T_{μ}g| |μ · n̂| ds ≤ C∥g∥_{L_{1}(Ω)},

where n̂ is the outward unit normal to Γ. The negative power of ϵ in (5.16) comes from the fact that the operator T_{μ} is smoothing in the direction of μ but not in other directions.

Applying Lemma 5.5, using (5.16) and (5.17), we thus obtain

∥(T − T_{n})T_{n}^{2}g∥_{L_{1}(Ω)} ≤ ∥(T − T_{n}) ∑_{(μ,ν)∈Δ_{ϵ}^{2}} ω_{μ}ω_{ν}T_{μ}T_{ν}g∥_{L_{1}(Ω)}

+ ∥(T − T_{n}) ∑_{(μ,ν)∈Δ^{2}\Δ_{ϵ}^{2}} ω_{μ}ω_{ν}T_{μ}T_{ν}g∥_{L_{1}(Ω)}

≤ C (\frac{1}{N^{4}} + \frac{1}{M} ) ∑_{(μ,ν)∈Δ_{ϵ}^{2}} ω_{μ}ω_{ν}∥T_{μ}T_{ν}g∥_{W_{1}^{1}(Ω)}

+ C ∑_{(μ,ν)∈Δ^{2}\Δ_{ϵ}^{2}} ω_{μ}ω_{ν}∥g∥_{L_{1}(Ω)}

≤ C ((\frac{1}{N^{4}} + \frac{1}{M} ) ϵ^{−a} + ∑_{(μ,ν)∈Δ^{2}\Δ_{ϵ}^{2}} ω_{μ}ω_{ν} ) ∥g∥_{L_{1}(Ω)}.

Choosing ϵ appropriately and using the fact that ∑_{(μ,ν)∈Δ^{2}\Δ_{ϵ}^{2}} ω_{μ}ω_{ν} → 0 as ϵ → 0, we find that

(5.18)  ∥(T − T_{n})T_{n}^{2}∥ → 0  as n̄ → ∞.

This proves (5.14), and hence that I − λT_{n} is one-to-one. Since T_{n} is not compact we may not conclude directly that I − λT_{n} is onto, but this can be shown by using (5.14) together with the fact that (5.15) splits T_{n}^{2} as T_{n}^{2} = A_{n} + B_{n}, where A_{n} : L_{1}(Ω) → W_{1}^{1}(Ω), i.e., A_{n} is compact, and ∥B_{n}∥ → 0 as n̄ → ∞. We refer to [8] for the details. □

We can now prove Theorem 3.1. Interpolating between (5.11) and the bound

∥(T − T_{n})g∥_{L_{1}(Ω)} ≤ C∥g∥_{L_{1}(Ω)},

obtained by stability, cf. (5.17), we get

(5.19)  ∥(T − T_{n})g∥_{L_{1}(Ω)} ≤ C \frac{1}{N^{4}} ∥g∥_{W_{1}^{1}(Ω)} + C \frac{1}{M^{2−ϵ}} ∥g∥_{(L_{1}(Ω),W_{1}^{2}(Ω))_{θ,q}},

for 0 ≤ 1 − ϵ/2 < θ ≤ 1. Applying this with g = λU + f, and recalling Theorem 4.1, we have in view of (5.1),

∥U − U_{n}∥_{L_{1}(Ω)} ≤ C (\frac{1}{N^{4}} + \frac{1}{M^{2−ϵ}} ) (∥U∥_{L_{1}(Ω)} + ∥f∥_{W_{1}^{1}(Ω)}) .

This completes the proof of Theorem 3.1.

6. An error estimate for the critical eigenvalue. The kernel of the compact integral operator T is symmetric and positive, see (4.1). Hence T is self-adjoint (on L_{2}(Ω)), and therefore has only real eigenvalues. Moreover, by the Krein-Rutman theory, its largest eigenvalue is positive and simple. Without loss of generality we may assume that M is even so that the set of quadrature points Δ is symmetric, i.e., μ ∈ Δ implies −μ ∈ Δ. Then it follows that T_{n} is self-adjoint and its eigenvalues are real.

Theorem 6.1. Assume that M is even and let κ and κ_{n} be the largest eigenvalues of T and T_{n}, respectively. For any ϵ > 0 there is C = C(ϵ, κ) such that, for n̄ sufficiently large, we have

|κ − κ_{n}| ≤ C (\frac{1}{N^{4}} + \frac{1}{M^{2−ϵ}} ) .

Proof. We shall show that

(6.1)  ∥T^{3} − T_{n}^{3}∥ → 0,  as n̄ → ∞,

i.e., T_{n}^{3} converges to T^{3} in the norm of bounded linear operators on L_{1}(Ω). It then follows that the family {T_{n}^{3}} is collectively compact, see [6, p. 251], and a result of Atkinson [4] yields

|κ^{3} − κ_{n}^{3}| ≤ C∥(T^{3} − T_{n}^{3})χ∥_{L_{1}(Ω)},

where χ is the eigenfunction of T associated with κ, normalized by ∥χ∥_{L_{1}(Ω)} = 1. Using the boundedness of T and T_{n}, the error bound (5.19), and the regularity estimate from Theorem 4.1, we thus have

\begin{aligned}|κ^{3} − κ_{n}^{3}| & ≤ C (\frac{1}{N^{4}} + \frac{1}{M^{2−ϵ}} ) ∥χ∥_{(L_{1}(Ω),W_{1}^{2}(Ω))_{θ,q}} \\ & ≤ C (\frac{1}{N^{4}} + \frac{1}{M^{2−ϵ}} ) ∥χ∥_{L_{1}(Ω)} = C (\frac{1}{N^{4}} + \frac{1}{M^{2−ϵ}} ) , \end{aligned}

which implies the desired result, since κ > 0.

It remains to prove (6.1). We first note that

T^{3} − T_{n}^{3} = T^{2}(T − T_{n}) + T (T − T_{n})T_{n} + (T − T_{n})T_{n}^{2}.

Here, by (5.18), the last term tends to zero in the operator norm. For the remaining terms we note that, with (·, ·) denoting the duality pairing between L_{1}(Ω) and L_{∞}(Ω), and using the symmetry of T and T_{n},

|(T (T − T_{n})f, g)| = |(f, (T − T_{n})T g)| ≤ ∥f∥_{L_{1}(Ω)}∥(T − T_{n})T g∥_{L_{∞}(Ω)}.

By interpolation, (5.12), stability in L_{∞} and Theorem 4.2 we get

∥(T − T_{n})T g∥_{L_{∞}(Ω)} ≤ C (\frac{1}{N^{2}} + \frac{1}{M} log M)^{θ} ∥T g∥_{(L_{∞}(Ω),W_{∞}^{1}(Ω))_{θ,q}}

≤ C (\frac{1}{N^{2}} + \frac{1}{M} log M)^{θ} ∥g∥_{L_{∞}(Ω)},

with 0 < θ < 1. Hence

∥T (T − T_{n})∥ ≤ C (\frac{1}{N^{2}} + \frac{1}{M} log M)^{θ} ,

and (6.1) follows. □

References

1. M. Asadzadeh, Analysis of a fully discrete scheme for neutron transport in two-dimensional geometry, SIAM J. Numer. Anal. 23 (1986), 543–561.

2. M. Asadzadeh, L_{p} and eigenvalue error estimates for two-dimensional neutron transport, SIAM J. Numer. Anal. 26 (1989), 66–87.

3. M. Asadzadeh, L_{2} error estimates for discrete ordinates method for three-dimensional neutron transport, Transport Theory and Statistical Physics 17 (1988), 1–24.

4. K. E. Atkinson, Convergence rates for approximate eigenvalues of compact integral operators, SIAM J. Numer. Anal. 12 (1975), 213–222.

5. J. Bergh and J. Löfström, Interpolation Spaces, Springer-Verlag, 1976.

6. F. Chatelin, Spectral Approximation of Linear Operators, Academic Press, 1983.

7. R. A. De Vore and L. R. Scott, Error bounds for Gaussian quadrature and weighted-L_{1} polynomial approximation, SIAM J. Numer. Anal. 21 (1984), 400–412.

8. C. Johnson and J. Pitkäranta, Convergence of a fully discrete scheme for two-dimensional neutron transport, SIAM J. Numer. Anal. 20 (1983), 951–966.

9. S. M. Nikolskii, Approximation of Functions of Several Variables and Imbedding Theorems, Springer-Verlag, 1975.

10. J. Pitkäranta, Estimates for the derivatives of solutions to weakly singular Fredholm integral equations, SIAM J. Math. Anal. 11 (1980), 952–968.

11. J. Pitkäranta and L. R. Scott, Error estimates for the combined spatial and angular approximations of the transport equation in slab geometry, SIAM J. Numer. Anal. 20 (1983), 922–950.

12. E. M. Stein, Singular Integrals and Differentiability Properties of Functions, Princeton University Press, 1970.

13. G. Szegő, Orthogonal Polynomials, AMS Colloquium Publications 23, American Mathematical Society, 1957.

Department of Mathematics, Chalmers University of Technology, S–412 96 Göteborg, Sweden

E-mail address: mohammad@math.chalmers.se, kumlin@math.chalmers.se, stig@math.chalmers.se