# Generalized Abel Inversion Using Extended Hat Functions Operational Matrix

Authors: Manoj P. Tripathi; Ram K. Pandey; Vipul K. Baranwal; Om P. Singh
Journal: International Journal of Analysis
 (2013)

---

## 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:I(y)=a(y)∫αy((rμ-1ε(r))/(yμ-rμ)γ)dr+b(y)∫yβ((rμ-1ε(r))/(rμ-yμ)γ)dr, α≤y≤β, 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 a(y)=0 or b(y)=0.


---

## Body

## 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(1)I(y)=2∫y1ε(r)rr2-y2dr,0≤y≤1,
where ε(r) and I(y) represent, respectively, the emissivity and measured intensity, as measured from outside the source [3].The analytical inversion formula for (1) is given as [4]
(2)ε(r)=-1π∫r11y2-r2dI(y)d(y)dy,0≤r≤1.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:
(3)I(y)=2∫0yε(r)rr2-y2dr,0≤y≤1.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:
(4)I(y)=a(y)∫αyrμ-1ε(r)(yμ-rμ)γdr+b(y)∫yβrμ-1ε(r)(rμ-yμ)γdr,α≤y≤β,0<γ<1,
where the coefficients a(y) and b(y) do not vanish simultaneously. But the numerical inversion is still needed for its application in physical models since the experimental data for the intensity I(y) 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 n×n 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)
a(y)=0,b(y)=2,β=1,μ=2;(ii)
a(y)=2,b(y)=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 n subintervals [ih,(i+1)h], i=0,1,2,…,n-1, of equal lengths h where h=1/n. The hat function’s family of first  (n+1) hat functions is defined as follows:
(5)ψ0(t)={h-th,0≤t<h,0,otherwise,(6)ψi(t)={t-(i-1)hh,(i-1)h≤t<ih,(i+1)h-th,ih≤t<(i+1)h,i=1,2,…,n-1,0,otherwise,(7)ψn(t)={t-(1-h)h,1-h≤t≤1,0,otherwise.
We modify these functions by adding characteristics functions χ[-τ,0) and χ(1,1+τ] to ψ0 and ψn, respectively, to yield a new class of extended hat functions ψ-i(t) defined over [-τ,1+τ] for τ>0, and these are given by
(8)ψ-0(t)=χ[-τ,0)(t)+ψ0(t),(9)ψ-i(t)=ψi(t),for1≤i≤n-1,(10)ψ-n(t)=ψn(t)+χ(1,1+τ](t).
Thus, the supports of ψ-0(t) and ψ-n(t) are extended to [-τ,h] and [1-h,1+τ], respectively. These extended hat functions ψ-j(t) are continuous, linearly independent and are in L2[-τ,1+τ]. As τ→0, obviously, the extended hat functions will converge to the traditional hat functions.A functionf∈L2[0,1] may be approximated in vector form as
(11)f(t)≃∑i=0i=nfiψ-i(t)=Fn+1TΨ-n+1(t)=Ψ-n+1T(t)Fn+1,
where
(12)Fn+1≜[f0,f1,f2,…,fn]T,Ψ-n+1(t)≜[ψ-0(t),ψ-1(t),ψ-2(t),…,ψ-n(t)]T.The important aspect of using extended hat functions in the approximation of functionf(t) lies in the fact that the coefficients fi in (11) are given by
(13)fi=f(ih),i=0,1,2…,n.
Taking α=0, β=1, and μ=2 and by change of variables, the Abel integral equation (4) reduces to
(14)I(y)=a(y)2∫0yε(r)(y-r)γdr+b(y)2∫y1ε(r)(r-y)γdr,0≤y≤1,
which may be written as
(15)I1(y)=a1(y)∫0yη(r)(y-r)γdr+b1(y)∫y1η(r)(r-y)γdr,0≤y≤1,
where I1(y)=I(y), η(r)=ε(r), a1(y)=a(y)/2, and b1(y)=b(y)/2.Instead of considering (15), we consider the more general equation of the form:
(16)I1(y)=a1(y)∫-τyη(r)(y-r)γdr+b1(y)∫y1+τη(r)(r-y)γdr,0≤y≤1.
Using (11), the functions I1(y) and η(r) may be approximated as
(17)I1(y)≃Fn+1TΨ-n+1(y),η(r)≃Cn+1TΨ-n+1(r).
Thus the problem of Abel inversion is reduced to finding the unknown matrix Cn+1. Substituting (17) into (16), we get
(18)Fn+1TΨ-n+1(y)=Cn+1T[a1(y)∫-τyΨ-n+1(r)(y-r)γdr+b1(y)∫y1+τΨ-n+1(r)(r-y)γdr],0≤y≤1.
The integrals in (18) involve, evaluating integrals of the type ∫y1+τ((Ψ-n+1(r))/(r-y)γ)dr and ∫-τy((ψ-i(r))/(y-r)γ)dr. Let
(19)ΦiL(y)=∫y1+τψ-i(r)(r-y)γdr,ΦiU(y)=∫-τyψ-i(r)(y-r)γdrfori=0,1,2…,n,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ΦiL(y)∈L2[-τ,1+τ] for i=0,1,2,…,n.Proof.
We prove the theorem fori=1,2…,n-1. The proofs for i=0 and i=n are skipped as they may be proved on the same pattern. Based on subdivision of interval [-τ,1+τ], we calculate ΦiL(y) by considering the following cases.(i)
When-τ≤y<(i-1)h, then
(20)ΦiL(y)=∫y1+τψ-i(r)(r-y)γdr=∫y(i-1)hψ-i(r)(r-y)γdr+∫(i-1)hihψ-i(r)(r-y)γdr+∫ih(i+1)hψ-i(r)(r-y)γdr+∫(i+1)h1+τψ-i(r)(r-y)γdr=∫(i-1)hihψ-i(r)(r-y)γdr+∫ih(i+1)hψ-i(r)(r-y)γdr(asthesupportofψ-i(t)liesin[(i-1)h,(i+1)h])=∫(i-1)hihr-(i-1)hh(r-y)γdr+∫ih(i+1)h(i+1)h-r(r-y)γdr.
Changing the variable, r-y=t, we get
(21)ΦiL(y)=1h(1-γ)(2-γ)[((i-1)h-y)2-γ-2(ih-y)2-γ+((i+1)h-y)2-γ].(ii)
When(i-1)h≤y<ih, then
(22)ΦiL(y)=∫yihψ-i(r)(r-y)γdr+∫ih(i+1)hψ-i(r)(r-y)γdr.
Adopting the same procedure as in (i), we get
(23)ΦiL(y)=1h(1-γ)(2-γ)×[((i+1)h-y)2-γ-2(ih-y)2-γ].(iii)
Whenih≤y<(i+1)h, then
(24)ΦiL(y)=∫y(i+1)hψ-i(r)(r-y)γdr=∫y(i+1)h(i+1)h-rh(r-y)γdr=1h(1-γ)(2-γ)[((i+1)h-y)2-γ].(iv)
When(i+1)h≤y≤1+τ, then ΦiL(y)=0. Hence(25)ΦiL(y)=1h(1-γ)(2-γ){((i-1)h-y)2-γ-2(ih-y)2-γ+((i+1)h-y)2-γ,-τ≤y<(i-1)h,((i+1)h-y)2-γ-2(ih-y)2-γ,(i-1)h≤y<ih,((i+1)h-y)2-γ,ih≤y<(i+1)h,0,otherwise.Thus, from (25), we see that ∥ΦiL(y)∥2<∞ and hence ΦiL(y)∈L2[-τ,1+τ] for τ>0 and bounded. This completes the proof.Therefore, from (11), we get
(26)ΦiL(y)≃∑j=0ncijψ-j(y),fori=0,1,2…,n.Theorem 2.
The coefficientscij in (26) are given by(i)
fori=0,1,2,…,n-1, j=0,1,2…,n,
(27)cij={β,forj=i,β[(i-j+1)2-γ-2(i-j)2-γ+(i-j-1)2-γ],forj<i,0,forj>i,(ii)
fori=n, j=0,1,2,…,n,
(28)cnj={β[(n-j-1)2-γ-(n-j)2-γ+(2-γ)(n-j)1-γ]+1(1-γ)[(1+τ-jh)1-γ-(1-jh)1-γ],forj<n,1(1-γ)τ1-γ,forj=n,
whereβ=h1-γ/((1-γ)(2-γ)).Proof.
(i) Whenj=i,cii=ΦiL(ih)=β follows from (25).
Whenj<i, which is equivalent to j≤i-1, we get
(29)cij=ΦiL(jh)=∫jh1+τψ-i(r)(r-jh)γdr=(∫(i-1)hih+∫ih(i+1)h)ψ-i(r)(r-jh)γdr,
since the support of  ψ-i(t) lies in [(i-1)h,(i+1)h]. Using (6) and (9), we get from change of variable
(30)cij=β[(i-j+1)2-γ-2(i-j)2-γ+(i-j-1)2-γ].
Similarly, for j>i, that is,  j≥i+1, cij=0 (it follows trivially from (19)).
(ii) Forj<n,
(31)cnj=ΦnL(jh)=(∫jh1+∫11+τ)ψ-n(r)(r-jh)γdr.
From (10), we have
(32)cnj=∫(n-1)hnh(r-(n-1)h)h(r-jh)γdr+∫11+τ1(r-jh)γdr,
so
(33)cnj=β[(n-j-1)2-γ-2(n-j)2-γ+(2-γ)(i-j)2-γ]+1(1-γ)[(1+τ-jh)1-γ-(1-jh)1-γ].
Similarly, cnn=ΦnL(1)=∫11+τ(ψ-n(r)/(r-1)γ)=(1/(1-γ))τ1-γ, thus, proving the theorem.Similar arguments prove the following theorem.Theorem 3.
The functionsΦiU(y)∈L2[-τ,1+τ] for i=0,1,2,…,n.From Theorem3 and (11), we have
(34)ΦiU(y)≃∑j=0ndijψ-j(y).The coefficientsdij’s are given as follows.(i)
Fori=0,1,2,…,n-1; j=0,1,2,…,n,
(35)dij={β,forj=i,β[(j-i+1)2-γ-2(j-i)2-γ+(j-i-1)2-γ],forj>i,0,forj<i,
and(ii)
fori=n; j=0,1,2…,n,
(36)dnj={1(1-γ)τ1-γ,forj=0,β[(j-1)2-γ-j2-γ+(2-γ)j1-γ]+1(1-γ)[(τ+jh)1-γ-(jh)1-γ],for1≤j≤n.Using (12) and (26), the following integral may be written as
(37)∫y1+τΨ-n+1(r)(r-y)γdr=[Φ0L(y),Φ1L(y),…,ΦnL(y)]T.Substituting the approximation ofΦiL(y) from (26) in (37), we get
(38)∫y1+τΨ-n+1(r)(r-y)γdr=Pn+1LΨ-n+1(y),
where Pn+1L is a (n+1)×(n+1) matrix whose (i+1,j+1)th entry is cij, given by (27) and (28) for i, j=0,1,2,…,n. The matrix Pn+1L is given as
(39)Pn+1L=h1-γ(1-γ)(2-γ)×[10000000ξ11000000ξ2ξ1100000ξ3ξ2ξ110000ξ4ξ3ξ2ξ11000⋯⋯⋯⋯⋯100ξn-1ξn-2⋯ξ3ξ2ξ110ζnζn-1⋯ζ4ζ3ζ2ζ1ζ0](n+1)×(n+1),
where
(40)ζ0=(2-γ)h1-γτ1-γ,(41)ζs=[(s-1)2-γ-s2-γ+(2-γ)s1-γ]+(2-γ)h1-γ[(τ+sh)1-γ-(sh)1-γ],s=1,2,3,…,n,(42)ξk=(k+1)2-γ-2k2-γ+(k-1)2-γ,k=1,2,3,…,n-1.
The matrix Pn+1L is called extended hat functions lower operational matrix for Abel’s inversion.Remark 4.
It is evident from (40) that when τ=0, then ζ0=0, and so the lower triangular matrix Pn+1L becomes   a singular matrix. In this case, the singularity of the matrix  Pn+1L makes it redundant for numerical computation since the invertibility of the matrix is required to obtain the solution. To make the matrix Pn+1L 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, Pn+1U, such that
(43)∫-τyΨ-n+1(r)(y-r)γdr=Pn+1UΨ-n+1(y),
where
(44)Pn+1U=h1-γ(1-γ)(2-γ)×[ζ0ζ1ζ2ζ3ζ4⋯ζn-1ζn01ξ1ξ2ξ3⋯ξn-2ξn-1001ξ1ξ2ξ3⋯ξn-20001ξ1ξ2⋯ξn-300001ξ1⋯ξn-4⋯⋯⋯⋯⋯⋯⋯⋯0⋯⋯⋯⋯01ξ10⋯⋯⋯⋯⋯01](n+1)×(n+1).
The various entriesζ0, ζs, and ξk are given by (40)–(42), for k=1,2,3,…,n-1 and s=1,2,3,…,n.
If we partition the matrixPn+1L in four blocks as Pn+1L=[AOBζ0], then  Pn+1U=[ζ0B′OAT], where B=[ζn⋯ζ2ζ1]1×n,B′=[ζ1ζ2⋯ζn]1×n, O is a n×1 null vector, and
(45)A=[100000ξ110000ξ2ξ11000ξ3ξ2ξ1100ξn-2⋯⋯ξ110ξn-1ξn-2⋯⋯ξ11]n×n.
Using (38) and (43), (18) may be written as
(46)Fn+1TΨ-n+1(y)=Cn+1T[a(y)Pn+1U+b(y)Pn+1L]Ψ-n+1(y).
Solving the above system of linear equations, we obtain
(47)Cn+1T=Fn+1T[a(y)Pn+1U+b(y)Pn+1L]-1.
Substituting the value of Cn+1T from (47) into (17), the approximate emissivity η(r) is given by
(48)η(r)=Fn+1T[a(y)Pn+1U+b(y)Pn+1L]-1Ψ-n+1(r).


## 3. Error Analysis

In this section, an error analysis of our proposed algorithm is given.  Letf∈L2[0,1], and then, using (11), it is approximated as
(49)f(t)≃∑i=0i=nfiψ-i(t)=∑i=0i=nf(ih)ψi(t).
The above approximation gives exact values at nodal points. We denote the right-hand side of (49) by f~(t), and then, for jh≤t<(j+1)h, j=0,1,…,n-1, we have
(50)f~(t)=∑i=0i=nf(ih)ψi(t)=f(jh)ψj(t)+f((j+1)h)ψj+1(t)=f(jh)((j+1)h-th)+f(jh+h)(t-jhh)(using(6))=(j+1)f(jh)+t(f(jh+h)-f(jh))h-jf(jh+h).
Expanding f(t) in Taylor series at t=jh, we obtain
(51)f(t)=∑r=0∞(t-jh)rr!f(r)(jh).
Thus, from (50) and (51), we get
(52)|f(t)-f~(t)|=|jh(f(jh+h)-f(jh))h-t(f(jh+h)-f(jh))h+∑r=1∞(t-jh)rr!f(r)(jh)|
as h→0, and we have
(53)|f(t)-f~(t)|=|jhf′(jh)-tf′(jh)+∑r=1∞(t-jh)rr!f(r)(jh)|=|∑r=2∞(t-jh)rr!f(r)(jh)|≤(t-jh)22!+O|(t-jh)3|≤n-22+O(n-3),
as (t-jh)<n-1, thus  proving  the following theorem.Theorem 5.
The absolute error|f(t)-f~(t)| associated with the approximation (49) is of the order O(n-2).


## 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 Δη(ri), average deviation σ also known as root mean square error (RMS). They are calculated using the following equations:
(54)Δη(ri)=|ηa(ri)-ηc(ri)|,(55)σn+1={1(n+1)∑i=0n[ηa(ri)-ηc(ri)]2}1/2={1(n+1)∑i=0n[Δη(ri)]2}1/2=∥Δη∥2,
where ηc(ri) is the emission coefficient calculated at point ri using  (48)  and ηa(ri) is the exact analytical emissivity at the corresponding point. Note that σ, henceforth, denoted by σn+1 (for computational convenience) is the discrete l2 norm of the absolute error Δη denoted by ∥Δη∥2. Note that the calculation of σn+1 in (55) is performed by taking different values of n. In all the test profiles, the exact and noisy intensity profiles are denoted by I1(y) and I1δ(y), respectively, where I1δ(y) is obtained by adding a noise δ to I1(y) such that I1δ(yi)=I1(yi)+δθi, where yi=ih, i=0,1,2,…,n, nh=1, and θi is the uniform random variable with values in [-1,1] such that Max0≤i≤n|I1δ(yi)-I1(yi)|≤δ.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δk as δ0=0, δ1=σn+1, and δ2=p% of μn+1, where we mean μn+1=(1/(n+1))∑i=0nI1(y). In each of the test problems given in this section we have taken positive parameter τ=0.0001 and p=0.1, except for Example 10, where τ=0 has been used. The absolute errors between exact and approximate emissivities, corresponding to different noises δk, k=0,1,2 have been denoted by E0(r),E1(r), and E2(r), respectively. In the text boxes of the figures, the notations E0(r), E1(r), and E2(r) have been used for E0(r), E1(r), and E2(r), 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ηcδ(r) (with δ noise) and ηc0(r) (without noise) are obtained with and without noise term in the intensity profile I1(y), and using (48) these are given by
(56)ηcδ(r)=Fn+1δT[a(y)Pn+1U+b(y)Pn+1L]-1Ψ-n+1(r),ηc0(r)=Fn+1T[a(y)Pn+1U+b(y)Pn+1L]-1Ψ-n+1(r),
where Fn+1δ and Fn+1 are known matrices, and they are obtained  from the following equations:
(57)I1δ(y)=I1(y)+δθi≃Fn+1δTΨ-n+1(y),I1(y)≃Fn+1TΨ-n+1(y).
Hence
(58)ηcδ(r)-ηc0(r)=(Fn+1δT-Fn+1T)[a(y)Pn+1U+b(y)Pn+1L]-1Ψ-n+1(r).
Writing Hn+1T=Fn+1δT-Fn+1T and replacing random noise δθi by its maximum value δ, we get
(59)∥ηcδ(r)-ηc0(r)∥=∥Hn+1T∥∥[a(y)Pn+1U+b(y)Pn+1L]-1∥∥Ψ-n+1(r)∥.
Let h(r)=ηcδ(r)-ηc0(r)=(Fn+1δT-Fn+1T)[a(y)Pn+1U+b(y)Pn+1L]-1Ψ-n+1(r), then h(r) 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 reductionh(r) for  n=100 at different values of δ,τ=0,γ=1/2.
r
h
(
r
) for δ=0.01
h
(
r
) for δ=0.001
h
(
r
) for δ=0.0001
0
0.019718
0.0019718
0.00019718
0.1
0.0040696
0.00040696
0.000040696
0.2
0.0035401
0.00035401
0.000035401
0.3
0.0033118
0.00033118
0.000033118
0.4
0.0032047
0.00032047
0.000032047
0.5
0.0031726
0.00031726
0.000031726
0.6
0.0032047
0.00032047
0.000032047
0.7
0.0033118
0.00033118
0.000033118
0.8
0.0035401
0.00035401
0.000035401
0.9
0.0040696
0.00040696
0.000040696
1
0.019718
0.0019718
0.00019718Figure 1
Noise reductionh(r) for n=100, δ=0.0001, τ=0, and γ=1/2.Table1 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 r=0.5, and the maximum reduction in noise is achieved at r=0.5 for all the three levels of noises δ=0.01, 0.001, and 0.0001 introduced in I1(y). 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:(60)43y3/2(1+y)+23y3/2[2(1-y)y3+y-y2]=a(y)∫0yη(r)(y-r)dr+b(y)∫y1η(r)(r-y)dr,0≤y≤1,
where a(y)=y+1,b(y)=y2, with the exact analytical solution η(r)=r.
The absolute errorsEk(r) have been calculated for n=1000 and are given in Table 2. The value of δ1 is 2.3829×10-15, for n=1000. As δ2=5.3347×10-4>1010δ1, the absolute error E2(r) is appreciably higher than E0(r) and E1(r). The Figure 2 compares the absolute errors E0(r) and E1(r) for noise δ1=5.9827×10-16,n=100.Table 2
The absolute errorsEk(r), at different nodal points r, for n=1000, in Example 6.
r
0.0
0.2
0.4
0.6
0.8
1.0
E
0
(
r
)
0
7.494
×
10
-
16
1.2212
×
10
-
15
1.8874
×
10
-
15
8.1046
×
10
-
15
1.5543
×
10
-
15
E
1
(
r
)
1.132
×
10
-
13
2.4591
×
10
-
14
8.3822
×
10
-
15
1.7764
×
10
-
14
2.931
×
10
-
14
3.5194
×
10
-
14
E
2
(
r
)
0.011106
0.0064263
0.0025148
0.01039
0.00029049
0.00014606Figure 2
Absolute errorsE0(r) and E1(r) for n=100, δ1=5.9827×10-16.Example 7.
In this example, we consider the following Abel’s integral equation [7, 20]:
(61)∫0yη(r)(y-r)dr=y11/2,0≤y≤1,withsolutionη(r)=21011[Γ(11/2)]22πΓ(11)r5.
The absolute errors corresponding to different noises are given in Table 3. The values of various parameters are given as:
σ
3001
=
9.4002
×
1
0
-
8
(
=
δ
1
,
n
=
3000
), σ2001=2.1116×10-7(=δ1,n=2000), σ1001=8.4162×10-7(=δ1,n=1000), μ3001=0.1540, μ2001=0.1540, and μ1001=0.1542. Taking p=0.1, the various values of respective δ2 are given in order as 1.5396×10-4, 1.5402×10-4, and 1.5419×10-4.
In Figure3, the exact and reconstructed emissivities (with δ2 noise) have been shown for n=50, and the two emissivities match very well even for higher noise δ2 introduced in the intensity profile. For n=100, Figures 4 and 5 show the absolute errors E0(r), E1(r) and E0(r), E2(r), respectively.Table 3
The absolute errorsEk(r), at different nodal points r, for n=3000,2000, and 1000.
n
r
=
0.0
r
=
0.2
r
=
0.4
r
=
0.6
r
=
0.8
r
=
1.0
E
0
(
r
)
3000
0
1.9712
×
10
-
9
1.5849
×
10
-
8
5.361
×
10
-
8
1.2724
×
10
-
7
2.4875
×
10
-
7
2000
0
4.4178
×
10
-
7
3.5563
×
10
-
8
1.2035
×
10
-
7
2.8575
×
10
-
7
5.5872
×
10
-
7
1000
0
1.7515
×
10
-
8
1.4137
×
10
-
7
4.7899
×
10
-
7
1.138
×
10
-
6
2.2262
×
10
-
6
E
1
(
r
)
3000
4.6073
×
10
-
6
1.5029
×
10
-
6
2.5454
×
10
-
6
5.3295
×
10
-
7
7.3491
×
10
-
7
2.0451
×
10
-
6
2000
1.0031
×
10
-
5
1.1159
×
10
-
6
4.5957
×
10
-
6
2.3615
×
10
-
7
4.4563
×
10
-
6
9.2247
×
10
-
7
1000
1.7061
×
10
-
5
7.7631
×
10
-
6
7.528
×
10
-
6
4.4903
×
10
-
6
1.5261
×
10
-
5
4.2806
×
10
-
6
E
2
(
r
)
3000
0.0061012
0.0011742
0.0011595
0.0039968
0.00064593
0.0029774
2000
0.0068177
0.0068177
0.00078732
0.00067568
0
.
0016351
4.2642
×
10
-
5
1000
0.0001341
0.00067286
0.0010734
0.00040813
0.00098541
0.00074723Figure 3
Exact and approximate emissivities with noiseδ2=1.6081×10-4 for n=50.Figure 4
Absolute errorsE0(r) and E1(r) for δ1=8.2435×10-5, n=100.Figure 5
Absolute errorsE0(r) and E2(r) for δ2=1.5732×10-4, n=100.Example 8.
In this example we consider the following Abel’s integral equation [22]:
(62)∫0yη(r)(y-r)dr=I1(y),0≤y≤1,
where
(63)I1(y)={43y3/2,0≤y<12,43y3/2-83(y-12)3/2,12≤y≤1.
The exact solution of the integral equation (62) is given by
(64)η(r)={r,0≤r<12,1-r,12≤r≤1.
In Table 4, the absolute errors for different noises have been shown. Various parameters used for Table 4 are σ1001=1.8576×10-15=δ1,μ1001=0.3446, and δ2=3.4462×10-4. Figure 6 shows the graph of exact and approximate emissivities η(r) (without noise) for n=50. Absolute errors E0(r) and E1(r), for n=50 and n=100, are shown in Figures 7 and 8, respectively.Table 4
The absolute errorsEk(r), at different nodal points r, for n=1000, in  Example 8.
r
0.0
0.2
0.4
0.6
0.8
1.0
E
0
(
r
)
0
4.4409
×
10
-
16
5.5511
×
10
-
17
2.7200
×
10
-
15
2.1649
×
10
-
15
3.3888
×
10
-
15
E
1
(
r
)
5.5359
×
10
-
14
7.8271
×
10
-
15
2.2093
×
10
-
14
1.1768
×
10
-
14
1.4433
×
10
-
14
1.5344
×
10
-
14
E
2
(
r
)
0.0042764
0.0018814
3.4439
×
10
-
5
0.0014343
0.0024313
0.00018694Figure 6
Exact and approximate solutions forn=50.Figure 7
Absolute errorsE0(r) and E1(r) for δ1=3.4218×10-16, n=50.Figure 8
Absolute errorsE0(r), E1(r) for δ1=5.5794×10-16, n=100.Example 9.
Consider the generalized Abel integral equation:(65)I1(y)=eysin(y)∫0yη(r)(y-r)1/3dr+e-ycos(y)∫y1η(r)(r-y)1/3dr,0<y<1,
where I1(y)=eysin(y)y5/3+e-ycos(y)(1-y)2/3(2+3y). The exact solution of (65) is η(r)=10r/9. In Figure 9, the comparison between E0(r) and E1(r) is shown, for n=100.Figure 9
Comparison betweenE0(r) and E1(r) for δ1=1.6911×10-15, n=100.Example 10.
For the pair [14, 15, 23]:
(66)η(r)=(1-r)2(1+12r),for0≤r≤1,I1(y)=[(38435y7/2-36815y5/2+403y3/2+2y)+16105(1-y)5/2(19+72y)],
with  a(y)=b(y)=1, and the various parameters are as follows:
δ
1
=
2.1111
×
1
0
-
7
(
n
=
3000
), δ1=4.7384×10-7(n=2000), δ1=1.8852×10-6(n=1000), and δ2=0.0036, for all the three chosen values   of n. 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 δ2=0.0036) whereas, in Figure 11, a comparison between E0(r) and E1(r) is shown for  δ1=1.8232×10-4, n=100.Table 5
The absolute errorsEk(r), at different nodal points r, for n=3000,2000, and 1000.
n
r
=
0.0
r
=
0.2
r
=
0.4
r
=
0.6
r
=
0.8
r
=
1.0
E
0
(
r
)
3000
5.4281
×
10
-
7
2.8998
×
10
-
7
1.5818
×
10
-
7
2.6121
×
10
-
8
1.0578
×
10
-
7
2.9684
×
10
-
7
2000
1.2088
×
10
-
6
6.5113
×
10
-
7
3.5535
×
10
-
7
5.887
×
10
-
8
2.3717
×
10
-
7
6.5839
×
10
-
7
1000
4.7364
×
10
-
6
2.5925
×
10
-
6
1.4164
e
×
10
-
6
2.3637
×
10
-
7
9.4123
×
10
-
7
2.559
×
10
-
6
E
1
(
r
)
3000
6.8261
×
10
-
6
2.1337
×
10
-
6
2.2043
×
10
-
6
1.4765
×
10
-
6
1.4642
×
10
-
6
3.421
×
10
-
6
2000
1.589
×
10
-
6
5.7191
×
10
-
6
2.8577
×
10
-
6
2.5629
×
10
-
6
4.4183
×
10
-
6
3.9074
×
10
-
8
1000
1.9156
×
10
-
5
4.4716
×
10
-
6
1.643
×
10
-
5
1.712
×
10
-
5
3.6022
×
10
-
6
3.6477
×
10
-
6
E
2
(
r
)
3000
0.01495
0.020628
0.021094
0.065933
0.011844
0.095356
2000
0.022427
0.010769
0.024946
0.014392
0.0068535
0.028831
1000
0.056845
0.033342
0.013319
0.018905
0.028086
0.045314Figure 10
Exact and approximate emissivities (without noise and with noiseδ2=0.0036) for n=50.Figure 11
Comparison betweenE0(r) and E1(r) for δ1=1.8232×10-4, n=100.Example 11.
Consider the generalized Abel integral equation (15) with γ=1/2, a(y)=(3/4)exp(y), b(y)=exp(2y)+(1/2π), for the pair η(r)=sinr and
(67)I1(y)=[exp(y)y13/2F2(1,54,74,-y24)+2π(exp(2y)+12π)(C(2-2yπ)siny+S(2-2yπ)cosy)],
where C(z) and S(z) in (67) are called Fresnel integrals. These are defined as
(68)C(z)=∫0zcos(πt22)dt,S(z)=∫0zsin(πt22)dt.
For, n=3000, different  absolute errors are given in Table 6. The various parameters for n=3000 are δ1=0.0173 and δ2=7.2830×10-4. In Figure 12, the exact and approximate emissivities (without noise) are shown for n=100.Table 6
The absolute errorsEk(r) for n=3000, for Example 11.
r
0.0
0.2
0.4
0.6
0.8
0.9
E
0
(
r
)
3.389
×
10
-
7
4.7548
×
10
-
7
7.3287
×
10
-
7
1.3451
×
10
-
6
3.7972
×
10
-
6
1.0734
×
10
-
5
E
1
(
r
)
0.38495
0.079121
0.17527
0.096373
0.40061
0.050512
E
2
(
r
)
0.0066376
0.0058214
0.0070449
0.018136
0.013987
0.0021626Figure 12
Exact and approximate emissivities (without noise) forn=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 t=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.



---

*Source: 652541-2013-05-19.xml*