Inversion of the Laplace transform from the real axis using an adaptive iterative method

In this paper 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 $f(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.


Introduction
Consider the Laplace transform : where We assume in (2) that f has compact support.This is not a restriction practically.Indeed, if lim t→∞ 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 ∈ L 1 (0, ∞), then Therefore, the contribution of the "tail" f b (t) of f , f b (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 (2) that f ∈ L 2 [0, ∞).One may also assume that f ∈ L 1 [0, ∞), or that |f (t)| ≤ c 1 e c 2 t , where c 1 , c 2 are positive constants.If the last assumption holds, then one may define the function g(t) := f (t)e −(c 2 +1)t .Then g(t) ∈ L 1 [0, ∞), and its Laplace transform G(p) = F (p + c 2 + 1) is known on the interval [c 2 + 1, c 2 + 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 L : X 0,b → L 2 [0, ∞) is compact.Therefore, the inversion of the Laplace transform ( 1) is an ill-posed problem (see [17], [20]).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 equation (1) stably: variational regularization, quasisolutions, iterative regularization (see e.g, [13], [17], [20], [21]).In this paper we propose an adaptive iterative method based on the Dynamical Systems Method (DSM) developed in [20], [21].Some methods have been developed earlier for the inversion of the Laplace transform (see [2], [5], [8], [12]).In many papers the data F (p) are assumed exact and given on the complex axis.In [16] 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 [5].There are several types of the Laplace inversion method compared in [5].The inversion formula for the Laplace transform is well known: 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 [6] and [15].Moreover, the ill-posedness of the Laplace transform inversion is not discussed in all the methods compared in [5].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 [1], [4], [7], [10], [16], [18], [19], [23], [24]).In [1] and [19] 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 [4] 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 [4] that if one uses some basis functions in X 0,b , the problem becomes extremely ill-conditioned if the number m of the basis functions exceeds 20.
In [10] a reproducing kernel method is used for the inversion of the Laplace transform.In the numerical experiments in [10] 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 [10], so this method may be not efficient in practice.A detailed description of the multiple precision technique can be found in [9] and [11].Moreover, the Laplace transform inversion with perturbed data is not discussed in [10].In [24] the authors develop an inversion formula, based on the eigenfunction expansion for the Laplace transform.The difficulties with this method are: a) the inversion formula is not applicable when the data are noisy, b) even for exact data the inversion formula is not suitable for numerical implementation.The Laplace transform as an operator from C 0k into L 2 , where [7].The finite difference method is used in [7] 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 [7].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 [7].
The novel points in our paper are: 1) the representation of the approximation solution (73) of the function f (t) which depends only on the kernel of the Laplace transform, 2) the adaptive iterative scheme (76) and adaptive stopping rule (87), which generate the regularization parameter, the discrete data F δ (p) and the number of terms in (73), needed for obtaining an approximation of the unknown function f (t).
We study the inversion problem using the pair of spaces (X 0,b , L 2 [0, d]), where X 0,b is defined in (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, e.g., with the double precision results given in paper [10].
The smoothness of the kernel allows one to use the compound Simpson's rule in approximating the Laplace transform.Our approach yields a representation (73) of the approximate inversion of the Laplace transform.The number of terms in approximation (73) and the regularization parameter are generated automatically by the proposed adaptive iterative method.Our iterative method is based on the iterative method proposed in [14].The adaptive stopping rule we propose here is based on the discrepancytype principle, established in [22].This stopping rule yields convergence of the approximation (73) to f (t) 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.

Description of the method
Let f ∈ X 0,b .Then equation (1) can be written as: Let us assume that the data F (p), the Laplace transform of f , are known only for 0 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 (5).Let be the inner product and norm in R m+1 , respectively, where w (m) j are the weights of the compound Simpson's rule (see [6, p.58]), i.e., where m is an even number.Then where and It follows from ( 5) and ( 10) that and where Lemma 2.1.Let w (m) j be defined in (8).
for any even number m.
Proof.From definition (8) one gets Lemma 2.1 is proved.
Lemma 2.2.The matrix Q (m) , defined in (14), is positive semidefinite and self-adjoint in R m+1 with respect to the inner product (7). and w (m) j are defined in (8).
We have Thus, Q (m) is self-adjoint with respect to inner product (7).We have where •, • X 0,b is defined in (11).This shows that H m is a Gram matrix.Therefore, This implies Thus, Q (m) is a positive semidefinite and self-adjoint matrix with respect to the inner product (7).
Lemma 2.3.Let T (m) be defined in (12).Then T (m) is self-adjoint and positive semidefinite operator in X 0,b with respect to inner product (11).
Proof.From definition (12) and inner product (11) we get Thus, T (m) is a self-adjoint operator with respect to inner product (11).Let us prove that T (m) is positive semidefinite.Using ( 12), ( 8), ( 7) and (11), one gets Let us approximate the unknown f (t) as follows: where p j are defined in (6), T a,m is defined in (34), and c (m) j are constants obtained by solving the linear algebraic system: where Q (m) is defined in (13), To prove the convergence of the approximate solution f (t), we use the following estimates, which are proved in [21], so their proofs are omitted.
Lemma 2.4.Let T (m) and Q (m) be defined in (12) and (13), respectively.Then, for a > 0, the following estimates hold: where I is the identity operator and a = const > 0.
Let us formulate an iterative method for obtaining the approximation solution of f (t) with the exact data F (p). Consider the following iterative scheme where L * is the adjoint of the operator L, i.e., k(p, t, z) is defined in (26), Lemma 2.5.Let T a be defined in (38), Lf = F , and f ⊥ N (L), where Proof.Since f ⊥ N (L), it follows from the spectral theorem that lim a→0 where E s is the resolution of the identity corresponding to L * L, and P is the orthogonal projector onto N (L).Lemma 2.5 is proved.
Theorem 2.6.Let Lf = F , and u n be defined in (35) Then Proof.By induction we get where T a is defined in (38), and Using the identities and we get Therefore, To prove relation (41) the following lemma is needed: Lemma 2.7.Let g(x) be a continuous function on (0, ∞), c > 0 and q ∈ (0, 1) be constants.If then Proof.Let where ω Take ǫ > 0 arbitrarily small.For sufficiently large fixed l(ǫ) one can choose n(ǫ) > l(ǫ), such that . This is possible because of (49).One has and if n(ǫ) is sufficiently large.Here we have used the relation Since ǫ > 0 is arbitrarily small, relation (50) follows.Lemma 2.7 is proved.
Lemma 2.8.Let T and T (m) be defined in (37) and (12), respectively.Then Proof.From definitions (37) and ( 12) we get where the following upper bound for the error of the compound Simpson's rule was used (see [6, p.58]): for f ∈ C (4) [x 0 , x 2l ], x 0 < x 2l , where and This implies Then where T and T (m) are defined in (37) and (12), respectively.
Lemma 2.9 leads to an adaptive iterative scheme: where q ∈ (0, 1), a n are defined in (39), T a,m is defined in (34), A m L is defined in (5), and p j are defined in (6).In the iterative scheme (61) we have used the finitedimensional operator T (m) approximating the operator T .Convergence of the iterative scheme (61) to the solution f of the equation Lf = F is established in the following lemma: Lemma 2.10.Let Lf = F and u n,mn be defined in (61).If m n are chosen by the rule where ⌈[x]⌉ is the smallest even number not less than x, then Proof.Consider the estimate where I 1 (n) := f − u n and I 2 (n) := u n − u n,mn .By Theorem 2.6, we get I 1 (n) → 0 as n → ∞.Let us prove that lim n→∞ I 2 (n) = 0. Let U n := u n − u n,mn .Then, from definitions (35) and (61), we get By induction we obtain where ω j are defined in (43).Using the identities one gets This together with the rule (63), estimate (32) and Lemma 2.8 yield Applying Lemma 2.5 and Lemma 2.7 with g(a) = a T −1 a f , we obtain lim n→∞ U n = 0. Lemma 2.10 is proved.

Noisy data
When the data F (p) are noisy, the approximate solution ( 27) is written as where the coefficients c (m,δ) j are obtained by solving the following linear algebraic system: are defined in (8), and p j are defined in (6).To get the approximation solution of the function f (t) with the noisy data F δ (p), we consider the following iterative scheme: where T a,m is defined in (34), a n are defined in (39), q ∈ (0, 1), F (m) δ is defined in (75), and m n are chosen by the rule (63).Let us assume that where δ j are random quantities generated from some statistical distributions, e.g., the uniform distribution on the interval [−δ, δ], and δ is the noise level of the data F (p).It follows from assumption (77), definition (8), Lemma 2.1 and the inner product (7) that Lemma 2.11.Let u n,mn and u δ n,mn be defined in (61) and (76), respectively.Then where a n are defined in (39).
Theorem 2.12.Suppose that conditions of Lemma 2.10 hold, and n δ satisfies the following conditions: Then Proof.Consider the estimate: This together with Lemma 2.11 yield Applying relations (83) in estimate (86), one gets relation (84).Theorem 2.12 is proved.
In the following subsection we propose a stopping rule which implies relations (83).

Stopping rule
In this subsection a stopping rule which yields relations (83) in Theorem 2.12 is given.We propose the stopping rule where • W m is defined in (7), (m) j and p j are defined in ( 8) and ( 6), respectively, and c (m,δ) j are obtained by solving linear algebraic system (74).
We observe that Thus, the sequence (88) can be written in the following form where • W m is defined in (7), and c (m,δ) solves the linear algebraic systems (74).It follows from estimates (78), ( 30) and (31) that This together with (91) yield or Lemma 2.13.The sequence (91) satisfies the following estimate: where a n are defined in (39). and Then estimate (94) can be rewritten as where the relation a n = qa n−1 was used.Let us prove estimate (95) by induction.For n = 0 we get Suppose estimate (95) is true for 0 ≤ n ≤ k.Then where the relation a k+1 = qa k was used.Lemma 2.13 is proved.
where G n,mn are defined in (91).Then there exist a unique integer n δ , satisfying the stopping rule (87) with C > √ d.
Proof.From Lemma 2.13 we get the estimate where a n are defined in (39).Therefore, where the relation lim n→∞ a n = 0 was used.This together with condition (101) yield the existence of the integer n δ .The uniqueness of the integer n δ follows from its definition.Lemma 2.14 is proved.
Lemma 2.15.Suppose conditions of Lemma 2.14 hold and n δ is chosen by the rule (87).Then Proof.From the stopping rule (87) and estimate (102) we get where C > √ d, ε ∈ (0, 1).This implies so, for ε ∈ (0, 1), and Proof.From the stopping rule (87) with the sequence G n defined in (91) one gets where c (m,δ) is obtained by solving linear algebraic system (74).This implies Thus, lim If F (m) = 0, then there exists a λ where is the resolution of the identity corresponding to the operator For a fixed number a > 0 we obtain λ 0 is a continuous operator, and Therefore, for the fixed number a > 0 we get for all sufficiently small δ > 0, where c 2 is a constant which does not depend on δ.Suppose lim δ→0 a n δ = 0. Then there exists a subsequence δ j → 0 as j → ∞, such that and where the rule (63) was used to obtain the parameters m n δ j .This together with (112) and (115) yield This contradicts relation (111).Thus, lim δ→0 a n δ = lim δ→0 a 0 q n δ = 0, i.e., lim δ→0 n δ = ∞.Lemma 2.16 is proved.
It follows from Lemma 2.15 and Lemma 2.16 that the stopping rule (87) yields the relations (83).We have proved the following theorem: Theorem 2.17.Suppose all the assumptions of Theorem 2.12 hold, m n are chosen by the rule (63), n δ is chosen by the rule (87) and G 1,m 1 > Cδ, where G n,mn are defined in (91), then (119)

The algorithm
Let us formulate the algorithm for obtaining the approximate solution f δ m : (1) The data F δ (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, a 0 > 0, q ∈ (0, 1), ε ∈ (0, 1), C > √ d, and set u δ 0,m 0 = 0, G 0 = 0, n = 1; (3) iterate, starting with n = 1, and stop when condition (126) ( see below) holds, (a) a n = a 0 q n , (b) choose m n by the rule (63), (c) construct the vector (d) construct the matrices H mn and D mn : where w (m) j are defined in ( 8), (e) solve the following linear algebraic systems: where (c (mn ,δ) of the approximate solution u δ n,mn (t) defined in (73) by the iterative formula: where Stop when for the first time the inequality holds, and get the approximation f δ (t) = u δ n δ ,mn δ (t) of the function f (t) by formula (124).
3 Numerical experiments 3.1 The parameters κ, a 0 , d From definition (39) and the rule (63) we conclude that m n → ∞ as a n → 0. Therefore, one needs to control the value of the parameter m n so that it will not grow too fast as a n decreases.The role of the parameter κ in (63) is to control the value of the parameter m n so that the value of the parameter m n will not be too large.Since for sufficiently small noise level δ, namely δ ∈ (10 −16 , 10 −6 ], the regularization parameter a n δ , obtained by the stopping rule (87), 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 a n = a 0 δ αn , where a 0 ∈ [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 60), ( 78), ( 79), ( 82), ( 94), (95), and (102) are replaced with the constant d − d 1 .If b = 1, i.e., 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: 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 δ (p j ), j = 0, 1, 2, . . ., m, over the interval [0, d] which are needed to get the approximation of the function f (t).

Experiments
To test the proposed method we consider some examples proposed in [1], [2], [3], [4], [5], [8], [10], [16], [18] and [24].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 (77) are obtained from the uniform probability density function over the interval [−δ, δ].In examples 1-12 we choose the value of the parameters as follows: a n = 0.1q n , 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 m n are not very large, namely m n ≤ 300.Therefore, the computation time for solving linear algebraic system (123) 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 (87) the following parameters are used: C = √ d + 0.01, ε = 0.99.In example 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: , t j = 0.01+0.1(j−1),j = 1, . . ., 100, (127) where f (t) is the exact solution and f δ mn δ (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 m n δ and a n δ generated by the proposed method are given in each experiment (see Tables 1-12).All the calculations are done in double precision generated by MATLAB.
• Example 1. (see [10]) , 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 [10].
The proposed method is stable with respect to the noise δ as shown in Table 1.
• Example 2. (see [4], [10] ) elsewhere, , p > 0.   The reconstruction of the function f 2 (t) is plotted in Figure 2. In [10] a high accuracy result is given by means of the multiple precision.But, as reported in [10], 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 [4].The result is comparable with the double precision results given in [10].For δ = 10 −6 and κ = 0.3 the value of the parameter m n δ is bounded by the constant 54.
• Example 3. (see [1], [4], [5], [18], [24]) 10e −(p+1)10 p + 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 [4].The mean absolute error M AE 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 [24].The value of the parameter m n δ is bounded by the constant 30 when the noise level δ = 10 −6 and κ = 0.3.
• Example 4. (see [4], [10]) 1 − e −0.5t , 0 ≤ t < 10, 0, elsewhere., p > 0.  4 gives the results of the stability of the proposed method with respect to the noise level δ.Moreover, the reconstruction of the function f 4 (t) obtained by the proposed method is better than the reconstruction of f 4 (t) shown in [4], and is comparable with the double precision reconstruction obtained in [10].In this example when δ = 10 −6 and κ = 0.3 the value of the parameter m n δ is bounded by the constant 109 as shown in Table 4.
• Example 5. (see [2], [4], [8])   This is an example of the damped sine function.In [2] and [8] 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 δ (p j ), j = 0, 1, 2, . . ., m, and get a satisfactory result which is comparable with the results given in [2] and [8] when the level noise δ = 10 −6 .The reconstruction of the exact solution f 5 (t) obtained by our method is better than this of the method given in [4].Moreover, our method yields stable solution with respect to the noise level δ as shown in Figure 5 and Table 5 show.In this example when κ = 0.3 the value of the parameter m n δ is bounded by 54 for the noise level δ = 10 −6 (see Table 5).
, p > 0.  Example 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 can not 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 [10].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.
• Example 7. (see [10]) elsewhere,   When the noise level δ = 10 −4 and δ = 10 −6 , we get numerical results which are comparable with the double precision results given in [10].Figure 7 and Table 7 show the stability of the proposed method for decreasing δ.
• Example 8. (see [3], [4])  The results of this example are similar to the results of Example 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 are relatively small.We get better reconstruction results than the results shown in [4].Our results are comparable with the results given in [3].• Example 9. (see [18]) As in Example 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   9.When the double precision is used, we get comparable results with the results shown in [18].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.
• Example 11. (see [5], [16])  Here the function f 11 (t) represents the class of periodic functions.It is mentioned in [16] 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 [5] and [16].The reconstruction of the function f 11 (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 m n δ is bounded by the constant 54 when the noise level δ = 10 −6 and κ = 0.3.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 [3] and [5].The stability of our method with respect to the noise level is shown in Table 12.Here the support of f 13 (t) is not compact.From the Laplace transform formula one gets  13 shows that the error decreases as the parameter b increases.The approximate solution obtained by the proposed method converges to the function f 13 (t) as b increases (see Figure 13).

Conclusion
We have tested the proposed algorithm on the wide class of examples considered in the literature.Using the rule (63) and the stopping rule (87), the number of terms in representation (73), the discrete data F δ (p j ), j = 0, 1, 2, . . ., m, and regularization parameter a n δ , which are used in computing the approximation f δ m (t) (see (73)) 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 seconds, and the proposed iterative scheme and the proposed adaptive stopping rule yield stable solution with respect to the noise level δ.The proposed method also works for f without compact support as shown in Example 13.Moreover, in the proposed method we only use a simple representation (73) which is based on the kernel of the Laplace transform integral, so it can be easily implemented numerically.

Figure 1 :
Figure 1: Example 1: the stability of the approximate solution

Figure 2 :
Figure 2: Example 2: the stability of the approximate solution

Figure 3 :
Figure 3: Example 3: the stability of the approximate solution

Figure 4 :
Figure 4: Example 4: the stability of the approximate solution

Figure 5 :
Figure 5: Example 5: the stability of the approximate solution

6 Figure 6 :
Figure 6: Example 6: the stability of the approximate solution

Figure 7 :
Figure 7: Example 7: the stability of the approximate solution

Figure 8 :
Figure 8: Example 8: the stability of the approximate solution

Figure 9 :
Figure 9: Example 9: the stability of the approximate solution

Figure 10 :
Figure 10: Example 10: the stability of the approximate solution

Figure 11 :
Figure 11: Example 11: the stability of the approximate solution

Figure 12 :
Figure 12: Example 12: the stability of the approximate solution

Figure 13 :
Figure 13: Example 13: the stability of the approximate solution