# Inversion of the Laplace Transform from the Real Axis Using an Adaptive Iterative Method

Authors: Sapto W. Indratno; Alexander G. Ramm
Journal: International Journal of Mathematics and Mathematical Sciences
 (2009)

---

## Abstract

A new method for inverting the Laplace transform from the real axis is formulated. This method is based on a quadrature formula. We assume that the unknown functionf(t) is continuous with (known) compact support. An adaptive iterative method and an
adaptive stopping rule, which yield the convergence of the approximate
solution to f(t), are proposed in this paper.


---

## Body

## 1. Introduction

Consider the Laplace transform(1.1)ℒf(p):=∫0∞e-ptf(t)dt=F(p),Rep>0,
where ℒ:X0,b→L2[0,∞),(1.2)X0,b:={f∈L2[0,∞)∣suppf⊂[0,b)},b>0.
We assume in (1.2) that f has compact support. This is not a restriction practically. Indeed, if limt→∞f(t)=0, then |f(t)|<δ for t>tδ, where δ>0 is an arbitrary small number. Therefore, one may assume that suppf⊂[0,tδ], and treat the values of f for t>tδ as noise. One may also note that if f∈L1(0,∞), then
(1.3)F(p):=∫0∞f(t)e-ptdt=∫0bf(t)e-ptdt+∫b∞f(t)e-ptdt:=F1(p)+F2(p),
and |F2(p)|≤e-bpδ, where ∫b∞|f(t)|dt≤δ. Therefore, the contribution of the “tail” fb(t) of f,
(1.4)fb(t):={0,t<b,f(t),t≥b,
can be considered as noise if b>0 is large and δ>0 is small. We assume in (1.2) that f∈L2[0,∞). One may also assume that f∈L1[0,∞), or that |f(t)|≤c1ec2t, where c1,c2 are positive constants. If the last assumption holds, then one may define the function g(t):=f(t)e-(c2+1)t. Then g(t)∈L1[0,∞), and its Laplace transform G(p)=F(p+c2+1) is known on the interval [c2+1,c2+1+b] of real axis if the Laplace transform F(p) of f(t) is known on the interval [0,b]. Therefore, our inversion methods are applicable to these more general classes of functions f as well.The operatorℒ:X0,b→L2[0,∞) is compact. Therefore, the inversion of the Laplace transform (1.1) is an ill-posed problem (see [1, 2]). Since the problem is ill-posed, a regularization method is needed to obtain a stable inversion of the Laplace transform. There are many methods to solve (1.1) stably: variational regularization, quasisolutions, and iterative regularization (see, e.g., [1–4]). In this paper we propose an adaptive iterative method based on the Dynamical Systems Method (DSM) developed in [2, 4]. Some methods have been developed earlier for the inversion of the Laplace transform (see [5–8]). In many papers the data F(p) are assumed exact and given on the complex axis. In [9] it is shown that the results of the inversion of the Laplace transform from the complex axis are more accurate than these of the inversion of the Laplace transform from the real axis. The reason is the ill-posedness of the Laplace transform inversion from the real axis. A survey regarding the methods of the Laplace transform inversion has been given in [6]. There are several types of the Laplace inversion method compared in [6]. The inversion formula for the Laplace transform that is well known(1.5)f(t)=12πi∫σ-i∞σ+i∞F(p)eptdp,σ>0,
is used in some of these methods, and then f(t) is computed by some quadrature formulas, and many of these formulas can be found in [10, 11]. Moreover, the ill-posedness of the Laplace transform inversion is not discussed in all the methods compared in [6]. The approximate f(t), obtained by these methods when the data are noisy, may differ significantly from f(t). There are some papers in which the inversion of the Laplace transform from the real axis was studied (see [9, 12–19]). In [12, 17] a method based on the Mellin transform is developed. In this method the Mellin transform of the data F(p) is calculated first and then inverted for f(t). In [13] a Fourier series method for the inversion of Laplace transform from the real axis is developed. The drawback of this method comes from the ill-conditioning of the discretized problem. It is shown in [13] that if one uses some basis functions in X0,b, the problem becomes extremely ill-conditioned if the number m of the basis functions exceeds 20. In [15] a reproducing kernel method is used for the inversion of the Laplace transform. In the numerical experiments in [15] the authors use double and multiple precision methods to obtain high accuracy inversion of the Laplace transform. The usage of the multiple precision increases the computation time significantly which is observed in [15], so this method may be not efficient in practice. A detailed description of the multiple precision technique can be found in [20, 21]. Moreover, the Laplace transform inversion with perturbed data is not discussed in [15]. In [19] the authors develop an inversion formula, based on the eigenfunction expansion for the Laplace transform. The difficulties with this method are (i) the inversion formula is not applicable when the data are noisy, (ii) even for exact data the inversion formula is not suitable for numerical implementation.The Laplace transform as an operator fromC0k into L2, where C0k={f(t)∈C[0,+∞)∣suppf⊂[0,k)}, k=const>0, L2:=L2[0,∞), is considered in [14]. The finite difference method is used in [14] to discretize the problem, where the size of the linear algebraic system obtained by this method is fixed at each iteration, so the computation time increases if one uses large linear algebraic systems. The method of choosing the size of the linear algebraic system is not given in [14]. Moreover, the inversion of the Laplace transform when the data F(p) is given only on a finite interval [0,d], d>0, is not discussed in [14].The novel points in our paper are the following:(1)
the representation of the approximation solution (2.69) of the function f(t) which depends only on the kernel of the Laplace transform,(2)
the adaptive iterative scheme (2.72) and adaptive stopping rule (2.83), which generate the regularization parameter, the discrete data Fδ(p), and the number of terms in (2.69), needed for obtaining an approximation of the unknown function f(t).We study the inversion problem using the pair of spaces(X0,b,L2[0,d]), where X0,b is defined in (1.2), develop an inversion method, which can be easily implemented numerically, and demonstrate in the numerical experiments that our method yields the results comparable in accuracy with the results, presented in the literature, for example, with the double precision results given in paper [15].The smoothness of the kernel allows one to use the compound Simpson's rule in approximating the Laplace transform. Our approach yields a representation (2.69) of the approximate inversion of the Laplace transform. A number of terms in approximation (2.69) and the regularization parameter are generated automatically by the proposed adaptive iterative method. Our iterative method is based on the iterative method proposed in [22]. The adaptive stopping rule we propose here is based on the discrepancy-type principle, established in [23, 24]. This stopping rule yields convergence of the approximation (2.69) to f(t) when the noise level δ→0.A detailed derivation of our inversion method is given in Section2. In Section 3 some results of the numerical experiments are reported. These results demonstrate the efficiency and stability of the proposed method.


## 2. Description of the Method

Letf∈X0,b. Then (1.1) can be written as(2.1)(ℒf)(p):=∫0be-ptf(t)dt=F(p),0≤p.
Let us assume that the data F(p), the Laplace transform of f, are known only for 0≤p≤d<∞. Consider the mapping ℒm:L2[0,b]→ℝm+1, where(2.2)(ℒmf)i:=∫0be-pitf(t)dt=F(pi),i=0,1,2,…,m,(2.3)pi:=ih,i=0,1,2,…,m,h:=dm,
and m is an even number which will be chosen later. Then the unknown function f(t) can be obtained from a finite-dimensional operator equation (2.2). Let(2.4)〈u,v〉Wm:=∑j=0mwj(m)ujvj,∥u∥Wm:=〈u,u〉Wm
be the inner product and norm in ℝm+1, respectively, where wj(m) are the weights of the compound Simpson's rule (see [10, page 58]), that is,(2.5)wj(m):={h3,j=0,m,4h3,j=2l-1,l=1,2,…,m2,2h3,j=2l,l=1,2,…,m-22,h=dm,
where m is an even number. Then(2.6)〈ℒmg,v〉Wm=∑j=0mwj(m)∫0be-pjtg(t)dtvj=∫0bg(t)∑j=0mwj(m)e-pjtvjdt=〈g,ℒm*v〉X0,b,
where(2.7)ℒm*v=∑j=0mwj(m)e-pjtvj,v:=(v0v1⋮vm)∈ℝm+1,(2.8)〈g,h〉X0,b:=∫0bg(t)h(t)dt.
It follows from (2.2) and (2.7) that(2.9)(ℒm*ℒmg)(t)=∑j=0mwj(m)e-pjt∫0be-pjzg(z)dz:=(T(m)g)(t),(2.10)ℒmℒm*v=(∫0be-p0t∑j=0mwj(m)e-pjtvjdt∫0be-p1t∑j=0mwj(m)e-pjtvjdt⋮∫0be-pmt∑j=0mwj(m)e-pjtvjdt):=Q(m)v,
where(2.11)(Q(m))ij:=wj(m)∫0be-(pi+pj)tdt=wj(m)1-e-b(pi+pj)pi+pj,i,j=0,1,2,…,m.From (2.7) we get Range[ℒm*]=span{wj(m)k(pj,·,0)}j=0m, whereLemma 2.1.
Letwj(m) be defined in (2.5). Then
(2.12)∑j=0mwj(m)=d,
for any even number m.Proof.
From definition (2.5) one gets
(2.13)∑j=0mwj(m)=w0(m)+wm(m)+∑j=1m/2w2j-1(m)+∑j=1(m-2)/2w2j(m)=2h3+∑j=1m/24h3+∑j=1(m-2)/22h3=2h3+2hm3+h(m-2)3=hm=dmm=d.
Lemma 2.1 is proved.Lemma 2.2.
The matrixQ(m), defined in (2.11), is positive semidefinite and self-adjoint in ℝm+1 with respect to the inner product (2.4).Proof.
Let(2.14)(Hm)ij:=∫0be-(pi+pj)tdt=1-e-b(pi+pj)pi+pj,(Dm)ij={wi(m),i=j,0,otherwise,wj(m) are defined in (2.5). Then 〈DmHmDmu,v〉ℝm+1=〈u,DmHmDmv〉ℝm+1, where
(2.15)〈u,v〉ℝm+1:=∑j=0mujvj,u,v∈ℝm+1.
We have
(2.16)〈Q(m)u,v〉Wm=∑j=0mwj(m)(Q(m)u)jvj=∑j=0m(DmHmDmu)jvj=〈DmHmDmu,v〉ℝm+1=〈u,DmHmDmv〉ℝm+1=∑j=0muj(DmHmDmv)j=∑j=0mujwj(m)(HmDmv)j=〈u,Q(m)v〉Wm.
Thus, Q(m) is self-adjoint with respect to inner product (2.4). We have
(2.17)(Hm)ij=∫0be-(pi+pj)tdt=∫0be-pite-pjtdt=〈ϕi,ϕj〉X0,b,ϕi(t):=e-pit,
where 〈·,·〉X0,b is defined in (2.8). This shows that Hm is a Gram matrix. Therefore,
(2.18)〈Hmu,u〉ℝm+1≥0,∀u∈ℝm+1.
This implies
(2.19)〈Q(m)u,u〉Wm=〈Q(m)u,Dmu〉ℝm+1=〈HmDmu,Dmu〉ℝm+1≥0.
Thus, Q(m) is a positive semidefinite and self-adjoint matrix with respect to the inner product (2.4).Lemma 2.3.
LetT(m) be defined in (2.9). Then T(m) is self-adjoint and positive semidefinite operator in X0,b with respect to inner product (2.8).Proof.
From definition (2.9) and inner product (2.8) we get
(2.20)〈T(m)g,h〉X0,b=∫0b∑j=0mwj(m)e-pjt∫0be-pjzg(z)dzh(t)dt=∫0bg(z)∑j=0mwj(m)e-pjz∫0be-pjth(t)dtdz=〈g,T(m)h〉X0,b.
Thus, T(m) is a self-adjoint operator with respect to inner product (2.8). Let us prove that T(m) is positive semidefinite. Using (2.9), (2.5), (2.4), and (2.8), one gets
(2.21)〈T(m)g,g〉X0,b=∫0b∑j=0mwj(m)e-pjt∫0be-pjzg(z)dzg(t)dt=∑j=0mwj(m)∫0be-pjzg(z)dz∫0be-pjtg(t)dt=∑j=0mwj(m)(∫0be-pjzg(z)dz)2≥0.
Lemma 2.3 is proved.(2.22)k(p,t,z):=e-p(t+z).
Let us approximate the unknown f(t) as follows:(2.23)f(t)≈∑j=0mcj(m)wj(m)e-pjt=Ta,m-1ℒm*F(m):=fm(t),
where pj are defined in (2.3), Ta,m is defined in (2.30), and cj(m) are constants obtained by solving the linear algebraic system:(2.24)(aI+Q(m))c(m)=F(m),
where Q(m) is defined in (2.10),(2.25)c(m):=(c0(m)c1(m)⋮cm(m)),F(m):=(F(p0)F(p1)⋮F(pm)).
To prove the convergence of the approximate solution f(t), we use the following estimates, which are proved in [4], so their proofs are omitted.Estimates (2.26) and (2.27) are used in proving inequality (2.88), while estimates (2.28) and (2.29) are used in the proof of Lemmas 2.9 and 2.10, respectively.Lemma 2.4.
LetT(m) and Q(m) be defined in (2.9) and (2.10), respectively. Then, for a>0, the following estimates hold:
(2.26)∥Qa,m-1ℒm∥≤12a,(2.27)a∥Qa,m-1∥≤1,(2.28)∥Ta,m-1∥≤1a,(2.29)∥Ta,m-1ℒm*∥≤12a,
where
(2.30)Qa,m:=Q(m)+aI,Ta,m:=T(m)+aI,I is the identity operator and a=const>0.Let us formulate an iterative method for obtaining the approximation solution off(t) with the exact data F(p). Consider the following iterative scheme:(2.31)un(t)=qun-1(t)+(1-q)Tan-1ℒ*F,u0(t)=0,
where ℒ* is the adjoint of the operator ℒ, that is,(2.32)(ℒ*g)(t)=∫0de-ptg(p)dp,(2.33)(Tf)(t):=(ℒ*ℒf)(t)=∫0b∫0dk(p,t,z)dpf(z)dz=∫0bf(z)t+z(1-e-d(t+z))dz,k(p,t,z) is defined in (2.22),(2.34)Ta:=aI+T,a>0,(2.35)an:=qan-1,a0>0,q∈(0,1).Lemma2.5 together with  Lemma 2.7 with g(a)=a∥Ta-1f∥ yieldsLemma 2.5.
LetTa be defined in (2.34), ℒf=F, and f⊥𝒩(ℒ), where 𝒩(ℒ) is the null space of ℒ. Then
(2.36)a∥Ta-1f∥→0asa→0.Proof.
Sincef⊥𝒩(ℒ), it follows from the spectral theorem that
(2.37)lima→0a2∥Ta-1f∥2=lima→0∫0∞a2(a+s)2d〈Esf,f〉=∥P𝒩(ℒ)f∥2=0,
where Es is the resolution of the identity corresponding to ℒ*ℒ, and P is the orthogonal projector onto 𝒩(ℒ).
Lemma2.5 is proved.Theorem 2.6.
Letℒf=F, and un be defined in (2.31). Then
(2.38)limn→∞∥f-un∥=0.Proof.
By induction we get(2.39)un=∑j=0n-1ωj(n)Taj+1-1ℒ*F,
where Ta is defined in (2.34), and
(2.40)ωj(n):=qn-j-1-qn-j.
Using the identities
(2.41)ℒf=F,Ta-1ℒ*ℒ=Ta-1(T+aI-aI)=I-aTa-1,∑j=0n-1ωj(n)=1-qn,
we get
(2.42)f-un=f-∑j=0n-1ωj(n)f+∑j=0n-1ωj(n)aj+1Taj+1-1f=qnf+∑j=0n-1ωj(n)aj+1Taj+1-1f.
Therefore,
(2.43)∥f-un∥≤qn∥f∥+∑j=0n-1ωj(n)aj+1∥Taj+1-1f∥.
To prove relation (2.38) the following lemma is needed.Lemma 2.7.
Letg(x) be a continuous function on (0,∞), c>0 and q∈(0,1) constants. If
(2.44)limx→0+g(x)=g(0):=g0,
then
(2.45)limn→∞∑j=0n-1(qn-j-1-qn-j)g(cqj+1)=g0.Proof.
Let(2.46)Fl(n)≔∑j=1l-1ωj(n)g(cqj),
where ωj(n)≔qn-j-qn+1-j. Then
(2.47)|Fn+1(n)-g0|≤|Fl(n)|+|∑j=lnωj(n)g(cqj)-g0|.
Take ϵ>0 arbitrarily small. For sufficiently large fixed l(ϵ) one can choose n(ϵ)>l(ϵ), such that
(2.48)|Fl(ϵ)(n)|≤ϵ2,∀n>n(ϵ),
because limn→∞qn=0. Fix l=l(ϵ) such that |g(cqj)-g0|≤ϵ/2 for j>l(ϵ). This is possible because of (2.44). One has
(2.49)|Fl(ϵ)(n)|≤ϵ2,n>n(ϵ)>l(ϵ),|∑j=l(ϵ)nωj(n)g(cqj)-g0|≤∑j=l(ϵ)nωj(n)|g(cqj+1)-g0|+|∑j=l(ϵ)nωj(n)-1||g0|≤ϵ2∑j=l(ϵ)nωj(n)+qn-l(ϵ)|g0|≤ϵ2+|g0|qn-l(ϵ)≤ϵ,
if n(ϵ) is sufficiently large. Here we have used the relation
(2.50)∑j=lnωj(n)=1-qn+1-l.
Since ϵ>0 is arbitrarily small, relation (2.45) follows.
Lemma2.7 is proved.(2.51)limn→∞∑j=0n-1ωj(n)aj+1∥Taj+1-1f∥=0.
This together with estimate (2.43) and condition q∈(0,1) yields relation (2.38). Theorem 2.6 is proved.Lemma2.9 leads to an adaptive iterative scheme:
(2.52)un,mn(t)=qun-1,mn-1+(1-q)Tan,mn-1ℒmn*F(mn),u0,m0(t)=0,
where q∈(0,1), an are defined in (2.35), Ta,m is defined in (2.30), Amℒ is defined in (2.2), andLemma 2.8.
LetT and T(m) be defined in (2.33) and (2.9), respectively. Then
(2.53)∥T-T(m)∥≤(2bd)554010m4.Proof.
From definitions (2.33) and (2.9) we get
(2.54)|(T-T(m))f(t)|≤∫0b|∫0dk(p,t,z)dp-∑j=0mwj(m)k(pj,t,z)||f(z)|dz≤∫0b|d5180m4maxp∈[0,d](t+z)4e-p(t+z)||f(z)|dz=∫0bd5180m4(t+z)4|f(z)|dz≤d5180m4(∫0b(t+z)8dz)1/2∥f∥X0,b=d5180m4[(t+b)9-t99]1/2∥f∥X0,b,
where the following upper bound for the error of the compound Simpson's rule was used (see [10, page 58]). For f∈C(4)[x0,x2l],x0<x2l,
(2.55)|∫x0x2lf(x)dx-h3[f0+4∑j=1lf2(j-1)+2∑j=1l-1f2j+fx2l]|≤Rl,
where
(2.56)fj:=f(xj),xj=x0+jh,j=0,1,2,…,2l,h=x2l-x02l,Rl=(x2l-x0)5180(2l)4|f(4)(ξ)|,x0<ξ<x2l.
This implies
(2.57)∥(T-T(m))f∥X0,b≤d5540m4[(2b)10-2b1010]1/2∥f∥X0,b≤(2bd)554010m4∥f∥X0,b,
so estimate (2.53) is obtained.
Lemma2.8 is proved.Lemma 2.9.
Let0<a<a0,
(2.58)m=κ(a0a)1/4,κ>0.
Then
(2.59)∥T-T(m)∥≤(2bd)554010a0κ4a,
where T and T(m) are defined in (2.33) and (2.9), respectively.Proof.
Inequality (2.59) follows from estimate (2.53) and formula (2.58).(2.60)F(m):=(F(p0)F(p1)⋮F(pm))∈ℝm+1,pj are   defined in (2.3). In the iterative scheme (2.52) we have used the finite-dimensional operator T(m) approximating the operator T. Convergence of the iterative scheme (2.52) to the solution f of the equation ℒf=F is established in the following lemma.Lemma 2.10.
Letℒf=F and un,mn be defined in (2.52). If mn are chosen by the rule
(2.61)mn=⌈[κ(a0an)1/4]⌉,an=qan-1,q∈(0,1),κ,a0>0,
where ⌈[x]⌉ is the smallest even number not less than x, then
(2.62)limn→∞∥f-un,mn∥=0.Proof.
Consider the estimate(2.63)∥f-un,mn∥≤∥f-un∥+∥un-un,mn∥:=I1(n)+I2(n),
where I1(n):=∥f-un∥ and I2(n):=∥un-un,mn∥. By  Theorem 2.6, we get I1(n)→0 as n→∞. Let us prove that limn→∞I2(n)=0. Let Un:=un-un,mn. Then, from definitions (2.31) and (2.52), we get
(2.64)Un=qUn-1+(1-q)(Tan-1ℒ*F-Tan,mn-1ℒmn*F(mn)),U0=0.
By induction we obtain
(2.65)Un=∑j=0n-1ωj(n)(Taj+1-1ℒ*F-Taj+1,mj+1-1(ℒmj+1)*F(mj+1)),
where ωj are defined in (2.40). Using the identities ℒf=F, ℒmf=F(m),
(2.66)Ta-1T=Ta-1(T+aI-aI)=I-aTa-1,Ta,m-1T(m)=Ta,m-1(T(m)+aI-aI)=I-aTa,m-1,Ta,m-1-Ta-1=Ta,m-1(T-T(m))Ta-1,
one gets
(2.67)Un=∑j=0n-1ωj(n)aj+1(Taj+1,mj+1-1-Taj+1-1)f=∑j=0n-1ωj(n)aj+1Taj+1,mj+1-1(T-T(mj+1))Taj+1-1f.
This together with the rule (2.61), estimate (2.28), and  Lemma 2.8 yields
(2.68)∥Un∥≤∑j=0n-1ωj(n)aj+1∥Taj+1,mj+1-1∥∥T-T(mj+1)∥∥Taj+1-1f∥≤(2bd)554010a0κ4∑j=0n-1ωj(n)aj+1∥Taj+1-1f∥.
Applying  Lemmas 2.5 and 2.7 with g(a)=a∥Ta-1f∥, we obtain limn→∞∥Un∥=0.
Lemma2.10 is proved.


### 2.1. Noisy Data

When the dataF(p) are noisy, the approximate solution (2.23) is written as(2.69)fmδ(t)=∑j=0mwj(m)cj(m,δ)e-pjt=Ta,m-1ℒm*Fδ(m),
where the coefficients cj(m,δ) are obtained by solving the following linear algebraic system:(2.70)Qa,mc(m,δ)=Fδ(m),Qa,m is defined in (2.30),(2.71)c(m,δ):=(c0(m,δ)c1(m,δ)⋮cm(m,δ)),Fδ(m):=(Fδ(p0)Fδ(p1)⋮Fδ(pm)),wj(m) are defined in (2.5), and pj are defined in (2.3).To get the approximation solution of the functionf(t) with the noisy data Fδ(p), we consider the following iterative scheme:(2.72)un,mnδ=qun-1,mn-1δ+(1-q)Tan,mn-1ℒmn*Fδ(mn),u0,m0δ=0,
where Ta,m is defined in (2.30), an are defined in (2.35), q∈(0,1), Fδ(m) is defined in (2.71), and mn are chosen by the rule (2.61). Let us assume that(2.73)Fδ(pj)=F(pj)+δj,0<|δj|≤δ,j=0,1,2,…,m,
where δj are random quantities generated from some statistical distributions, for example, the uniform distribution on the interval [-δ,δ], and δ is the noise level of the data F(p). It follows from assumption (2.73), definition (2.5),  Lemma 2.1, and the inner product (2.4) that(2.74)∥Fδ(m)-F(m)∥Wm2=∑j=0mwj(m)δj2≤δ2∑j=0mwj(m)=δ2d.Lemma 2.11.
Letun,mn and un,mnδ be defined in (2.52) and (2.72), respectively. Then
(2.75)∥un,mn-un,mnδ∥≤dδ2an(1-qn),q∈(0,1),
where an are defined in (2.35).Proof.
LetUnδ:=un,mn-un,mnδ. Then, from definitions (2.52) and (2.72),
(2.76)Unδ=qUn-1δ+(1-q)Tan,mn-1ℒmn*(F(mn)-Fδ(mn)),U0δ=0.
By induction we obtain
(2.77)Unδ=∑j=0n-1ωj(n)Taj+1,mj+1-1(ℒmj+1)*(F(mj+1)-Fδ(mj+1)),
where ωj(n) are defined in (2.40). Using estimates (2.74) and inequality (2.29), one gets
(2.78)∥Unδ∥≤d∑j=0n-1ωj(n)δ2aj+1≤dδ2an∑j=0mωj(n)=dδ2an(1-qn),
where ωj are defined in (2.40).
Lemma2.11 is proved.Theorem 2.12.
Suppose that conditions of  Lemma2.10 hold, and nδ satisfies the following conditions:
(2.79)limδ→0nδ=∞,limδ→0δanδ=0.
Then
(2.80)limδ→0∥f-unδ,mnδδ∥=0.Proof.
Consider the estimate(2.81)∥f-unδ,mnδδ∥≤∥f-unδ,mnδ∥+∥unδ,mnδ-unδ,mnδδ∥.
This together with  Lemma 2.11 yields
(2.82)∥f-unδ,mnδδ∥≤∥f-unδ,mnδ∥+dδ2anδ(1-qn).
Applying relations (2.79) in estimate (2.82), one gets relation (2.80).
Theorem2.12 is proved.In the following subsection we propose a stopping rule which implies relations (2.79).


### 2.2. Stopping Rule

In this subsection a stopping rule which yields relations (2.79) in  Theorem 2.12 is given. We propose the stopping rule(2.83)Gnδ,mnδ≤Cδε<Gn,mn,1≤n<nδ,C>d,ε∈(0,1),
where(2.84)Gn,mn=qGn-1,mn-1+(1-q)∥ℒmnz(mn,δ)-Fδ(mn)∥Wmn,G0,m0=0,∥·∥Wm is defined in (2.4),(2.85)z(m,δ):=∑j=0mcj(m,δ)wj(m)e-pjt,wj(m) and pj are defined in (2.5) and (2.3), respectively, and cj(m,δ) are obtained by solving linear algebraic system (2.70).We observe that(2.86)ℒmnz(mn,δ)-Fδ(mn)=Q(mn)c(mn,δ)-Fδ(mn)=Q(mn)(anI+Q(mn))-1Fδ(mn)-Fδ(mn)=(Q(mn)+anI-anI)(anI+Q(mn))-1Fδ(mn)-Fδ(mn)=-an(anI+Q(mn))-1Fδ(mn)=-anc(mn,δ).
Thus, the sequence (2.84) can be written in the following form:(2.87)Gn,mn=qGn-1,mn-1+(1-q)an∥c(mn,δ)∥Wmn,G0,m0=0,
where ∥·∥Wm is defined in (2.4), and c(m,δ) solves the linear algebraic system (2.70).It follows from estimates (2.74), (2.26), and (2.27) that(2.88)an∥c(mn,δ)∥Wmn=an∥(anI+Q(mn))-1Fδ(mn)∥Wmn≤an∥(anI+Q(mn))-1(Fδ(mn)-F(mn))∥Wmn+an∥(anI+Q(mn))-1F(mn)∥Wmn≤∥Fδ(mn)-F(mn)∥Wmn+an∥(anI+Q(mn))-1ℒmnf∥Wmn≤δd+an∥f∥X0,b.
This together with (2.87) yields(2.89)Gn,mn≤qGn-1,mn-1+(1-q)(δd+an∥f∥X0,b),
or(2.90)Gn,mn-δd≤q(Gn-1,mn-1-δd)+(1-q)an∥f∥X0,b.Lemma 2.13.
The sequence (2.87) satisfies the following estimate:
(2.91)Gn,mn-δd≤(1-q)an∥f∥X0,b1-q,
where an are defined in (2.35).Proof.
Define(2.92)Ψn:=Gn,mn-δd,ψn:=(1-q)an∥f∥X0,b.
Then estimate (2.90) can be rewritten as
(2.93)Ψn≤qΨn-1+qψn-1,
where the relation an=qan-1 was used. Let us prove estimate (2.91) by induction. For n=0 we get
(2.94)Ψ0=-δd≤(1-q)a0∥f∥X0,b1-q.
Suppose that estimate (2.91) is true for 0≤n≤k. Then
(2.95)Ψk+1≤qΨk+qψk≤q1-qψk+qψk=q1-qψk=q1-qψkψk+1ψk+1=q1-qakak+1ψk+1=11-qψk+1,
where the relation ak+1=qak was used.
Lemma2.13 is proved.Lemma 2.14.
Suppose(2.96)G1,m1>δd,
where Gn,mn are defined in (2.87). Then there exists a unique integer nδ, satisfying the stopping rule (2.83) with C>d.Proof.
From  Lemma2.13 we get the estimate
(2.97)Gn,mn≤δd+(1-q)an∥f∥X0,b1-q,
where an are defined in (2.35). Therefore,
(2.98)lim supn→∞Gn,mn≤δd,
where the relation limn→∞an=0 was used. This together with condition (2.96) yields the existence of the integer nδ. The uniqueness of the integer nδ follows from its definition.
Lemma2.14 is proved.Lemma 2.15.
Suppose that conditions of  Lemma2.14 hold and nδ is chosen by the rule (2.83), then
(2.99)limδ→0δanδ=0.Proof.
From the stopping rule (2.83) and estimate (2.97) we get
(2.100)Cδε≤Gnδ-1,mnδ-1≤δd+(1-q)anδ-1∥f∥X0,b1-q,
where C>d,ε∈(0,1). This implies
(2.101)δ(Cδε-1-d)anδ-1≤(1-q)∥f∥X0,b1-q,
so, for ε∈(0,1), and anδ=qanδ-1, one gets
(2.102)limδ→0δanδ=limδ→0δqanδ-1≤limδ→0(1-q)δ1-ε∥f∥X0,b(q-q)(C-δ1-εd)=0.
Lemma2.15 is proved.Lemma 2.16.
Consider the stopping rule (2.83), where the parameters mn are chosen by rule (2.61). If nδ is chosen by the rule (2.83) then
(2.103)limδ→0nδ=∞.Proof.
From the stopping rule (2.83) with the sequence Gn defined in (2.87) one gets
(2.104)qCδε+(1-q)anδ∥c(mnδ,δ)∥Wmnδ≤qGnδ-1,mnδ-1+(1-q)anδ∥c(mnδ,δ)∥Wmnδ=Gnδ,mnδ<Cδε,
where c(m,δ) is obtained by solving linear algebraic system (2.70). This implies
(2.105)0<anδ∥c(mnδ,δ)∥Wmnδ≤Cδε.
Thus,
(2.106)limδ→0anδ∥c(mnδ,δ)∥Wmnδ=0.
If F(m)≠0, then there exists a λ0(m)>0 such that
(2.107)Eλ0(m)(m)F(m)≠0,〈Eλ0(m)F(m),F(m)〉Wm:=ξ(m)>0,
where Es(m) is the resolution of the identity corresponding to the operator Q(m):=ℒmℒm*. Let
(2.108)hm(δ,α):=α2∥Qm,α-1Fδ(m)∥Wm2,Qm,a:=aI+Q(m).
For a fixed number a>0 we obtain
(2.109)hm(δ,a)=a2∥Qm,a-1Fδ(m)∥Wm2=∫0∞a2(a+s)2d〈Es(m)Fδ(m),Fδ(m)〉Wm≥∫0λ0(m)a2(a+s)2d〈Es(m)Fδ(m),Fδ(m)〉Wm≥a2(a+λ0)2∫0λ0(m)d〈Es(m)Fδ(m),Fδ(m)〉Wm=a2∥Eλ0(m)(m)Fδ(m)∥Wm2(a+λ0(m))2.
Since Eλ0(m) is a continuous operator, and ∥F(m)-Fδ(m)∥Wm<dδ, it follows from (2.107) that
(2.110)limδ→0〈Eλ0(m)Fδ(m),Fδ(m)〉Wm=〈Eλ0(m)F(m),F(m)〉Wm>0.
Therefore, for the fixed number a>0 we get
(2.111)hm(δ,a)≥c2>0
for all sufficiently small δ>0, where c2 is a constant which does not depend on δ. Suppose limδ→0anδ≠0. Then there exists a subsequence δj→0 as j→∞, such that
(2.112)anδj≥c1>0,0<mnδj=⌈[κ(a0anδj)1/4]⌉≤⌈[κ(a0c1)1/4]⌉:=c3<∞,κ,a0>0,
where rule (2.61) was used to obtain the parameters mnδj. This together with (2.107) and (2.111) yields
(2.113)limj→∞hmnδj(δj,anδj)≥limj→∞anδj2∥Eλ0(mnδj)(mnδj)Fδj(mnδj)∥Wmnδj2(anδj+λ0(mnδj))2≥lim infj→∞c12∥Eλ0(mnδj)(mnδj)F(mnδj)∥Wmnδj2(c1+λ0(mnδj))2>0.
This contradicts relation (2.106). Thus, limδ→0anδ=limδ→0a0qnδ=0, that is, limδ→0nδ=∞.
Lemma2.16 is proved.It follows from  Lemmas2.15 and 2.16 that the stopping rule (2.83) yields the relations (2.79). We have proved the following theorem.Theorem 2.17.
Suppose all the assumptions of  Theorem2.12 hold, mn are chosen by rule (2.61), nδ is chosen by rule (2.83), and G1,m1>Cδ, where Gn,mn are defined in (2.87), then
(2.114)limδ→0∥f-unδ,mnδδ∥=0.


### 2.3. The Algorithm

Let us formulate the algorithm for obtaining the approximate solutionfmδ:(1)
The dataFδ(p) on the interval [0,d], d>0, the support of the function f(t), and the noise level δ.(2)
Initialization: choose the parametersκ>0, a0>0, q∈(0,1), ε∈(0,1), C>d, and set u0,m0δ=0, G0=0, n=1.(3)
Iterate, starting withn=1, and stop when condition (2.120) (see below) holds,
(a)
an=a0qn,(b)
choosemn by the rule (2.61),(c)
construct the vectorFδ(mn):
(2.115)(Fδ(mn))l=Fδ(pl),pl=lh,h=dmn,l=0,1,…,m,(d)
construct the matricesHmn and Dmn:
(2.116)(Hmn)ij:=∫0be-(pi+pj)tdt=1-e-b(pi+pj)pi+pj,i,j=1,2,3,…,mn,(Dmn)ij={wi(mn),i=j,0,otherwise,
where wj(m) are defined in (2.5),(e)
solve the following linear algebraic system:(2.117)(anI+HmnDmn)c(mn,δ)=Fδ(mn),
where (c(mn,δ))i=ci(mn,δ),(f)
update the coefficientcj(mn,δ) of the approximate solution un,mnδ(t) defined in (2.69) by the iterative formula:
(2.118)un,mnδ(t)=qun-1,mn-1δ(t)+(1-q)∑j=1mnc(mn,δ)wj(mn)e-pjt,
where
(2.119)u0,m0δ(t)=0.Stop when for the first time the inequality(2.120)Gn,mn=qGn-1,mn-1+an∥c(mn,δ)∥Wmn≤Cδε
holds, and get the approximation fδ(t)=unδ,mnδδ(t) of the function f(t) by formula (2.118).


## 2.1. Noisy Data

When the dataF(p) are noisy, the approximate solution (2.23) is written as(2.69)fmδ(t)=∑j=0mwj(m)cj(m,δ)e-pjt=Ta,m-1ℒm*Fδ(m),
where the coefficients cj(m,δ) are obtained by solving the following linear algebraic system:(2.70)Qa,mc(m,δ)=Fδ(m),Qa,m is defined in (2.30),(2.71)c(m,δ):=(c0(m,δ)c1(m,δ)⋮cm(m,δ)),Fδ(m):=(Fδ(p0)Fδ(p1)⋮Fδ(pm)),wj(m) are defined in (2.5), and pj are defined in (2.3).To get the approximation solution of the functionf(t) with the noisy data Fδ(p), we consider the following iterative scheme:(2.72)un,mnδ=qun-1,mn-1δ+(1-q)Tan,mn-1ℒmn*Fδ(mn),u0,m0δ=0,
where Ta,m is defined in (2.30), an are defined in (2.35), q∈(0,1), Fδ(m) is defined in (2.71), and mn are chosen by the rule (2.61). Let us assume that(2.73)Fδ(pj)=F(pj)+δj,0<|δj|≤δ,j=0,1,2,…,m,
where δj are random quantities generated from some statistical distributions, for example, the uniform distribution on the interval [-δ,δ], and δ is the noise level of the data F(p). It follows from assumption (2.73), definition (2.5),  Lemma 2.1, and the inner product (2.4) that(2.74)∥Fδ(m)-F(m)∥Wm2=∑j=0mwj(m)δj2≤δ2∑j=0mwj(m)=δ2d.Lemma 2.11.
Letun,mn and un,mnδ be defined in (2.52) and (2.72), respectively. Then
(2.75)∥un,mn-un,mnδ∥≤dδ2an(1-qn),q∈(0,1),
where an are defined in (2.35).Proof.
LetUnδ:=un,mn-un,mnδ. Then, from definitions (2.52) and (2.72),
(2.76)Unδ=qUn-1δ+(1-q)Tan,mn-1ℒmn*(F(mn)-Fδ(mn)),U0δ=0.
By induction we obtain
(2.77)Unδ=∑j=0n-1ωj(n)Taj+1,mj+1-1(ℒmj+1)*(F(mj+1)-Fδ(mj+1)),
where ωj(n) are defined in (2.40). Using estimates (2.74) and inequality (2.29), one gets
(2.78)∥Unδ∥≤d∑j=0n-1ωj(n)δ2aj+1≤dδ2an∑j=0mωj(n)=dδ2an(1-qn),
where ωj are defined in (2.40).
Lemma2.11 is proved.Theorem 2.12.
Suppose that conditions of  Lemma2.10 hold, and nδ satisfies the following conditions:
(2.79)limδ→0nδ=∞,limδ→0δanδ=0.
Then
(2.80)limδ→0∥f-unδ,mnδδ∥=0.Proof.
Consider the estimate(2.81)∥f-unδ,mnδδ∥≤∥f-unδ,mnδ∥+∥unδ,mnδ-unδ,mnδδ∥.
This together with  Lemma 2.11 yields
(2.82)∥f-unδ,mnδδ∥≤∥f-unδ,mnδ∥+dδ2anδ(1-qn).
Applying relations (2.79) in estimate (2.82), one gets relation (2.80).
Theorem2.12 is proved.In the following subsection we propose a stopping rule which implies relations (2.79).


## 2.2. Stopping Rule

In this subsection a stopping rule which yields relations (2.79) in  Theorem 2.12 is given. We propose the stopping rule(2.83)Gnδ,mnδ≤Cδε<Gn,mn,1≤n<nδ,C>d,ε∈(0,1),
where(2.84)Gn,mn=qGn-1,mn-1+(1-q)∥ℒmnz(mn,δ)-Fδ(mn)∥Wmn,G0,m0=0,∥·∥Wm is defined in (2.4),(2.85)z(m,δ):=∑j=0mcj(m,δ)wj(m)e-pjt,wj(m) and pj are defined in (2.5) and (2.3), respectively, and cj(m,δ) are obtained by solving linear algebraic system (2.70).We observe that(2.86)ℒmnz(mn,δ)-Fδ(mn)=Q(mn)c(mn,δ)-Fδ(mn)=Q(mn)(anI+Q(mn))-1Fδ(mn)-Fδ(mn)=(Q(mn)+anI-anI)(anI+Q(mn))-1Fδ(mn)-Fδ(mn)=-an(anI+Q(mn))-1Fδ(mn)=-anc(mn,δ).
Thus, the sequence (2.84) can be written in the following form:(2.87)Gn,mn=qGn-1,mn-1+(1-q)an∥c(mn,δ)∥Wmn,G0,m0=0,
where ∥·∥Wm is defined in (2.4), and c(m,δ) solves the linear algebraic system (2.70).It follows from estimates (2.74), (2.26), and (2.27) that(2.88)an∥c(mn,δ)∥Wmn=an∥(anI+Q(mn))-1Fδ(mn)∥Wmn≤an∥(anI+Q(mn))-1(Fδ(mn)-F(mn))∥Wmn+an∥(anI+Q(mn))-1F(mn)∥Wmn≤∥Fδ(mn)-F(mn)∥Wmn+an∥(anI+Q(mn))-1ℒmnf∥Wmn≤δd+an∥f∥X0,b.
This together with (2.87) yields(2.89)Gn,mn≤qGn-1,mn-1+(1-q)(δd+an∥f∥X0,b),
or(2.90)Gn,mn-δd≤q(Gn-1,mn-1-δd)+(1-q)an∥f∥X0,b.Lemma 2.13.
The sequence (2.87) satisfies the following estimate:
(2.91)Gn,mn-δd≤(1-q)an∥f∥X0,b1-q,
where an are defined in (2.35).Proof.
Define(2.92)Ψn:=Gn,mn-δd,ψn:=(1-q)an∥f∥X0,b.
Then estimate (2.90) can be rewritten as
(2.93)Ψn≤qΨn-1+qψn-1,
where the relation an=qan-1 was used. Let us prove estimate (2.91) by induction. For n=0 we get
(2.94)Ψ0=-δd≤(1-q)a0∥f∥X0,b1-q.
Suppose that estimate (2.91) is true for 0≤n≤k. Then
(2.95)Ψk+1≤qΨk+qψk≤q1-qψk+qψk=q1-qψk=q1-qψkψk+1ψk+1=q1-qakak+1ψk+1=11-qψk+1,
where the relation ak+1=qak was used.
Lemma2.13 is proved.Lemma 2.14.
Suppose(2.96)G1,m1>δd,
where Gn,mn are defined in (2.87). Then there exists a unique integer nδ, satisfying the stopping rule (2.83) with C>d.Proof.
From  Lemma2.13 we get the estimate
(2.97)Gn,mn≤δd+(1-q)an∥f∥X0,b1-q,
where an are defined in (2.35). Therefore,
(2.98)lim supn→∞Gn,mn≤δd,
where the relation limn→∞an=0 was used. This together with condition (2.96) yields the existence of the integer nδ. The uniqueness of the integer nδ follows from its definition.
Lemma2.14 is proved.Lemma 2.15.
Suppose that conditions of  Lemma2.14 hold and nδ is chosen by the rule (2.83), then
(2.99)limδ→0δanδ=0.Proof.
From the stopping rule (2.83) and estimate (2.97) we get
(2.100)Cδε≤Gnδ-1,mnδ-1≤δd+(1-q)anδ-1∥f∥X0,b1-q,
where C>d,ε∈(0,1). This implies
(2.101)δ(Cδε-1-d)anδ-1≤(1-q)∥f∥X0,b1-q,
so, for ε∈(0,1), and anδ=qanδ-1, one gets
(2.102)limδ→0δanδ=limδ→0δqanδ-1≤limδ→0(1-q)δ1-ε∥f∥X0,b(q-q)(C-δ1-εd)=0.
Lemma2.15 is proved.Lemma 2.16.
Consider the stopping rule (2.83), where the parameters mn are chosen by rule (2.61). If nδ is chosen by the rule (2.83) then
(2.103)limδ→0nδ=∞.Proof.
From the stopping rule (2.83) with the sequence Gn defined in (2.87) one gets
(2.104)qCδε+(1-q)anδ∥c(mnδ,δ)∥Wmnδ≤qGnδ-1,mnδ-1+(1-q)anδ∥c(mnδ,δ)∥Wmnδ=Gnδ,mnδ<Cδε,
where c(m,δ) is obtained by solving linear algebraic system (2.70). This implies
(2.105)0<anδ∥c(mnδ,δ)∥Wmnδ≤Cδε.
Thus,
(2.106)limδ→0anδ∥c(mnδ,δ)∥Wmnδ=0.
If F(m)≠0, then there exists a λ0(m)>0 such that
(2.107)Eλ0(m)(m)F(m)≠0,〈Eλ0(m)F(m),F(m)〉Wm:=ξ(m)>0,
where Es(m) is the resolution of the identity corresponding to the operator Q(m):=ℒmℒm*. Let
(2.108)hm(δ,α):=α2∥Qm,α-1Fδ(m)∥Wm2,Qm,a:=aI+Q(m).
For a fixed number a>0 we obtain
(2.109)hm(δ,a)=a2∥Qm,a-1Fδ(m)∥Wm2=∫0∞a2(a+s)2d〈Es(m)Fδ(m),Fδ(m)〉Wm≥∫0λ0(m)a2(a+s)2d〈Es(m)Fδ(m),Fδ(m)〉Wm≥a2(a+λ0)2∫0λ0(m)d〈Es(m)Fδ(m),Fδ(m)〉Wm=a2∥Eλ0(m)(m)Fδ(m)∥Wm2(a+λ0(m))2.
Since Eλ0(m) is a continuous operator, and ∥F(m)-Fδ(m)∥Wm<dδ, it follows from (2.107) that
(2.110)limδ→0〈Eλ0(m)Fδ(m),Fδ(m)〉Wm=〈Eλ0(m)F(m),F(m)〉Wm>0.
Therefore, for the fixed number a>0 we get
(2.111)hm(δ,a)≥c2>0
for all sufficiently small δ>0, where c2 is a constant which does not depend on δ. Suppose limδ→0anδ≠0. Then there exists a subsequence δj→0 as j→∞, such that
(2.112)anδj≥c1>0,0<mnδj=⌈[κ(a0anδj)1/4]⌉≤⌈[κ(a0c1)1/4]⌉:=c3<∞,κ,a0>0,
where rule (2.61) was used to obtain the parameters mnδj. This together with (2.107) and (2.111) yields
(2.113)limj→∞hmnδj(δj,anδj)≥limj→∞anδj2∥Eλ0(mnδj)(mnδj)Fδj(mnδj)∥Wmnδj2(anδj+λ0(mnδj))2≥lim infj→∞c12∥Eλ0(mnδj)(mnδj)F(mnδj)∥Wmnδj2(c1+λ0(mnδj))2>0.
This contradicts relation (2.106). Thus, limδ→0anδ=limδ→0a0qnδ=0, that is, limδ→0nδ=∞.
Lemma2.16 is proved.It follows from  Lemmas2.15 and 2.16 that the stopping rule (2.83) yields the relations (2.79). We have proved the following theorem.Theorem 2.17.
Suppose all the assumptions of  Theorem2.12 hold, mn are chosen by rule (2.61), nδ is chosen by rule (2.83), and G1,m1>Cδ, where Gn,mn are defined in (2.87), then
(2.114)limδ→0∥f-unδ,mnδδ∥=0.


## 2.3. The Algorithm

Let us formulate the algorithm for obtaining the approximate solutionfmδ:(1)
The dataFδ(p) on the interval [0,d], d>0, the support of the function f(t), and the noise level δ.(2)
Initialization: choose the parametersκ>0, a0>0, q∈(0,1), ε∈(0,1), C>d, and set u0,m0δ=0, G0=0, n=1.(3)
Iterate, starting withn=1, and stop when condition (2.120) (see below) holds,
(a)
an=a0qn,(b)
choosemn by the rule (2.61),(c)
construct the vectorFδ(mn):
(2.115)(Fδ(mn))l=Fδ(pl),pl=lh,h=dmn,l=0,1,…,m,(d)
construct the matricesHmn and Dmn:
(2.116)(Hmn)ij:=∫0be-(pi+pj)tdt=1-e-b(pi+pj)pi+pj,i,j=1,2,3,…,mn,(Dmn)ij={wi(mn),i=j,0,otherwise,
where wj(m) are defined in (2.5),(e)
solve the following linear algebraic system:(2.117)(anI+HmnDmn)c(mn,δ)=Fδ(mn),
where (c(mn,δ))i=ci(mn,δ),(f)
update the coefficientcj(mn,δ) of the approximate solution un,mnδ(t) defined in (2.69) by the iterative formula:
(2.118)un,mnδ(t)=qun-1,mn-1δ(t)+(1-q)∑j=1mnc(mn,δ)wj(mn)e-pjt,
where
(2.119)u0,m0δ(t)=0.Stop when for the first time the inequality(2.120)Gn,mn=qGn-1,mn-1+an∥c(mn,δ)∥Wmn≤Cδε
holds, and get the approximation fδ(t)=unδ,mnδδ(t) of the function f(t) by formula (2.118).


## 3. Numerical Experiments

### 3.1. The Parametersκ,a0,d

From definition (2.35) and the rule (2.61) we conclude that mn→∞ as an→0. Therefore, one needs to control the value of the parameter mn so that it will not grow too fast as an decreases. The role of the parameter κ in (2.61) is to control the value of the parameter mn so that the value of the parameter mn will not be too large. Since for sufficiently small noise level δ, namely, δ∈(10-16,10-6], the regularization parameter anδ, obtained by the stopping rule (2.83), is at most O(10-9), we suggest to choose κ in the interval (0,1]. For the noise level δ∈(10-6,10-2] one can choose κ∈(1,3]. To reduce the number of iterations we suggest to choose the geometric sequence an=a0δαn, where a0∈[0.1,0.2] and α∈[0.5,0.9]. One may assume without loss of generality that b=1, because a scaling transformation reduces the integral over (0,b) to the integral over (0,1). We have assumed that the data F(p) are defined on the interval J:=[0,d]. In the case the interval J=[d1,d], 0<d1<d, the constant d in estimates (2.59), (2.74), (2.75), (2.78), (2.90), (2.91), and (2.97) are replaced with the constant d-d1. If b=1, that is, f(t)=0 for t>1, then one has to take d not too large. Indeed, if f(t)=0 for t>1, then an integration by parts yields F(p)=[f(0)-e-pf(1)]/p+O(1/p2),p→∞. If the data are noisy, and the noise level is δ, then the data becomes indistinguishable from noise for p=O(1/δ). Therefore it is useless to keep the data Fδ(p) for d>O(1/δ). In practice one may get a satisfactory accuracy of inversion by the method, proposed in this paper, when one uses the data with d∈[1,20] when δ≤10-2. In all the numerical examples we have used d=5. Given the interval [0,d], the proposed method generates automatically the discrete data Fδ(pj), j=0,1,2,…,m, over the interval [0,d] which are needed to get the approximation of the function f(t).


### 3.2. Experiments

To test the proposed method we consider some examples proposed in [5–7, 9, 12, 13, 15, 16, 19, 25]. To illustrate the numerical stability of the proposed method with respect to the noise, we use the noisy data Fδ(p) with various noise levels δ=10-2,δ=10-4, and δ=10-6. The random quantities δj in (2.73) are obtained from the uniform probability density function over the interval [-δ,δ]. In Examples 3.1–3.12 we choose the value of the parameters as follows: an=0.1qn, q=δ1/2, and d=5. The parameter κ=1 is used for the noise levels δ=10-2 and δ=10-4. When δ=10-6 we choose κ=0.3 so that the value of the parameters mn are not very large, namely, mn≤300. Therefore, the computation time for solving linear algebraic system (2.117) can be reduced significantly. We assume that the support of the function f(t) is in the interval [0,b] with b=10. In the stopping rule (2.83) the following parameters are used: C=d+0.01, ε=0.99. In Example 3.13 the function f(t)=e-t is used to test the applicability of the proposed method to functions without compact support. The results are given in Table 13 and Figure 13.For a comparison with the exact solutions we use the mean absolute error(3.1)MAE:=[∑j=1100(f(ti)-fmnδδ(ti))2100]1/2,tj=0.01+0.1(j-1),j=1,…,100,
where f(t) is the exact solution and fmnδδ(t) is the approximate solution. The computation time (CPU time) for obtaining the approximation of f(t), the number of iterations (Iter.), and the parameters mnδ and anδ generated by the proposed method is given in each experiment (see Tables 1–12). All the calculations are done in double precision generated by MATLAB.Table 1
Example3.1.
δMAEmnδIter.CPU time (second)anδ1.00×10-29.62×10-23033.13×10-22.00×10-31.00×10-45.99×10-23246.25×10-22.00×10-71.00×10-64.74×10-25453.28×10-12.00×10-10Table 2
Example3.2.
δMAEmnδIter.CPU time (seconds)anδ1.00×10-21.09×10-13023.13×10-22.00×10-31.00×10-48.47×10-23236.25×10-22.00×10-61.00×10-67.41×10-25454.38×10-12.00×10-12Table 3
Example3.3.
δMAEmnδIter.CPU time (seconds)anδ1.00×10-22.42×10-23023.13×10-22.00×10-31.00×10-41.08×10-33033.13×10-22.00×10-61.00×10-64.02×10-43044.69×10-22.00×10-9Table 4
Example3.4.
δMAEmnδIter.CPU time (seconds)anδ1.00×10-21.59×10-23023.13×10-22.00×10-31.00×10-48.26×10-43039.400×10-22.00×10-61.00×10-61.24×10-43041.250×10-12.00×10-9Table 5
Example3.5.
δMAEmnδIter.CPU time (seconds)anδ1.00×10-24.26×10-23036.300×10-22.00×10-31.00×10-41.25×10-23039.38×10-22.00×10-61.00×10-61.86×10-35443.13×10-22.00×10-9Table 6
Example3.6.
δMAEmnδIter.CPU time (seconds)anδ1.00×10-24.19×10-23024.700×10-22.00×10-31.00×10-41.64×10-23239.38×10-22.00×10-61.00×10-61.22×10-25443.13×10-22.00×10-9Table 7
Example3.7.
δMAEmnδIter.CPU time (seconds)anδ1.00×10-21.52×10-23024.600×10-22.00×10-31.00×10-42.60×10-33039.38×10-22.00×10-61.00×10-62.02×10-33043.13×10-22.00×10-9Table 8
Example3.8.
δMAEmnδIter.CPU time (seconds)anδ1.00×10-22.74×10-23021.100×10-22.00×10-31.00×10-43.58×10-33033.13×10-22.00×10-61.00×10-65.04×10-43044.69×10-22.00×10-9Table 9
Example3.9.
δMAEmnδIter.CPU time (seconds)anδ1.00×10-22.07×10-13036.25×10-22.00×10-61.00×10-47.14×10-23243.44×10-12.00×10-91.00×10-62.56×10-25453.75×10-12.00×10-12Table 10
Example3.10.
δMAEmnδIter.CPU time (seconds)anδ1.00×10-22.09×10-13033.13×10-22.00×10-61.00×10-41.35×10-23249.38×10-22.00×10-91.00×10-63.00×10-35442.66×10-12.00×10-9Table 11
Example3.11.
δMAEmnδIter.CPU time (seconds)anδ1.00×10-22.47×10-13039.38×10-22.00×10-61.00×10-44.91×10-23242.50×10-12.00×10-91.00×10-62.46×10-25454.38×10-12.00×10-12Table 12
Example3.12.
δMAEmnδIter.CPU time (seconds)anδ1.00×10-21.37×1009639.38×10-22.00×10-61.00×10-45.98×10-110042.66×10-12.00×10-91.00×10-62.24×10-130053.44×10-12.00×10-12Table 13
Example3.13.
bMAEmδIter.CPU time (seconds)51.487×10-2243.125×10-282.183×10-4243.125×10-2204.517×10-9243.125×10-2301.205×10-13243.125×10-2Example 3.1 (see [15]).
We have(3.2)f1(t)={1,12≤t≤32,0,otherwise,F1(p)={1,p=0,e-p/2-e-3p/2p,p>0.
The reconstruction of the exact solution for different values of the noise level δ is shown in Figure 1. When the noise level δ=10-6, our result is comparable with the double precision results shown in [15]. The proposed method is stable with respect to the noise δ as shown in Table 1.Figure 1
Example3.1: the stability of the approximate solution.Example 3.2 (see [13, 15]).
We have(3.3)f2(t)={12,t=1,1,1<t<10,0,elsewhere,F2(p)={9,p=0,e-p-e-10pp,p>0.
The reconstruction of the function f2(t) is plotted in Figure 2. In [15] a high accuracy result is given by means of the multiple precision. But, as reported in [15], to get such high accuracy results, it takes 7 hours. From Table 2 and Figure 2 we can see that the proposed method yields stable solution with respect to the noise level δ. The reconstruction of the exact solution obtained by the proposed method is better than the reconstruction shown in [13]. The result is comparable with the double precision results given in [15]. For δ=10-6 and κ=0.3 the value of the parameter mnδ is bounded by the constant 54.Figure 2
Example3.2: the stability of the approximate solution.Example 3.3 (see [6, 12, 13, 16, 19]).
We have(3.4)f3(t)={te-t,0≤t<10,0,otherwise,F3(p)=1-e-(p+1)10(p+1)2-10e-(p+1)10p+1.
We get an excellent agreement between the approximate solution and the exact solution when the noise level δ=10-4 and 10-6 as shown in Figure 3. The results obtained by the proposed method are better than the results given in [13]. The mean absolute error MAE decreases as the noise level decreases which shows the stability of the proposed method. Our results are more stable with respect to the noise δ than the results presented in [19]. The value of the parameter mnδ is bounded by the constant 30 when the noise level δ=10-6 and κ=0.3.Figure 3
Example3.3: the stability of the approximate solution.Figure 4
Example3.4: the stability of the approximate solution.Example 3.4 (see [13, 15]).
We have(3.5)f4(t)={1-e-0.5t,0≤t<10,0,elsewhere,F4(p)={8+2e-5,p=0,1-e-10pp-1-e-(p+1/2)10p+0.5,p>0.
As in our Example 3.3 when the noise δ=10-4 and 10-6 are used, we get a satisfactory agreement between the approximate solution and the exact solution. Table 4 gives the results of the stability of the proposed method with respect to the noise level δ. Moreover, the reconstruction of the function f4(t) obtained by the proposed method is better than the reconstruction of f4(t) shown in [13], and is comparable with the double precision reconstruction obtained in [15].
In this example whenδ=10-6 and κ=0.3 the value of the parameter mnδ is bounded by the constant 109 as shown in Table 4.Example 3.5 (see [5, 7, 13]).
We have(3.6)f5(t)=23e-t/2sin(t3/2)F5(p)=1-cos(103/2)e-10(p+0.5)[(p+0.5)2+3/4]-2(p+0.5)e-10(p+0.5)sin(103/2)3[(p+0.5)2+3/4].
This is an example of the damped sine function. In [5, 7] the knowledge of the exact data F(p) in the complex plane is required to get the approximate solution. Here we only use the knowledge of the discrete perturbed data Fδ(pj), j=0,1,2,…,m, and get a satisfactory result which is comparable with the results given in [5, 7] when the level noise δ=10-6. The reconstruction of the exact solution f5(t) obtained by our method is better than this of the method given in [13]. Moreover, our method yields stable solution with respect to the noise level δ as shown in Figure 5 and Table 5. In this example when κ=0.3 the value of the parameter mnδ is bounded by 54 for the noise level δ=10-6 (see Table 5).Figure 5
Example3.5: the stability of the approximate solution.Example 3.6 (see [15]).
We have(3.7)f6(t)={t,0≤t<1,32-t2,1≤t<3,0,elsewhere,F6(p)={32,p=0,1-e-p(1+p)p2+e-3p+e2p(2p-1)2p2,p>0.
Example 3.6 represents a class of piecewise continuous functions. From Figure 6 the value of the exact solution at the points where the function is not differentiable cannot be well approximated for the given levels of noise by the proposed method. When the noise level δ=10-6, our result is comparable with the results given in [15]. Table 6 reports the stability of the proposed method with respect to the noise δ. It is shown in Table 6 that the value of the parameter m generated by the proposed adaptive stopping rule is bounded by the constant 54 for the noise level δ=10-6 and κ=0.3 which gives a relatively small computation time.Figure 6
Example3.6: the stability of the approximate solution.Example 3.7 (see [15]).
We have(3.8)f7(t)={-te-t-e-t+1,0≤t<1,1-2e-1,1≤t<10,0,elsewhere,F7(p)={3e-1+9(1-2e),p=0,e-1-p(e1+p-e(1+p)2+p(3+2p)/p(p+1)2)+(e-2)e-1-p-10p((e10p-ep)/p),p>0.
When the noise level δ=10-4 and δ=10-6, we get numerical results which are comparable with the double precision results given in [15]. Figure 7 and Table 7 show the stability of the proposed method for decreasing δ.Figure 7
Example3.7: the stability of the approximate solution.Example 3.8 (see [13, 25]).
We have(3.9)f8(t)={4t2e-2t,0≤t<10,0,elsewhere,F8(p)=8+4e-10(2+p)[-2-20(2+p)-100(2-p)2](2+p)3.
The results of this example are similar to the results of Example 3.3. The exact solution can be well reconstructed by the approximate solution obtained by our method at the levels noise δ=10-4 and δ=10-6 (see Figure 8). Table 8 shows that the MAE decreases as the noise level decreases which shows the stability of the proposed method with respect to the noise. In all the levels of noise δ the computation time of the proposed method in obtaining the approximate solution is relatively small. We get better reconstruction results than the results shown in [13]. Our results are comparable with the results given in [25].Figure 8
Example3.8: the stability of the approximate solution.Example 3.9 (see [16]).
We have(3.10)f9(t)={5-t,0≤t<5,0,elsewhere,F9(p)={252,p=0,e-5p+5p-1p2,p>0.
As in Example 3.6 the error of the approximate solution at the point where the function is not differentiable dominates the error of the approximation. The reconstruction of the exact solution can be seen in Figure 9. The detailed results are presented in Table 9. When the double precision is used, we get comparable results with the results shown in [16].Figure 9
Example3.9: the stability of the approximate solution.Example 3.10 (see [6]).
We have(3.11)f10(t)={t,0≤t<10,0,elsewhere,F10(p)={50,p=0,1-e-10pp2-10e-10pp,p>0.
Table 10 shows the stability of the solution obtained by our method with respect to the noise level δ. We get an excellent agreement between the exact solution and the approximate solution for all the noise levels δ as shown in Figure 10.Figure 10
Example3.10: the stability of the approximate solution.Example 3.11 (see [6, 9]).
We have(3.12)f11(t)={sin(t),0≤t<10,0,elsewhere,F11(p)=1-e-10p(psin(10)+cos(10))1+p2.
Here the function f11(t) represents the class of periodic functions. It  is mentioned in [9] that oscillating function can be found with acceptable accuracy only for relatively small values of t. In this example the best approximation is obtained when the noise level δ=10-6 which is comparable with the results given in [6, 9]. The reconstruction of the function f11(t) for various levels of the noise δ are given in Figure 11. The stability of the proposed method with respect to the noise δ is shown in Table 11. In this example the parameter mnδ is bounded by the constant 54 when the noise level δ=10-6 and κ=0.3.Figure 11
Example3.11: the stability of the approximate solution.Example 3.12 (see [6, 25]).
We have(3.13)f12(t)={tcos(t),0≤t<10,0,elsewhere,F12(p)=(p2-1)-e-10p(-1+p2+10p+10p3)cos(10)(1+p2)2+e-10p(2p+10+10p2)sin(10)(1+p2)2.
Here we take an increasing function which oscillates as the variable t increases over the interval [0,10). A poor approximation is obtained when the noise level δ=10-2. Figure 12 shows that the exact solution can be approximated very well when the noise level δ=10-6. The results of our method are comparable with these of the methods given in [6, 25]. The stability of our method with respect to the noise level is shown in Table 12.Figure 12
Example3.12: the stability of the approximate solution.Example3.13: the stability of the approximate solution.
(a)(b)(c)(d)Example 3.13.
We have(3.14)f13(t)=e-t,F13(p)=11+p.
Here the support of f13(t) is not compact. From the Laplace transform formula one gets
(3.15)F13(p)=∫0∞e-te-ptdt=∫0be-(1+p)tdt+∫b∞e-(1+p)tdt=∫0bf13(t)e-ptdt+e-(1+p)b1+p:=I1+I2,
where δ(b):=e-b. Therefore, I2 can be considered as noise of the data F13(p), that is,
(3.16)F13δ(p):=F13(p)-δ(b),
where δ(b):=e-b. In this example the following parameters are used: d=2, κ=10-1 for δ=e-5 and κ=10-5 for δ=10-8,10-20 and 10-30. Table 13 shows that the error decreases as the parameter b increases. The approximate solution obtained by the proposed method converges to the function f13(t) as b increases (see Figure 13).


## 3.1. The Parametersκ,a0,d

From definition (2.35) and the rule (2.61) we conclude that mn→∞ as an→0. Therefore, one needs to control the value of the parameter mn so that it will not grow too fast as an decreases. The role of the parameter κ in (2.61) is to control the value of the parameter mn so that the value of the parameter mn will not be too large. Since for sufficiently small noise level δ, namely, δ∈(10-16,10-6], the regularization parameter anδ, obtained by the stopping rule (2.83), is at most O(10-9), we suggest to choose κ in the interval (0,1]. For the noise level δ∈(10-6,10-2] one can choose κ∈(1,3]. To reduce the number of iterations we suggest to choose the geometric sequence an=a0δαn, where a0∈[0.1,0.2] and α∈[0.5,0.9]. One may assume without loss of generality that b=1, because a scaling transformation reduces the integral over (0,b) to the integral over (0,1). We have assumed that the data F(p) are defined on the interval J:=[0,d]. In the case the interval J=[d1,d], 0<d1<d, the constant d in estimates (2.59), (2.74), (2.75), (2.78), (2.90), (2.91), and (2.97) are replaced with the constant d-d1. If b=1, that is, f(t)=0 for t>1, then one has to take d not too large. Indeed, if f(t)=0 for t>1, then an integration by parts yields F(p)=[f(0)-e-pf(1)]/p+O(1/p2),p→∞. If the data are noisy, and the noise level is δ, then the data becomes indistinguishable from noise for p=O(1/δ). Therefore it is useless to keep the data Fδ(p) for d>O(1/δ). In practice one may get a satisfactory accuracy of inversion by the method, proposed in this paper, when one uses the data with d∈[1,20] when δ≤10-2. In all the numerical examples we have used d=5. Given the interval [0,d], the proposed method generates automatically the discrete data Fδ(pj), j=0,1,2,…,m, over the interval [0,d] which are needed to get the approximation of the function f(t).


## 3.2. Experiments

To test the proposed method we consider some examples proposed in [5–7, 9, 12, 13, 15, 16, 19, 25]. To illustrate the numerical stability of the proposed method with respect to the noise, we use the noisy data Fδ(p) with various noise levels δ=10-2,δ=10-4, and δ=10-6. The random quantities δj in (2.73) are obtained from the uniform probability density function over the interval [-δ,δ]. In Examples 3.1–3.12 we choose the value of the parameters as follows: an=0.1qn, q=δ1/2, and d=5. The parameter κ=1 is used for the noise levels δ=10-2 and δ=10-4. When δ=10-6 we choose κ=0.3 so that the value of the parameters mn are not very large, namely, mn≤300. Therefore, the computation time for solving linear algebraic system (2.117) can be reduced significantly. We assume that the support of the function f(t) is in the interval [0,b] with b=10. In the stopping rule (2.83) the following parameters are used: C=d+0.01, ε=0.99. In Example 3.13 the function f(t)=e-t is used to test the applicability of the proposed method to functions without compact support. The results are given in Table 13 and Figure 13.For a comparison with the exact solutions we use the mean absolute error(3.1)MAE:=[∑j=1100(f(ti)-fmnδδ(ti))2100]1/2,tj=0.01+0.1(j-1),j=1,…,100,
where f(t) is the exact solution and fmnδδ(t) is the approximate solution. The computation time (CPU time) for obtaining the approximation of f(t), the number of iterations (Iter.), and the parameters mnδ and anδ generated by the proposed method is given in each experiment (see Tables 1–12). All the calculations are done in double precision generated by MATLAB.Table 1
Example3.1.
δMAEmnδIter.CPU time (second)anδ1.00×10-29.62×10-23033.13×10-22.00×10-31.00×10-45.99×10-23246.25×10-22.00×10-71.00×10-64.74×10-25453.28×10-12.00×10-10Table 2
Example3.2.
δMAEmnδIter.CPU time (seconds)anδ1.00×10-21.09×10-13023.13×10-22.00×10-31.00×10-48.47×10-23236.25×10-22.00×10-61.00×10-67.41×10-25454.38×10-12.00×10-12Table 3
Example3.3.
δMAEmnδIter.CPU time (seconds)anδ1.00×10-22.42×10-23023.13×10-22.00×10-31.00×10-41.08×10-33033.13×10-22.00×10-61.00×10-64.02×10-43044.69×10-22.00×10-9Table 4
Example3.4.
δMAEmnδIter.CPU time (seconds)anδ1.00×10-21.59×10-23023.13×10-22.00×10-31.00×10-48.26×10-43039.400×10-22.00×10-61.00×10-61.24×10-43041.250×10-12.00×10-9Table 5
Example3.5.
δMAEmnδIter.CPU time (seconds)anδ1.00×10-24.26×10-23036.300×10-22.00×10-31.00×10-41.25×10-23039.38×10-22.00×10-61.00×10-61.86×10-35443.13×10-22.00×10-9Table 6
Example3.6.
δMAEmnδIter.CPU time (seconds)anδ1.00×10-24.19×10-23024.700×10-22.00×10-31.00×10-41.64×10-23239.38×10-22.00×10-61.00×10-61.22×10-25443.13×10-22.00×10-9Table 7
Example3.7.
δMAEmnδIter.CPU time (seconds)anδ1.00×10-21.52×10-23024.600×10-22.00×10-31.00×10-42.60×10-33039.38×10-22.00×10-61.00×10-62.02×10-33043.13×10-22.00×10-9Table 8
Example3.8.
δMAEmnδIter.CPU time (seconds)anδ1.00×10-22.74×10-23021.100×10-22.00×10-31.00×10-43.58×10-33033.13×10-22.00×10-61.00×10-65.04×10-43044.69×10-22.00×10-9Table 9
Example3.9.
δMAEmnδIter.CPU time (seconds)anδ1.00×10-22.07×10-13036.25×10-22.00×10-61.00×10-47.14×10-23243.44×10-12.00×10-91.00×10-62.56×10-25453.75×10-12.00×10-12Table 10
Example3.10.
δMAEmnδIter.CPU time (seconds)anδ1.00×10-22.09×10-13033.13×10-22.00×10-61.00×10-41.35×10-23249.38×10-22.00×10-91.00×10-63.00×10-35442.66×10-12.00×10-9Table 11
Example3.11.
δMAEmnδIter.CPU time (seconds)anδ1.00×10-22.47×10-13039.38×10-22.00×10-61.00×10-44.91×10-23242.50×10-12.00×10-91.00×10-62.46×10-25454.38×10-12.00×10-12Table 12
Example3.12.
δMAEmnδIter.CPU time (seconds)anδ1.00×10-21.37×1009639.38×10-22.00×10-61.00×10-45.98×10-110042.66×10-12.00×10-91.00×10-62.24×10-130053.44×10-12.00×10-12Table 13
Example3.13.
bMAEmδIter.CPU time (seconds)51.487×10-2243.125×10-282.183×10-4243.125×10-2204.517×10-9243.125×10-2301.205×10-13243.125×10-2Example 3.1 (see [15]).
We have(3.2)f1(t)={1,12≤t≤32,0,otherwise,F1(p)={1,p=0,e-p/2-e-3p/2p,p>0.
The reconstruction of the exact solution for different values of the noise level δ is shown in Figure 1. When the noise level δ=10-6, our result is comparable with the double precision results shown in [15]. The proposed method is stable with respect to the noise δ as shown in Table 1.Figure 1
Example3.1: the stability of the approximate solution.Example 3.2 (see [13, 15]).
We have(3.3)f2(t)={12,t=1,1,1<t<10,0,elsewhere,F2(p)={9,p=0,e-p-e-10pp,p>0.
The reconstruction of the function f2(t) is plotted in Figure 2. In [15] a high accuracy result is given by means of the multiple precision. But, as reported in [15], to get such high accuracy results, it takes 7 hours. From Table 2 and Figure 2 we can see that the proposed method yields stable solution with respect to the noise level δ. The reconstruction of the exact solution obtained by the proposed method is better than the reconstruction shown in [13]. The result is comparable with the double precision results given in [15]. For δ=10-6 and κ=0.3 the value of the parameter mnδ is bounded by the constant 54.Figure 2
Example3.2: the stability of the approximate solution.Example 3.3 (see [6, 12, 13, 16, 19]).
We have(3.4)f3(t)={te-t,0≤t<10,0,otherwise,F3(p)=1-e-(p+1)10(p+1)2-10e-(p+1)10p+1.
We get an excellent agreement between the approximate solution and the exact solution when the noise level δ=10-4 and 10-6 as shown in Figure 3. The results obtained by the proposed method are better than the results given in [13]. The mean absolute error MAE decreases as the noise level decreases which shows the stability of the proposed method. Our results are more stable with respect to the noise δ than the results presented in [19]. The value of the parameter mnδ is bounded by the constant 30 when the noise level δ=10-6 and κ=0.3.Figure 3
Example3.3: the stability of the approximate solution.Figure 4
Example3.4: the stability of the approximate solution.Example 3.4 (see [13, 15]).
We have(3.5)f4(t)={1-e-0.5t,0≤t<10,0,elsewhere,F4(p)={8+2e-5,p=0,1-e-10pp-1-e-(p+1/2)10p+0.5,p>0.
As in our Example 3.3 when the noise δ=10-4 and 10-6 are used, we get a satisfactory agreement between the approximate solution and the exact solution. Table 4 gives the results of the stability of the proposed method with respect to the noise level δ. Moreover, the reconstruction of the function f4(t) obtained by the proposed method is better than the reconstruction of f4(t) shown in [13], and is comparable with the double precision reconstruction obtained in [15].
In this example whenδ=10-6 and κ=0.3 the value of the parameter mnδ is bounded by the constant 109 as shown in Table 4.Example 3.5 (see [5, 7, 13]).
We have(3.6)f5(t)=23e-t/2sin(t3/2)F5(p)=1-cos(103/2)e-10(p+0.5)[(p+0.5)2+3/4]-2(p+0.5)e-10(p+0.5)sin(103/2)3[(p+0.5)2+3/4].
This is an example of the damped sine function. In [5, 7] the knowledge of the exact data F(p) in the complex plane is required to get the approximate solution. Here we only use the knowledge of the discrete perturbed data Fδ(pj), j=0,1,2,…,m, and get a satisfactory result which is comparable with the results given in [5, 7] when the level noise δ=10-6. The reconstruction of the exact solution f5(t) obtained by our method is better than this of the method given in [13]. Moreover, our method yields stable solution with respect to the noise level δ as shown in Figure 5 and Table 5. In this example when κ=0.3 the value of the parameter mnδ is bounded by 54 for the noise level δ=10-6 (see Table 5).Figure 5
Example3.5: the stability of the approximate solution.Example 3.6 (see [15]).
We have(3.7)f6(t)={t,0≤t<1,32-t2,1≤t<3,0,elsewhere,F6(p)={32,p=0,1-e-p(1+p)p2+e-3p+e2p(2p-1)2p2,p>0.
Example 3.6 represents a class of piecewise continuous functions. From Figure 6 the value of the exact solution at the points where the function is not differentiable cannot be well approximated for the given levels of noise by the proposed method. When the noise level δ=10-6, our result is comparable with the results given in [15]. Table 6 reports the stability of the proposed method with respect to the noise δ. It is shown in Table 6 that the value of the parameter m generated by the proposed adaptive stopping rule is bounded by the constant 54 for the noise level δ=10-6 and κ=0.3 which gives a relatively small computation time.Figure 6
Example3.6: the stability of the approximate solution.Example 3.7 (see [15]).
We have(3.8)f7(t)={-te-t-e-t+1,0≤t<1,1-2e-1,1≤t<10,0,elsewhere,F7(p)={3e-1+9(1-2e),p=0,e-1-p(e1+p-e(1+p)2+p(3+2p)/p(p+1)2)+(e-2)e-1-p-10p((e10p-ep)/p),p>0.
When the noise level δ=10-4 and δ=10-6, we get numerical results which are comparable with the double precision results given in [15]. Figure 7 and Table 7 show the stability of the proposed method for decreasing δ.Figure 7
Example3.7: the stability of the approximate solution.Example 3.8 (see [13, 25]).
We have(3.9)f8(t)={4t2e-2t,0≤t<10,0,elsewhere,F8(p)=8+4e-10(2+p)[-2-20(2+p)-100(2-p)2](2+p)3.
The results of this example are similar to the results of Example 3.3. The exact solution can be well reconstructed by the approximate solution obtained by our method at the levels noise δ=10-4 and δ=10-6 (see Figure 8). Table 8 shows that the MAE decreases as the noise level decreases which shows the stability of the proposed method with respect to the noise. In all the levels of noise δ the computation time of the proposed method in obtaining the approximate solution is relatively small. We get better reconstruction results than the results shown in [13]. Our results are comparable with the results given in [25].Figure 8
Example3.8: the stability of the approximate solution.Example 3.9 (see [16]).
We have(3.10)f9(t)={5-t,0≤t<5,0,elsewhere,F9(p)={252,p=0,e-5p+5p-1p2,p>0.
As in Example 3.6 the error of the approximate solution at the point where the function is not differentiable dominates the error of the approximation. The reconstruction of the exact solution can be seen in Figure 9. The detailed results are presented in Table 9. When the double precision is used, we get comparable results with the results shown in [16].Figure 9
Example3.9: the stability of the approximate solution.Example 3.10 (see [6]).
We have(3.11)f10(t)={t,0≤t<10,0,elsewhere,F10(p)={50,p=0,1-e-10pp2-10e-10pp,p>0.
Table 10 shows the stability of the solution obtained by our method with respect to the noise level δ. We get an excellent agreement between the exact solution and the approximate solution for all the noise levels δ as shown in Figure 10.Figure 10
Example3.10: the stability of the approximate solution.Example 3.11 (see [6, 9]).
We have(3.12)f11(t)={sin(t),0≤t<10,0,elsewhere,F11(p)=1-e-10p(psin(10)+cos(10))1+p2.
Here the function f11(t) represents the class of periodic functions. It  is mentioned in [9] that oscillating function can be found with acceptable accuracy only for relatively small values of t. In this example the best approximation is obtained when the noise level δ=10-6 which is comparable with the results given in [6, 9]. The reconstruction of the function f11(t) for various levels of the noise δ are given in Figure 11. The stability of the proposed method with respect to the noise δ is shown in Table 11. In this example the parameter mnδ is bounded by the constant 54 when the noise level δ=10-6 and κ=0.3.Figure 11
Example3.11: the stability of the approximate solution.Example 3.12 (see [6, 25]).
We have(3.13)f12(t)={tcos(t),0≤t<10,0,elsewhere,F12(p)=(p2-1)-e-10p(-1+p2+10p+10p3)cos(10)(1+p2)2+e-10p(2p+10+10p2)sin(10)(1+p2)2.
Here we take an increasing function which oscillates as the variable t increases over the interval [0,10). A poor approximation is obtained when the noise level δ=10-2. Figure 12 shows that the exact solution can be approximated very well when the noise level δ=10-6. The results of our method are comparable with these of the methods given in [6, 25]. The stability of our method with respect to the noise level is shown in Table 12.Figure 12
Example3.12: the stability of the approximate solution.Example3.13: the stability of the approximate solution.
(a)(b)(c)(d)Example 3.13.
We have(3.14)f13(t)=e-t,F13(p)=11+p.
Here the support of f13(t) is not compact. From the Laplace transform formula one gets
(3.15)F13(p)=∫0∞e-te-ptdt=∫0be-(1+p)tdt+∫b∞e-(1+p)tdt=∫0bf13(t)e-ptdt+e-(1+p)b1+p:=I1+I2,
where δ(b):=e-b. Therefore, I2 can be considered as noise of the data F13(p), that is,
(3.16)F13δ(p):=F13(p)-δ(b),
where δ(b):=e-b. In this example the following parameters are used: d=2, κ=10-1 for δ=e-5 and κ=10-5 for δ=10-8,10-20 and 10-30. Table 13 shows that the error decreases as the parameter b increases. The approximate solution obtained by the proposed method converges to the function f13(t) as b increases (see Figure 13).


## 4. Conclusion

We have tested the proposed algorithm on the wide class of examples considered in the literature. Using the rule (2.61) and the stopping rule (2.83), the number of terms in representation (2.69), the discrete data Fδ(pj), j=0,1,2,…,m, and regularization parameter anδ, which are used in computing the approximation fmδ(t) (see (2.69)) of the unknown function f(t), are obtained automatically. Our numerical experiments show that the computation time (CPU time) for approximating the function f(t) is small, namely, CPU time ≤1 second, and the proposed iterative scheme and the proposed adaptive stopping rule yield stable solution given the noisy data with the noise level δ. The proposed method also works for f with non-compact support as shown in Example 3.13. Moreover, in the proposed method we only use a simple representation (2.69) which is based on the kernel of the Laplace transform integral, so it can be easily implemented numerically.



---

*Source: 898195-2010-01-05.xml*