Generalized Abel Inversion Using Extended Hat Functions Operational Matrix

Abstract

Abel type integral equations play a vital role in the study of compressible flows around axially symmetric bodies. The relationship between emissivity and the measured intensity, as measured from the outside cylindrically symmetric, optically thin extended radiation source, is given by this equation as well. The aim of the present paper is to propose a stable algorithm for the numerical inversion of the following generalized Abel integral equation: 𝐼(𝑦)=𝑎(𝑦)∫^𝑦_𝛼((𝑟^𝜇−1𝜀(𝑟))/(𝑦^𝜇−𝑟^𝜇)^𝛾)𝑑𝑟+𝑏(𝑦)∫^𝛽_𝑦((𝑟^𝜇−1𝜀(𝑟))/(𝑟^𝜇−𝑦^𝜇)^𝛾)𝑑𝑟, 𝛼≤𝑦≤𝛽, 0<𝛾<1, using our newly constructed extended hat functions operational matrix of integration, and give an error analysis of the algorithm. The earlier numerical inversions available for the above equation assumed either 𝑎(𝑦)=0 or 𝑏(𝑦)=0.

1. Introduction

Abel integral equation [1] occurs in many branches of science and technology, such as plasma diagnostics and flame studies, where the most common problem of deduction of radial distributions of some important physical quantity from measurement of line-of-sight projected values is encountered. For a cylindrically symmetric, optically thin plasma source, the relation between radial distribution of the emission coefficient and the intensity measured from outside of the radial source is described by Abel transform. The challenging task of reconstruction of emission coefficient from its projection is known as Abel inversion. The earliest application, due to Mach [2], arose in the study of compressible flows around axially symmetric bodies.

The Abel integral equation is given by 𝐼𝑦=2¹_𝑦𝜀(𝑟)𝑟𝑟²−𝑦²𝑑𝑟,0≤𝑦≤1, (1) where 𝜀(𝑟) and 𝐼(𝑦) represent, respectively, the emissivity and measured intensity, as measured from outside the source [3].

The analytical inversion formula for (1) is given as [4] 𝜀(𝑟)1=−𝜋¹_𝑟1𝑦²−𝑟²𝑑𝐼𝑦𝑑𝑦𝑑𝑦,0≤𝑟≤1. (2)

There are several analytic and numerical inversion formulae available in the literature [1, 5–20].

Singh et al. [19] constructed an operational matrix of integration based on orthonormal Bernstein polynomials and used it to propose a stable algorithm to invert the following form of Abel integral equation: 𝐼𝑦=2^𝑦₀𝜀(𝑟)𝑟𝑟²−𝑦²𝑑𝑟,0≤𝑦≤1. (3)

In 2010, Singh et al. [20] constructed yet another operational matrix of integration based on orthonormal Bernstein polynomials and used it to propose an algorithm to invert the Abel integral equation (1).

In 2008, Chakrabarti [21] employed a direct function theoretic method to determine the closed form solution of the following generalized Abel integral equation: 𝐼𝑦=𝑎𝑦^𝑦_𝛼𝑟^𝜇−1𝜀(𝑟)𝑦^𝜇−𝑟^𝜇𝛾𝑑𝑟+𝑏𝑦^𝛽_𝑦𝑟^𝜇−1𝜀(𝑟)𝑟^𝜇−𝑦^𝜇𝛾𝑑𝑟,𝛼≤𝑦≤𝛽,0<𝛾<1, (4) where the coefficients 𝑎(𝑦) and 𝑏(𝑦) do not vanish simultaneously. But the numerical inversion is still needed for its application in physical models since the experimental data for the intensity 𝐼(𝑦) is available only at a discrete set of points, and it may also be distorted by the noise.

This motivated us to look for a stable algorithm which can be used for numerical inversion of the Abel integral equation (4) obtained by joining the two integrals (1) and (3). In this paper, we construct extended hat functions operational matrix of integration to invert the generalized Abel integral equation (4). Using hat functions for approximation of emissivity and intensity profiles has an edge over the earlier works of Singh et al. [19, 20], where they have used orthonormal Bernstein polynomials to approximate those physical quantities in the sense that a general formula for 𝑛×𝑛 operational matrix of integration is obtained in the earlier case whereas no such formula is available for the latter case. In Sections 3 and 4, we give the error estimate and the stability analysis followed by numerical examples to illustrate the efficiency and stability of the proposed algorithm.

The above two forms (1) and (3) of Abel integral equations are obtained by taking 𝛾=1/2 and(i)𝑎(𝑦)=0,𝑏(𝑦)=2,𝛽=1,𝜇=2;(ii)𝑎(𝑦)=2,𝑏(𝑦)=0,𝛼=0,𝜇=1, respectively, in (4).

Mostly for 𝜇=1,2 and 𝛾=1/2 the generalized Abel integral equation models the physical problems but the integral equation for 𝜇=2 can be reduced to the case 𝜇=1, by change of variables. So we restrict ourselves to 𝜇=1 only.

2. Extended Hat Functions and Their Operational Matrices for Abel Inversion

Hat functions are defined on the domain [0,1]. These are continuous functions with shape of hats, when plotted on two-dimensional planes. The interval [0,1] is divided into 𝑛 subintervals [𝑖ℎ,(𝑖+1)ℎ], 𝑖=0,1,2,…,𝑛−1, of equal lengths ℎ where ℎ=1/𝑛. The hat function’s family of first (𝑛+1) hat functions is defined as follows: 𝜓₀(𝑡)=ℎ−𝑡ℎ𝜓,0≤𝑡<ℎ,0,otherwise,_𝑖(𝑡)=𝑡−(𝑖−1)ℎℎ,()()𝑖−1ℎ≤𝑡<𝑖ℎ,𝑖+1ℎ−𝑡ℎ()𝜓,𝑖ℎ≤𝑡<𝑖+1ℎ,𝑖=1,2,…,𝑛−1,0,otherwise,_𝑛(𝑡)=𝑡−(1−ℎ)ℎ,1−ℎ≤𝑡≤1,0,otherwise. (5)(6)(7) We modify these functions by adding characteristics functions 𝜒_[−𝜏,0) and 𝜒_(1,1+𝜏] to 𝜓₀ and 𝜓_𝑛, respectively, to yield a new class of extended hat functions 𝜓_𝑖(𝑡) defined over [−𝜏,1+𝜏] for 𝜏>0, and these are given by 𝜓₀(𝑡)=𝜒_[−𝜏,0)(𝑡)+𝜓₀(𝑡),𝜓_𝑖(𝑡)=𝜓_𝑖(𝑡),for1≤𝑖≤𝑛−1,𝜓_𝑛(𝑡)=𝜓_𝑛(𝑡)+𝜒_(1,1+𝜏](𝑡). (8)(9)(10) Thus, the supports of 𝜓₀(𝑡) and 𝜓_𝑛(𝑡) are extended to [−𝜏,ℎ] and [1−ℎ,1+𝜏], respectively. These extended hat functions 𝜓_𝑗(𝑡) are continuous, linearly independent and are in 𝐿²[−𝜏,1+𝜏]. As 𝜏→0, obviously, the extended hat functions will converge to the traditional hat functions.

A function 𝑓∈𝐿²[0,1] may be approximated in vector form as 𝑓(𝑡)≃^𝑖=𝑛_𝑖=0𝑓_𝑖𝜓_𝑖(𝑡)=𝐹^𝑇_𝑛+1Ψ_𝑛+1(𝑡)=Ψ^𝑇_𝑛+1(𝑡)𝐹_𝑛+1, (11) where 𝐹_𝑛+1≜𝑓₀,𝑓₁,𝑓₂,…,𝑓_𝑛^𝑇,Ψ_𝑛+1(𝑡)≜𝜓₀(𝑡),𝜓₁(𝑡),𝜓₂(𝑡),…,𝜓_𝑛(𝑡)^𝑇. (12)

The important aspect of using extended hat functions in the approximation of function 𝑓(𝑡) lies in the fact that the coefficients 𝑓_𝑖 in (11) are given by 𝑓_𝑖=𝑓(𝑖ℎ),𝑖=0,1,2…,𝑛. (13) Taking 𝛼=0, 𝛽=1, and 𝜇=2 and by change of variables, the Abel integral equation (4) reduces to 𝐼√𝑦=𝑎√𝑦2^𝑦₀𝜀√𝑟𝑦−𝑟^𝛾𝑏√𝑑𝑟+𝑦2¹_𝑦𝜀√𝑟𝑟−𝑦^𝛾𝑑𝑟,0≤𝑦≤1, (14) which may be written as 𝐼₁𝑦=𝑎₁𝑦^𝑦₀𝜂(𝑟)𝑦−𝑟^𝛾𝑑𝑟+𝑏₁𝑦¹_𝑦𝜂(𝑟)𝑟−𝑦^𝛾𝑑𝑟,0≤𝑦≤1, (15) where 𝐼₁(𝑦)=𝐼(√𝑦), 𝜂(𝑟)=𝜀(√𝑟), 𝑎₁(𝑦)=𝑎(√𝑦)/2, and 𝑏₁(𝑦)=𝑏(√𝑦)/2.

Instead of considering (15), we consider the more general equation of the form: 𝐼₁𝑦=𝑎₁𝑦^𝑦_−𝜏𝜂(𝑟)𝑦−𝑟^𝛾𝑑𝑟+𝑏₁𝑦_𝑦^1+𝜏𝜂(𝑟)𝑟−𝑦^𝛾𝑑𝑟,0≤𝑦≤1. (16) Using (11), the functions 𝐼₁(𝑦) and 𝜂(𝑟) may be approximated as 𝐼₁𝑦≃𝐹^𝑇_𝑛+1Ψ_𝑛+1𝑦,𝜂(𝑟)≃𝐶^𝑇_𝑛+1Ψ_𝑛+1(𝑟). (17) Thus the problem of Abel inversion is reduced to finding the unknown matrix 𝐶_𝑛+1. Substituting (17) into (16), we get 𝐹^𝑇_𝑛+1Ψ_𝑛+1𝑦=𝐶^𝑇_𝑛+1𝑎₁𝑦^𝑦_−𝜏Ψ_𝑛+1(𝑟)𝑦−𝑟^𝛾𝑑𝑟+𝑏₁𝑦_𝑦^1+𝜏Ψ_𝑛+1(𝑟)𝑟−𝑦^𝛾𝑑𝑟,0≤𝑦≤1. (18) The integrals in (18) involve, evaluating integrals of the type ∫_𝑦^1+𝜏((Ψ_𝑛+1(𝑟))/(𝑟−𝑦)^𝛾)𝑑𝑟 and ∫^𝑦_−𝜏((𝜓_𝑖(𝑟))/(𝑦−𝑟)^𝛾)𝑑𝑟. Let Φ^𝐿_𝑖𝑦=_𝑦^1+𝜏𝜓_𝑖(𝑟)𝑟−𝑦^𝛾𝑑𝑟,Φ^𝑈_𝑖𝑦=^𝑦_−𝜏𝜓_𝑖(𝑟)𝑦−𝑟^𝛾𝑑𝑟for𝑖=0,1,2…,𝑛, (19)

and compute the two operational matrices of integration to evaluate these integrals. The scheme of derivation of these two operational matrices is based on the following theorems.

Theorem 1. The functions Φ^𝐿_𝑖(𝑦)∈𝐿²[−𝜏,1+𝜏] for 𝑖=0,1,2,…,𝑛.

Proof. We prove the theorem for 𝑖=1,2…,𝑛−1. The proofs for 𝑖=0 and 𝑖=𝑛 are skipped as they may be proved on the same pattern. Based on subdivision of interval [−𝜏,1+𝜏], we calculate Φ^𝐿_𝑖(𝑦) by considering the following cases. (i)When −𝜏≤𝑦<(𝑖−1)ℎ, then Φ^𝐿_𝑖𝑦=_𝑦^1+𝜏𝜓_𝑖(𝑟)𝑟−𝑦^𝛾=𝑑𝑟_{𝑦(𝑖−1)ℎ}𝜓_𝑖(𝑟)𝑟−𝑦^𝛾𝑑𝑟+_()ℎ^𝑖ℎ_𝑖−1𝜓_𝑖(𝑟)𝑟−𝑦^𝛾+𝑑𝑟_{(𝑖+1)ℎ𝑖ℎ}𝜓_𝑖(𝑟)𝑟−𝑦^𝛾𝑑𝑟+_()ℎ^1+𝜏_𝑖+1𝜓_𝑖(𝑟)𝑟−𝑦^𝛾=𝑑𝑟_{()ℎ𝑖ℎ𝑖−1}𝜓_𝑖(𝑟)𝑟−𝑦^𝛾𝑑𝑟+^(𝑖+1)ℎ_𝑖ℎ𝜓_𝑖(𝑟)𝑟−𝑦^𝛾𝑑𝑟asthesupportof𝜓_𝑖(𝑡))()ℎ]=liesin[(𝑖−1ℎ,𝑖+1_{()ℎ𝑖ℎ𝑖−1}𝑟−(𝑖−1)ℎℎ𝑟−𝑦^𝛾𝑑𝑟+^(𝑖+1)ℎ_𝑖ℎ(𝑖+1)ℎ−𝑟𝑟−𝑦^𝛾𝑑𝑟. (20) Changing the variable, 𝑟−𝑦=𝑡, we get Φ^𝐿_𝑖𝑦=1ℎ)1−𝛾2−𝛾(𝑖−1ℎ−𝑦^2−𝛾−2𝑖ℎ−𝑦^2−𝛾+(𝑖+1)ℎ−𝑦^2−𝛾. (21) (ii) When (𝑖−1)ℎ≤𝑦<𝑖ℎ, then Φ^𝐿_𝑖𝑦=_𝑦^𝑖ℎ𝜓_𝑖(𝑟)𝑟−𝑦^𝛾𝑑𝑟+^(𝑖+1)ℎ_𝑖ℎ𝜓_𝑖(𝑟)𝑟−𝑦^𝛾𝑑𝑟. (22) Adopting the same procedure as in (i), we get Φ^𝐿_𝑖𝑦=1ℎ×()1−𝛾2−𝛾𝑖+1ℎ−𝑦^2−𝛾−2𝑖ℎ−𝑦^2−𝛾. (23) (iii) When 𝑖ℎ≤𝑦<(𝑖+1)ℎ, then Φ^𝐿_𝑖𝑦=_𝑦^(𝑖+1)ℎ𝜓_𝑖(𝑟)𝑟−𝑦^𝛾𝑑𝑟=_𝑦^(𝑖+1)ℎ(𝑖+1)ℎ−𝑟ℎ𝑟−𝑦^𝛾=1𝑑𝑟ℎ)1−𝛾2−𝛾(𝑖+1ℎ−𝑦^2−𝛾. (24) (iv) When (𝑖+1)ℎ≤𝑦≤1+𝜏, then Φ^𝐿_𝑖(𝑦)=0. HenceΦ^𝐿_𝑖𝑦=1ℎ()1−𝛾2−𝛾𝑖−1ℎ−𝑦^2−𝛾−2𝑖ℎ−𝑦^2−𝛾+(𝑖+1)ℎ−𝑦^2−𝛾()(),−𝜏≤𝑦<𝑖−1ℎ,𝑖+1ℎ−𝑦^2−𝛾−2𝑖ℎ−𝑦^2−𝛾,()()𝑖−1ℎ≤𝑦<𝑖ℎ,𝑖+1ℎ−𝑦^2−𝛾,𝑖ℎ≤𝑦<(𝑖+1)ℎ,0,otherwise. (25)Thus, from (25), we see that ‖Φ^𝐿_𝑖(𝑦)‖₂<∞ and hence Φ^𝐿_𝑖(𝑦)∈𝐿²[−𝜏,1+𝜏] for 𝜏>0 and bounded. This completes the proof.

Therefore, from (11), we get Φ^𝐿_𝑖𝑦≃^𝑛_𝑗=0𝑐_𝑖𝑗𝜓_𝑗𝑦,for𝑖=0,1,2…,𝑛. (26)

Theorem 2. The coefficients 𝑐_𝑖𝑗 in (26) are given by (i) for 𝑖=0,1,2,…,𝑛−1, 𝑗=0,1,2…,𝑛, 𝑐_𝑖𝑗=𝛽𝛽,𝑓𝑜𝑟𝑗=𝑖,𝑖−𝑗+1^2−𝛾−2𝑖−𝑗^2−𝛾+𝑖−𝑗−1^2−𝛾,𝑓𝑜𝑟𝑗<𝑖,0,𝑓𝑜𝑟𝑗>𝑖, (27) (ii) for 𝑖=𝑛, 𝑗=0,1,2,…,𝑛, 𝑐_𝑛𝑗=𝛽𝑛−𝑗−1^2−𝛾−𝑛−𝑗^2−𝛾+2−𝛾𝑛−𝑗^1−𝛾+11−𝛾1+𝜏−𝑗ℎ^1−𝛾−1−𝑗ℎ^1−𝛾1,𝑓𝑜𝑟𝑗<𝑛,𝜏1−𝛾^1−𝛾,𝑓𝑜𝑟𝑗=𝑛, (28)
where 𝛽=ℎ^1−𝛾/((1−𝛾)(2−𝛾)).

Proof.  (i) When 𝑗=𝑖,𝑐_𝑖𝑖=Φ^𝐿_𝑖(𝑖ℎ)=𝛽 follows from (25).
When 𝑗<𝑖, which is equivalent to 𝑗≤𝑖−1, we get 𝑐_𝑖𝑗=Φ^𝐿_𝑖𝑗ℎ=^1+𝜏_𝑗ℎ𝜓_𝑖(𝑟)𝑟−𝑗ℎ^𝛾=𝑑𝑟_{()ℎ𝑖ℎ𝑖−1}+^(𝑖+1)ℎ_𝑖ℎ𝜓_𝑖(𝑟)𝑟−𝑗ℎ^𝛾𝑑𝑟, (29) since the support of 𝜓_𝑖(𝑡) lies in [(𝑖−1)ℎ,(𝑖+1)ℎ]. Using (6) and (9), we get from change of variable 𝑐_𝑖𝑗=𝛽𝑖−𝑗+1^2−𝛾−2𝑖−𝑗^2−𝛾+𝑖−𝑗−1^2−𝛾. (30) Similarly, for 𝑗>𝑖, that is, 𝑗≥𝑖+1, 𝑐_𝑖𝑗=0 (it follows trivially from (19)).
 (ii) For 𝑗<𝑛, 𝑐_𝑛𝑗=Φ^𝐿_𝑛𝑗ℎ=¹_𝑗ℎ+₁^1+𝜏𝜓_𝑛(𝑟)𝑟−𝑗ℎ^𝛾𝑑𝑟. (31) From (10), we have 𝑐_𝑛𝑗=_()ℎ^𝑛ℎ_𝑛−1(𝑟−(𝑛−1)ℎ)ℎ𝑟−𝑗ℎ^𝛾𝑑𝑟+₁^1+𝜏1𝑟−𝑗ℎ^𝛾𝑑𝑟, (32) so 𝑐_𝑛𝑗=𝛽𝑛−𝑗−1^2−𝛾−2𝑛−𝑗^2−𝛾+2−𝛾𝑖−𝑗^2−𝛾+11−𝛾1+𝜏−𝑗ℎ^1−𝛾−1−𝑗ℎ^1−𝛾. (33) Similarly, 𝑐_𝑛𝑛=Φ^𝐿_𝑛(1)=∫₁^1+𝜏(𝜓_𝑛(𝑟)/(𝑟−1)^𝛾)=(1/(1−𝛾))𝜏^1−𝛾, thus, proving the theorem.

Similar arguments prove the following theorem.

Theorem 3. The functions Φ^𝑈_𝑖(𝑦)∈𝐿²[−𝜏,1+𝜏] for 𝑖=0,1,2,…,𝑛.

From Theorem 3 and (11), we have Φ^𝑈_𝑖𝑦≃^𝑛_𝑗=0𝑑_𝑖𝑗𝜓_𝑗𝑦. (34)

The coefficients 𝑑_𝑖𝑗’s are given as follows. (i) For 𝑖=0,1,2,…,𝑛−1; 𝑗=0,1,2,…,𝑛, 𝑑_𝑖𝑗=𝛽𝛽,for𝑗=𝑖,𝑗−𝑖+1^2−𝛾−2𝑗−𝑖^2−𝛾+𝑗−𝑖−1^2−𝛾,for𝑗>𝑖,0,for𝑗<𝑖, (35) and (ii) for 𝑖=𝑛; 𝑗=0,1,2…,𝑛, 𝑑_𝑛𝑗=1𝜏1−𝛾^1−𝛾𝛽,for𝑗=0,𝑗−1^2−𝛾−𝑗^2−𝛾+2−𝛾𝑗^1−𝛾+11−𝛾𝜏+𝑗ℎ^1−𝛾−𝑗ℎ^1−𝛾,for1≤𝑗≤𝑛. (36)

Using (12) and (26), the following integral may be written as _𝑦^1+𝜏Ψ_𝑛+1(𝑟)𝑟−𝑦^𝛾𝑑𝑟=Φ^𝐿₀𝑦,Φ^𝐿₁𝑦,…,Φ^𝐿_𝑛𝑦^𝑇. (37)

Substituting the approximation of Φ^𝐿_𝑖(𝑦) from (26) in (37), we get _𝑦^1+𝜏Ψ_𝑛+1(𝑟)𝑟−𝑦^𝛾𝑑𝑟=𝑃^𝐿_𝑛+1Ψ_𝑛+1𝑦, (38) where 𝑃^𝐿_𝑛+1 is a (𝑛+1)×(𝑛+1) matrix whose (𝑖+1,𝑗+1)th entry is 𝑐_𝑖𝑗, given by (27) and (28) for 𝑖, 𝑗=0,1,2,…,𝑛. The matrix 𝑃^𝐿_𝑛+1 is given as 𝑃^𝐿_𝑛+1=ℎ^1−𝛾×𝜉1−𝛾2−𝛾10000000₁𝜉1000000₂𝜉₁𝜉100000₃𝜉₂𝜉₁𝜉10000₄𝜉₃𝜉₂𝜉₁𝜉1000⋯⋯⋯⋯⋯100_𝑛−1𝜉_𝑛−2⋯𝜉₃𝜉₂𝜉₁𝜁10_𝑛𝜁_𝑛−1⋯𝜁₄𝜁₃𝜁₂𝜁₁𝜁_{0(𝑛+1)×(𝑛+1)}, (39) where 𝜁₀=2−𝛾ℎ^1−𝛾𝜏^1−𝛾,𝜁_𝑠=(𝑠−1)^2−𝛾−𝑠^2−𝛾+2−𝛾𝑠^1−𝛾+2−𝛾ℎ^1−𝛾(𝜏+𝑠ℎ)^1−𝛾−(𝑠ℎ)^1−𝛾𝜉,𝑠=1,2,3,…,𝑛,_𝑘=(𝑘+1)^2−𝛾−2𝑘^2−𝛾+(𝑘−1)^2−𝛾,𝑘=1,2,3,…,𝑛−1. (40)(41)(42) The matrix 𝑃^𝐿_𝑛+1 is called extended hat functions lower operational matrix for Abel’s inversion.

Remark 4. It is evident from (40) that when 𝜏=0, then 𝜁₀=0, and so the lower triangular matrix 𝑃^𝐿_𝑛+1 becomes a singular matrix. In this case, the singularity of the matrix 𝑃^𝐿_𝑛+1 makes it redundant for numerical computation since the invertibility of the matrix is required to obtain the solution. To make the matrix 𝑃^𝐿_𝑛+1 invertible, we introduced a positive parameter 𝜏 and extend the traditional hat function to the interval [−𝜏,1+𝜏].
Similarly, using (19) and (34)–(36), we construct the extended hat functions upper operational matrix for Abel’s inversion, 𝑃^𝑈_𝑛+1, such that ^𝑦_−𝜏Ψ_𝑛+1(𝑟)𝑦−𝑟^𝛾𝑑𝑟=𝑃^𝑈_𝑛+1Ψ_𝑛+1𝑦, (43) where 𝑃^𝑈_𝑛+1=ℎ^1−𝛾×𝜁1−𝛾2−𝛾₀𝜁₁𝜁₂𝜁₃𝜁₄⋯𝜁_𝑛−1𝜁_𝑛01𝜉₁𝜉₂𝜉₃⋯𝜉_𝑛−2𝜉_𝑛−1001𝜉₁𝜉₂𝜉₃⋯𝜉_𝑛−20001𝜉₁𝜉₂⋯𝜉_𝑛−300001𝜉₁⋯𝜉_𝑛−4⋯⋯⋯⋯⋯⋯⋯⋯0⋯⋯⋯⋯01𝜉₁0⋯⋯⋯⋯⋯01_{(𝑛+1)×(𝑛+1)}. (44)
The various entries 𝜁₀, 𝜁_𝑠, and 𝜉_𝑘 are given by (40)–(42), for 𝑘=1,2,3,…,𝑛−1 and 𝑠=1,2,3,…,𝑛.
If we partition the matrix 𝑃^𝐿_𝑛+1 in four blocks as 𝑃^𝐿_𝑛+1=^𝐴𝑂_𝐵𝜁0, then 𝑃^𝑈_𝑛+1=^𝜁0𝐵′_𝑂𝐴𝑇, where 𝐵=𝜁_𝑛⋯𝜁₂𝜁_11×𝑛,𝐵^′=𝜁₁𝜁₂⋯𝜁_𝑛1×𝑛, 𝑂 is a 𝑛×1 null vector, and 𝜉𝐴=100000₁𝜉10000₂𝜉₁𝜉1000₃𝜉₂𝜉₁𝜉100_𝑛−2⋯⋯𝜉₁𝜉10_𝑛−1𝜉_𝑛−2⋯⋯𝜉₁1_𝑛×𝑛. (45) Using (38) and (43), (18) may be written as 𝐹^𝑇_𝑛+1Ψ_𝑛+1𝑦=𝐶^𝑇_𝑛+1𝑎𝑦𝑃^𝑈_𝑛+1+𝑏𝑦𝑃^𝐿_𝑛+1Ψ_𝑛+1𝑦. (46) Solving the above system of linear equations, we obtain 𝐶^𝑇_𝑛+1=𝐹^𝑇_𝑛+1𝑎𝑦𝑃^𝑈_𝑛+1+𝑏𝑦𝑃^𝐿_𝑛+1⁻¹. (47) Substituting the value of 𝐶^𝑇_𝑛+1 from (47) into (17), the approximate emissivity 𝜂(𝑟) is given by 𝜂(𝑟)=𝐹^𝑇_𝑛+1𝑎𝑦𝑃^𝑈_𝑛+1+𝑏𝑦𝑃^𝐿_𝑛+1⁻¹Ψ_𝑛+1(𝑟). (48)

3. Error Analysis

In this section, an error analysis of our proposed algorithm is given. Let 𝑓∈𝐿²[0,1], and then, using (11), it is approximated as 𝑓(𝑡)≃^𝑖=𝑛_𝑖=0𝑓_𝑖𝜓_𝑖(𝑡)=^𝑖=𝑛_𝑖=0𝑓(𝑖ℎ)𝜓_𝑖(𝑡). (49) The above approximation gives exact values at nodal points. We denote the right-hand side of (49) by 𝑓(𝑡), and then, for 𝑗ℎ≤𝑡<(𝑗+1)ℎ, 𝑗=0,1,…,𝑛−1, we have 𝑓(𝑡)=^𝑖=𝑛_𝑖=0𝑓(𝑖ℎ)𝜓_𝑖(𝑡)=𝑓𝑗ℎ𝜓_𝑗(𝑡)+𝑓𝑗+1ℎ𝜓_𝑗+1(𝑡)=𝑓𝑗ℎ𝑗+1ℎ−𝑡ℎ+𝑓𝑗ℎ+ℎ𝑡−𝑗ℎℎ(6)=𝑓𝑓using𝑗+1𝑗ℎ+𝑡𝑗ℎ+ℎ−𝑓𝑗ℎℎ−𝑗𝑓𝑗ℎ+ℎ. (50) Expanding 𝑓(𝑡) in Taylor series at 𝑡=𝑗ℎ, we obtain 𝑓(𝑡)=^∞_𝑟=0𝑡−𝑗ℎ^𝑟𝑓𝑟!^(𝑟)𝑗ℎ. (51) Thus, from (50) and (51), we get 𝑓(𝑡)−𝑓(𝑡=𝑓)𝑗ℎ𝑗ℎ+ℎ−𝑓𝑗ℎℎ𝑓−𝑡𝑗ℎ+ℎ−𝑓𝑗ℎℎ+^∞_𝑟=1𝑡−𝑗ℎ^𝑟𝑓𝑟!^(𝑟)𝑗ℎ (52) as ℎ→0, and we have 𝑓(𝑡)−𝑓(𝑡)=𝑗ℎ𝑓^′𝑗ℎ−𝑡𝑓^′𝑗ℎ+^∞_𝑟=1𝑡−𝑗ℎ^𝑟𝑓𝑟!^(𝑟)=𝑗ℎ_∞𝑟=2𝑡−𝑗ℎ^𝑟𝑓𝑟!^(𝑟)≤𝑗ℎ𝑡−𝑗ℎ²2!+𝑂𝑡−𝑗ℎ³≤𝑛⁻²2+𝑂𝑛⁻³, (53) as (𝑡−𝑗ℎ)<𝑛⁻¹, thus proving the following theorem.

Theorem 5. The absolute error |𝑓(𝑡)−𝑓(𝑡)| associated with the approximation (49) is of the order 𝑂(𝑛⁻²).

4. Numerical Results and Stability Analysis

In this section, we discuss the implementation of our proposed algorithm and investigate its accuracy and stability by applying it on test functions with known analytical Abel inverse.

For, it is always desirable to test the behaviour of a numerical inversion method using simulated data for which the exact results are known, and thus making a comparison between inverted results and theoretical data is possible. We have tested our algorithm on several well-known test profiles that are commonly encountered in experimental data and widely used by researchers [7, 15, 20]. The accuracy of the proposed algorithm is demonstrated by calculating the parameters of absolute error Δ𝜂(𝑟_𝑖), average deviation 𝜎 also known as root mean square error (RMS). They are calculated using the following equations: Δ𝜂𝑟_𝑖=𝜂_𝑎𝑟_𝑖−𝜂_𝑐𝑟_𝑖,𝜎_𝑛+1=1(𝑛+1)^𝑛_𝑖=0𝜂_𝑎𝑟_𝑖−𝜂_𝑐𝑟_𝑖^21/2=1(𝑛+1)^𝑛_𝑖=0Δ𝜂𝑟_𝑖^21/2=Δ𝜂₂, (54)(55) where 𝜂_𝑐(𝑟_𝑖) is the emission coefficient calculated at point 𝑟_𝑖 using (48) and 𝜂_𝑎(𝑟_𝑖) is the exact analytical emissivity at the corresponding point. Note that 𝜎, henceforth, denoted by 𝜎_𝑛+1 (for computational convenience) is the discrete 𝑙² norm of the absolute error Δ𝜂 denoted by ‖Δ𝜂‖₂. Note that the calculation of 𝜎_𝑛+1 in (55) is performed by taking different values of 𝑛. In all the test profiles, the exact and noisy intensity profiles are denoted by 𝐼₁(𝑦) and 𝐼^𝛿₁(𝑦), respectively, where 𝐼^𝛿₁(𝑦) is obtained by adding a noise 𝛿 to 𝐼₁(𝑦) such that 𝐼^𝛿₁(𝑦_𝑖)=𝐼₁(𝑦_𝑖)+𝛿𝜃_𝑖, where 𝑦_𝑖=𝑖ℎ, 𝑖=0,1,2,…,𝑛, 𝑛ℎ=1, and 𝜃_𝑖 is the uniform random variable with values in [−1,1] such that Max_{0≤𝑖≤𝑛}|𝐼^𝛿₁(𝑦_𝑖)−𝐼₁(𝑦_𝑖)|≤𝛿.

The following test problems are solved with and without noise to illustrate the efficiency and stability of our method by choosing three different values of the noises 𝛿_𝑘 as 𝛿₀=0, 𝛿₁=𝜎_𝑛+1, and 𝛿₂=𝑝% of 𝜇_𝑛+1, where we mean 𝜇_𝑛+1=(1/(𝑛+1))∑^𝑛_𝑖=0𝐼₁(𝑦). In each of the test problems given in this section we have taken positive parameter 𝜏=0.0001 and 𝑝=0.1, except for Example 10, where 𝜏=0 has been used. The absolute errors between exact and approximate emissivities, corresponding to different noises 𝛿_𝑘, 𝑘=0,1,2 have been denoted by 𝐸₀(𝑟),𝐸₁(𝑟), and 𝐸₂(𝑟), respectively. In the text boxes of the figures, the notations 𝐸0(𝑟), 𝐸1(𝑟), and 𝐸2(𝑟) have been used for 𝐸₀(𝑟), 𝐸₁(𝑟), and 𝐸₂(𝑟), respectively.

Though the stability of the algorithm is illustrated by various numerical experiments performed in this section, we analyze it also mathematically as follows.

The reconstructed emissivities 𝜂_𝑐𝛿(𝑟) (with 𝛿 noise) and 𝜂_𝑐0(𝑟) (without noise) are obtained with and without noise term in the intensity profile 𝐼₁(𝑦), and using (48) these are given by 𝜂_𝑐𝛿(𝑟)=𝐹^𝛿𝑇_𝑛+1𝑎𝑦𝑃^𝑈_𝑛+1+𝑏𝑦𝑃^𝐿_𝑛+1⁻¹Ψ_𝑛+1(𝑟),𝜂_𝑐0(𝑟)=𝐹^𝑇_𝑛+1𝑎𝑦𝑃^𝑈_𝑛+1+𝑏𝑦𝑃^𝐿_𝑛+1⁻¹Ψ_𝑛+1(𝑟), (56) where 𝐹^𝛿_𝑛+1 and 𝐹_𝑛+1 are known matrices, and they are obtained from the following equations: 𝐼^𝛿₁𝑦=𝐼₁𝑦+𝛿𝜃_𝑖≃𝐹^𝛿𝑇_𝑛+1Ψ_𝑛+1𝑦,𝐼₁𝑦≃𝐹^𝑇_𝑛+1Ψ_𝑛+1𝑦. (57) Hence 𝜂_𝑐𝛿(𝑟)−𝜂_𝑐0(𝑟)=𝐹^𝛿𝑇_𝑛+1−𝐹^𝑇_𝑛+1𝑎𝑦𝑃^𝑈_𝑛+1+𝑏𝑦𝑃^𝐿_𝑛+1⁻¹Ψ_𝑛+1(𝑟). (58) Writing 𝐻^𝑇_𝑛+1=𝐹^𝛿𝑇_𝑛+1−𝐹^𝑇_𝑛+1 and replacing random noise 𝛿𝜃_𝑖 by its maximum value 𝛿, we get 𝜂_𝑐𝛿(𝑟)−𝜂_𝑐0(𝑟=𝐻_𝑇𝑛+1)𝑎𝑦𝑃^𝑈_𝑛+1+𝑏𝑦𝑃^𝐿_𝑛+1⁻¹Ψ_𝑛+1(𝑟). (59) Let ℎ(𝑟)=𝜂_𝑐𝛿(𝑟)−𝜂_𝑐0(𝑟)=(𝐹^𝛿𝑇_𝑛+1−𝐹^𝑇_𝑛+1)[𝑎(𝑦)𝑃^𝑈_𝑛+1+𝑏(𝑦)𝑃^𝐿_𝑛+1]⁻¹Ψ_𝑛+1(𝑟), then ℎ(𝑟) reflects the noise reduction capability of the algorithm and its values at various points, and its graph is shown in Table 1 and Figure 1, respectively.

Table 1 

Noise reduction ℎ(𝑟) for 𝑛=100 at different values of 𝛿,𝜏=0,𝛾=1/2.

Figure 1 

Noise reduction ℎ(𝑟) for 𝑛=100, 𝛿=0.0001, 𝜏=0, and 𝛾=1/2.

Table 1 demonstrates the noise filtering capability of the algorithm for three different noise outputs. From Table 1 and Figure 1 we see that noise reduction is symmetric about the point 𝑟=0.5, and the maximum reduction in noise is achieved at 𝑟=0.5 for all the three levels of noises 𝛿=0.01, 0.001, and 0.0001 introduced in 𝐼₁(𝑦). The general behaviour of the noise reduction is the same irrespective of the value of 𝛿. In the interval [0.02,0.98] the algorithm is stable, whereas the noise filtering capability decreases continuously and then jumps symmetrically in [0,0.02)∪(0.98,1].

Example 6. Consider the generalized Abel integral equation: 43𝑦^3/2+21+𝑦3𝑦^3/221−𝑦𝑦³+𝑦−𝑦²=𝑎𝑦^𝑦₀𝜂(𝑟)𝑦𝑦−𝑟𝑑𝑟+𝑏¹_𝑦𝜂(𝑟)𝑟−𝑦𝑑𝑟,0≤𝑦≤1, (60) where 𝑎(𝑦)=𝑦+1,𝑏(𝑦)=𝑦², with the exact analytical solution 𝜂(𝑟)=𝑟.
The absolute errors 𝐸_𝑘(𝑟) have been calculated for 𝑛=1000 and are given in Table 2. The value of 𝛿₁ is 2.3829×10⁻¹⁵, for 𝑛=1000. As 𝛿₂=5.3347×10⁻⁴>10¹⁰𝛿₁, the absolute error 𝐸₂(𝑟) is appreciably higher than 𝐸₀(𝑟) and 𝐸₁(𝑟). The Figure 2 compares the absolute errors 𝐸₀(𝑟) and 𝐸₁(𝑟) for noise 𝛿₁=5.9827×10⁻¹⁶,𝑛=100.

Table 2 

The absolute errors 𝐸_𝑘(𝑟), at different nodal points 𝑟, for 𝑛=1000, in Example 6.

Figure 2 

Absolute errors 𝐸0(𝑟) and 𝐸1(𝑟) for 𝑛=100, 𝛿1=5.9827×10−16.

Example 7. In this example, we consider the following Abel’s integral equation [7, 20]: ^𝑦₀𝜂(𝑟)𝑦−𝑟𝑑𝑟=𝑦^11/2(𝑟)=2,0≤𝑦≤1,withsolution𝜂¹⁰11[Γ(11/2)]²()𝑟2𝜋Γ11⁵. (61) The absolute errors corresponding to different noises are given in Table 3. The values of various parameters are given as:
𝜎₃₀₀₁=9.4002×10⁻⁸(=𝛿₁,𝑛=3000), 𝜎₂₀₀₁=2.1116×10⁻⁷(=𝛿₁,𝑛=2000), 𝜎₁₀₀₁=8.4162×10⁻⁷(=𝛿₁,𝑛=1000), 𝜇₃₀₀₁=0.1540, 𝜇₂₀₀₁=0.1540, and 𝜇₁₀₀₁=0.1542. Taking 𝑝=0.1, the various values of respective 𝛿₂ are given in order as 1.5396×10⁻⁴, 1.5402×10⁻⁴, and 1.5419×10⁻⁴.
In Figure 3, the exact and reconstructed emissivities (with 𝛿₂ noise) have been shown for 𝑛=50, and the two emissivities match very well even for higher noise 𝛿₂ introduced in the intensity profile. For 𝑛=100, Figures 4 and 5 show the absolute errors 𝐸₀(𝑟), 𝐸₁(𝑟) and 𝐸₀(𝑟), 𝐸₂(𝑟), respectively.

Table 3 

The absolute errors 𝐸_𝑘(𝑟), at different nodal points 𝑟, for 𝑛=3000,2000, and 1000.

Figure 3 

Exact and approximate emissivities with noise 𝛿2=1.6081×10−4 for 𝑛=50.

Figure 4 

Absolute errors 𝐸0(𝑟) and 𝐸1(𝑟) for 𝛿1=8.2435×10−5, 𝑛=100.

Figure 5 

Absolute errors 𝐸0(𝑟) and 𝐸2(𝑟) for 𝛿2=1.5732×10−4, 𝑛=100.

Example 8. In this example we consider the following Abel’s integral equation [22]: ^𝑦₀𝜂(𝑟)𝑦−𝑟𝑑𝑟=𝐼₁𝑦,0≤𝑦≤1, (62) where 𝐼₁𝑦=43𝑦^3/21,0≤𝑦<2,43𝑦^3/2−831𝑦−2^3/2,12≤𝑦≤1. (63) The exact solution of the integral equation (62) is given by 𝜂(𝑟)=1𝑟,0≤𝑟<2,11−𝑟,2≤𝑟≤1. (64) In Table 4, the absolute errors for different noises have been shown. Various parameters used for Table 4 are 𝜎₁₀₀₁=1.8576×10⁻¹⁵=𝛿₁,𝜇₁₀₀₁=0.3446, and 𝛿₂=3.4462×10⁻⁴. Figure 6 shows the graph of exact and approximate emissivities 𝜂(𝑟) (without noise) for 𝑛=50. Absolute errors 𝐸₀(𝑟) and 𝐸₁(𝑟), for 𝑛=50 and 𝑛=100, are shown in Figures 7 and 8, respectively.

Table 4 

The absolute errors 𝐸_𝑘(𝑟), at different nodal points 𝑟, for 𝑛=1000, in Example 8.

Figure 6 

Exact and approximate solutions for 𝑛=50.

Figure 7 

Absolute errors 𝐸0(𝑟) and 𝐸1(𝑟) for 𝛿1=3.4218×10−16, 𝑛=50.

Figure 8 

Absolute errors 𝐸0(𝑟), 𝐸1(𝑟) for 𝛿1=5.5794×10−16, 𝑛=100.

Example 9. Consider the generalized Abel integral equation: 𝐼₁𝑦=𝑒^𝑦sin𝑦^𝑦₀𝜂(𝑟)𝑦−𝑟^1/3𝑑𝑟+𝑒^−𝑦cos𝑦¹_𝑦𝜂(𝑟)𝑟−𝑦^1/3𝑑𝑟,0<𝑦<1, (65) where 𝐼₁(𝑦)=𝑒^𝑦sin(𝑦)𝑦^5/3+𝑒^−𝑦cos(𝑦)(1−𝑦)^2/3(2+3𝑦). The exact solution of (65) is 𝜂(𝑟)=10𝑟/9. In Figure 9, the comparison between 𝐸₀(𝑟) and 𝐸₁(𝑟) is shown, for 𝑛=100.

Figure 9 

Comparison between 𝐸0(𝑟) and 𝐸1(𝑟) for 𝛿1=1.6911×10−15, 𝑛=100.

Example 10. For the pair [14, 15, 23]: 𝜂(𝑟)=(1−𝑟)²()𝐼1+12𝑟,for0≤𝑟≤1,₁𝑦=384𝑦35^7/2−368𝑦15^5/2+403𝑦^3/2+2√𝑦+161051−𝑦^5/219+72𝑦, (66) with 𝑎(𝑦)=𝑏(𝑦)=1, and the various parameters are as follows:
𝛿₁=2.1111×10⁻⁷(𝑛=3000), 𝛿₁=4.7384×10⁻⁷(𝑛=2000), 𝛿₁=1.8852×10⁻⁶(𝑛=1000), and 𝛿₂=0.0036, for all the three chosen values of 𝑛. The absolute errors corresponding to different noises are given in Table 5. Figure 10 shows the exact and approximate emissivities (without noise and with noise 𝛿₂=0.0036) whereas, in Figure 11, a comparison between 𝐸₀(𝑟) and 𝐸₁(𝑟) is shown for 𝛿₁=1.8232×10⁻⁴, 𝑛=100.

Table 5 

The absolute errors 𝐸_𝑘(𝑟), at different nodal points 𝑟, for 𝑛=3000,2000, and 1000.

Figure 10 

Exact and approximate emissivities (without noise and with noise 𝛿2=0.0036) for 𝑛=50.

Figure 11 

Comparison between 𝐸0(𝑟) and 𝐸1(𝑟) for 𝛿1=1.8232×10−4, 𝑛=100.

Example 11. Consider the generalized Abel integral equation (15) with 𝛾=1/2, 𝑎(𝑦)=(3/4)exp(𝑦), 𝑏(𝑦)=exp(2𝑦)+(1/√2𝜋), for the pair 𝜂(𝑟)=sin𝑟 and 𝐼₁𝑦=exp𝑦𝑦₁^3/2𝐹₂51,4,74𝑦,−²4+√+12𝜋exp2𝑦√𝐶2𝜋2−2𝑦√𝜋sin𝑦+𝑆2−2𝑦√𝜋,cos𝑦 (67) where 𝐶(𝑧) and 𝑆(𝑧) in (67) are called Fresnel integrals. These are defined as 𝐶(𝑧)=^𝑧₀cos𝜋𝑡²2𝑑𝑡,𝑆(𝑧)=^𝑧₀sin𝜋𝑡²2𝑑𝑡. (68) For, 𝑛=3000, different absolute errors are given in Table 6. The various parameters for 𝑛=3000 are 𝛿₁=0.0173 and 𝛿₂=7.2830×10⁻⁴. In Figure 12, the exact and approximate emissivities (without noise) are shown for 𝑛=100.

Table 6 

The absolute errors 𝐸_𝑘(𝑟) for 𝑛=3000, for Example 11.

Figure 12 

Exact and approximate emissivities (without noise) for 𝑛=100.

5. Conclusions

We have constructed operational matrices of integration based on extended hat functions and used them to propose a stable algorithm for the numerical inversion of the generalized Abel integral equation. The earlier numerical inversions were restricted to a part of the general Abel’s integral equation. We have extended the hat function beyond their domain [0,1] to avoid the singularity of the matrix at 𝑡=0,1. These operational matrices are given by the general formulae (39) and (44), thus making the evaluation of these matrices of any order extremely easy whereas in case of Bernstein operational matrices no such formula was available [19, 20]. The stability with respect to the data is restored and good accuracy is obtained even for high noise levels in the data. An error analysis and stability analysis are also given.

Acknowledgments

The authors are grateful to the learned reviewer for his valuable suggestions which have led to the improvement of the paper in the present form. Also, the first author acknowledges the financial support from UGC New-Delhi, India, under Faculty Improvement Program (FIP), whereas the second and the third authors acknowledge the financial support from UGC and CSIR New-Delhi, India, respectively, under JRF schemes.