On an inverse source problem for the full radiative transfer equation with incomplete data
OON AN INVERSE SOURCE PROBLEM FOR THE FULL RADIATIVETRANSFER EQUATION WITH INCOMPLETE DATA ∗ ALEXEY V. SMIRNOV † , MICHAEL V. KLIBANOV † , AND
LOC H. NGUYEN † Abstract.
A new numerical method to solve an inverse source problem for the radiative transferequation involving the absorption and scattering terms, with incomplete data, is proposed. Norestrictive assumption on those absorption and scattering coefficients is imposed. The original inversesource problem is reduced to boundary value problem for a system of coupled partial differentialequations of the first order. The unknown source function is not a part of this system. Next, wewrite this system in the fully discrete form of finite differences. That discrete problem is solved viathe quasi-reversibility method. We prove the existence and uniqueness of the regularized solution.Especially, we prove the convergence of regularized solutions to the exact one as the noise level inthe data tends to zero via a new discrete Carleman estimate. Numerical simulations demonstrategood performance of this method even when the data is highly noisy.
Key words. radiative transfer equation, absorption term, scattering term, inverse source prob-lem, discrete Carleman estimate, quasi-reversibility method
AMS subject classifications.
1. Introduction.
The stationary radiative transfer equation (RTE) is com-monly used in optics, tomography, astrophysics, atmospheric science and remotesensing to describe the propagation of the radiation field in media with absorbing,emitting and scattering radiation. A significant number of studies is dedicated tothe recovery of the parameters of the observed objects from the measured data; i.e.,to the solutions of the inverse source problems (ISOPs) [1, 24] and coefficient inverseproblems (CIPs) [2, 33]. A number of inverse problems may be formulated, dependingon the object’s parameters of one’s interest.The first reconstruction formula for the problem of the attenuated tomographywas obtained by Novikov [31]. We also refer to [3, 12, 29] for reconstruction formulaeand as well as to [10, 29] for numerical results for the attenuated tomography withcomplete data and with the scattering phase function K ≡
0. Uniqueness and sta-bility results for similar ISOPs with complete data were obtained in [3, 34]. It wasassumed in [3] that | K | is sufficiently small. The assumption of [34] is that functions σ and K belong to certain dense sets of some function spaces. The scattering phasefunction K is involved in RTE as the kernel of a certain integral operator, the atten-uation coefficient is σ = µ a + µ s , where µ a and µ s are the absorption and scatteringcoefficients respectively, see Section 2.In this paper, we propose a new numerical approach for the ISOP with limitedangle data for the stationary RTE and prove its convergence. This is the first publica-tion, in which a rigorously derived numerical method for the ISOP for the RTE doesnot use any restrictive assumptions neither on µ a , nor on µ s , nor on K , except thesmoothness and the requirement that functions µ a and µ s are compactly supported.Also, for the first time, a discrete Carleman estimate is applied here for the conver-gence analysis of an inverse problem. We note that discrete Carleman estimates are ∗ Submitted to the editors DATE.
Funding:
This work was supported by US Army Research Laboratory and US Army ResearchOffice grant W911NF-19-1-0044. In addition, the effort of Nguyen was supported by research fundsno. FRG 111172 provided by The University of North Carolina at Charlotte. † Department of Mathematics and Statistics, University of North Carolina Charlotte, Charlotte,NC, 28223 ([email protected], [email protected] (corresponding author), [email protected]).1 a r X i v : . [ m a t h . NA ] A p r ALEXEY V. SMIRNOV, MICHAEL V. KLIBANOV AND LOC H. NGUYEN very rare, unlike the continuous ones. In addition, we prove the Lipschitz stabilityand uniqueness for our statement of the ISOP.Our method is based on the solution of an overdetermined boundary value prob-lem for a linear system of coupled integro-differential equations, in which the unknownsource function is not present. The solution of this problem directly yields the solu-tion of the desired ISOP. A similar idea was recently used in [23]. However, unlikethe current paper, a quite restrictive condition σ ≡ K ≡ L ( a, b ) , ( a, b ) ⊂ R , which was recently introduced in [20]. Thisbasis has proven to be effective for numerical studies [21, 22, 23]. We use a truncatedFourier series with respect to this basis. We estimate an optimal number of terms ofthis series numerically and assume that this approximation still satisfies the RTE, i.e.we work with an approximate mathematical model, also, see Remark 4 at the end ofSection 5.We solve the above mentioned overdetermined boundary value problem by thequasi-reversibility method (QRM), which is known to be effective to solve overdeter-mined boundary value problems. We consider a fully discrete form of our system,which is similar to what we use in the numerical tests. Next, we establish a newdiscrete Carleman estimate and use it to prove uniqueness and existence of the reg-ularized solution for the QRM in the fully discrete form, in which partial derivativeswith respect to spatial variables are written via finite differences. This Carleman es-timate is also used to establish the convergence rate of regularized solutions. Finally,we conduct numerical testing for several different regimes of absorption and scatteringto show the method’s potential for solving problems in real-world tomography.The QRM was originally introduced by Lattes and Lions in 1969 [25]. We alsorefer to, e.g. [7, 8, 12, 23] for this method. The second author has shown in the surveypaper [19] that as long as a proper Carleman estimate for an ill-posed problem for alinear PDE is available, the convergent QRM can be constructed for this problem.For brevity, we consider in this paper only the 2D case. The considerations inthe 3D case are similar. We state both forward and inverse problems in Section 2.In Section 3 we derive the above mentioned over-determined boundary value problemfor a system of coupled partial differential equations of the first order. To solve thisproblem, we apply the QRM by stating a Minimization Problem. In Section 4 weintroduce the fully discrete version of the quasi-reversibility method to solve thatproblem. Next, we derive a new discrete Carleman estimate. This estimate is usedin Section 5 to prove the existence and uniqueness of the minimizer of the QRM andalso to establish the convergence rate of the minimizers to the exact solution as thelevel of noise in the measured data tends to zero. Section 6 is devoted to numericalstudies. Everywhere below we work only with real-valued functions. NVERSE PROBLEM FOR RADIATIVE TRANSFER EQUATION (a) The source/detector configuration of theproblem in the case when the source locatedat x α with | x α | < R. (b) The source/detector configuration of theproblem in the case when the source locatedat x α with | x α | > R. Fig. 2.1:
A schematic diagram of measurements for the 2D case. L ( x , x α ) is a straightline, connecting the detector x with the source x α .
2. Statements of Forward and Inverse Problems.
Let x = ( x, y ) denotean arbitrary point in R . Let a, b, d and R be the positive numbers, where 1 < a < b and d ≥ R . Define the rectangular domain Ω ⊂ R (Figure 2.1) as(2.1) Ω = { ( x, y ) : − R < x < R, a < y < b } . Let Γ d be the line with external sources(2.2) Γ d = { x α = ( α,
0) : α ∈ [ − d, d ] } . Let u ( x , α ) denotes the steady-state radiance at the point x generated by the externalsource located at x α = ( α, ∈ Γ d . Then, the function u ( x , α ) satisfies the followingradiative transfer equation, see, e.g. [11](2.3) ν ( x , α ) · ∇ x u ( x , α ) + ( µ a ( x ) + µ s ( x )) u ( x , α )= µ s ( x ) (cid:90) Γ d K ( x , α, β ) u ( x , β ) dβ + f ( x ) for all x ∈ Ω . In the equation above, the function f ( x ) ∈ L ( R ) is called the source function whilethe functions µ a ( x ), µ s ( x ) ∈ C (cid:0) R (cid:1) denote the absorption and scattering coefficientsrespectively. We assume that(2.4) µ a ( x ) = µ s ( x ) = f ( x ) = 0 for all x ∈ R \ Ω . The function K ( x , α, β ) ∈ C ( R × [ − d, d ] ) represents the so-called “scattering phasefunction”. Scattering phase function is the probability density of a particle scatteringfrom ν ( x , β )-direction into ν ( x , α )-direction. As the probability density, K ( x , α, β )possesses the following properties, discussed in detail in [11](2.5) K ( x , α, β ) ≥ , (cid:90) Γ d (cid:90) Γ d K ( x , α, β ) dαdβ = 1 . ALEXEY V. SMIRNOV, MICHAEL V. KLIBANOV AND LOC H. NGUYEN
Finally, ν ( x , α ) is the R − vector, showing the direction of particles propagating fromthe external source located at x α = ( α,
0) to x , (2.6) ν ( x , α ) = (cid:18) x − α | x − x α | , y | x − x α | (cid:19) , x ∈ [ − R, R ] , y ∈ [ a, b ] , α ∈ [ − d, d ] . For a fixed α , let ∂ Ω + = { x ∈ ∂ Ω : ν ( x , α ) · n ( x ) ≤ } , where n ( x ) is the unit outward normal vector at ∂ Ω at point x . Assuming that allfunctions in equation (2.3), except u ( x , α ) , are known in Ω, we formulate the followingforward problem. Problem
For each α ∈ [ − d, d ] , find the function u ( x , α ) , satisfying equation (2.3) in the domain Ω as well as the following boundary condition (2.7) u ( x , α ) = 0 for all x ∈ ∂ Ω + . In Appendix we prove existence and uniqueness of the solution of the boundary valueproblem (2.3), (2.7) and; moreover, discuss a numerical method to solve it. Conversely,assume now that the function f ( x ) is unknown and the information of u ( x , α ) on ∂ Ωis known. The main goal of this paper is to numerically solve the following inversesource problem:
Problem
Assume that equation (2.3) and con-ditions (2.4) , (2.5) hold. Also, let the vector ν ( x , α ) in (2.3) has the form (2.6) .Reconstruct the function f ( x ) , x ∈ Ω , given the following boundary data (2.8) F ( x , α ) = u ( x , α ) , for all x ∈ ∂ Ω , α ∈ [ − d, d ] , where u ( x , α ) is the solution of Problem 2.1 and (2.9) F ( x , α ) = 0 for x ∈ ∂ Ω + . Remark
In the particular case when µ a ( x ) ≡ µ s ( x , α ) ≡ , this InverseSource Problem is exactly the problem of X-ray tomography with incomplete data,which was considered in [23]. However, the main focus of this paper is to develop anumerical method for this problem allowing the presence of µ a , µ s , K . Especially, notechnical condition is imposed on these interesting terms.
3. Numerical Method for the Inverse Source Problem.3.1. An orthonormal basis in L ( − d, d ) . First, we recall a special orthonor-mal basis in the space L ( − d, d ), which was introduced in [20]. For α ∈ [ − d, d ] con-sider the set of linearly independent functions (cid:8) α n − e α (cid:9) ∞ n =1 . These functions form acomplete set in L ( − d, d ). Applying the classical Gram-Schmidt orthonormalizationprocedure to this set, we obtain the orthonormal basis { Ψ n ( α ) } ∞ n =1 in L ( − d, d ). Thisbasis has the following properties [20]:1. The functions Ψ n ∈ C [ − d, d ] and Ψ (cid:48) n ( α ) is not identically 0, for all n = 1 , , . . . NVERSE PROBLEM FOR RADIATIVE TRANSFER EQUATION a nn = 1 and a mn = 0 for all m, n = 1 , , . . . such that n < m, where a mn = (cid:90) d − d Ψ (cid:48) n ( α ) , Ψ m ( α ) dα == (cid:26) m = m = n, m > n.m > n. Item 2 implies that the matrix(3.1) M N = ( a mn ) Nm,n =1 is invertible for all N = 1 , , . . . . Hence, the function u ( x , α ) can be written as the following Fourier series converg-ing in L ( − d, d ) for every point x ∈ Ω u ( x , α ) = ∞ (cid:88) n =1 u n ( x )Ψ n ( α ) for all α ∈ [ − d, d ]where u n ( x ) = (cid:90) d − d u ( x , α )Ψ n ( α ) dα. We approximate the function u ( x , α ) via the truncated Fourier series, and the samefor u α ( x , α ) , u ( x , α ) ≈ N (cid:88) n =1 u n ( x )Ψ n ( α ) , x ∈ Ω , α ∈ [ − d, d ] , (3.2) u α ( x , α ) ≈ N (cid:88) n =1 u n ( x )Ψ (cid:48) n ( α ) , x ∈ Ω , α ∈ [ − d, d ] . (3.3)where N ≥ α as in (3.3). These assumptions form our approximatemathematical model mentioned in Section 1. Just as in thefirst step of the above mentioned BK method [9], we eliminate the unknown sourcefunction f ( x ) from equation (2.3) via the differentiation of that equation with respectto the parameter α from which f ( x ) does not depend. We obtain(3.4) ν ( x, y, α ) · ∇ u α − y | x − x α | u x + ( x − α ) y | x − x α | u y + ( µ a + µ s ) ( x ) u α − µ s ( x ) (cid:90) Γ d K α ( x , α, β ) u ( x , β ) dβ = 0for all x = ( x, y ) ∈ Ω . Multiplying equation (3.4) by | x − x α | /y , we obtain u y,α + x − αy u x,α + y | x − x α | u x + ( x − α ) | x − x α | u y ++ | x − x α | y (cid:20) ( µ a + µ s ) ( x ) u α − µ s ( x ) (cid:90) Γ d K α ( x , α, β ) u ( x , β ) dβ (cid:21) = 0 . (3.5) ALEXEY V. SMIRNOV, MICHAEL V. KLIBANOV AND LOC H. NGUYEN
Substituting representations (3.2) and (3.3) into equation (3.5), multiplying theresulting equation by functions Ψ m ( α ) , for each m ∈ { , , . . . , N } we obtain(3.6) N (cid:88) n =1 ∂u n ∂y Ψ (cid:48) n ( α )Ψ m ( α ) + x − αy N (cid:88) n =1 ∂u n ∂x Ψ (cid:48) n ( α )Ψ m ( α )+ y | x − x α | N (cid:88) n =1 ∂u n ∂x Ψ n ( α )Ψ m ( α ) + ( x − α ) | x − x α | N (cid:88) n =1 ∂u n ∂y Ψ n ( α )Ψ m ( α )+ | x − x α | y Ψ m ( α ) N (cid:88) n =1 [( µ a + µ s ) ( x ) u n Ψ (cid:48) n ( α )] − | x − x α | y Ψ m ( α ) N (cid:88) n =1 (cid:20) µ s ( x ) (cid:90) Γ d K α ( x , α, β ) u n ( x )Ψ n ( β ) dβ (cid:21) = 0 . Integrate equation (3.6) with respect to α ∈ ( − d, d ) . Recalling the definition of thematrix M N in (3.1), we obtain(3.7) M N U y = A U y + B U x + C U, U ( x ) = ( u ( x ) , . . . , u N ( x )) T , Here A , B and C are N × N matrices with the following entries:( A ) mn = (cid:90) Γ d ( x − α ) | x − x α | Ψ n ( α )Ψ m ( α ) dα, (3.8) ( B ) mn = (cid:90) Γ d (cid:20) x − αy Ψ (cid:48) n ( α )Ψ m ( α ) + y | x − x α | Ψ n ( α )Ψ m ( α ) (cid:21) dα, (3.9) ( C ) mn = (cid:90) Γ d | x − x α | y ( µ a + µ s ) ( x ) Ψ (cid:48) n ( α )Ψ m ( α ) dα − (cid:90) Γ d | x − x α | y µ s ( x ) (cid:18)(cid:90) Γ d K α ( x , α, β )Ψ n ( β ) dβ (cid:19) Ψ m ( α ) dα. (3.10)Everywhere below the norm of a matrix is the square root of the sum of square normsof its entries. Since in the definition of the domain Ω the number a >
1, the followingestimates follow from (3.8)-(3.10):max x ∈ Ω (cid:12)(cid:12)(cid:12)(cid:12) A ( x ) (cid:12)(cid:12)(cid:12)(cid:12) ≤ C a , max x ∈ Ω (cid:12)(cid:12)(cid:12)(cid:12) B ( x ) (cid:12)(cid:12)(cid:12)(cid:12) ≤ C a , max x ∈ Ω (cid:12)(cid:12)(cid:12)(cid:12) C ( x ) (cid:12)(cid:12)(cid:12)(cid:12) ≤ C a , where the number C = C ( R, d ) > a = a ( N, R, d ) > a > a the matrix ˜ A = M N (Id − M − N A ) is invertible. Everywhere below we assume withoutfurther mentioning that a > a . Denote A := ˜ A − B , A := ˜ A − C . Therefore, equation (3.7) is equivalent to(3.11) U y − A U x − A U = 0; A = A ( x, y ) , A = A ( x, y ) , ( x, y ) ∈ Ω . Using (2.8) and (2.9) we complement equation (3.11) with the following Dirichletboundary condition(3.12) U = F ( x, y ) , for ( x, y ) ∈ ∂ Ω . NVERSE PROBLEM FOR RADIATIVE TRANSFER EQUATION U ( x ) = ( u ( x ) , . . . , u N ( x )) T of theboundary value problem (3.11)–(3.12) directly yields the desired numerical solutionto Problem 2.2 via the substitution of (3.2) in (2.3). The problem (3.11)–(3.12) isan overdetermined one. Indeed, although equations (3.11) are of the first order, theboundary condition (3.12) is given on the entire boundary ∂ Ω . To find an approximatesolution to this problem, we use the QRM, which, in general works properly foroverdetermined problems for PDEs. Thus, we consider the following minimizationproblem for the Tikhonov-like functional J (cid:15) with the regularization parameter (cid:15) ∈ (0 ,
1) :(3.13) J (cid:15) ( U ) = (cid:90) Ω | U y − A U x − A U | dxdy + (cid:15) (cid:107) U (cid:107) H (Ω) . When we say below that a vector function belongs to a Hilbert space, we mean thateach of its components belongs to this space and its norm is the square root of thesum of norms in that space of its components.
Problem
Minimize the functional J (cid:15) on the setof N -dimensional vector valued functions U ∈ H (Ω) satisfying boundary condition(3.12).
4. The Fully Discrete Form of the QRM.
To solve Problem 3.1, we write U x , U y in the functional J (cid:15) ( U ) in its finite difference version and minimize it withrespect to values of the vector function U at grid points. Hence, we formulate theQRM in this section in the fully discrete form of finite differences. We prove existenceand uniqueness of the minimizer and establish convergence rate of minimizers to theexact solution, which is also written via finite differences. Consider the followinguniform 2-dimensional grid points on Ω whose x and y coordinates are given by − R = x < x < · · · < x M x = R, x i +1 − x i = h x , ∀ i ∈ { , , . . . , M x − } , (4.1) a = y < y < · · · < y M y = b, y j +1 − y j = h y , ∀ j ∈ { , , . . . , M y − } . (4.2)Denote h = ( h x , h y ) . We define the discrete set Ω h asΩ h = { ( x, y ) : { ( x i , y j ) } , i ∈ { , . . . , M x − } , j ∈ { , . . . , M y − }} ,∂ Ω h = { ( x, y ) : { ( x i , y j ) } for i = 0 , M x , j = 0 , M y } , Ω h = Ω h ∪ ∂ Ω h . For any N − D matrix Q ( x, y ) ∈ C (Ω) we introduce the following notations Q hi,j = Q ( x i , y j ) , i ∈ { , . . . , M x − } , j ∈ { , . . . , M y − } , (cid:101) Q hi,j = Q ( x i , y j ) , i ∈ { , . . . , M x } , j ∈ { , . . . , M y } , Q h = (cid:8) Q hi,j (cid:9) M x − ,M y − i,j =1 is an ( M x − × ( M y −
1) matrix , (cid:101) Q h = (cid:8) Q hi,j (cid:9) M x ,M y i,j =0 is an ( M x + 1) × ( M y + 1) matrix . (4.3) ALEXEY V. SMIRNOV, MICHAEL V. KLIBANOV AND LOC H. NGUYEN
Note that the matrix Q h , in contrast to (cid:101) Q h , does not include boundary terms ofthe form Q h ,j = Q ( − R, y j ) , Q hi, = Q ( x i , a ) , Q hM x ,j = Q ( R, y j ) , Q hi,M y = Q ( x i , b ) . Recall the forward finite difference formulae for the vector function Q h :( U h ) (cid:48) x = (cid:8) ( U hi,j ) (cid:48) x (cid:9) M x − ,M y − i,j =0 , ( U hi,j ) (cid:48) x = U hi +1 ,j − U hi,j h x , i ∈ { , . . . , M x − } , (4.4) ( U h ) (cid:48) y = (cid:8) ( U hi,j ) (cid:48) y (cid:9) M x − ,M y − i,j =0 , ( U hi,j ) (cid:48) y = U hi,j +1 − U hi,j h y , j ∈ { , . . . , M y − } . (4.5)Hence, we obtain the following finite difference analog of (3.11)–(3.12) L h (cid:0) U h (cid:1) = ( U h ) (cid:48) y − A h ( U h ) (cid:48) x + A h U h = 0; A h = A , A h = A in Ω h , (4.6) (cid:101) U h = F h on ∂ Ω h . (4.7)where the boundary matrix F h is defined using the values of the matrix F ( x, y ) onthe grid (4.1), (4.2). We define the following discrete functional spaces for matrices Q h , (cid:101) Q h : L ,h (Ω h ) = (cid:110) Q h : (cid:107) Q h (cid:107) L ,h (Ω h ) = h y h x M y − (cid:88) j =1 M x − (cid:88) i =1 [ Q hi,j ] < ∞ (cid:111) , and H ,h (Ω h ) = (cid:110) Q h : (cid:107) Q h (cid:107) H ,h (Ω h ) = h y h x M x − (cid:88) i =1 M y − (cid:88) j =1 ([( Q hi,j ) (cid:48) x ] + [( Q hi,j ) (cid:48) y ] + [ Q hi,j ] ) < ∞ (cid:111) . We define the inner products in these spaces in the obvious manner and denote themas (cid:0) · , · (cid:1) and (cid:2) · , · (cid:3) for L ,h (Ω h ) and H ,h (Ω h ) respectively. Remark
Here and everywhere below if a matrix Q h is defined as in (4.3),then (cid:101) Q h denotes the matrix Q h , complemented by boundary conditions at ∂ Ω h . Remark
Below we fix the number h ∈ (0 , and restrict h x from the belowas h x ∈ [ h , . However, we do not restrict from the below h y > by a positiveconstant. Then it follows from (4.5) that there exists a constant B h > dependingonly on h such that if Q h ,j = Q hM x ,j = 0; j = 1 , ..., M y − , then (4.8) (cid:107) ( Q h ) (cid:48) x (cid:107) L ,h (Ω h ) ≤ B h (cid:107) Q h (cid:107) L ,h (Ω h ) , ∀ h ∈ [ h , , ∀ Q h : (cid:101) Q h ∈ H ,h (Ω h ) . The fully discrete QRM applied to problem (4.6)-(4.7) leads to the followingdiscrete version of the above Minimization Problem:
Problem
Minimize the functional (4.9) J h(cid:15) ( (cid:101) U h ) = (cid:107) ( U h ) (cid:48) y − A h ( U h ) (cid:48) x − A h U h (cid:107) L ,h (Ω h ) + (cid:15) (cid:107) U h (cid:107) H ,h (Ω h ) on the set of matrices (cid:101) U h , satisfying the boundary condition (4.7) . The minimizer of J h(cid:15) ( (cid:101) U h ) satisfying boundary condition (4.7) is called the regularizedsolution of the problem (4.6)–(4.7). NVERSE PROBLEM FOR RADIATIVE TRANSFER EQUATION We now derive a discrete Carleman es-timate for the finite difference version of the differential operator d/dy . Consider auniform partition of the interval ( a, b ) ⊂ R of the real line into M subintervals withthe grid step size h y ,(4.10) a = y < y < . . . < y M = b, y j +1 − y j = h y , j ∈ { , , . . . , M − } . Following the book [32], for any discrete function u defined on this grid denote u j = u ( y j ) and define both its forward u (cid:48) j and backward u (cid:48) j finite difference derivatives,which are the finite difference analogs of the differential operator d/dy, as(4.11) u (cid:48) j = ( u j +1 − u j ) h y , j ∈ { , . . . , M − } , u (cid:48) j = ( u j − u j − ) h y , ∀ j ∈ { , . . . , M } . Lemma
For any discrete function w , defined on the grid (4.10) the followinginequality holds: − h y M − (cid:88) j =1 w j w (cid:48) j ≥ − ( w M − w ) . Proof.
Using the summation by parts formula for the discrete function w [32],we obtain h y M − (cid:88) j =1 w j w (cid:48) j = ( w M − w w ) − h y M (cid:88) i =1 w i w i (cid:48) . Next, h y M (cid:88) i =1 w i w i (cid:48) = h y M − (cid:88) j =0 w j +1 w (cid:48) j = h y M − (cid:88) j =0 ( w j + w (cid:48) j h y ) w (cid:48) j = h y ( w + hw (cid:48) ) w (cid:48) + h y M − (cid:88) j =1 ( w j + w (cid:48) j h y ) w (cid:48) j = w ( w − w ) + h y M − (cid:88) j =1 ( w j + w (cid:48) j h y ) w (cid:48) j Combining all equalities written above, we obtain h y M − (cid:88) j =1 w j w (cid:48) j = (cid:0) w M − w w (cid:1) − ( w − w w ) − h y M − (cid:88) j =1 (cid:0) w j + w (cid:48) j h y (cid:1) w (cid:48) j = (cid:0) w M − w (cid:1) − h y M − (cid:88) j =1 w j w (cid:48) j − h y M − (cid:88) j =1 ( w (cid:48) j ) . Hence, − h y M − (cid:88) j =1 w j w (cid:48) j = − ( w M − w ) + h y M − (cid:88) j =1 ( w (cid:48) j ) ≥ − (cid:0) w M − w (cid:1) . (cid:3) ALEXEY V. SMIRNOV, MICHAEL V. KLIBANOV AND LOC H. NGUYEN
Theorem
For any positive number λ > ,the following discrete Carleman estimate holds for any discrete function u , defined onthe grid (4.10) h y M y − (cid:88) j =1 e λy j (cid:0) u (cid:48) j (cid:1) ≥ h y M y − (cid:88) j =1 (cid:18) − e − λh y h y (cid:19) e λy j u j + 2 e − λh y (cid:18) − e − λh y h y (cid:19) ( e λy u − e λy M u M ) . Proof . For each j , we define(4.12) w j = e λy j u j , u j = e − λy j w j . Hence, according to (4.11), the forward difference derivative of the function u at y j is u (cid:48) j = e − λ ( y j + h y ) w j +1 − e − λy j w j h y = e − λy j (cid:0) e − λh y w j +1 − w j (cid:1) h y = e − λy j (cid:18) e − λh y w j +1 − w j h y + e − λh y w j − w j h y (cid:19) = e − λy j (cid:18) w (cid:48) j e − λh y − − e − λh y h y w j (cid:19) . Hence, we have for each j = 1 , ..., M y − e λy j ( u (cid:48) j ) = (cid:18) w (cid:48) j e − λh y − − e − λh y h y w j (cid:19) ≥ (cid:18) − e − λh y h y (cid:19) ( w j ) − e − λh y (cid:0) − e − λh y (cid:1) h y w (cid:48) j w j . As a result, h y M y − (cid:88) j =1 e λy j ( u (cid:48) j ) ≥ h y M y − (cid:88) j =1 (cid:18) − e − λh y h y (cid:19) w j − e − λh y (1 − e − λh y ) h y h y M − (cid:88) j =1 w (cid:48) j w j . Applying Lemma 4.1 to the second term in the right hand side, we obtain h y M y − (cid:88) j =1 e λy j ( u (cid:48) j ) ≥ h y M y − (cid:88) j =1 (cid:18) − e − λh y h y (cid:19) w j + 2 e − λh y (cid:18) − e − λh y h y (cid:19) ( w − w M ) . The statement of Theorem 4.2 follows from this estimate and (4.12). (cid:3)
Lemma
Let u be a discrete function, defined on the grid (4.10) , such that u M = 0 . Then for any two numbers λ, h y > such that λh y < the followinginequality holds (4.13) h y M y − (cid:88) j =1 e λy j (cid:0) u (cid:48) j (cid:1) ≥ λ h y M y − (cid:88) j =1 e λy j u j . NVERSE PROBLEM FOR RADIATIVE TRANSFER EQUATION Proof.
By Taylor formula e − λh y = 1 − λh y + e − ξ λh y ) = 1 − λh y (cid:18) − e − ξ λh y (cid:19) , where ξ ∈ (0 , λh y ) is a certain number. Hence, 1 − e − λh y ≥ λh y / . Hence, (cid:18) − e − λh y h y (cid:19) ≥ λ . Therefore, using Theorem 4.2, we obtain (4.13). (cid:3)
Remark
This lemma is a discrete analog of the Carleman estimate in [23,Lemma 4.1] for the continuous case of the operator d/dy.
5. Convergence Analysis.5.1. Existence of the solution of the Discrete Minimization Problem.
Theorem
Assume that there exists a matrix G h ∈ H ,h (Ω h ) such that (cid:101) G h | ∂ Ω h = F h . Then for each (cid:15) > , there exists unique minimizer U h min ,(cid:15) ∈ H ,h (Ω h ) of the functional (4.9) satisfying boundary condition (4.7). Proof.
Let H ,h (Ω h ) be the subspace of the space H ,h (Ω h ) consisting on suchmatrices (cid:101) Q h ∈ H ,h (Ω h ) that (cid:101) Q h | ∂ Ω h = 0 . Recalling notation (4.6) for the operator L h , we rewrite the functional J h(cid:15) ( (cid:101) U h ) in the following form: J h(cid:15) ( W h ) = (cid:13)(cid:13) L h (cid:0) W h (cid:1) + L h (cid:0) G h (cid:1)(cid:13)(cid:13) L ,h (Ω h ) + (cid:15) (cid:13)(cid:13) W h + G h (cid:13)(cid:13) H ,h (Ω h ) , (5.1) (cid:102) W h | ∂ Ω = ( (cid:101) U h − (cid:101) G h ) | ∂ Ω = 0 . (5.2)Thus, in order to work with zero boundary condition in (5.1)-(5.2), we consider thefunction W h = U h − G h instead of U h . Let W h min ,(cid:15) with (cid:102) W h min ,(cid:15) ∈ H ,h (Ω h ) be any minimizer of functional (5.1). Bythe variational principle the following identity holds for all (cid:101) P h ∈ H ,h (Ω h ) :(5.3) (cid:0) L h (cid:0) W h min ,(cid:15) (cid:1) , L h (cid:0) P h (cid:1)(cid:1) + (cid:15) (cid:2) W h min ,(cid:15) , P h (cid:3) = − (cid:0) L h (cid:0) G h (cid:1) , L h (cid:0) P h (cid:1)(cid:1) − (cid:15) (cid:2) G h , P h (cid:3) . The left hand side of the identity (5.3) generates a new scalar product {· , ·} in thesubspace H ,h (Ω h ) . Consider the corresponding norm {·} , (5.4) (cid:8) Q h (cid:9) = (cid:13)(cid:13) L h (cid:0) Q h (cid:1)(cid:13)(cid:13) L ,h (Ω h ) + (cid:15) (cid:13)(cid:13) Q h (cid:13)(cid:13) H ,h (Ω h ) , ∀ (cid:101) Q h ∈ H ,h (Ω h ) . Obviously, there exists a certain constant C = C (cid:0) L h , h, Ω h , (cid:15) (cid:1) > , which dependsonly on listed parameters such that(5.5) (cid:15) (cid:13)(cid:13) Q h (cid:13)(cid:13) H ,h (Ω h ) ≤ (cid:8) Q h (cid:9) ≤ C (cid:13)(cid:13) Q h (cid:13)(cid:13) H ,h (Ω h ) , ∀ (cid:101) Q h ∈ H ,h (Ω h ) , Below C denotes different positive numbers depending on the same parameters.Hence, norms (cid:8) Q h (cid:9) and (cid:13)(cid:13) Q h (cid:13)(cid:13) H ,h (Ω h ) are equivalent for (cid:101) Q h ∈ H ,h (Ω h ) . There-fore, (5.3) is equivalent with(5.6) (cid:110) (cid:102) W h min ,µ , P h (cid:111) = − (cid:0) L h (cid:0) G h (cid:1) , L h (cid:0) P h (cid:1)(cid:1) − (cid:15) (cid:2) G h , P h (cid:3) , ∀ (cid:101) P h ∈ H ,h (Ω h ) . ALEXEY V. SMIRNOV, MICHAEL V. KLIBANOV AND LOC H. NGUYEN
Using the Cauchy-Schwarz inequality, (5.4) and (5.5), we obtain (cid:12)(cid:12) − (cid:0) L h (cid:0) G h (cid:1) , L h (cid:0) P h (cid:1)(cid:1) − (cid:15) (cid:2) G h , P h (cid:3)(cid:12)(cid:12) ≤ C (cid:8) G h (cid:9) (cid:8) P h (cid:9) , for all (cid:101) P h ∈ H ,h (Ω h ) . Hence, the right hand side of (5.6) can be considered as a bounded linear functionalmapping the space H ,h (Ω h ) in R . Since the regular norm in H ,h (Ω h ) is equiva-lent with the norm generated by new scalar product {· , ·} , then Riesz representationtheorem implies that there exists unique matrix (cid:101) Φ h ∈ H ,h (Ω h ) such that (cid:110) (cid:102) W h min ,(cid:15) , (cid:101) P h (cid:111) = (cid:110) (cid:101) Φ h , (cid:101) P h (cid:111) , for all (cid:101) P h ∈ H ,h (cid:0) Ω h (cid:1) . Therefore, (cid:102) W h min ,(cid:15) = (cid:101) Φ h and W h min ,(cid:15) = Φ h . Finally, the matrix U h min ,(cid:15) = W h min ,(cid:15) + G h is the unique minimizer claimed by this theorem. (cid:3) The minimizer (cid:101) U h min ,(cid:15) is called the regularized solution of problem (4.6),(4.7). In this section, we establish the convergence rate of regularized solutions to theexact one when the noise in the data tends to zero. In addition, we establish Lipschitzstability estimate and uniqueness for the problem (4.6), (4.7).Let a matrix P h ∈ L ,h (cid:0) Ω h (cid:1) . Denote(5.7) (cid:13)(cid:13) P h e λy (cid:13)(cid:13) L ,h (Ω h ) = h y h x M y − (cid:88) j =1 M x − (cid:88) i =1 (cid:0) P hi,j (cid:1) e λy j . Hence, by Lemma 4.3 for all λh y ∈ (0 ,
1) and for all (cid:101) P h ∈ H ,h (cid:0) Ω h (cid:1) (5.8) (cid:13)(cid:13)(cid:13)(cid:0) P h (cid:1) (cid:48) y e λy (cid:13)(cid:13)(cid:13) L ,h (Ω h ) ≥ λ (cid:13)(cid:13) P h e λy (cid:13)(cid:13) L ,h (Ω h ) . Let U ∗ h ∈ H ,h (cid:0) Ω h (cid:1) be the exact solution of problem (4.6), (4.7) with the exactboundary data F ∗ ,h . We assume that there exists an exact matrix G ∗ h such that(5.9) G ∗ h ∈ H ,h (Ω h ) , (cid:101) G ∗ ,h | ∂ Ω h = F ∗ ,h . As to the boundary data F h , we assume, as in Theorem 2, that there exists a matrix G h ∈ H ,h (Ω h ) such that (cid:101) G h | ∂ Ω h = F h . In addition, we assume that G h is givenwith a noise of the level δ ∈ (0 , , i.e.(5.10) (cid:107) G ∗ ,h − G h (cid:107) H ,h (Ω h ) ≤ δ. Our main goal now is to estimate the difference between U h min ,(cid:15) and U ∗ h via δ and (cid:15). Lemma
There exists a number C = C ( N, d, R, h , a, b, L h ) > and a suf-ficiently small number h y = h y ( N, d, R, h , a, b, L h ) ∈ (0 , , both depending only onlisted parameters, such that for h x ∈ [ h , , h y ∈ (cid:0) , h y (cid:3) the following estimate isvalid (5.11) (cid:13)(cid:13) L h (cid:0) Q h (cid:1)(cid:13)(cid:13) L ,h (Ω h ) ≥ C (cid:13)(cid:13) Q h (cid:13)(cid:13) L ,h (Ω h ) , for all (cid:101) Q h ∈ H ,h (Ω h ) . NVERSE PROBLEM FOR RADIATIVE TRANSFER EQUATION Proof.
Below C > L h in (4.6) as well as (5.4), (5.5) andthe Cauchy-Schwarz inequality, we obtain (cid:13)(cid:13) L h (cid:0) Q h (cid:1)(cid:13)(cid:13) L ,h (Ω h ) = (cid:13)(cid:13) L h (cid:0) Q h (cid:1) e λy e − λy (cid:13)(cid:13) L ,h (Ω h ) ≥ e − λb (cid:13)(cid:13) L h (cid:0) Q h (cid:1) e λy (cid:13)(cid:13) L ,h (Ω h ) ≥ e − λb (cid:20)(cid:13)(cid:13)(cid:13)(cid:0) Q h (cid:1) (cid:48) y e λy (cid:13)(cid:13)(cid:13) L ,h (Ω h ) − C (cid:13)(cid:13) Q h (cid:13)(cid:13) L ,h (Ω h ) (cid:21) . Choose h y ∈ (0 ,
1) so small that 1 / (cid:0) h y (cid:1) > C and let h y ∈ (cid:0) , h y (cid:1) . Set λ =1 / (2 h y ) . Then λh y < / < λ / / (cid:16)
32 ( h y ) (cid:17) > C . Hence, by (5.8)and (5.11) it follows from the above inequality (cid:13)(cid:13) L h (cid:0) Q h (cid:1)(cid:13)(cid:13) L ,h (Ω h ) ≥ e − λb (cid:20) λ (cid:13)(cid:13) Q h e λy (cid:13)(cid:13) L ,h (Ω h ) − C (cid:13)(cid:13) Q h e λy (cid:13)(cid:13) L ,h (Ω h ) (cid:21) ≥ e − λb (cid:16) C (cid:13)(cid:13) Q h e λy (cid:13)(cid:13) L ,h (Ω h ) − C (cid:13)(cid:13) Q h e λy (cid:13)(cid:13) L ,h (Ω h ) (cid:17) = 12 e − λb C (cid:13)(cid:13) Q h e λy (cid:13)(cid:13) L ,h (Ω h ) . This estimate immediately implies (5.10) with a new constant C > . (cid:3) Theorem
Assume that condi-tions of Theorem 5.1 as well as (5.10) and (5.11) hold. Let U h min ,(cid:15) ∈ H ,h (Ω h ) bethe unique minimizer of the functional (4.9) that satisfies boundary condition (4.7)(see Theorem 5.1). Suppose that h x ∈ [ h , and h y ∈ (cid:0) , h y (cid:3) , where the number h y is defined in Lemma 5.2. Then for any (cid:15) > the following convergence rate ofregularized solutions holds (5.12) (cid:107) U h min ,(cid:15) − U ∗ h (cid:107) L ,h (Ω h ) ≤ C ( δ + √ (cid:15) (cid:107) U ∗ h (cid:107) H ,h (Ω h ) ) . Proof.
Define the matrix (cid:102) W h min ,(cid:15) ∈ H ,h (Ω h ) as in Theorem 5.1, i.e. (cid:102) W h min ,(cid:15) = (cid:101) U h min ,(cid:15) − (cid:101) G h . Similarly define (cid:102) W ∗ ,h = (cid:101) U ∗ ,h − (cid:101) G ∗ ,h ∈ H ,h (Ω h ) . Then (5.3) is valid for W h min ,(cid:15) . As to W ∗ ,h , (4.6) and (4.7) imply that the following analog of (5.3) is validfor all (cid:101) P h ∈ H ,h (Ω h ) : (cid:0) L h (cid:0) W ∗ ,h (cid:1) , L h (cid:0) P h (cid:1)(cid:1) + (cid:15) (cid:2) W ∗ ,h , P h (cid:3) = − (cid:0) L h (cid:0) G ∗ h (cid:1) , L h (cid:0) P h (cid:1)(cid:1) + (cid:15) (cid:2) W ∗ ,h , P h (cid:3) . Denote (cid:101) V h = (cid:102) W h min ,(cid:15) − W ∗ ,h ∈ H ,h (Ω h ) , Y h = G h − G ∗ h and subtract (5.14) from(5.3). We obtain for all (cid:101) P h ∈ H ,h (Ω h )(5.13) (cid:0) L h (cid:0) V h (cid:1) , L h (cid:0) P h (cid:1)(cid:1) + (cid:15) (cid:2) V h , P h (cid:3) = − (cid:0) L h (cid:0) Y h (cid:1) , L h (cid:0) P h (cid:1)(cid:1) − (cid:15) (cid:2) W ∗ ,h , P h (cid:3) . Set in (5.13) P h = V h and use the Cauchy-Schwarz inequality. We obtain (cid:13)(cid:13) L h (cid:0) V h (cid:1)(cid:13)(cid:13) L ,h (Ω h ) ≤ (cid:13)(cid:13) L h (cid:0) Y h (cid:1)(cid:13)(cid:13) L ,h (Ω h ) + (cid:15) (cid:13)(cid:13) W ∗ ,h (cid:13)(cid:13) H ,h (Ω h ) ≤ C δ + (cid:15) (cid:13)(cid:13) W ∗ ,h (cid:13)(cid:13) H ,h (Ω h ) . ALEXEY V. SMIRNOV, MICHAEL V. KLIBANOV AND LOC H. NGUYEN
Next, using (5.11) and (5.13), we obtain(5.14) (cid:13)(cid:13) V h (cid:13)(cid:13) L ,h (Ω h ) ≤ C (cid:16) δ + (cid:15) (cid:13)(cid:13) W ∗ ,h (cid:13)(cid:13) H ,h (Ω h ) (cid:17) . The target estimate (5.12) follows immediately from (5.14). (cid:3)
Theorem
Suppose that there exist twomatrices G h , G h ∈ H ,h (Ω h ) such that (cid:101) G h | ∂ Ω h = F h and (cid:101) G h | ∂ Ω h = F h , where F h and F h are two different boundary conditions in (4.7). Suppose that there existsolutions (cid:101) U h ∈ H ,h (cid:0) Ω h (cid:1) and (cid:101) U h ∈ H ,h (cid:0) Ω h (cid:1) of boundary value problem (4.6)-(4.7)with boundary conditions F h and F h respectively. Assume that h x ∈ [ h , and h y ∈ (cid:0) , h y (cid:3) , where the number h y is defined in Lemma 5.2. Then the followingLipschitz stability estimate is valid (5.15) (cid:13)(cid:13) U h − U h (cid:13)(cid:13) L ,h (Ω h ) ≤ C (cid:13)(cid:13) G h − G h (cid:13)(cid:13) L ,h (Ω h ) . Next, suppose that F h = F h , but the existence of the function (cid:101) G h is not assumed.Then U h = U h , where U h , U h ∈ H ,h (cid:0) Ω h (cid:1) are two possible solution of boundaryvalue problem (4.6), (4.7). Proof . Since (cid:101) U h and (cid:101) U h are two exact solutions of problem (4.6), (4.7) with twodifferent boundary conditions, then by (5.12)(5.16) (cid:0) L h (cid:0) U hi (cid:1) , L h (cid:0) P h (cid:1)(cid:1) = − (cid:0) L h (cid:0) G hi (cid:1) , L h (cid:0) P h (cid:1)(cid:1) , ∀ P h ∈ H ,h (cid:0) Ω h (cid:1) , i = 1 , , where (cid:102) W hi = (cid:101) U hi − (cid:101) G hi . Setting S h = W h − W h , X h = G h − G h and then setting P h = S h , we obtain from (5.16) and (4.8) (cid:13)(cid:13) L h (cid:0) S h (cid:1)(cid:13)(cid:13) L ,h (Ω h ) ≤ (cid:13)(cid:13) L h (cid:0) X h (cid:1)(cid:13)(cid:13) L ,h (Ω h ) ≤ C (cid:13)(cid:13) X h (cid:13)(cid:13) L ,h (Ω h ) . Hence, by (5.11) (cid:13)(cid:13) S h (cid:13)(cid:13) L ,h (Ω h ) ≤ C (cid:13)(cid:13) X h (cid:13)(cid:13) L ,h (Ω h ) . Therefore,(5.17) (cid:13)(cid:13)(cid:0) U h − U h (cid:1) − (cid:0) G h − G h (cid:1)(cid:13)(cid:13) L ,h (Ω h ) ≤ C (cid:13)(cid:13) G h − G h (cid:13)(cid:13) L ,h (Ω h ) . Thus, (5.15) follows from (5.17) and the triangle inequality. As to the uniquenesspart, since F h = F h , then we extend the boundary condition (cid:0) F h − F h (cid:1) = 0 in thedomain Ω h as G h − G h ≡ . Hence, (5.15) implies that U h − U h ≡ . (cid:3) Remark
Due to the ill-posedness of Problem 2.2, we cannot prove conver-gence of our solutions to the correct one as N → ∞ . We note that the truncatedFourier series is used both quite often and quite successfully in numerical methods formany inverse problems. Although convergences at N → ∞ are not proven in manycases, numerical results are usually good ones, see, e.g. [10] for the attenuated tomog-raphy with complete data, [14] for the 2D version of the Gelfand-Levitan method, [26]for the inverse problem of finding initial condition of heat equation, [30] for the inversesource problem for the Helmholtz equation, and [18, 21, 22] for the convexification.
6. Numerical Implementation.
In this section, we describe the numericalimplementation of the minimization procedure for the functional J h(cid:15) . While invertingthe matrix M N of (3.7) is convenient for the convergence analysis, we have discoveredthat it is better in real computations not to invert while still considering a problem NVERSE PROBLEM FOR RADIATIVE TRANSFER EQUATION J h(cid:15) ,(cid:15) ( (cid:101) U h ) = (cid:107) ( M N − A h )( U h ) (cid:48) y − B h ( U h ) (cid:48) x − C h U h (cid:107) L ,h (Ω h ) ++ (cid:15) (cid:107) U h (cid:107) L ,h (Ω h ) + (cid:15) (cid:107)∇ U h (cid:107) L ,h (Ω h ) , where A h , B h , C h are operators (3.8)-(3.10), with the domain Ω h . Moreover, in con-trast to the original functional, we use in our computations two regularization pa-rameters (cid:15) and (cid:15) , instead of just one parameter (cid:15) . This yields better reconstructionresults. The regularization parameters (cid:15) , (cid:15) in our numerical tests were found by atrial and error procedure. They were the same for all the tests we have conducted.To minimize the functional J h(cid:15) ,(cid:15) ( (cid:101) U h ) , we first have to simulate the boundarydata via solving Problem 2.1. We discuss the solution of this forward problem inAppendix. Using this solution, we generate the noisy data as, see (2.8), (2.9) and(5.10): F ( x , α ) = (cid:26) u compδ ( x , α ) = u comp ( x , α )(1 + δ (2rand( x ) − , x ∈ ∂ Ω \ ∂ Ω + , , x ∈ ∂ Ω + . where u comp ( x , α ) is the boundary data computed via the solution of the forwardproblem, δ > · ) is the function that generates uniformlydistributed random numbers in the interval [0 , δ = 0 . h x = h y ,M x = M y . We rewrite the functional J h(cid:15) ,(cid:15) in the following discrete form J h(cid:15) ,(cid:15) ( (cid:101) U h ) = h x M x − (cid:88) i,j =1 (cid:32) ( M N − A h ) ij U hi,j +1 − U hi,j h x − ( B h ) ij U hi +1 ,j − U hi,j h x − ( C h ) ij U hi,j (cid:33) + (cid:15) h x M x − (cid:88) i,j =1 (cid:0) U hi,j (cid:1) + (cid:15) h x M x − (cid:88) i,j =1 (cid:32) | U hi,j +1 − U hi,j | h x + | U hi +1 ,j − U hi,j | h x (cid:33) . Denote u m ( x i , y j ) = u i,jm . Since U := ( u ( x, y ) , . . . , u N ( x, y )) T , then J h(cid:15) ,(cid:15) ( (cid:101) U h ) = h x M x − (cid:88) i,j =1 N (cid:88) m =1 (cid:18) ( M N − A h ) ij u i,j +1 m − u i,jm h x − ( B h ) ij u i +1 ,jm − u i,jm h x − ( C h ) ij u i,jm (cid:19) + (cid:15) h x M x − (cid:88) i,j =1 N (cid:88) m =1 (cid:0) u i,jm ) (cid:1) + (cid:15) h x M x − (cid:88) i,j =1 N (cid:88) m =1 (cid:18) | u i,j +1 m − u i,jm | h x + | u i +1 ,jm − u i,jm | h x (cid:19) . Introduce the “lined-up” versions of the matrices (cid:101) U h , M N − A h , B h , C h . The( M x + 1) N dimensional vector U , (6.1) U m = u m ( x i , y j ) 1 ≤ i, j ≤ M x + 1 , ≤ m ≤ N, and the ( M x + 1) N × ( M x + 1) N dimensional matrices A h , B h , C h , correspondingto M N − A h , B h , C h , where(6.2) m = ( i − M x + 1) N + ( j − N + m. ALEXEY V. SMIRNOV, MICHAEL V. KLIBANOV AND LOC H. NGUYEN
We introduce the map { , . . . , M x + 1 } × { , . . . , M x + 1 } × { , . . . , N } → { , . . . , ( M x + 1) N } that sends ( i, j, m ) to m as in (6.2) is onto and one-to-one. The functional J h(cid:15) ,(cid:15) ( (cid:101) U h )is rewritten in terms of the lined-up vector U as(6.3) J h(cid:15) ,(cid:15) ( U ) = h x (cid:0) |LU| + (cid:15) |U| + (cid:15) |D x U| + (cid:15) |D x U| (cid:1) , where D x and D y are the matrices that provide the finite difference analogs of thepartial derivatives of U with respect to x and y , defined similarly to (4.4),(4.5). L isthe ( M x + 1) N × ( M x + 1) N matrix defined as follows. For each(6.4) m = ( i − M x + 1) N + ( j − N + m, ≤ i, j ≤ M x , ≤ m ≤ N L mn = ( − ( A h ) mn + ( B h ) mn ) /h x − ( C h ) mn , if m corresponds to ( i, j, n ) in thesense of (6.4) for each n ∈ { , . . . , N } ,2. L mn = ( A h ) mn /h x , if m corresponds to ( i, j + 1 , n ) in the sense of (6.4) foreach n ∈ { , . . . , N } ,3. L mn = − ( B h ) mn /h x , if m corresponds to ( i + 1 , j, n ) in the sense of (6.4) foreach n ∈ { , . . . , N } ,4. L mn = 0 otherwise.Next, we consider the “lined-up” version of the boundary condition (4.7). Let D be the ( M x + 1) N × ( M x + 1) N diagonal matrix with m th diagonal entries takingvalue 1 while the others equal 0. This Dirichlet boundary constraint of the vector U becomes DU = ˜ F . Here, the vector ˜ F is the “lined-up” vector of the data F N in thesame manner when we defined U , see (6.1). We solve the Inverse Source Problem bycomputing the vector U , subject to constraint DU = ˜ F , such that(6.5) L µ U = (cid:0) L T L + (cid:15) Id + (cid:15) D Tx D x + (cid:15) D Ty D y (cid:1) U = (cid:126) , which is equivalent to the minimization of the functional (6.3).The knowledge of U yields the knowledge of (cid:101) U h via (6.1). Denote the resultobtained by the procedure of this section as (cid:101) U comp = ( u comp ( x, y ) , . . . , u compN ( x, y )) T .Using this vector function, we calculate function u comp ( x, y, α ) according via (3.2).Next, the reconstructed function f comp ( x, y ) is determined as the averaged over α value of the source function f comp ( x, y, α ) calculated via the substitution of u comp ( x, y, α ) in (2.3).These arguments lead to the Algorithm 6.1 for solving Problem 2.2 Algorithm 6.1
The procedure to solve Problem 2.2 Choose a number N . Construct the functions Ψ m , 1 ≤ m ≤ N, in Section 3.1and compute the matrix M N as in (3.1). Calculate the boundary data F h for the vector valued function (cid:101) U h on ∂ Ω h viasolving Problem 2.1 with f = f true . Find an approximate solution of (4.6)-(4.7) by the quasi-reversibility method. Having (cid:101) U h , calculate u comp ( x , α ) for x ∈ Ω h via (3.2). Compute the reconstructed function f comp by (2.3). NVERSE PROBLEM FOR RADIATIVE TRANSFER EQUATION We test our method using numerical simulations for dif-ferent types of absorption and scattering coefficients. Test 1 demonstrates the stabilityof the solution in the “no scattering” model, which corresponds to the µ s ( x , α ) ≡ (cid:15) , (cid:15) , d, N, which we use in subsequent tests. The optimal values for thoselisted parameters are: (cid:15) = 0 . , (cid:15) = 0 . , d = 5 , N = 12.In Test 2 we consider the “uniform scattering” version of the RTE for the objectwith a non-smooth boundary, where µ s ( x ) (cid:54)≡ , K ( x , α, β ) = 1 / (2 d ) . Test 3 demonstrates the performance of our method for both sophisticated formof the absorption coefficient µ a ( x ) and “strongly forward-peaked scattering” case,i.e. the case when the scattered particles move in the direction preferentially closeto the one they were moving. That corresponds to the case when µ s ( x ) (cid:54)≡ K ( x , α, β ) (cid:54) = 0 whenever | α − β | ≤ C , for some constant C >
0. We choose our scat-tering phase function K ( x , α, β ) to be a 2-dimensional Henyey-Greenstein function,which, possesses the property mentioned above and is convenient to describe stronglyforward-peaked scattering, according to [11].All numerical tests were provided on the domain Ω h defined in Section 4 on theuniform 100 ×
100 grid (4.1),(4.2) with a = 1 , b = 3 , d = 5, R = 1 , see (2.1), (2.2). Weuse the uniform grid for α as well − d = α < α < · · · < α M α = d, where M α = 50 . Similarly with [23], we apply a 2-step post-processing procedure. Define m = max x ∈ Ω h ( f compδ ( x )). Then, the first step of the procedure removes undesirableartifacts by setting ˜ f compδ ( x ) = (cid:26) f compδ ( x ) if f compδ ( x ) > . m f compδ ( x ) function out. For every gridpoint x of Ω h , the value of ˜ f compδ ( x ) at that point is replaced by the mean valueˆ f compδ ( x ) of the function over neighboring grid points. For brevity, we use f compδ ( x )notation instead of ˆ f compδ ( x ) below. In Tests 1-3 we display the exact functions f true ( x, y ) and computed functions f comp ( x, y ) , f compδ ( x, y ) before and after the useof the post processing procedure. Test 1. Circle-shaped smooth inclusion, no scattering
The function f true is a smoothed circle of the radius r = 0 .
3, centered at(0 , h . f true is depicted on the Figure 6.1(a). µ s ( x ) ≡ , µ a ( x ) = (cid:26) . x + y ) < . , Test 2. X-shaped non-smooth inclusion, uniform scattering
The function f true is depicted on the Figure 6.2 (a). We define absorptionand scattering coefficients as µ a ( x ) = (cid:26) . x + y ) < . , µ s ( x , α ) = (cid:26) .
01 if ( x + y ) < . , ALEXEY V. SMIRNOV, MICHAEL V. KLIBANOV AND LOC H. NGUYEN(a) The function f true (b) The reconstructed func-tion f comp , no noise (c) The reconstructed func-tion f comp δ , noise 90%(d) The incomplete tomo-graphic data, no noise (e) Post-processed f comp , nonoise (f) Post-processed f comp δ ,noise level 90% Fig. 6.1:
Test 1. The true and reconstructed source functions for circle-shaped smoothinclusion without scattering. and the constant scattering phase function K ( x , α, β ) = 1 / (2 d ) . The recon-struction is displayed in Figure 6.2.3.
Test 3. Y-shaped inclusion, strongly forward-peaked scattering
The function f true depicted on Figure 6.3(a) is smooth and compactly sup-ported in the domain. We define absorption and scattering coefficients as µ a ( x ) = .
15 if ( x, y ) ∈ supp( f true ) , . x + y ) < . , µ s ( x ) = (cid:26) .
01 if ( x + y ) < . , K ( x , α, β ) = H ( α, β ) = 12 d (cid:20) − g g − g cos ( α − β ) (cid:21) , where H ( α, β ) is the 2-dimensional Henyey-Greenstein function. In (6.6), g = g ( x, y ) is the Henyey-Greenstein factor, it is the smoothed version of thefunction ˜ g ( x, y ) = (cid:26) . x + y ) < . , . . The numerical solution for Test 3 is depicted on Figure 6.3.
NVERSE PROBLEM FOR RADIATIVE TRANSFER EQUATION (a) The function f true (b) The reconstructed func-tion f comp , no noise (c) The reconstructed func-tion f comp δ , noise 30%(d) The incomplete tomo-graphic data, no noise (e) Post-processed f comp , nonoise (f) Post-processed f comp δ ,noise level 30% Fig. 6.2:
Test 2. The true and reconstructed source functions for for X-shaped non-smooth inclusion, for the case of uniform scattering (a) The function f true (b) The reconstructed func-tion f comp δ , noise 60% (c) Post-processed f comp δ ,noise level 60%(d) Post-processed f comp , nonoise (e) The absorption coefficient (f) The scattering coefficient Fig. 6.3:
Test 3. The true and reconstructed source functions for Y-shaped inclusionfor the case of strongly forward-peaked scattering ALEXEY V. SMIRNOV, MICHAEL V. KLIBANOV AND LOC H. NGUYEN
Remark
In Tests 1 and 2, Figures 6.1d and 6.2d are given only for theillustration purpose. They represent the incomplete Radon transform data, obtainedby the well-known Radon transform of the function f true , in which case µ a = µ s = 0 . The case when the data is available for all x and x α such that the set of lines L ( x , x α ) contains all possible lines intersecting domain Ω is considered to be the tomographicinverse problem with complete data. In contrast to this, we study the case when, thepoint source x α runs along the straight line as shown in Figure 2.1. In that scenario,the data in our problem is said to be incomplete and Figures 6.1d and 6.2d illustratethis. In Test 3 we omit the image representing incomplete Radon data, for the reasonthat it does not differ qualitatively from Figures 6.1d and 6.2d, while the distributionsof the absorption and scattering coefficients differ significantly from the previous twotests.
7. Appendix.
Theorem
Consider a rectangular domain Ω , defined in (2.1). Assume that functions f, µ a , µ s ∈ C (cid:0) R (cid:1) , f ∈ L (cid:0) R (cid:1) and the function K ∈ C (cid:16) R × [ − d, d ] (cid:17) . Also,assume that (2.4) holds. Then there exists a unique solution u ( x , α ) of Problem2.1 in the domain Ω such that u ( x, y, α ) ∈ C ( R × [0 , b ] × [ − d, d ]) . Furthermore, u ( x, y, α ) = 0 for y ∈ (0 , a ) as well as for sufficiently large | x | . Proof.
Let L ( x , x α ) be the segment of the straight line connecting points x and x α . Then for any appropriate function φ ( x ) (cid:90) L ( x , x α ) φ ( ξ ) ds ξ = | x − x α | y (cid:90) y φ (cid:18) α + w ( x − α ) y , w (cid:19) dw, where ds ξ is the arc length. It follows from (2.1)-(2.4) that u ( x , α ) = 0 for y ∈ (0 , a ) . Therefore, the Forward Problem in the domain Ω is equivalent to(7.1) u ( x , α ) = χ − ( x , α ) | x − x α | y y (cid:90) a χ ( z ) µ s ( z ) d (cid:90) − d K ( z , α, β ) u ( z , β ) dβdw + χ − ( x , α ) | x − x α | y y (cid:90) a ( f χ )( z ) dw, and(7.2) z ( w, x, α ) = (cid:18) α + w ( x − α ) y , w (cid:19) , ln χ ( x , α ) = | x − x α | y (cid:90) y ( µ a + µ s )( z ) dw. Estimate from the below the absolute value of the first argument in z ( w, x, α ) in (7.2).By (2.1), (2.4) and and (7.2) the left hand side of equation (7.1) is not zero only if y ∈ ( a, b ) . Since in z ( w, x, α ) we have w ∈ ( a, y ) , α ∈ ( − d, d ) , then(7.3) (cid:12)(cid:12)(cid:12)(cid:12) α + w ( x − α ) y (cid:12)(cid:12)(cid:12)(cid:12) ≥ wy | x − α | − | α | ≥ ab | x | − (cid:16) ab (cid:17) | α | ≥ ab | x | − (cid:16) ab (cid:17) d. Suppose that | x | ≥ X and X is so large that 1 − X (cid:0) ba (cid:1) d > , a b X > R.
Then(7.3) implies that | α + w ( x − α ) /y | > R. Hence, by (2.1) and (7.2) the right hand of
NVERSE PROBLEM FOR RADIATIVE TRANSFER EQUATION | x | ≥ X. Let A ( u ) : C (cid:0) R × [ − d, d ] (cid:1) → C (( | x | ≤ X ) × [ − d, d ] × [ − d, d ])be the operator in the right hand side of (7.1). Then (7.1) can be considered asthe equation u = A ( u ) with the Volterra-like integral operator, where the “Volterraproperty” is due to the integration with respect to y . Therefore, the latter equationcan be solved iteratively, as it is usually done for the Volterra integral equations. It isobvious from the above discussion that all iterates A ( u n ) ( x , α ) = 0 for | x | ≥ X. Thus,the solution of equation (7.1) in the space C (cid:0) R × [ − d, d ] (cid:1) exists and is unique. The C − smoothness of this solution with respect to x, y, α obviously follows from the wellknown convergence estimate for the iterates of a Volterra integral equation. Due tothe above mentioned equivalence, this implies uniqueness and existence of the solutionof the Forward Problem (cid:3) The numerical solution of the Forward Problem was performed via the iterativesolution of the Volterra-like integral equation (7.1) for ( x , α ) ∈ ( | x | ≤ X ) × [ a, b ] × [ − d, d ]. REFERENCES[1]
M.A. Anastasio, J. Zhang, D. Modgil and P. J. La Rivi`ere , Application of inversesource concepts to photoacoustic tomography , Inverse Problems, 23 (2007), pp. S21–S35.[2]
A.B. Bakushinsky, M.Yu. Kokurin and A. Smirnova,
Iterative Methods for Ill-Posed Prob-lems: An Introduction , de Gruyter, 2011.[3]
G. Bal and A. Tamasan , Inverse source problems in transport equations , SIAM J. Math.Anal., 39, 57-76, 2007.[4]
L. Baudouin, M. de Buhan and S. Ervedoza , Convergent algorithm based on Carlemanestimates for the recovery of a potential in the wave equation , SIAM J. Numer. Anal. 55(2017), pp. 1578–1613.[5]
L. Beilina and M.V. Klibanov , Approximate Global Convergence and Adaptivity for Coef-ficient Inverse Problems,
Springer, New York, 2012.[6]
M. Bellassoued and M. Yamamoto , Carleman Estimates and Applications to InverseProblems for Hyperbolic Systems (Japan: Springer), 2017.[7]
L. Bourgeois and J. Dard´e , A duality-based method of quasi-reversibility to solve the Cauchyproblem in the presence of noisy data, Inverse Problems, 26 (2010), 095016.[8]
L. Bourgeois, D. Ponomarev, and J. Dard´e , An inverse obstacle problem for the waveequation in a finite time domain , Inverse Problems and Imaging, 13 (2019) 377-400.[9]
A. L. Bukhgeim and M. V. Klibanov , Uniqueness in the large of a class of multidimensionalinverse problems , Soviet Math. Doklady, 17 (1981), pp. 244–247.[10]
J.P. Guillement and R.G. Novikov,
Inversion of weighted Radon transforms via finiteFourier series weight approximation , Inverse Problems in Science and Engineering, 22,787-802, 2013.[11]
J. Heino, S. Arridge, J. Sikora and E. Somersalo , Anisotropic effects in highly scatteringmedia , Physical Review E, 68, 03198, 2003.[12]
V. Isakov , Inverse Problems for Partial Differential Equations , Third Edition, Springer, NewYork, 2017.[13]
M. Jiang, T. Zhou, J. Cheng, W. Cong and G. Wang , Image reconstruction for biolumines-cence tomography from partial measurement , Opt. Express, 15 (2007), pp. 11095–11116.[14]
S. I. Kabanikhin, A. D. Satybaev, and M. A. Shishlenin , Direct Methods of Solving InverseHyperbolic Problems , VSP, Utrecht, 2005.[15]
M. V. Klibanov , Inverse problems and Carleman estimates , Inverse Problems, 8 (1992),pp. 575–596.[16]
M. V. Klibanov and A. Timonov , Carleman Estimates for Coefficient Inverse Problems andNumerical Applications , Inverse and Ill-Posed Problems Series, VSP, Utrecht, 2004.[17]
M. V. Klibanov , Carleman estimates for global uniqueness, stability and numerical methodsfor coefficient inverse problems , J. Inverse and Ill-Posed Problems, 21 (2013), pp. 477–560.[18]
M. V. Klibanov and N. T. Th`anh , Recovering of dielectric constants of explosives via aglobally strictly convex cost functional , SIAM J. Appl. Math., 75 (2015), pp. 518–537. ALEXEY V. SMIRNOV, MICHAEL V. KLIBANOV AND LOC H. NGUYEN[19]
M. V. Klibanov , Carleman estimates for the regularization of ill-posed Cauchy problems ,Applied Numerical Mathematics, 94 (2015), pp. 46–74.[20]
M. V. Klibanov , Convexification of restricted Dirichlet to Neumann map , J. Inverse andIll-Posed Problems, 25 (2017), pp. 669–685.[21]
M. V. Klibanov, A. E. Kolesov, A. Sullivan, and L. Nguyen , A new version of the con-vexification method for a 1-D coefficient inverse problem with experimental data , InverseProblems, 34 (2018), 115014.[22]
M. V. Klibanov, J. Li, and W. Zhang , Electrical impedance tomography with restricteddirichlet-to-neumann map data , Inverse Problems, 35 (2019), 35005.[23]
M. V. Klibanov and L. H. Nguyen , PDE-based numerical method for a limited angle X-raytomography , Inverse Problems, 35 (2019), 045009.[24]
A. D. Klose, V. Ntziachristos, and A. H. Hielscher , The inverse source problem basedon the radiative transfer equation in optical molecular imaging , Journal of ComputationalPhysics, 202 (2005), pp. 323 – 345.[25] R.
Lattes and J. L. Lions , The Method of Quasireversibility: Applications to Partial Differ-ential Equations , Elsevier, New York, 1969.[26]
Q. Li and L. H. Nguyen , Recovering the initial condition of parabolic equations from lateralCauchy data via the quasi-reversibility method, , preprint, arXiv:1902.07637, (2019).[27]
A. K. Louis , Incomplete data problems in x-ray computerized tomography , Numerische Math-ematik, 48 (1986), pp. 251–262.[28]
F. Natterer , The Mathematics of Computerized Tomography , SIAM, Philadelphia, 2001.[29]
F. Natterer , Inversion of the Attenuated Radon Transform ,Inverse Problems, 17, 113-119,2001.[30]
L. H. Nguyen, Q. Li, and M. V. Klibanov , A convergent numerical method for a multi-frequency inverse source problem in inhomogenous media , preprint, arXiv:1901.10047,(2019).[31]
R.G. Novikov , An inversion formula for the attenuated X-ray transformation , Ark. Mat., 40,145-167, 2002.[32]
A.A. Samarskii , The Theory of Difference Schemes , Marcel Dekker, Inc., New York.[33]
A.J. Silva Neto and M.N. ¨Ozis¸ik , An inverse problem of simultaneous estimation of radi-ation phase function, albedo and optical thickness , Journal of Quantitative Spectroscopyand Radiative Transfer, 53 (1995), pp. 397 – 409.[34]