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

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 function 𝑓(𝑡) is continuous with (known) compact support. An adaptive iterative method and an adaptive stopping rule, which yield the convergence of the approximate solution to 𝑓(𝑡), are proposed in this paper.

1. Introduction

Consider the Laplace transform

ℒ𝑓(𝑝)∶=^∞₀𝑒^−𝑝𝑡𝑓(𝑡)𝑑𝑡=𝐹(𝑝),Re𝑝>0, (1.1) where ℒ∶𝑋_0,𝑏→𝐿²[0,∞),

𝑋_0,𝑏∶={𝑓∈𝐿²[0,∞)∣supp𝑓⊂[0,𝑏)},𝑏>0. (1.2) We assume in (1.2) that 𝑓 has compact support. This is not a restriction practically. Indeed, if lim_𝑡→∞𝑓(𝑡)=0, then |𝑓(𝑡)|<𝛿 for 𝑡>𝑡_𝛿, where 𝛿>0 is an arbitrary small number. Therefore, one may assume that supp𝑓⊂[0,𝑡_𝛿], and treat the values of 𝑓 for 𝑡>𝑡_𝛿 as noise. One may also note that if 𝑓∈𝐿¹(0,∞), then 𝐹(𝑝)∶=^∞₀𝑓(𝑡)𝑒^−𝑝𝑡𝑑𝑡=^𝑏₀𝑓(𝑡)𝑒^−𝑝𝑡𝑑𝑡+^∞_𝑏𝑓(𝑡)𝑒^−𝑝𝑡𝑑𝑡∶=𝐹₁(𝑝)+𝐹₂(𝑝), (1.3) and |𝐹₂(𝑝)|≤𝑒^−𝑏𝑝𝛿, where ∫^∞_𝑏|𝑓(𝑡)|𝑑𝑡≤𝛿. Therefore, the contribution of the “tail” 𝑓_𝑏(𝑡) of 𝑓, 𝑓_𝑏(𝑡)∶=0,𝑡<𝑏,𝑓(𝑡),𝑡≥𝑏, (1.4) can be considered as noise if 𝑏>0 is large and 𝛿>0 is small. We assume in (1.2) that 𝑓∈𝐿²[0,∞). One may also assume that 𝑓∈𝐿¹[0,∞), or that |𝑓(𝑡)|≤𝑐₁𝑒^𝑐2𝑡, where 𝑐₁,𝑐₂ are positive constants. If the last assumption holds, then one may define the function 𝑔(𝑡)∶=𝑓(𝑡)𝑒^{−(𝑐2+1)𝑡}. Then 𝑔(𝑡)∈𝐿¹[0,∞), and its Laplace transform 𝐺(𝑝)=𝐹(𝑝+𝑐₂+1) is known on the interval [𝑐₂+1,𝑐₂+1+𝑏] of real axis if the Laplace transform 𝐹(𝑝) of 𝑓(𝑡) is known on the interval [0,𝑏]. Therefore, our inversion methods are applicable to these more general classes of functions 𝑓 as well.

The operator ℒ∶𝑋_0,𝑏→𝐿²[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 𝐹(𝑝) 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𝑓(𝑡)=2𝜋𝑖^𝜎+𝑖∞_{𝜎−𝑖∞}𝐹(𝑝)𝑒^𝑝𝑡𝑑𝑝,𝜎>0, (1.5) is used in some of these methods, and then 𝑓(𝑡) 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 𝑓(𝑡), obtained by these methods when the data are noisy, may differ significantly from 𝑓(𝑡). 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 𝐹(𝑝) is calculated first and then inverted for 𝑓(𝑡). 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 𝑋_0,𝑏, the problem becomes extremely ill-conditioned if the number 𝑚 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 from 𝐶_0𝑘 into 𝐿², where 𝐶_0𝑘={𝑓(𝑡)∈𝐶[0,+∞)∣supp𝑓⊂[0,𝑘)}, 𝑘=const>0, 𝐿²∶=𝐿²[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 𝐹(𝑝) is given only on a finite interval [0,𝑑], 𝑑>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 𝑓(𝑡) 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 𝐹_𝛿(𝑝), and the number of terms in (2.69), needed for obtaining an approximation of the unknown function 𝑓(𝑡).

We study the inversion problem using the pair of spaces (𝑋_0,𝑏,𝐿²[0,𝑑]), where 𝑋_0,𝑏 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 𝑓(𝑡) when the noise level 𝛿→0.

A detailed derivation of our inversion method is given in Section 2. 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

Let 𝑓∈𝑋_0,𝑏. Then (1.1) can be written as

(ℒ𝑓)(𝑝)∶=^𝑏₀𝑒^−𝑝𝑡𝑓(𝑡)𝑑𝑡=𝐹(𝑝),0≤𝑝. (2.1) Let us assume that the data 𝐹(𝑝), the Laplace transform of 𝑓, are known only for 0≤𝑝≤𝑑<∞. Consider the mapping ℒ_𝑚∶𝐿²[0,𝑏]→ℝ^𝑚+1, where

(ℒ_𝑚𝑓)_𝑖∶=^𝑏₀𝑒^{−𝑝𝑖𝑡}𝑓(𝑡)𝑑𝑡=𝐹(𝑝_𝑖𝑝),𝑖=0,1,2,…,𝑚,_𝑖𝑑∶=𝑖ℎ,𝑖=0,1,2,…,𝑚,ℎ∶=𝑚, (2.2)(2.3) and 𝑚 is an even number which will be chosen later. Then the unknown function 𝑓(𝑡) can be obtained from a finite-dimensional operator equation (2.2). Let

⟨𝑢,𝑣⟩_𝑊^𝑚∶=^𝑚_𝑗=0𝑤_𝑗^(𝑚)𝑢_𝑗𝑣_𝑗,‖𝑢‖_𝑊^𝑚∶=⟨𝑢,𝑢⟩_𝑊^𝑚 (2.4) be the inner product and norm in ℝ^𝑚+1, respectively, where 𝑤_𝑗^(𝑚) are the weights of the compound Simpson's rule (see [10, page 58]), that is,

𝑤_𝑗^(𝑚)ℎ∶=3,𝑗=0,𝑚,4ℎ3𝑚,𝑗=2𝑙−1,𝑙=1,2,…,2,2ℎ3,𝑗=2𝑙,𝑙=1,2,…,𝑚−22,𝑑ℎ=𝑚, (2.5) where 𝑚 is an even number. Then

⟨ℒ_𝑚𝑔,𝑣⟩_𝑊𝑚=^𝑚_𝑗=0𝑤_𝑗^(𝑚)𝑏₀𝑒^{−𝑝𝑗𝑡}𝑔(𝑡)𝑑𝑡𝑣_𝑗=^𝑏₀𝑔(𝑡)^𝑚_𝑗=0𝑤_𝑗^(𝑚)𝑒^{−𝑝𝑗𝑡}𝑣_𝑗𝑑𝑡=⟨𝑔,ℒ^∗_𝑚𝑣⟩_𝑋0,𝑏, (2.6) where

ℒ^∗_𝑚𝑣=^𝑚_𝑗=0𝑤_𝑗^(𝑚)𝑒^{−𝑝𝑗𝑡}𝑣_𝑗𝑣,𝑣∶=⁰𝑣₁⋮𝑣_𝑚∈ℝ^𝑚+1,⟨𝑔,ℎ⟩_𝑋0,𝑏∶=^𝑏₀𝑔(𝑡)ℎ(𝑡)𝑑𝑡. (2.7)(2.8) It follows from (2.2) and (2.7) that

(ℒ^∗_𝑚ℒ_𝑚𝑔)(𝑡)=^𝑚_𝑗=0𝑤_𝑗^(𝑚)𝑒^{−𝑝𝑗𝑡𝑏}₀𝑒^{−𝑝𝑗𝑧}𝑔(𝑧)𝑑𝑧∶=(𝑇^(𝑚)ℒ𝑔)(𝑡),_𝑚ℒ^∗_𝑚𝑣=^𝑏0𝑒^{−𝑝0𝑡𝑚}_𝑗=0𝑤_𝑗^(𝑚)𝑒^{−𝑝𝑗𝑡}𝑣_𝑗𝑑𝑡_𝑏0𝑒^{−𝑝1𝑡𝑚}_𝑗=0𝑤_𝑗^(𝑚)𝑒^{−𝑝𝑗𝑡}𝑣_𝑗⋮𝑑𝑡_𝑏0𝑒^{−𝑝𝑚𝑡𝑚}_𝑗=0𝑤_𝑗^(𝑚)𝑒^{−𝑝𝑗𝑡}𝑣_𝑗𝑑𝑡∶=𝑄^(𝑚)𝑣, (2.9)(2.10) where

(𝑄^(𝑚))_𝑖𝑗∶=𝑤_𝑗^(𝑚)𝑏₀𝑒^{−(𝑝𝑖+𝑝𝑗)𝑡}𝑑𝑡=𝑤_𝑗^(𝑚)1−𝑒^{−𝑏(𝑝𝑖+𝑝𝑗)}𝑝_𝑖+𝑝_𝑗,𝑖,𝑗=0,1,2,…,𝑚. (2.11)

From (2.7) we get Range[ℒ^∗_𝑚]=𝑠𝑝𝑎𝑛{𝑤_𝑗^(𝑚)𝑘(𝑝_𝑗,⋅,0)}^𝑚_𝑗=0, where

Lemma 2.1. Let 𝑤_𝑗^(𝑚) be defined in (2.5). Then ^𝑚_𝑗=0𝑤_𝑗^(𝑚)=𝑑, (2.12) for any even number 𝑚.

Proof. From definition (2.5) one gets ^𝑚_𝑗=0𝑤_𝑗^(𝑚)=𝑤₀^(𝑚)+𝑤_𝑚^(𝑚)+^𝑚/2_𝑗=1𝑤^(𝑚)_2𝑗−1+^(𝑚−2)/2_𝑗=1𝑤^(𝑚)_2𝑗=2ℎ3+^𝑚/2_𝑗=14ℎ3+^(𝑚−2)/2_𝑗=12ℎ3=2ℎ3+2ℎ𝑚3+ℎ(𝑚−2)3𝑑=ℎ𝑚=𝑚𝑚=𝑑. (2.13) Lemma 2.1 is proved.

Lemma 2.2. The matrix 𝑄^(𝑚), defined in (2.11), is positive semidefinite and self-adjoint in ℝ^𝑚+1 with respect to the inner product (2.4).

Proof. Let (𝐻_𝑚)_𝑖𝑗∶=^𝑏₀𝑒^{−(𝑝𝑖+𝑝𝑗)𝑡}𝑑𝑡=1−𝑒^{−𝑏(𝑝𝑖+𝑝𝑗)}𝑝_𝑖+𝑝_𝑗,(𝐷_𝑚)_𝑖𝑗=𝑤_𝑖^(𝑚),𝑖=𝑗,0,otherwise, (2.14)𝑤_𝑗^(𝑚) are defined in (2.5). Then ⟨𝐷_𝑚𝐻_𝑚𝐷_𝑚𝑢,𝑣⟩_ℝ^𝑚+1=⟨𝑢,𝐷_𝑚𝐻_𝑚𝐷_𝑚𝑣⟩_ℝ^𝑚+1, where ⟨𝑢,𝑣⟩_ℝ^𝑚+1∶=^𝑚_𝑗=0𝑢_𝑗𝑣_𝑗,𝑢,𝑣∈ℝ^𝑚+1. (2.15) We have ⟨𝑄^(𝑚)𝑢,𝑣⟩_𝑊^𝑚=^𝑚_𝑗=0𝑤_𝑗^(𝑚)(𝑄^(𝑚)𝑢)_𝑗𝑣_𝑗=^𝑚_𝑗=0(𝐷_𝑚𝐻_𝑚𝐷_𝑚𝑢)_𝑗𝑣_𝑗=⟨𝐷_𝑚𝐻_𝑚𝐷_𝑚𝑢,𝑣⟩_ℝ𝑚+1=⟨𝑢,𝐷_𝑚𝐻_𝑚𝐷_𝑚𝑣⟩_ℝ𝑚+1=^𝑚_𝑗=0𝑢_𝑗(𝐷_𝑚𝐻_𝑚𝐷_𝑚𝑣)_𝑗=^𝑚_𝑗=0𝑢_𝑗𝑤_𝑗^(𝑚)(𝐻_𝑚𝐷_𝑚𝑣)_𝑗=⟨𝑢,𝑄^(𝑚)𝑣⟩_𝑊^𝑚. (2.16) Thus, 𝑄^(𝑚) is self-adjoint with respect to inner product (2.4). We have (𝐻_𝑚)_𝑖𝑗=^𝑏₀𝑒^{−(𝑝𝑖+𝑝𝑗)𝑡}𝑑𝑡=^𝑏₀𝑒^{−𝑝𝑖𝑡}𝑒^{−𝑝𝑗𝑡}𝑑𝑡=⟨𝜙_𝑖,𝜙_𝑗⟩_𝑋0,𝑏,𝜙_𝑖(𝑡)∶=𝑒^{−𝑝𝑖𝑡}, (2.17) where ⟨⋅,⋅⟩_𝑋0,𝑏 is defined in (2.8). This shows that 𝐻_𝑚 is a Gram matrix. Therefore, ⟨𝐻_𝑚𝑢,𝑢⟩_ℝ𝑚+1≥0,∀𝑢∈ℝ^𝑚+1. (2.18) This implies ⟨𝑄^(𝑚)𝑢,𝑢⟩_𝑊^𝑚=⟨𝑄^(𝑚)𝑢,𝐷_𝑚𝑢⟩_ℝ𝑚+1=⟨𝐻_𝑚𝐷_𝑚𝑢,𝐷_𝑚𝑢⟩_ℝ𝑚+1≥0. (2.19) Thus, 𝑄^(𝑚) is a positive semidefinite and self-adjoint matrix with respect to the inner product (2.4).

Lemma 2.3. Let 𝑇^(𝑚) be defined in (2.9). Then 𝑇^(𝑚) is self-adjoint and positive semidefinite operator in 𝑋_0,𝑏 with respect to inner product (2.8).

Proof. From definition (2.9) and inner product (2.8) we get ⟨𝑇^(𝑚)𝑔,ℎ⟩_𝑋0,𝑏=^𝑏₀^𝑚_𝑗=0𝑤_𝑗^(𝑚)𝑒^{−𝑝𝑗𝑡𝑏}₀𝑒^{−𝑝𝑗𝑧}=𝑔(𝑧)𝑑𝑧ℎ(𝑡)𝑑𝑡_𝑏0𝑔(𝑧)^𝑚_𝑗=0𝑤_𝑗^(𝑚)𝑒^{−𝑝𝑗𝑧𝑏}₀𝑒^{−𝑝𝑗𝑡}ℎ(𝑡)𝑑𝑡𝑑𝑧=⟨𝑔,𝑇^(𝑚)ℎ⟩_𝑋0,𝑏. (2.20) Thus, 𝑇^(𝑚) is a self-adjoint operator with respect to inner product (2.8). Let us prove that 𝑇^(𝑚) is positive semidefinite. Using (2.9), (2.5), (2.4), and (2.8), one gets ⟨𝑇^(𝑚)𝑔,𝑔⟩_𝑋0,𝑏=^𝑏₀^𝑚_𝑗=0𝑤_𝑗^(𝑚)𝑒^{−𝑝𝑗𝑡𝑏}₀𝑒^{−𝑝𝑗𝑧}=𝑔(𝑧)𝑑𝑧𝑔(𝑡)𝑑𝑡_𝑚𝑗=0𝑤_𝑗^(𝑚)𝑏₀𝑒^{−𝑝𝑗𝑧}𝑔(𝑧)𝑑𝑧^𝑏₀𝑒^{−𝑝𝑗𝑡}=𝑔(𝑡)𝑑𝑡_𝑚𝑗=0𝑤_𝑗^(𝑚)(^𝑏₀𝑒^{−𝑝𝑗𝑧}𝑔(𝑧)𝑑𝑧)²≥0. (2.21) Lemma 2.3 is proved.

𝑘(𝑝,𝑡,𝑧)∶=𝑒^{−𝑝(𝑡+𝑧)}. (2.22) Let us approximate the unknown 𝑓(𝑡) as follows:

𝑓(𝑡)≈^𝑚_𝑗=0𝑐_𝑗^(𝑚)𝑤_𝑗^(𝑚)𝑒^{−𝑝𝑗𝑡}=𝑇⁻¹_𝑎,𝑚ℒ^∗_𝑚𝐹^(𝑚)∶=𝑓_𝑚(𝑡), (2.23) where 𝑝_𝑗 are defined in (2.3), 𝑇_𝑎,𝑚 is defined in (2.30), and 𝑐_𝑗^(𝑚) are constants obtained by solving the linear algebraic system:

(𝑎𝐼+𝑄^(𝑚))𝑐^(𝑚)=𝐹^(𝑚), (2.24) where 𝑄^(𝑚) is defined in (2.10),

𝑐^(𝑚)𝑐∶=^0(𝑚)𝑐₁^(𝑚)⋮𝑐_𝑚^(𝑚),𝐹^(𝑚)∶=𝐹(𝑝₀)𝐹(𝑝₁)⋮𝐹(𝑝_𝑚). (2.25) To prove the convergence of the approximate solution 𝑓(𝑡), 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. Let 𝑇^(𝑚) and 𝑄^(𝑚) be defined in (2.9) and (2.10), respectively. Then, for 𝑎>0, the following estimates hold: ‖𝑄⁻¹_𝑎,𝑚ℒ_𝑚1‖≤2√𝑎,𝑎‖𝑄⁻¹_𝑎,𝑚‖≤1,‖𝑇⁻¹_𝑎,𝑚1‖≤𝑎,‖𝑇⁻¹_𝑎,𝑚ℒ^∗_𝑚1‖≤2√𝑎, (2.26)(2.27)(2.28)(2.29) where 𝑄_𝑎,𝑚∶=𝑄^(𝑚)+𝑎𝐼,𝑇_𝑎,𝑚∶=𝑇^(𝑚)+𝑎𝐼, (2.30)𝐼 is the identity operator and 𝑎=𝑐𝑜𝑛𝑠𝑡>0.

Let us formulate an iterative method for obtaining the approximation solution of 𝑓(𝑡) with the exact data 𝐹(𝑝). Consider the following iterative scheme:

𝑢_𝑛(𝑡)=𝑞𝑢_𝑛−1(𝑡)+(1−𝑞)𝑇_𝑎⁻¹_𝑛ℒ^∗𝐹,𝑢₀(𝑡)=0, (2.31) where ℒ^∗ is the adjoint of the operator ℒ, that is,

(ℒ^∗𝑔)(𝑡)=^𝑑₀𝑒^−𝑝𝑡𝑔(𝑝)𝑑𝑝,(𝑇𝑓)(𝑡)∶=(ℒ^∗ℒ𝑓)(𝑡)=^𝑏₀^𝑑₀𝑘(𝑝,𝑡,𝑧)𝑑𝑝𝑓(𝑧)𝑑𝑧=^𝑏₀𝑓(𝑧)𝑡+𝑧(1−𝑒^{−𝑑(𝑡+𝑧)})𝑑𝑧, (2.32)(2.33)𝑘(𝑝,𝑡,𝑧) is defined in (2.22),

𝑇_𝑎𝑎∶=𝑎𝐼+𝑇,𝑎>0,_𝑛∶=𝑞𝑎_𝑛−1,𝑎₀>0,𝑞∈(0,1). (2.34)(2.35)

Lemma 2.5 together with Lemma 2.7 with 𝑔(𝑎)=𝑎‖𝑇_𝑎⁻¹𝑓‖ yields

Lemma 2.5. Let 𝑇_𝑎 be defined in (2.34), ℒ𝑓=𝐹, and 𝑓⟂𝒩(ℒ), where 𝒩(ℒ) is the null space of ℒ. Then 𝑎‖𝑇_𝑎⁻¹𝑓‖→0as𝑎→0. (2.36)

Proof. Since 𝑓⟂𝒩(ℒ), it follows from the spectral theorem that lim_𝑎→0𝑎²‖𝑇_𝑎⁻¹𝑓‖²=lim_𝑎→0^∞₀𝑎²(𝑎+𝑠)²𝑑⟨𝐸_𝑠𝑓,𝑓⟩=‖𝑃_𝒩(ℒ)𝑓‖²=0, (2.37) where 𝐸_𝑠 is the resolution of the identity corresponding to ℒ^∗ℒ, and 𝑃 is the orthogonal projector onto 𝒩(ℒ).
Lemma 2.5 is proved.

Theorem 2.6. Let ℒ𝑓=𝐹, and 𝑢_𝑛 be defined in (2.31). Then lim_𝑛→∞‖𝑓−𝑢_𝑛‖=0. (2.38)

Proof. By induction we get 𝑢_𝑛=^𝑛−1_𝑗=0𝜔_𝑗^(𝑛)𝑇_𝑎⁻¹_𝑗+1ℒ^∗𝐹, (2.39) where 𝑇_𝑎 is defined in (2.34), and 𝜔_𝑗^(𝑛)∶=𝑞^{𝑛−𝑗−1}−𝑞^𝑛−𝑗. (2.40) Using the identities 𝑇ℒ𝑓=𝐹,_𝑎−1ℒ^∗ℒ=𝑇_𝑎⁻¹(𝑇+𝑎𝐼−𝑎𝐼)=𝐼−𝑎𝑇_𝑎⁻¹,_{𝑛−1𝑗=0}𝜔_𝑗^(𝑛)=1−𝑞^𝑛, (2.41) we get 𝑓−𝑢_𝑛=𝑓−^𝑛−1_𝑗=0𝜔_𝑗^(𝑛)𝑓+^𝑛−1_𝑗=0𝜔_𝑗^(𝑛)𝑎_𝑗+1𝑇_𝑎⁻¹_𝑗+1𝑓=𝑞^𝑛𝑓+^𝑛−1_𝑗=0𝜔_𝑗^(𝑛)𝑎_𝑗+1𝑇_𝑎⁻¹_𝑗+1𝑓. (2.42) Therefore, ‖𝑓−𝑢_𝑛‖≤𝑞^𝑛‖𝑓‖+^𝑛−1_𝑗=0𝜔_𝑗^(𝑛)𝑎_𝑗+1‖𝑇_𝑎⁻¹_𝑗+1𝑓‖. (2.43) To prove relation (2.38) the following lemma is needed.

Lemma 2.7. Let 𝑔(𝑥) be a continuous function on (0,∞), 𝑐>0 and 𝑞∈(0,1) constants. If lim_𝑥→0+𝑔(𝑥)=𝑔(0)∶=𝑔₀, (2.44) then lim_𝑛→∞^𝑛−1_𝑗=0(𝑞^{𝑛−𝑗−1}−𝑞^𝑛−𝑗)𝑔(𝑐𝑞^𝑗+1)=𝑔₀. (2.45)

Proof. Let 𝐹_𝑙(𝑛)∶=^𝑙−1_𝑗=1𝜔_𝑗^(𝑛)𝑔(𝑐𝑞^𝑗), (2.46) where 𝜔_𝑗^(𝑛)∶=𝑞^𝑛−𝑗−𝑞^{𝑛+1−𝑗}. Then |𝐹_𝑛+1(𝑛)−𝑔₀|≤|𝐹_𝑙(𝑛)|+|^𝑛_𝑗=𝑙𝜔_𝑗^(𝑛)𝑔(𝑐𝑞^𝑗)−𝑔₀|. (2.47) Take 𝜖>0 arbitrarily small. For sufficiently large fixed 𝑙(𝜖) one can choose 𝑛(𝜖)>𝑙(𝜖), such that |𝐹_𝑙(𝜖)𝜖(𝑛)|≤2,∀𝑛>𝑛(𝜖), (2.48) because lim_𝑛→∞𝑞^𝑛=0. Fix 𝑙=𝑙(𝜖) such that |𝑔(𝑐𝑞^𝑗)−𝑔₀|≤𝜖/2 for 𝑗>𝑙(𝜖). This is possible because of (2.44). One has |𝐹_𝑙(𝜖)𝜖(𝑛)|≤2|,𝑛>𝑛(𝜖)>𝑙(𝜖),_{𝑛𝑗=𝑙(𝜖)}𝜔_𝑗^(𝑛)𝑔(𝑐𝑞^𝑗)−𝑔₀|≤^𝑛_{𝑗=𝑙(𝜖)}𝜔_𝑗^(𝑛)|𝑔(𝑐𝑞^𝑗+1)−𝑔₀|+|^𝑛_{𝑗=𝑙(𝜖)}𝜔_𝑗^(𝑛)−1||𝑔₀|≤𝜖2^𝑛_{𝑗=𝑙(𝜖)}𝜔_𝑗^(𝑛)+𝑞^{𝑛−𝑙(𝜖)}|𝑔₀|≤𝜖2+|𝑔₀|𝑞^{𝑛−𝑙(𝜖)}≤𝜖, (2.49) if 𝑛(𝜖) is sufficiently large. Here we have used the relation ^𝑛_𝑗=𝑙𝜔_𝑗^(𝑛)=1−𝑞^{𝑛+1−𝑙}. (2.50) Since 𝜖>0 is arbitrarily small, relation (2.45) follows.
Lemma 2.7 is proved.

lim_𝑛→∞^𝑛−1_𝑗=0𝜔_𝑗^(𝑛)𝑎_𝑗+1‖𝑇_𝑎⁻¹_𝑗+1𝑓‖=0. (2.51) This together with estimate (2.43) and condition 𝑞∈(0,1) yields relation (2.38). Theorem 2.6 is proved.

Lemma 2.9 leads to an adaptive iterative scheme: 𝑢_{𝑛,𝑚𝑛}(𝑡)=𝑞𝑢_{𝑛−1,𝑚𝑛−1}+(1−𝑞)𝑇_𝑎⁻¹_{𝑛,𝑚𝑛}ℒ^∗_𝑚𝑛𝐹^(𝑚𝑛),𝑢_0,𝑚0(𝑡)=0, (2.52) where 𝑞∈(0,1), 𝑎_𝑛 are defined in (2.35), 𝑇_𝑎,𝑚 is defined in (2.30), 𝐴_𝑚ℒ is defined in (2.2), and

Lemma 2.8. Let 𝑇 and 𝑇^(𝑚) be defined in (2.33) and (2.9), respectively. Then ‖𝑇−𝑇^(𝑚)‖≤(2𝑏𝑑)⁵540√10𝑚⁴. (2.53)

Proof. From definitions (2.33) and (2.9) we get |(𝑇−𝑇^(𝑚))𝑓(𝑡)|≤^𝑏₀|^𝑑₀𝑘(𝑝,𝑡,𝑧)𝑑𝑝−^𝑚_𝑗=0𝑤_𝑗^(𝑚)𝑘(𝑝_𝑗≤,𝑡,𝑧)||𝑓(𝑧)|𝑑𝑧_𝑏0|𝑑⁵180𝑚⁴max_{𝑝∈[0,𝑑]}(𝑡+𝑧)⁴𝑒^{−𝑝(𝑡+𝑧)}=||𝑓(𝑧)|𝑑𝑧_𝑏0𝑑⁵180𝑚⁴(𝑡+𝑧)⁴𝑑|𝑓(𝑧)|𝑑𝑧≤⁵180𝑚⁴(^𝑏₀(𝑡+𝑧)⁸𝑑𝑧)^1/2‖𝑓‖_𝑋0,𝑏=𝑑⁵180𝑚⁴[(𝑡+𝑏)⁹−𝑡⁹9]^1/2‖𝑓‖_𝑋0,𝑏, (2.54) where the following upper bound for the error of the compound Simpson's rule was used (see [10, page 58]). For 𝑓∈𝐶⁽⁴⁾[𝑥₀,𝑥_2𝑙],𝑥₀<𝑥_2𝑙, |^𝑥2𝑙_𝑥0ℎ𝑓(𝑥)𝑑𝑥−3[𝑓₀+4^𝑙_𝑗=1𝑓_2(𝑗−1)+2^𝑙−1_𝑗=1𝑓_2𝑗+𝑓_𝑥2𝑙]|≤𝑅_𝑙, (2.55) where 𝑓_𝑗∶=𝑓(𝑥_𝑗),𝑥_𝑗=𝑥₀𝑥+𝑗ℎ,𝑗=0,1,2,…,2𝑙,ℎ=^2𝑙−𝑥₀,𝑅2𝑙_𝑙=(𝑥_2𝑙−𝑥₀)⁵180(2𝑙)⁴|𝑓⁽⁴⁾(𝜉)|,𝑥₀<𝜉<𝑥_2𝑙. (2.56) This implies ‖(𝑇−𝑇^(𝑚))𝑓‖_𝑋0,𝑏≤𝑑⁵540𝑚⁴[(2𝑏)¹⁰−2𝑏¹⁰]10^1/2‖𝑓‖_𝑋0,𝑏≤(2𝑏𝑑)⁵540√10𝑚⁴‖𝑓‖_𝑋0,𝑏, (2.57) so estimate (2.53) is obtained.
Lemma 2.8 is proved.

Lemma 2.9. Let 0<𝑎<𝑎₀, 𝑎𝑚=𝜅(⁰𝑎)^1/4,𝜅>0. (2.58) Then ‖𝑇−𝑇^(𝑚)‖≤(2𝑏𝑑)⁵540√10𝑎₀𝜅⁴𝑎, (2.59) where 𝑇 and 𝑇^(𝑚) are defined in (2.33) and (2.9), respectively.

Proof. Inequality (2.59) follows from estimate (2.53) and formula (2.58).

𝐹^(𝑚)∶=𝐹(𝑝₀)𝐹(𝑝₁)⋮𝐹(𝑝_𝑚)∈ℝ^𝑚+1, (2.60)𝑝_𝑗 are defined in (2.3). In the iterative scheme (2.52) we have used the finite-dimensional operator 𝑇^(𝑚) approximating the operator 𝑇. Convergence of the iterative scheme (2.52) to the solution 𝑓 of the equation ℒ𝑓=𝐹 is established in the following lemma.

Lemma 2.10. Let ℒ𝑓=𝐹 and 𝑢_{𝑛,𝑚𝑛} be defined in (2.52). If 𝑚_𝑛 are chosen by the rule 𝑚_𝑛𝑎=⌈[𝜅(⁰𝑎_𝑛)^1/4]⌉,𝑎_𝑛=𝑞𝑎_𝑛−1,𝑞∈(0,1),𝜅,𝑎₀>0, (2.61) where ⌈[𝑥]⌉ is the smallest even number not less than 𝑥, then lim_𝑛→∞‖𝑓−𝑢_{𝑛,𝑚𝑛}‖=0. (2.62)

Proof. Consider the estimate ‖𝑓−𝑢_{𝑛,𝑚𝑛}‖≤‖𝑓−𝑢_𝑛‖+‖𝑢_𝑛−𝑢_{𝑛,𝑚𝑛}‖∶=𝐼₁(𝑛)+𝐼₂(𝑛), (2.63) where 𝐼₁(𝑛)∶=‖𝑓−𝑢_𝑛‖ and 𝐼₂(𝑛)∶=‖𝑢_𝑛−𝑢_{𝑛,𝑚𝑛}‖. By Theorem 2.6, we get 𝐼₁(𝑛)→0 as 𝑛→∞. Let us prove that lim_𝑛→∞𝐼₂(𝑛)=0. Let 𝑈_𝑛∶=𝑢_𝑛−𝑢_{𝑛,𝑚𝑛}. Then, from definitions (2.31) and (2.52), we get 𝑈_𝑛=𝑞𝑈_𝑛−1+(1−𝑞)(𝑇_𝑎⁻¹_𝑛ℒ^∗𝐹−𝑇_𝑎⁻¹_{𝑛,𝑚𝑛}ℒ^∗_𝑚𝑛𝐹^(𝑚𝑛)),𝑈₀=0. (2.64) By induction we obtain 𝑈_𝑛=^𝑛−1_𝑗=0𝜔_𝑗^(𝑛)(𝑇_𝑎⁻¹_𝑗+1ℒ^∗𝐹−𝑇_𝑎⁻¹_{𝑗+1,𝑚𝑗+1}(ℒ_𝑚𝑗+1)^∗𝐹^(𝑚𝑗+1)), (2.65) where 𝜔_𝑗 are defined in (2.40). Using the identities ℒ𝑓=𝐹, ℒ_𝑚𝑓=𝐹^(𝑚), 𝑇_𝑎⁻¹𝑇=𝑇_𝑎⁻¹(𝑇+𝑎𝐼−𝑎𝐼)=𝐼−𝑎𝑇_𝑎⁻¹,𝑇⁻¹_𝑎,𝑚𝑇^(𝑚)=𝑇⁻¹_𝑎,𝑚(𝑇^(𝑚)+𝑎𝐼−𝑎𝐼)=𝐼−𝑎𝑇⁻¹_𝑎,𝑚,𝑇⁻¹_𝑎,𝑚−𝑇_𝑎⁻¹=𝑇⁻¹_𝑎,𝑚(𝑇−𝑇^(𝑚))𝑇_𝑎⁻¹, (2.66) one gets 𝑈_𝑛=^𝑛−1_𝑗=0𝜔_𝑗^(𝑛)𝑎_𝑗+1(𝑇_𝑎⁻¹_{𝑗+1,𝑚𝑗+1}−𝑇_𝑎⁻¹_𝑗+1=)𝑓_{𝑛−1𝑗=0}𝜔_𝑗^(𝑛)𝑎_𝑗+1𝑇_𝑎⁻¹_{𝑗+1,𝑚𝑗+1}(𝑇−𝑇^(𝑚𝑗+1))𝑇_𝑎⁻¹_𝑗+1𝑓. (2.67) This together with the rule (2.61), estimate (2.28), and Lemma 2.8 yields ‖𝑈_𝑛‖≤^𝑛−1_𝑗=0𝜔_𝑗^(𝑛)𝑎_𝑗+1‖𝑇_𝑎⁻¹_{𝑗+1,𝑚𝑗+1}‖‖𝑇−𝑇^(𝑚𝑗+1)‖‖𝑇_𝑎⁻¹_𝑗+1≤𝑓‖(2𝑏𝑑)⁵540√10𝑎₀𝜅^4𝑛−1_𝑗=0𝜔_𝑗^(𝑛)𝑎_𝑗+1‖𝑇_𝑎⁻¹_𝑗+1𝑓‖. (2.68) Applying Lemmas 2.5 and 2.7 with 𝑔(𝑎)=𝑎‖𝑇_𝑎⁻¹𝑓‖, we obtain lim_𝑛→∞‖𝑈_𝑛‖=0.
Lemma 2.10 is proved.

2.1. Noisy Data

When the data 𝐹(𝑝) are noisy, the approximate solution (2.23) is written as

𝑓^𝛿_𝑚(𝑡)=^𝑚_𝑗=0𝑤_𝑗^(𝑚)𝑐_𝑗^(𝑚,𝛿)𝑒^{−𝑝𝑗𝑡}=𝑇⁻¹_𝑎,𝑚ℒ^∗_𝑚𝐹_𝛿^(𝑚), (2.69) where the coefficients 𝑐_𝑗^(𝑚,𝛿) are obtained by solving the following linear algebraic system:

𝑄_𝑎,𝑚𝑐^(𝑚,𝛿)=𝐹_𝛿^(𝑚), (2.70)𝑄_𝑎,𝑚 is defined in (2.30),

𝑐^(𝑚,𝛿)𝑐∶=^0(𝑚,𝛿)𝑐₁^(𝑚,𝛿)⋮𝑐_𝑚^(𝑚,𝛿),𝐹_𝛿^(𝑚)𝐹∶=^𝛿(𝑝₀)𝐹_𝛿(𝑝₁)⋮𝐹_𝛿(𝑝_𝑚), (2.71)𝑤_𝑗^(𝑚) are defined in (2.5), and 𝑝_𝑗 are defined in (2.3).

To get the approximation solution of the function 𝑓(𝑡) with the noisy data 𝐹_𝛿(𝑝), we consider the following iterative scheme:

𝑢^𝛿_{𝑛,𝑚𝑛}=𝑞𝑢^𝛿_{𝑛−1,𝑚𝑛−1}+(1−𝑞)𝑇_𝑎⁻¹_{𝑛,𝑚𝑛}ℒ^∗_𝑚𝑛𝐹^(𝑚𝑛)_𝛿,𝑢^𝛿_0,𝑚0=0, (2.72) where 𝑇_𝑎,𝑚 is defined in (2.30), 𝑎_𝑛 are defined in (2.35), 𝑞∈(0,1), 𝐹_𝛿^(𝑚) is defined in (2.71), and 𝑚_𝑛 are chosen by the rule (2.61). Let us assume that

𝐹_𝛿(𝑝_𝑗)=𝐹(𝑝_𝑗)+𝛿_𝑗,0<|𝛿_𝑗|≤𝛿,𝑗=0,1,2,…,𝑚, (2.73) where 𝛿_𝑗 are random quantities generated from some statistical distributions, for example, the uniform distribution on the interval [−𝛿,𝛿], and 𝛿 is the noise level of the data 𝐹(𝑝). It follows from assumption (2.73), definition (2.5), Lemma 2.1, and the inner product (2.4) that

‖𝐹_𝛿^(𝑚)−𝐹^(𝑚)‖²_𝑊𝑚=^𝑚_𝑗=0𝑤_𝑗^(𝑚)𝛿²_𝑗≤𝛿^2𝑚_𝑗=0𝑤_𝑗^(𝑚)=𝛿²𝑑. (2.74)

Lemma 2.11. Let 𝑢_{𝑛,𝑚𝑛} and 𝑢^𝛿_{𝑛,𝑚𝑛} be defined in (2.52) and (2.72), respectively. Then ‖𝑢_{𝑛,𝑚𝑛}−𝑢^𝛿_{𝑛,𝑚𝑛}√‖≤𝑑𝛿2√𝑎_𝑛(1−𝑞^𝑛),𝑞∈(0,1), (2.75) where 𝑎_𝑛 are defined in (2.35).

Proof. Let 𝑈^𝛿_𝑛∶=𝑢_{𝑛,𝑚𝑛}−𝑢^𝛿_{𝑛,𝑚𝑛}. Then, from definitions (2.52) and (2.72), 𝑈^𝛿_𝑛=𝑞𝑈^𝛿_𝑛−1+(1−𝑞)𝑇_𝑎⁻¹_{𝑛,𝑚𝑛}ℒ^∗_𝑚𝑛(𝐹^(𝑚𝑛)−𝐹^(𝑚𝑛)_𝛿),𝑈^𝛿₀=0. (2.76) By induction we obtain 𝑈^𝛿_𝑛=^𝑛−1_𝑗=0𝜔_𝑗^(𝑛)𝑇_𝑎⁻¹_{𝑗+1,𝑚𝑗+1}(ℒ_𝑚𝑗+1)^∗(𝐹^(𝑚𝑗+1)−𝐹^(𝑚𝑗+1)_𝛿), (2.77) where 𝜔_𝑗^(𝑛) are defined in (2.40). Using estimates (2.74) and inequality (2.29), one gets ‖𝑈^𝛿_𝑛‖≤√𝑑^𝑛−1_𝑗=0𝜔_𝑗^(𝑛)𝛿2𝑎_𝑗+1≤√𝑑𝛿2√𝑎_𝑛^𝑚_𝑗=0𝜔_𝑗^(𝑛)=√𝑑𝛿2√𝑎_𝑛(1−𝑞^𝑛), (2.78) where 𝜔_𝑗 are defined in (2.40).
Lemma 2.11 is proved.

Theorem 2.12. Suppose that conditions of Lemma 2.10 hold, and 𝑛_𝛿 satisfies the following conditions: lim_𝛿→0𝑛_𝛿=∞,lim_𝛿→0𝛿𝑎_𝑛𝛿=0. (2.79) Then lim_𝛿→0‖𝑓−𝑢^𝛿_{𝑛𝛿,𝑚𝑛𝛿}‖=0. (2.80)

Proof. Consider the estimate ‖𝑓−𝑢^𝛿_{𝑛𝛿,𝑚𝑛𝛿}‖≤‖𝑓−𝑢_{𝑛𝛿,𝑚𝑛𝛿}‖+‖𝑢_{𝑛𝛿,𝑚𝑛𝛿}−𝑢^𝛿_{𝑛𝛿,𝑚𝑛𝛿}‖. (2.81) This together with Lemma 2.11 yields ‖𝑓−𝑢^𝛿_{𝑛𝛿,𝑚𝑛𝛿}‖≤‖𝑓−𝑢_{𝑛𝛿,𝑚𝑛𝛿}√‖+𝑑𝛿2𝑎_𝑛𝛿(1−𝑞^𝑛). (2.82) Applying relations (2.79) in estimate (2.82), one gets relation (2.80).
Theorem 2.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

𝐺_{𝑛𝛿,𝑚𝑛𝛿}≤𝐶𝛿^𝜀<𝐺_{𝑛,𝑚𝑛},1≤𝑛<𝑛_𝛿,𝐶>√𝑑,𝜀∈(0,1), (2.83) where

𝐺_{𝑛,𝑚𝑛}=𝑞𝐺_{𝑛−1,𝑚𝑛−1}+(1−𝑞)‖ℒ_𝑚𝑛𝑧^{(𝑚𝑛,𝛿)}−𝐹^(𝑚𝑛)_𝛿‖_𝑊𝑚𝑛,𝐺_0,𝑚0=0, (2.84)‖⋅‖_𝑊^𝑚 is defined in (2.4),

𝑧^(𝑚,𝛿)∶=^𝑚_𝑗=0𝑐_𝑗^(𝑚,𝛿)𝑤_𝑗^(𝑚)𝑒^{−𝑝𝑗𝑡}, (2.85)𝑤_𝑗^(𝑚) and 𝑝_𝑗 are defined in (2.5) and (2.3), respectively, and 𝑐_𝑗^(𝑚,𝛿) are obtained by solving linear algebraic system (2.70).

We observe that

ℒ_𝑚𝑛𝑧^{(𝑚𝑛,𝛿)}−𝐹^(𝑚𝑛)_𝛿=𝑄^(𝑚𝑛)𝑐^{(𝑚𝑛,𝛿)}−𝐹^(𝑚𝑛)_𝛿=𝑄^(𝑚𝑛)(𝑎_𝑛𝐼+𝑄^(𝑚𝑛))⁻¹𝐹^(𝑚𝑛)_𝛿−𝐹^(𝑚𝑛)_𝛿=(𝑄^(𝑚𝑛)+𝑎_𝑛𝐼−𝑎_𝑛𝐼)(𝑎_𝑛𝐼+𝑄^(𝑚𝑛))⁻¹𝐹^(𝑚𝑛)_𝛿−𝐹^(𝑚𝑛)_𝛿=−𝑎_𝑛(𝑎_𝑛𝐼+𝑄^(𝑚𝑛))⁻¹𝐹^(𝑚𝑛)_𝛿=−𝑎_𝑛𝑐^{(𝑚𝑛,𝛿)}. (2.86) Thus, the sequence (2.84) can be written in the following form:

𝐺_{𝑛,𝑚𝑛}=𝑞𝐺_{𝑛−1,𝑚𝑛−1}+(1−𝑞)𝑎_𝑛‖𝑐^{(𝑚𝑛,𝛿)}‖_𝑊^𝑚_𝑛,𝐺_0,𝑚0=0, (2.87) where ‖⋅‖_𝑊^𝑚 is defined in (2.4), and 𝑐^(𝑚,𝛿) solves the linear algebraic system (2.70).

It follows from estimates (2.74), (2.26), and (2.27) that

𝑎_𝑛‖𝑐^{(𝑚𝑛,𝛿)}‖_𝑊^𝑚_𝑛=𝑎_𝑛‖(𝑎_𝑛𝐼+𝑄^(𝑚𝑛))⁻¹𝐹^(𝑚𝑛)_𝛿‖_𝑊𝑚𝑛≤𝑎_𝑛‖(𝑎_𝑛𝐼+𝑄^(𝑚𝑛))⁻¹(𝐹^(𝑚𝑛)_𝛿−𝐹^(𝑚𝑛))‖_𝑊𝑚𝑛+𝑎_𝑛‖(𝑎_𝑛𝐼+𝑄^(𝑚𝑛))⁻¹𝐹^(𝑚𝑛)‖_𝑊𝑚𝑛≤‖𝐹^(𝑚𝑛)_𝛿−𝐹^(𝑚𝑛)‖_𝑊𝑚𝑛+𝑎_𝑛‖(𝑎_𝑛𝐼+𝑄^(𝑚𝑛))⁻¹ℒ_𝑚𝑛𝑓‖_𝑊𝑚𝑛≤𝛿√𝑑+√𝑎_𝑛‖𝑓‖_𝑋0,𝑏. (2.88) This together with (2.87) yields

𝐺_{𝑛,𝑚𝑛}≤𝑞𝐺_{𝑛−1,𝑚𝑛−1}+(1−𝑞)(𝛿√𝑑+√𝑎_𝑛‖𝑓‖_𝑋0,𝑏), (2.89) or

𝐺_{𝑛,𝑚𝑛}−𝛿√𝑑≤𝑞(𝐺_{𝑛−1,𝑚𝑛−1}−𝛿√𝑑)+(1−𝑞)√𝑎_𝑛‖𝑓‖_𝑋0,𝑏. (2.90)

Lemma 2.13. The sequence (2.87) satisfies the following estimate: 𝐺_{𝑛,𝑚𝑛}−𝛿√√𝑑≤(1−𝑞)𝑎_𝑛‖𝑓‖_𝑋0,𝑏1−√𝑞, (2.91) where 𝑎_𝑛 are defined in (2.35).

Proof. Define Ψ_𝑛∶=𝐺_{𝑛,𝑚𝑛}−𝛿√𝑑,𝜓_𝑛∶=(1−𝑞)√𝑎_𝑛‖𝑓‖_𝑋0,𝑏. (2.92) Then estimate (2.90) can be rewritten as Ψ_𝑛≤𝑞Ψ_𝑛−1+√𝑞𝜓_𝑛−1, (2.93) where the relation 𝑎_𝑛=𝑞𝑎_𝑛−1 was used. Let us prove estimate (2.91) by induction. For 𝑛=0 we get Ψ₀=−𝛿√√𝑑≤(1−𝑞)𝑎₀‖𝑓‖_𝑋0,𝑏1−√𝑞. (2.94) Suppose that estimate (2.91) is true for 0≤𝑛≤𝑘. Then Ψ_𝑘+1≤𝑞Ψ_𝑘+√𝑞𝜓_𝑘≤𝑞1−√𝑞𝜓_𝑘+√𝑞𝜓_𝑘=√𝑞1−√𝑞𝜓_𝑘=√𝑞1−√𝑞𝜓_𝑘𝜓_𝑘+1𝜓_𝑘+1=√𝑞1−√𝑞√𝑎_𝑘√𝑎_𝑘+1𝜓_𝑘+1=11−√𝑞𝜓_𝑘+1, (2.95) where the relation 𝑎_𝑘+1=𝑞𝑎_𝑘 was used.
Lemma 2.13 is proved.

Lemma 2.14. Suppose 𝐺_1,𝑚1>𝛿√𝑑, (2.96) where 𝐺_{𝑛,𝑚𝑛} are defined in (2.87). Then there exists a unique integer 𝑛_𝛿, satisfying the stopping rule (2.83) with 𝐶>√𝑑.

Proof. From Lemma 2.13 we get the estimate 𝐺_{𝑛,𝑚𝑛}≤𝛿√√𝑑+(1−𝑞)𝑎_𝑛‖𝑓‖_𝑋0,𝑏1−√𝑞, (2.97) where 𝑎_𝑛 are defined in (2.35). Therefore, limsup_𝑛→∞𝐺_{𝑛,𝑚𝑛}≤𝛿√𝑑, (2.98) where the relation lim_𝑛→∞𝑎_𝑛=0 was used. This together with condition (2.96) yields the existence of the integer 𝑛_𝛿. The uniqueness of the integer 𝑛_𝛿 follows from its definition.
Lemma 2.14 is proved.

Lemma 2.15. Suppose that conditions of Lemma 2.14 hold and 𝑛_𝛿 is chosen by the rule (2.83), then lim_𝛿→0𝛿𝑎_𝑛𝛿=0. (2.99)

Proof. From the stopping rule (2.83) and estimate (2.97) we get 𝐶𝛿^𝜀≤𝐺_{𝑛𝛿−1,𝑚𝑛𝛿−1}≤𝛿√𝑑+(1−𝑞)𝑎_𝑛𝛿−1‖𝑓‖_𝑋0,𝑏1−√𝑞, (2.100) where 𝐶>√𝑑,𝜀∈(0,1). This implies 𝛿(𝐶𝛿^𝜀−1−√𝑑)𝑎_𝑛𝛿−1≤(1−𝑞)‖𝑓‖_𝑋0,𝑏1−√𝑞, (2.101) so, for 𝜀∈(0,1), and 𝑎_𝑛𝛿=𝑞𝑎_𝑛𝛿−1, one gets lim_𝛿→0𝛿𝑎_𝑛𝛿=lim_𝛿→0𝛿√𝑞𝑎_𝑛𝛿−1≤lim_𝛿→0(1−𝑞)𝛿^1−𝜀‖𝑓‖_𝑋0,𝑏(√𝑞−𝑞)(𝐶−𝛿^1−𝜀√𝑑)=0. (2.102)
Lemma 2.15 is proved.

Lemma 2.16. Consider the stopping rule (2.83), where the parameters 𝑚_𝑛 are chosen by rule (2.61). If 𝑛_𝛿 is chosen by the rule (2.83) then lim_𝛿→0𝑛_𝛿=∞. (2.103)

Proof. From the stopping rule (2.83) with the sequence 𝐺_𝑛 defined in (2.87) one gets 𝑞𝐶𝛿^𝜀+(1−𝑞)𝑎_𝑛𝛿‖𝑐^{(𝑚𝑛𝛿,𝛿)}‖_𝑊^𝑚𝑛_𝛿≤𝑞𝐺_{𝑛𝛿−1,𝑚𝑛𝛿−1}+(1−𝑞)𝑎_𝑛𝛿‖𝑐^{(𝑚𝑛𝛿,𝛿)}‖_𝑊^𝑚𝑛_𝛿=𝐺_{𝑛𝛿,𝑚𝑛𝛿}<𝐶𝛿^𝜀, (2.104) where 𝑐^(𝑚,𝛿) is obtained by solving linear algebraic system (2.70). This implies 0<𝑎_𝑛𝛿‖𝑐^{(𝑚𝑛𝛿,𝛿)}‖_𝑊^𝑚𝑛_𝛿≤𝐶𝛿^𝜀. (2.105) Thus, lim_𝛿→0𝑎_𝑛𝛿‖𝑐^{(𝑚𝑛𝛿,𝛿)}‖_𝑊^𝑚𝑛_𝛿=0. (2.106) If 𝐹^(𝑚)≠0, then there exists a 𝜆₀^(𝑚)>0 such that 𝐸_𝜆^(𝑚)_0(𝑚)𝐹^(𝑚)≠0,⟨𝐸_𝜆^(𝑚)₀𝐹^(𝑚),𝐹^(𝑚)⟩_𝑊𝑚∶=𝜉^(𝑚)>0, (2.107) where 𝐸_𝑠^(𝑚) is the resolution of the identity corresponding to the operator 𝑄^(𝑚)∶=ℒ_𝑚ℒ^∗_𝑚. Let ℎ_𝑚(𝛿,𝛼)∶=𝛼²‖𝑄⁻¹_𝑚,𝛼𝐹_𝛿^(𝑚)‖²_𝑊𝑚,𝑄_𝑚,𝑎∶=𝑎𝐼+𝑄^(𝑚). (2.108) For a fixed number 𝑎>0 we obtain ℎ_𝑚(𝛿,𝑎)=𝑎²‖𝑄⁻¹_𝑚,𝑎𝐹_𝛿^(𝑚)‖²_𝑊𝑚=^∞₀𝑎²(𝑎+𝑠)²𝑑⟨𝐸_𝑠^(𝑚)𝐹_𝛿^(𝑚),𝐹_𝛿^(𝑚)⟩_𝑊𝑚≥^𝜆0(𝑚)₀𝑎²(𝑎+𝑠)²𝑑⟨𝐸_𝑠^(𝑚)𝐹_𝛿^(𝑚),𝐹_𝛿^(𝑚)⟩_𝑊𝑚≥𝑎²(𝑎+𝜆₀)^{2𝜆0(𝑚)0}𝑑⟨𝐸_𝑠^(𝑚)𝐹_𝛿^(𝑚),𝐹_𝛿^(𝑚)⟩_𝑊𝑚=𝑎²‖𝐸_𝜆^(𝑚)_0(𝑚)𝐹_𝛿^(𝑚)‖²_𝑊𝑚(𝑎+𝜆₀^(𝑚))². (2.109) Since 𝐸_𝜆^(𝑚)₀ is a continuous operator, and ‖𝐹^(𝑚)−𝐹_𝛿^(𝑚)‖_𝑊𝑚<√𝑑𝛿, it follows from (2.107) that lim_𝛿→0⟨𝐸_𝜆^(𝑚)₀𝐹_𝛿^(𝑚),𝐹_𝛿^(𝑚)⟩_𝑊𝑚=⟨𝐸_𝜆^(𝑚)₀𝐹^(𝑚),𝐹^(𝑚)⟩_𝑊𝑚>0. (2.110) Therefore, for the fixed number 𝑎>0 we get ℎ_𝑚(𝛿,𝑎)≥𝑐₂>0 (2.111) for all sufficiently small 𝛿>0, where 𝑐₂ is a constant which does not depend on 𝛿. Suppose lim_𝛿→0𝑎_𝑛𝛿≠0. Then there exists a subsequence 𝛿_𝑗→0 as 𝑗→∞, such that 𝑎_𝑛𝛿𝑗≥𝑐₁>0,0<𝑚_𝑛𝛿𝑗𝑎=⌈[𝜅(⁰𝑎_𝑛𝛿𝑗)^1/4𝑎]⌉≤⌈[𝜅(⁰𝑐₁)^1/4]⌉∶=𝑐₃<∞,𝜅,𝑎₀>0, (2.112) where rule (2.61) was used to obtain the parameters 𝑚_𝑛𝛿𝑗. This together with (2.107) and (2.111) yields lim_𝑗→∞ℎ_{𝑚𝑛𝛿𝑗}(𝛿_𝑗,𝑎_𝑛𝛿𝑗)≥lim_𝑗→∞𝑎²_𝑛𝛿𝑗‖𝐸^{(𝑚𝑛𝛿𝑗)}_{𝜆𝑛𝛿𝑗)0(𝑚}𝐹^{(𝑚𝑛𝛿𝑗)}_𝛿𝑗‖²_{𝑊𝑚𝑛𝛿𝑗}(𝑎_𝑛𝛿𝑗+𝜆^{(𝑚𝑛𝛿𝑗)}₀)²≥liminf_𝑗→∞𝑐²₁‖𝐸^{(𝑚𝑛𝛿𝑗)}_{𝜆𝑛𝛿𝑗)0(𝑚}𝐹^{(𝑚𝑛𝛿𝑗)}‖²_{𝑊𝑚𝑛𝛿𝑗}(𝑐₁+𝜆^{(𝑚𝑛𝛿𝑗)}₀)²>0. (2.113) This contradicts relation (2.106). Thus, lim_𝛿→0𝑎_𝑛𝛿=lim_𝛿→0𝑎₀𝑞^𝑛𝛿=0, that is, lim_𝛿→0𝑛_𝛿=∞.
Lemma 2.16 is proved.

It follows from Lemmas 2.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 Theorem 2.12 hold, 𝑚_𝑛 are chosen by rule (2.61), 𝑛_𝛿 is chosen by rule (2.83), and 𝐺_1,𝑚1>𝐶𝛿, where 𝐺_{𝑛,𝑚𝑛} are defined in (2.87), then lim_𝛿→0‖𝑓−𝑢^𝛿_{𝑛𝛿,𝑚𝑛𝛿}‖=0. (2.114)

2.3. The Algorithm

Let us formulate the algorithm for obtaining the approximate solution 𝑓^𝛿_𝑚:

(1)The data 𝐹_𝛿(𝑝) on the interval [0,𝑑], 𝑑>0, the support of the function 𝑓(𝑡), and the noise level 𝛿. (2)Initialization: choose the parameters 𝜅>0, 𝑎₀>0, 𝑞∈(0,1), 𝜀∈(0,1), 𝐶>√𝑑, and set 𝑢^𝛿_0,𝑚0=0, 𝐺₀=0, 𝑛=1. (3)Iterate, starting with 𝑛=1, and stop when condition (2.120) (see below) holds, (a)𝑎_𝑛=𝑎₀𝑞^𝑛, (b)choose 𝑚_𝑛 by the rule (2.61), (c)construct the vector 𝐹^(𝑚𝑛)_𝛿: (𝐹^(𝑚𝑛)_𝛿)_𝑙=𝐹_𝛿(𝑝_𝑙),𝑝_𝑙𝑑=𝑙ℎ,ℎ=𝑚_𝑛,𝑙=0,1,…,𝑚, (2.115)(d)construct the matrices 𝐻_𝑚𝑛 and 𝐷_𝑚𝑛: (𝐻_𝑚𝑛)_𝑖𝑗∶=^𝑏₀𝑒^{−(𝑝𝑖+𝑝𝑗)𝑡}𝑑𝑡=1−𝑒^{−𝑏(𝑝𝑖+𝑝𝑗)}𝑝_𝑖+𝑝_𝑗,𝑖,𝑗=1,2,3,…,𝑚_𝑛,(𝐷_𝑚𝑛)_𝑖𝑗=𝑤^(𝑚𝑛)_𝑖,𝑖=𝑗,0,otherwise, (2.116) where 𝑤_𝑗^(𝑚) are defined in (2.5), (e)solve the following linear algebraic system: (𝑎_𝑛𝐼+𝐻_𝑚𝑛𝐷_𝑚𝑛)𝑐^{(𝑚𝑛,𝛿)}=𝐹^(𝑚𝑛)_𝛿, (2.117) where (𝑐^{(𝑚𝑛,𝛿)})_𝑖=𝑐^(𝑚𝑛_𝑖^,𝛿), (f)update the coefficient 𝑐^(𝑚𝑛_𝑗^,𝛿) of the approximate solution 𝑢^𝛿_{𝑛,𝑚𝑛}(𝑡) defined in (2.69) by the iterative formula: 𝑢^𝛿_{𝑛,𝑚𝑛}(𝑡)=𝑞𝑢^𝛿_{𝑛−1,𝑚𝑛−1}(𝑡)+(1−𝑞)^𝑚𝑛_𝑗=1𝑐^{(𝑚𝑛,𝛿)}𝑤^(𝑚𝑛)_𝑗𝑒^{−𝑝𝑗𝑡}, (2.118) where 𝑢^𝛿_0,𝑚0(𝑡)=0. (2.119)

Stop when for the first time the inequality

𝐺_{𝑛,𝑚𝑛}=𝑞𝐺_{𝑛−1,𝑚𝑛−1}+𝑎_𝑛‖𝑐^{(𝑚𝑛,𝛿)}‖_𝑊^𝑚_𝑛≤𝐶𝛿^𝜀 (2.120) holds, and get the approximation 𝑓^𝛿(𝑡)=𝑢^𝛿_{𝑛𝛿,𝑚𝑛𝛿}(𝑡) of the function 𝑓(𝑡) by formula (2.118).

3. Numerical Experiments

3.1. The Parameters 𝜅,𝑎₀,𝑑

From definition (2.35) and the rule (2.61) we conclude that 𝑚_𝑛→∞ as 𝑎_𝑛→0. Therefore, one needs to control the value of the parameter 𝑚_𝑛 so that it will not grow too fast as 𝑎_𝑛 decreases. The role of the parameter 𝜅 in (2.61) is to control the value of the parameter 𝑚_𝑛 so that the value of the parameter 𝑚_𝑛 will not be too large. Since for sufficiently small noise level 𝛿, namely, 𝛿∈(10⁻¹⁶,10⁻⁶], the regularization parameter 𝑎_𝑛𝛿, obtained by the stopping rule (2.83), is at most 𝑂(10⁻⁹), we suggest to choose 𝜅 in the interval (0,1]. For the noise level 𝛿∈(10⁻⁶,10⁻²] one can choose 𝜅∈(1,3]. To reduce the number of iterations we suggest to choose the geometric sequence 𝑎_𝑛=𝑎₀𝛿^𝛼𝑛, where 𝑎₀∈[0.1,0.2] and 𝛼∈[0.5,0.9]. One may assume without loss of generality that 𝑏=1, because a scaling transformation reduces the integral over (0,𝑏) to the integral over (0,1). We have assumed that the data 𝐹(𝑝) are defined on the interval 𝐽∶=[0,𝑑]. In the case the interval 𝐽=[𝑑₁,𝑑], 0<𝑑₁<𝑑, the constant 𝑑 in estimates (2.59), (2.74), (2.75), (2.78), (2.90), (2.91), and (2.97) are replaced with the constant 𝑑−𝑑₁. If 𝑏=1, that is, 𝑓(𝑡)=0 for 𝑡>1, then one has to take 𝑑 not too large. Indeed, if 𝑓(𝑡)=0 for 𝑡>1, then an integration by parts yields 𝐹(𝑝)=[𝑓(0)−𝑒^−𝑝𝑓(1)]/𝑝+𝑂(1/𝑝²),𝑝→∞. If the data are noisy, and the noise level is 𝛿, then the data becomes indistinguishable from noise for 𝑝=𝑂(1/𝛿). Therefore it is useless to keep the data 𝐹_𝛿(𝑝) for 𝑑>𝑂(1/𝛿). In practice one may get a satisfactory accuracy of inversion by the method, proposed in this paper, when one uses the data with 𝑑∈[1,20] when 𝛿≤10⁻². In all the numerical examples we have used 𝑑=5. Given the interval [0,𝑑], the proposed method generates automatically the discrete data 𝐹_𝛿(𝑝_𝑗), 𝑗=0,1,2,…,𝑚, over the interval [0,𝑑] which are needed to get the approximation of the function 𝑓(𝑡).

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 𝐹_𝛿(𝑝) with various noise levels 𝛿=10⁻²,𝛿=10⁻⁴, and 𝛿=10⁻⁶. The random quantities 𝛿_𝑗 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: 𝑎_𝑛=0.1𝑞^𝑛, 𝑞=𝛿^1/2, and 𝑑=5. The parameter 𝜅=1 is used for the noise levels 𝛿=10⁻² and 𝛿=10⁻⁴. When 𝛿=10⁻⁶ we choose 𝜅=0.3 so that the value of the parameters 𝑚_𝑛 are not very large, namely, 𝑚_𝑛≤300. Therefore, the computation time for solving linear algebraic system (2.117) can be reduced significantly. We assume that the support of the function 𝑓(𝑡) is in the interval [0,𝑏] with 𝑏=10. In the stopping rule (2.83) the following parameters are used: 𝐶=√𝑑+0.01, 𝜀=0.99. In Example 3.13 the function 𝑓(𝑡)=𝑒^−𝑡 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

∑MAE∶=[^100𝑗=1(𝑓(𝑡_𝑖)−𝑓^𝛿_𝑚𝑛𝛿(𝑡_𝑖))²]100^1/2,𝑡_𝑗=0.01+0.1(𝑗−1),𝑗=1,…,100, (3.1) where 𝑓(𝑡) is the exact solution and 𝑓^𝛿_𝑚𝑛𝛿(𝑡) is the approximate solution. The computation time (CPU time) for obtaining the approximation of 𝑓(𝑡), the number of iterations (Iter.), and the parameters 𝑚_𝑛𝛿 and 𝑎_𝑛𝛿 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 

Example 3.1.

Table 2 

Example 3.2.

Table 3 

Example 3.3.

Table 4 

Example 3.4.

Table 5 

Example 3.5.

Table 6 

Example 3.6.

Table 7 

Example 3.7.

Table 8 

Example 3.8.

Table 9 

Example 3.9.

Table 10 

Example 3.10.

Table 11 

Example 3.11.

Table 12 

Example 3.12.

Table 13 

Example 3.13.

Example 3.1 (see [15]). We have 𝑓₁1(𝑡)=1,23≤𝑡≤2,𝐹0,otherwise,₁(𝑝)=𝑒1,𝑝=0,_−𝑝/2−𝑒^−3𝑝/2𝑝,𝑝>0. (3.2) The reconstruction of the exact solution for different values of the noise level 𝛿 is shown in Figure 1. When the noise level 𝛿=10⁻⁶, 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 

Example 3.1: the stability of the approximate solution.

Example 3.2 (see [13, 15]). We have 𝑓₂1(𝑡)=2𝐹,𝑡=1,1,1<𝑡<10,0,elsewhere,₂(𝑝)=𝑒9,𝑝=0,_−𝑝−𝑒^−10𝑝𝑝,𝑝>0. (3.3) The reconstruction of the function 𝑓₂(𝑡) 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⁻⁶ and 𝜅=0.3 the value of the parameter 𝑚_𝑛𝛿 is bounded by the constant 54.

Figure 2 

Example 3.2: the stability of the approximate solution.

Example 3.3 (see [6, 12, 13, 16, 19]). We have 𝑓₃(𝑡)=𝑡𝑒^−𝑡𝐹,0≤𝑡<10,0,otherwise,₃(𝑝)=1−𝑒^{−(𝑝+1)10}(𝑝+1)²−10𝑒^{−(𝑝+1)10}.𝑝+1 (3.4) We get an excellent agreement between the approximate solution and the exact solution when the noise level 𝛿=10⁻⁴ and 10⁻⁶ 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 𝑚_𝑛𝛿 is bounded by the constant 30 when the noise level 𝛿=10⁻⁶ and 𝜅=0.3.

Figure 3 

Example 3.3: the stability of the approximate solution.

Figure 4 

Example 3.4: the stability of the approximate solution.

Example 3.4 (see [13, 15]). We have 𝑓₄(𝑡)=1−𝑒^−0.5𝑡𝐹,0≤𝑡<10,0,elsewhere,₄(𝑝)=8+2𝑒⁻⁵,𝑝=0,1−𝑒^−10𝑝𝑝−1−𝑒^{−(𝑝+1/2)10}𝑝+0.5,𝑝>0. (3.5) As in our Example 3.3 when the noise 𝛿=10⁻⁴ and 10⁻⁶ 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 𝑓₄(𝑡) obtained by the proposed method is better than the reconstruction of 𝑓₄(𝑡) shown in [13], and is comparable with the double precision reconstruction obtained in [15].
In this example when 𝛿=10⁻⁶ and 𝜅=0.3 the value of the parameter 𝑚_𝑛𝛿 is bounded by the constant 109 as shown in Table 4.

Example 3.5 (see [5, 7, 13]). We have 𝑓₅2(𝑡)=√3𝑒^−𝑡/2sin(𝑡√𝐹3/2)₅√(𝑝)=1−cos(103/2)𝑒^{−10(𝑝+0.5)}[(𝑝+0.5)²−+3/4]2(𝑝+0.5)𝑒^{−10(𝑝+0.5)}sin(10√3/2)√3[(𝑝+0.5)².+3/4] (3.6) This is an example of the damped sine function. In [5, 7] the knowledge of the exact data 𝐹(𝑝) in the complex plane is required to get the approximate solution. Here we only use the knowledge of the discrete perturbed data 𝐹_𝛿(𝑝_𝑗), 𝑗=0,1,2,…,𝑚, and get a satisfactory result which is comparable with the results given in [5, 7] when the level noise 𝛿=10⁻⁶. The reconstruction of the exact solution 𝑓₅(𝑡) 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 𝑚_𝑛𝛿 is bounded by 54 for the noise level 𝛿=10⁻⁶ (see Table 5).

Figure 5 

Example 3.5: the stability of the approximate solution.

Example 3.6 (see [15]). We have 𝑓₆3(𝑡)=𝑡,0≤𝑡<1,2−𝑡2𝐹,1≤𝑡<3,0,elsewhere,₆3(𝑝)=2,𝑝=0,1−𝑒^−𝑝(1+𝑝)𝑝²+𝑒^−3𝑝+𝑒^2𝑝(2𝑝−1)2𝑝²,𝑝>0. (3.7) 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⁻⁶, 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 𝑚 generated by the proposed adaptive stopping rule is bounded by the constant 54 for the noise level 𝛿=10⁻⁶ and 𝜅=0.3 which gives a relatively small computation time.

Figure 6 

Example 3.6: the stability of the approximate solution.

Example 3.7 (see [15]). We have 𝑓₇(𝑡)=−𝑡𝑒^−𝑡−𝑒^−𝑡+1,0≤𝑡<1,1−2𝑒⁻¹𝐹,1≤𝑡<10,0,elsewhere,₇3(𝑝)=𝑒2−1+9(1−𝑒𝑒),𝑝=0,_{−1−𝑝(𝑒1+𝑝−𝑒(1+𝑝)2+𝑝(3+2𝑝)/𝑝(𝑝+1)2)}+(𝑒−2)𝑒^{−1−𝑝−10𝑝((𝑒10𝑝−𝑒𝑝)/𝑝)},𝑝>0. (3.8) When the noise level 𝛿=10⁻⁴ and 𝛿=10⁻⁶, 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 

Example 3.7: the stability of the approximate solution.

Example 3.8 (see [13, 25]). We have 𝑓₈(𝑡)=4𝑡²𝑒^−2𝑡𝐹,0≤𝑡<10,0,elsewhere,₈(𝑝)=8+4𝑒^{−10(2+𝑝)}[−2−20(2+𝑝)−100(2−𝑝)²](2+𝑝)³. (3.9) 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⁻⁴ and 𝛿=10⁻⁶ (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 

Example 3.8: the stability of the approximate solution.

Example 3.9 (see [16]). We have 𝑓₉𝐹(𝑡)=5−𝑡,0≤𝑡<5,0,elsewhere,₉(𝑝)=252𝑒,𝑝=0,_−5𝑝+5𝑝−1𝑝²,𝑝>0. (3.10) 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 

Example 3.9: the stability of the approximate solution.

Example 3.10 (see [6]). We have 𝑓₁₀𝐹(𝑡)=𝑡,0≤𝑡<10,0,elsewhere,₁₀(𝑝)=50,𝑝=0,1−𝑒^−10𝑝𝑝²−10𝑒^−10𝑝𝑝,𝑝>0. (3.11) 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 

Example 3.10: the stability of the approximate solution.

Example 3.11 (see [6, 9]). We have 𝑓₁₁𝐹(𝑡)=sin(𝑡),0≤𝑡<10,0,elsewhere,₁₁(𝑝)=1−𝑒^−10𝑝(𝑝sin(10)+cos(10))1+𝑝². (3.12) Here the function 𝑓₁₁(𝑡) 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 𝑡. In this example the best approximation is obtained when the noise level 𝛿=10⁻⁶ which is comparable with the results given in [6, 9]. The reconstruction of the function 𝑓₁₁(𝑡) 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 𝑚_𝑛𝛿 is bounded by the constant 54 when the noise level 𝛿=10⁻⁶ and 𝜅=0.3.

Figure 11 

Example 3.11: the stability of the approximate solution.

Example 3.12 (see [6, 25]). We have 𝑓₁₂𝐹(𝑡)=𝑡cos(𝑡),0≤𝑡<10,0,elsewhere,₁₂(𝑝)=(𝑝²−1)−𝑒^−10𝑝(−1+𝑝²+10𝑝+10𝑝³)cos(10)(1+𝑝²)²+𝑒^−10𝑝(2𝑝+10+10𝑝²)sin(10)(1+𝑝²)². (3.13) Here we take an increasing function which oscillates as the variable 𝑡 increases over the interval [0,10). A poor approximation is obtained when the noise level 𝛿=10⁻². Figure 12 shows that the exact solution can be approximated very well when the noise level 𝛿=10⁻⁶. 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 

Example 3.12: the stability of the approximate solution.

(a)

(b)

(c)

(d)

Figure 13 

Example 3.13: the stability of the approximate solution.

Example 3.13. We have 𝑓₁₃(𝑡)=𝑒^−𝑡,𝐹₁₃1(𝑝)=.1+𝑝 (3.14) Here the support of 𝑓₁₃(𝑡) is not compact. From the Laplace transform formula one gets 𝐹₁₃(𝑝)=^∞₀𝑒^−𝑡𝑒^−𝑝𝑡𝑑𝑡=^𝑏₀𝑒^{−(1+𝑝)𝑡}𝑑𝑡+^∞_𝑏𝑒^{−(1+𝑝)𝑡}=𝑑𝑡_𝑏0𝑓₁₃(𝑡)𝑒^−𝑝𝑡𝑒𝑑𝑡+^{−(1+𝑝)𝑏}1+𝑝∶=𝐼₁+𝐼₂, (3.15) where 𝛿(𝑏)∶=𝑒^−𝑏. Therefore, 𝐼₂ can be considered as noise of the data 𝐹₁₃(𝑝), that is, 𝐹^𝛿₁₃(𝑝)∶=𝐹₁₃(𝑝)−𝛿(𝑏), (3.16) where 𝛿(𝑏)∶=𝑒^−𝑏. In this example the following parameters are used: 𝑑=2, 𝜅=10⁻¹ for 𝛿=𝑒⁻⁵ and 𝜅=10⁻⁵ for 𝛿=10⁻⁸,10⁻²⁰ and 10⁻³⁰. Table 13 shows that the error decreases as the parameter 𝑏 increases. The approximate solution obtained by the proposed method converges to the function 𝑓₁₃(𝑡) as 𝑏 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 𝐹_𝛿(𝑝_𝑗), 𝑗=0,1,2,…,𝑚, and regularization parameter 𝑎_𝑛𝛿, which are used in computing the approximation 𝑓^𝛿_𝑚(𝑡) (see (2.69)) of the unknown function 𝑓(𝑡), are obtained automatically. Our numerical experiments show that the computation time (CPU time) for approximating the function 𝑓(𝑡) 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 𝑓 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.