# A Low Memory Solver for Integral Equations of Chandrasekhar Type in the Radiative Transfer Problems

Authors: Mohammed Yusuf Waziri; Wah June Leong; Malik Abu Hassan; Mansor Monsi
Journal: Mathematical Problems in Engineering
 (2011)

---

## Abstract

The problems of radiative transfer give rise to interesting integral equations
that must be faced with efficient numerical solver. Very often the integral equations are discretized to large-scale nonlinear equations and solved by Newton's-like methods. Generally, these kind of methods require the computation and storage
of the Jacobian matrix or its approximation. In this paper, we present a new approach that was based on approximating the Jacobian inverse into a diagonal matrix by means of variational technique. Numerical results on well-known benchmarks integral equations involved in the radiative transfer authenticate the reliability and efficiency of the approach. The fact that the proposed method can solve the integral equations without function derivative and matrix storage can be considered as a clear advantage over some other variants of Newton's method.


---

## Body

## 1. Introduction

The study of Chandrasekhar’s integral equation involved in the radiative transfer problem has been a foremost subject of much investigations and was first formulated by Chandrasekhar [1] in 1960. It arose originally in connection with scattering through a homogeneous semi-infinite plane atmosphere and since it has been used to model diverse forms of scattering via the H-function of Chandrasekhar [2], defined by
(1.1)H(x)=1+H(x)∫01xx+tψ(t)H(t)dt.
Chandrasekhar H-function plays a crucial role in radiative transfer and transport theory [3, 4]. Since then, there have been diverse solvers of (1.1). It is well known that the numerical solution of Chandrasekhar integral equation is difficult to be obtained [5], and thus it is convenient to have a reliable and efficient solver. The problem of finding approximate solution of such integral equations is still popular today, and various methods of solving these integral equations have been established [5–7]. The common approach for the approximate solution of (1.1) is at first discretizing (1.1) by a vector x¯∈Rn, then replacing the integrals by quadrature sums and the derivatives by difference quotients involving only the components of x¯∈Rn (see [8], e.g.). By doing so, (1.1) becomes a problem of finding the solution of system of n nonlinear equations with n unknowns
(1.2)F(x¯)=0,
where F:Rn→Rn is a nonlinear mapping. Often, the mapping F is assumed to satisfy the following assumptions:(A1)
there existsx¯*∈Rn s.t. F(x¯*)=0,(A2)
F is a continuously differentiable mapping in an open neighborhood of x¯*,(A3)
F′(x¯*) is invertible.The famous iterative method for solving (1.2) is the classical Newton’s method, where the Newtonian iteration is given by
(1.3)xk+1=xk-(F′(xk))-1F(xk),k=0,1,2,….
The convergence rate for the Newton’s method is quadratic from any initial point x0 in the neighborhood of x¯* [9]. However, an iteration of (1.3) turns to be expensive, because it requires to compute and store the Jacobian matrix, as well as solving Newton’s system which is a linear system in each iteration. The major difficulty of Newton’s type method is the matrix storage requirements especially when handling large systems of nonlinear equations [5, 6, 9]. There are quite a number of revised Newton’s type methods being introduced, which include fixed Newton’s and quasi-Newton’s, to diminish the weakness of (1.3). Fixed Newton method [10] for the determination of solution x* is given by
(1.4)xk+1=xk-(F′(x0))-1F(xk),k=0,1,2,….
The method avoids computation and storing the Jacobian in each iteration (except at k=0). However, it still requires solving the systems of n linear equations and may consume more CPU time as the system’s dimension increases [10].A quasi-Newton’s method is another variant of Newton’s type methods, and it replaces the Jacobian or its inverse with an approximation which can be updated at each iteration [11], and its updating scheme is given by
(1.5)xk+1=xk-Bk-1F(xk),
where the matrix Bk is the approximation of the Jacobian at xk. The main idea behind quasi-Newton’s method is to eliminate the evaluation cost of the Jacobian matrix, in which if function evaluations are very expensive, the cost of finding a solution by quasi-Newton’s methods could be much smaller than some other Newton’s-like methods [7, 12, 13]. Various Jacobian approximations matrices such as the Broyden’s method [11, 14] are proposed. However, the most critical part of such solvers is that they need the storage of full matrix of the approximate Jacobian, which can be a very expensive task as the dimension of systems increases [15]. In this paper, we propose an alternative approximation to the Jacobian inverse into a diagonal matrix by means of variational techniques. It is worth mentioning that the suggested method can be applied to solve (1.2) without the cost of computing or storing the true Jacobian. Hence, it can reduce computational cost, storage requirement, processing time (CPU time) and also eliminates the need for solving n linear equations in each iteration. The proposed method works efficiently, and the results so far are very encouraging. This paper is arranged as follows; we present our proposed method in Section 2; numerical results are reported in Section 3; finally conclusion is given in Section 4.


## 2. ChandrasekharH-Equation

In this section, we present the detailed process of discretizing the Chandrasekhar-type integral equations in the radiative transfer problem. Chandrasekhar and Breen [16] compute H-equation as the solution of the nonlinear integral equation
(2.1)H(x)-cx2H(x)∫01H(y)x+ydy=1,
where c∈[0,1] and H:[0,1]→ℝ is an unknown continuous function. From (2.1), we obtain
(2.2)H(x)[1-c2∫01xH(y)x+ydy]=1.
Let us partition [0,1] into n subinterval, 0<x1<⋯<xj=j/n<⋯<1. Denote Hk as H(xk), then the evaluation of (2.1) at every xi yields the equation
(2.3)Hi[1-c2∫01xiH(y)xi+ydy]=1,i=1,2,…,n.
After multiplying both sides of (2.2) by [1-(c/2)∫01(xH(y)/(x+y))dy]-1 and performing some algebra, we arrive at (2.4), which is known as the Chandrasekhar H-equation [15]
(2.4)F(H)(x)=H(x)-(1-c2∫01xH(y)dyx+y)-1=0.
If (2.4) is discretized by using the midpoint quadrature formula
(2.5)∫01f(t)dt=1n∑j=0nf(tj),
for tj=(j-0.5)h, 0≤j≤1, i=2,…,n, h=1/n, c∈(0,1), then we have the following:
(2.6)Fi=xi-(1-c2n∑j=1ntixjti+tj)-1.
Function (2.6) is called the discretized Chandrasekhar H-equation which can be solved by some iterative methods.Nevertheless, the most difficult part in solving (2.6) arises dramatically as c approaches 1, since its Jacobian is singular at c=1. Due to this disadvantage, we aim to derive a method that hopefully will not be affected by this difficulty.


## 3. Derivation of the Method (LMSI)

Firstly, note that by the mean value theorem, we have(3.1)F′¯(xk)(xk+1-xk)=F(xk+1)-F(xk),
where F′¯(xk)=∫01F′(xk+θ(xk+1-xk))dθ.Let us denoteΔxk=xk+1-xk and ΔFk=F(xk+1)-F(xk), then (3.1) becomes
(3.2)F′¯(xk)Δxk=ΔFk.
Equation (3.2) is always regarded as the secant equation. Alternatively, we can rearrange (3.2) to obtain
(3.3)Δxk=(F′¯(xk))-1ΔFk.
Here, we propose to use a diagonal matrix, say D, to approximate (F′¯(xk))-1, that is,
(3.4)(F′¯(xk))-1≈Dk.
Let us consider an updating scheme for D, in which we should update D by adding a correction M which is also a diagonal matrix at every iteration
(3.5)Dk+1=Dk+Mk.
In order to incorporate correct information of the Jacobian inverse into the updating matrix, Dk+1, we require that Dk+1 satisfies the secant equation (3.2), that is,
(3.6)Δxk=(Dk+Mk)ΔFk.
However, since it is difficult to have a diagonal matrix that satisfies the secant equation, in particular, because Jacobian approximations are not usually done in element wise, we consider the use of the weak secant condition [17] instead,
(3.7)ΔFkTΔxk=ΔFkT(Dk+Mk)ΔFk.
To encourage good condition number as well as numerical stability in the approximation, we attempt to control the growth error of the correction by minimizing its magnitude under some norms (here, we consider the Frobenuis norm), such that (3.7) holds. To this end, we consider the following problem:
(3.8)min⁡12‖Dk+1-Dk‖F2s.t.ΔFkTDk+1ΔFk=ΔFkTΔxk,
where ∥·∥F is the Frobenuis norm. If we let Dk+1-Dk=Mk=diag(β1,β2,…,βn) and ΔFk=(ΔFk(1),ΔFk(2),…,ΔFk(n)), the above problem can be expressed as follows:
(3.9)min⁡12(β12+β22+⋯+βn2)s.t.∑i=1nΔFk(i)2βi-ΔFkTΔxk+ΔFkTDkΔFk=0.
Since the objective function and the constraint are convex, we will have unique solution for (3.9). The solution can be obtained by considering the Lagrangian function of problem (3.9)
(3.10)L(βi,λ)=12(β12+β22+⋯+βn2)+λ(∑i=1nΔFk(i)2βi-ΔFkTΔxk+ΔFkTDkΔFk),
where λ is the corresponding Lagrangian multiplier.Taking the partial derivatives of (3.10) with respect to each βi and λ, respectively, and setting them equal to zero, we have
(3.11)∂L∂βi=βi+λ(ΔFk(i))2=0,i=0,1,2,…,n,(3.12)∂L∂λ=∑i=1n(ΔFk(i))2βi-ΔFkTΔxk+ΔFkTDkΔFk=0.
Premultiplying both sides of (3.11) by ΔFk(i)2 and summing them all yield
(3.13)∑i=1n(ΔFk(i))2βi+λ∑i=1n(ΔFk(i))4=0.
It follows from (3.13) that
(3.14)∑i=1n(ΔFk(i))2βi=-λ∑i=1n(ΔFk(i))4.
Invoking the constraint (3.12), we have
(3.15)∑i=1n(ΔFk(i))2βi=ΔFkTΔxk-ΔFkTDkΔFk.
Equating (3.14) with (3.15) gives
(3.16)λ=-ΔFkTΔxk-ΔFkTDkΔFk∑i=1n(ΔFk(i))4.
Substituting (3.16) into (3.14) and after some simplifications, we obtain
(3.17)βi=(ΔFkTΔxk-ΔFkTDkΔFk)∑i=1n(ΔFk(i))4(ΔFk(i))2,i=1,2,…,n.
Denoting Gk=diag((ΔFk(1))2,(ΔFk(2))2,…,(ΔFk(n))2) and ∑i=1n(ΔFk(i))4=Tr(Gk2) where Tr is the trace operation, we obtain, therefore,
(3.18)Mk=(ΔFkTΔxk-ΔFkTDkΔFk)Tr⁡(Gk2)Gk.
Finally, the proposed updating formula for the approximation of the Jacobian inverse is given as follows:
(3.19)Dk+1=Dk+(ΔFkTΔxk-ΔFkTDkΔFk)Tr⁡(Gk2)Gk.
To safeguard possibly very small ΔFk and Tr(Gk2), we require that ∥ΔFk∥≥ϵ1 for some chosen small ϵ1>0. Else, we will skip the update by setting Dk+1=Dk.Now, we can describe the algorithm for our proposed method (LMSI) as follows.Algorithm LMSI
Steps are the following.
Step 1.
Givenx0 and D0, set k=0.
Step 2.
ComputeF(xk) and xk+1=xk-DkF(xk).
Step 3.
If∥Δxk∥2+∥F(xk)∥2≤10-4, stop. Else, go to Step 4.
Step 4.
If∥ΔFk∥2≥ϵ1 where ϵ1=10-4, compute Dk+1, if not, Dk+1=Dk. Set k:=k+1 and go to Step 2.


## 4. Local Convergence Results

In this section, we will give some convergence properties of LMSI method. Before we proceed further, we will make the following standard assumptions on nonlinear systemsF.Assumption 4.1.
We have the following.(i)
F is differentiable in an open-convex set E in ℜn.(ii)
There existsx*∈E such that F(x*)=0, and F′(x) is continuous for all x.(iii)
F′(x) satisfies Lipschitz condition of order one, that is, there exists a positive constant μ such that
(4.1)‖F′(x)-F′(y)‖≤μ‖x-y‖,
for all x,y∈ℜn.(iv)
There exists constantsc1≤c2 such that c1∥ω∥2≤ωTF′(x)ω≤c2∥ω∥2 for all x∈E and ω∈ℜn.We will also need the following result which is a special case of a more general theorem of [15].Theorem 4.2.
Assume that Assumption4.1 holds. If there exists KB>0, δ>0, and δ1>0, such that for x0∈B(δ) and the matrix-valued function B(x) satisfies ∥I-B(x)F′(x*)∥=ρ(x)<δ1 for all x∈B(δ), then the iteration
(4.2)xk+1=xk-B(xk)F(xk)
converges linearly to x*.For the proof of Theorem4.2, see [15].Using Assumption4.1 and Theorem 4.2, one has the following result.Theorem 4.3.
Assume that Assumption4.1 holds. There exist β>0, δ>0, α>0, and γ>0, such that if x0∈E and D0 satisfies ∥I-D0F′(x*)∥F<δ for all xk∈E, then for iteration
(4.3)xk+1=xk-DkF(xk),Dk defined by (3.19),
(4.4)‖I-DkF′(x*)‖F<δk,
holds for some constant δk>0, k≥0.Proof.
Since∥Dk+1∥F=∥Dk+Mk∥F, it follows that
(4.5)‖Dk+1‖F≤‖Dk‖F+‖Mk‖F.
For k=0 and assuming D0=I, we have
(4.6)|M0(i)|=|ΔF0TΔx0-ΔF0TD0ΔF0Tr⁡(G2)ΔF0(i)2|≤|ΔF0TΔx0-ΔF0TD0ΔF0|Tr⁡(G02)ΔF0(max⁡)2,
where (ΔF0(max))2 is the largest element among (ΔF0(i)), i=1,2,…,n. After multiplying (4.6) by (ΔF0(max))2/(ΔF0(max))2 and substituting Tr(G02)=∑i=1n(ΔF0(i))4, we have
(4.7)|M0(i)|≤|ΔF0TΔx0-ΔF0TD0ΔF0|(ΔF0(max⁡))2∑i=1n(ΔF(i))4(ΔF0(max⁡))4.
Since (ΔF0(max))4/∑i=1n(ΔF0(i))4≤1, then (4.7) turns into
(4.8)|M0(i)|≤|ΔF0TF′(x)ΔF0-ΔF0TD0ΔF0|(ΔF0(max⁡))2.
From Assumption 4.1(iv) and D0=I, (4.8) becomes
(4.9)|M0(i)|≤|c-1|(ΔF0TΔF0)(ΔF0(max⁡))2,
where c=max{|c1|,|c2|}.
SinceΔ(F0(i))2≤(ΔF0(max))2 for i=1,…,n, it follows that(4.10)|M0(i)|≤|c-1|n(ΔF0(max⁡))2(ΔF0(max⁡))2.
Hence, we obtain
(4.11)‖M0‖F≤n3/2|c-1|.
Suppose that α=n3/2|c-1|, then
(4.12)‖M0‖F≤α.
From the fact that ∥D0∥F=n, it follows that
(4.13)‖D1‖F≤β,
where β=n+α>0.
Therefore, if we assume that∥I-D0F′(x*)∥F<δ, then(4.14)‖I-D1F′(x*)‖F=‖I-(D0+M0)F′(x*)‖F≤‖I-D0F′(x*)‖F+‖M0F′(x*)‖F≤‖I-D0F′(x*)‖F+‖M0‖F‖F′(x*)‖F,
hence ∥I-D1F′(x*)∥F<δ+αϕ=δ1.
And hence, by induction,∥I-DkF′(x*)∥F<δk for all k.


## 5. Numerical Results

In this section, we compare the performance of LMSI method with that of the Newton’s method (NM), fixed Newton’s method (FN), and Broyden’s method (BM). We apply the algorithms to the well-known benchmarks integral equations involved in radiative transfer. The comparison is based upon the following criterion: number of iterations, CPU time in seconds, and storage requirement. The computations are done in MATLAB 7.0 using double-precision computer. The stopping criterion used is(5.1)‖Δxk‖+‖F(xk)‖≤10-4.
The starting point x0 is given by (1,1,…1)T.The symbol “−” is used to indicate a failure due to the following:(1)
The number of iteration is at least 200, but no point ofxk satisfying (5.1) is obtained,(2)
CPU time in second reaches 200,(3)
insufficient memory to initial the run.The numerical results of the methods when solving ChandrasekharH-Equation in different parameter are reported in Table 1. The first column of the table contains the parameter of problem. Generally, with our choice of c, the corresponding Jacobian is not diagonally dominate; however, when c→1, the Jacobian is nearly singular. From Table 1, it was shown that only LMSI is able to solve problems where n>2000. This is due to the fact that LMSI requires very low-storage requirement in building the approximation of the Jacobian inverse. Indeed, the size of the updating matrix increases in O(n) as the dimension of the system increases, as opposed to NM, FN, and BM methods that increase in O(n2).Table 1
Results of ChandrasekharH-equation (number of iteration/CPU time).
nNMFNBMLMSI2004/6.91347/5.16736/4.71125/0.03125004/14.09457/11.04188/6.98135/0.062410004/64.58139/27.29748/9.74165/0.1716c=0.920004/95.61809/118.20948/11.16905/1.64105000———5/2.355610000———5/2.820920000———5/3.79152006/8.18348/6.42827/5.19466/0.07805006/18.24578/14.61327/7.61376/0.082410006/83.05699/30.05428/11.04326/0.2184c=0.9920006/121.53099/153.03518/15.97244/1.31045000———4/1.497610000———5/1.791020000———5/2.03712008/15.1792—9/5.03885/0.33125008/36.0330—10/7.14125/0.602410008/134.9045—10/12.06446/0.7803c=0.999920008/175.0521—10/14.90646/1.25375000———6/1.513210000———6/1.580820000———6/1.6301Moreover, we observe that LMSI method has a 100% of success rate (convergence to the solution) when compared with NM method having 57%, FN method with 39% and BM with 71%, respectively. In addition, it is worth mentioning that the result of LMSI in solving problem 1 whenc=0.9999 shows that the method could be a good solver even when the Jacobian is nearly singular. Figures 1, 2, and 3 reveal that the CPU time of LMSI method increases linearly as the dimension of the systems increases, whereas for NM, FN, and BM, the rates grow exponentially. This also suggests that our solver is a good alternative when the dimension of the problem is very high.Figure 1
Comparison of NM, FN, BM, and LMSI methods whenc=0.9, in terms of CPU time.Figure 2
Comparison of NM, FN, BM, and LMSI methods whenc=0.99, in terms of CPU time.Figure 3
Comparison of NM, FN, BM, and LMSI methods whenc=0.9999, in terms of CPU time.


## 6. Conclusion

In this paper, we present a low memory solver for integral equation of Chandrasekhar type in the radiative transfer problems. Our approach is based on approximating the Jacobian inverse into a diagonal matrix. The fact that the LMSI method can solve the discretized integral equations without computing and storing the Jacobian makes clear the advantage over NM and FN methods. It is also worth mentioning that the method is capable of significantly reducing the execution time (CPU time), as compared to NM, FN, and BM methods while maintaining good accuracy of the numerical solution to some extend. Another fact that makes the LMSI method appealing is that throughout the numerical experiments it never fails to converge. Finally, we conclude that our method (LMSI) is a good alternative to Newton-type methods for solving large-scale nonlinear equations with nearly singular Jacobian.



---

*Source: 467017-2011-11-10.xml*