Mathematical Models and Methods in Applied Sciences c❢World Scientific Publishing Company arXiv:math-ph/9911038v1 26 Nov 1999 CONTINUOUS ANALOG OF THE GAUSS-NEWTON METHOD RUBEN G. AIRAPETYAN ∗, ALEXANDER G. RAMM †, and ALEXANDRA B. SMIRNOVA Department of Mathematics, Kansas State University, Manhattan, Kansas 66506-2602, U.S.A. A Continuous Analog of discrete Gauss-Newton Method (CAGNM) for numerical solution of nonlinear problems is suggested. In order to avoid the ill-posed inversion of the Fréchet derivative operator some regularization function is introduced. For the CAGNM a convergence theorem is proved. The proposed method is illustrated by a numerical example in which a nonlinear inverse problem of gravimetry is considered. Based on the results of the numerical experiments practical recommendations for the choice of the regularization function are given. Keywords: Continuous Gauss-Newton method; iterative scheme; Fréchet derivative; regularization. 1. Introduction Let H1 and H2 be real Hilbert spaces and ϕ : H1 → H2 a nonlinear operator. Let us consider the equation: ϕ(x) = 0. (1.1) We assume that the following condition on ϕ holds. Condition A: Problem (1.1) has a solution x̂, not necessarily unique. In the well-known Newton’s method (5 ) one constructs a sequence {xn } for n = 0, 1, . . . which converges to a solution (in general non unique) of equation (1.1). The first term x0 is an initial approximation point and the other terms are constructed by means of the following iterative process: xk+1 = xk − ϕ′ (xk )−1 ϕ(xk ), (1.2) where ϕ′ (x) is the Fréchet derivative of the operator ϕ. Recall that ϕ′ (x) is a linear operator from H1 to H2 . The usual necessary condition for the realization of the Newton method is the bounded invertibility of ϕ′ (xk ), that is, the existence ∗E-mail: airapet@math.ksu.edu ramm@math.ksu.edu ‡ E-mail: smirn@math.ksu.edu † E-mail: ‡ of a bounded linear operator [ϕ′ (xk )]−1 for all k. Actually in order to provide the convergence of Newton iterations one needs bounded invertibility of ϕ′ in a ball B(x̂, R) := {x : x ∈ H1 , ||x − x̂|| ≤ R}. However this condition does not hold in many important applications. In order to avoid this restriction several modifications of the Newton method have been developed. In this paper we consider the GaussNewton procedure for equation (1.1) (see e.g. (5 )): xk+1 = xk − [ϕ′∗ (xk )ϕ′ (xk )]−1 ϕ′∗ (xk )ϕ(xk ), x0 = x0 . (1.3) If the operator ϕ′∗ (x)ϕ′ (x) is not boundedly invertible one needs some regularization procedure. In order to construct such a procedure one can introduce a sequence of positive numbers αk , αk → 0, and replace iterative method (1.3) by the following one (5,2 ): xk+1 = xk − [ϕ′∗ (xk )ϕ′ (xk ) + αk I]−1 [ϕ′∗ (xk )ϕ(xk ) + αk (xk − x0 )], (1.4) where I is the identity operator. The methods constructed above can be also considered as discrete analogs of some continuous methods (called sometimes continuation methods). In (3 ) the following Cauchy problem has been considered as a continuous analog of (1.2): ẋ(t) = −ϕ′ (x(t))ϕ(x(t)), x(0) = x0 , ẋ(t) := dx . dt (1.5) A solution to problem (1.1) can be obtained as limit of the function x(t) for t → ∞. If one solves this Cauchy problem by means of Euler’s method with a stepsize τ = 1 one gets (1.2) with xk = x(k). Continuous analogs of iterative methods have several advantages over the discrete ones. Convergence theorems for continuous methods usually can be obtained easier. If a convergence theorem is proved for a continuous method, that is, for the Cauchy problem for a differential equation, for instance (1.5), one can construct various finite difference schemes for the solution of this Cauchy problem. These difference schemes give discrete methods for the solution of equation (1.1). For instance the methods of Euler and Runge-Kutta can be used. More detailed information about the applications and modifications of continuous Newton methods can be found in (3,8,1 ). The aim of this paper is to construct a continuous analog of iterative scheme (1.4), to prove a convergence theorem for this continuous analog of (1.4), and to test the method numerically by applying it to a practically interesting nonlinear inverse problem of gravimetry. The paper is organized as follows. In section 2 a continuous analog of method (1.4) is described and a convergence theorem for this method is formulated. In section 3 this convergence theorem is proved. In section 4 an inverse gravimetry problem is considered and the proposed method is numerically tested. In our numerical experiments comparison of different regularization functions is done. Based on the results of the numerical experiments some recommendations are given for the choice of the regularization function. 2. Continuous Gauss-Newton Method and Convergence Theorem In order to describe convergence rates we introduce the following Definition 2.1. A positive function α(t) ∈ C 1 [0, ∞) is said to be a convergence rate function if α(t) decreases monotonically to zero as t → ∞, α(t) ∈ C 1 [0, ∞) and ln α(t) is concave, that is, α̇(t)/α(t) is monotonically increasing. Remark 2.2. The number α(0) can be chosen sufficiently large and simultaneously the number |α̇(0)/α(0)| can be sufficiently small. Here and below the over dot denotes the derivative with respect to time ẋ := dx/dt. For example, one can choose α(t) = b/(t + a), where a and b are positive constants such that a and b/a are sufficiently large. A continuous analog of iterative process (1.4) is the following Cauchy problem: ẋ(t) = −[ϕ′∗ (x(t))ϕ′ (x(t)) + α(t)I]−1 [ϕ′∗ (x(t))ϕ(x(t)) + α(t)(x(t) − x0 )], (2.1) x(0) = x0 . Denote by Ran(L) the range of the linear operator L. The convergence of the continuous analog of Gauss-Newton method (CAGNM) is established by the following theorem, in which (and throughout this paper) the norms ||ϕ′ (x)|| and ||ϕ′′ (x)|| are the norms of linear and bilinear operators from H1 to H2 and from H1 × H1 to H2 respectively. Theorem 2.3. Let α(t) be a convergence rate function. Assume that there exists a positive number R for which Condition A and the following conditions hold: (i) The Fréchet derivatives ϕ′ (x) and ϕ′′ (x) exist in the ball B(x̂, R) and satisfy the following inequalities: ||ϕ′ (x)|| ≤ N1 , where α(0) N1 N2 ||ϕ′′ (x)|| ≤ N2 ∀x ∈ B(x̂, R), (2.2)   α̇(0) ≤ R. 1 − 2N1 N2 ||v|| + α(0) (ii) x0 ∈ B(x̂, R) ∩ [x̂ + Ran(ϕ′∗ (x̂)ϕ′ (x̂))]. (2.3) (iii) For some v, such that x̂ − x0 = ϕ′∗ (x̂)ϕ′ (x̂)v, the following inequalities hold: 1 − 2N1 N2 ||v|| +  α̇(0) 1 − 2N1 N2 ||v|| + α(0) α̇(0) > 0, α(0) 2 − 2N1 N2 ||v|| > 0. (2.4) (2.5) Then the following conclusions hold: (i) The solution x = x(t) of problem (2.1) exists, and x(t) ∈ B(x̂, R) for t ∈ [0, ∞), (ii) ||x(t) − x̂|| = O(α(t)) for t → ∞. Remark 2.4. Condition (ii) in Theorem 2.3 gives some restriction on the choice of an initial approximation point. It is not easy to verify this condition algorithmically. However some kind of this condition is necessary if one works with equation (1.1) with the operator ϕ′∗ (x)ϕ′ (x), which is not boundedly invertible. If the operator ϕ′∗ ϕ′ is injective but is not boundedly invertible, then the image of the linear selfadjoint operator ϕ′∗ ϕ′ is dense in B(x̂, R) and consequently the set of the suitable initial approximation points satisfying condition (ii) is also dense in B(x̂, R). As our numerical results show (see section 4) the proposed method is practically efficient. 3. Proof of Theorem 2.3 The main part of the proof is to show that the solution to problem (2.1) does not leave the ball B(x̂, R) (Lemma 3.3). In order to prove it, let us assume that there exists such a point t1 ∈ [0, ∞) that x(t) intersects the boundary of B(x̂, R) for the first time at t = t1 . Hence x(t) belongs to the interior of the B(x̂, R) for t ∈ [0, t1 ) and ||x(t1 ) − x̂|| = R. Let us introduce an auxiliary function w(t) := ||x(t) − x̂||/α(t). (3.1) First, in Lemma 3.1, we derive a nonlinear differential inequality for w(t). From this differential inequality we get the estimate which shows that for all t ∈ [0, t1 ] the points of the integral curve of problem (2.1) belong to the interior of the ball B(x̂, R). This contradiction proves that the integral curve of the solution does not leave the above ball, and consequently problem (2.1) has the global solution for t ∈ [0, ∞). Also we show the boundedness of the w(t) and this implies, by formula (3.1), strong convergence of x(t) to x̂ for t → ∞. Lemma 3.1. If the assumptions of Theorem 2.3 hold then the differential inequality dw ≤ C1 w2 − C2 w + C3 , (3.2) dt is valid for t ∈ [0, t1 ], where C1 = N1 N2 , 2 C2 = 1 − 2N1 N2 ||v|| + α̇(0) , α(0) C3 = ||v||. (3.3) Proof. The Gâteaux derivative ϕ′′ (x, ξ1 , ξ2 ) is a bilinear operator such that ϕ′ (x+ξ1 )ξ2 −ϕ′ (x)ξ2 := ϕ′′ (x, ξ1 , ξ2 )+ηξ2 , and ||η||·||ξ1 ||−1 → 0 for ξ1 → 0, ||ξ1 || > 0. Let us define operators K, G : B(x̂, R) × H1 × H1 → H2 by the formulas: K(x, ξ1 , ξ2 ) = Z1 Z1 0 and G(x, ξ1 , ξ2 ) = ϕ′′ (x + stξ1 , ξ1 , ξ2 )tdtds (3.4) ϕ′′ (x + tξ1 , ξ1 , ξ2 )dt. (3.5) 0 Z1 0 Then from (2.2) we get |K(x, ξ1 , ξ2 )| ≤ N2 ||ξ1 || · ||ξ2 ||, 2 |G(x, ξ1 , ξ2 )| ≤ N2 ||ξ1 || · ||ξ2 ||. (3.6) The following formulas will be used: ϕ(x̂) − ϕ(x) = ϕ′ (x)(x̂ − x) + K(x, x̂ − x, x̂ − x), (3.7) (ϕ′ (x̂) − ϕ′ (x))ξ = G(x, x̂ − x, ξ), (3.8) and where K and G are defined by (3.4) and (3.5) respectively. Let us derive formulas (3.7) and (3.8). One has ϕ(x̂) − ϕ(x) = Z1 d ϕ(x + t(x̂ − x))dt = dt Z1 ϕ′ (x + t(x̂ − x))(x̂ − x)dt 0 0 = ϕ′ (x)(x̂ − x) + Z1 [ϕ′ (x + t(x̂ − x)) − ϕ′ (x)](x̂ − x)dt 0 ′ = ϕ (x)(x̂ − x) + Z1 0 tdt Z1 ds[ϕ′′ (x + st(x̂ − x))](x̂ − x)(x̂ − x) 0 and ′ ′ [ϕ (x̂) − ϕ (x)]ξ = Z1 0 d ′ ϕ (x + t(x̂ − x))ξdt = dt Z1 ϕ′′ (x + t(x̂ − x), x̂ − x, ξ)dt. 0 Since x̂ solves (1.1), one can rewrite equation (2.1) as dx = −[ϕ′∗ (x)ϕ′ (x) + αI]−1 [ϕ′∗ (x)(ϕ(x) − ϕ(x̂)) + α(x − x0 )]. dt From the condition (iii) of Theorem 2.3 and from (3.7) we get dx = −[ϕ′∗ (x)ϕ′ (x) + αI]−1 [−ϕ′∗ (x)(ϕ′ (x)(x̂ − x) + K(x, x̂ − x, x̂ − x))+ dt α(x − x̂) + αϕ′∗ (x̂)ϕ′ (x̂)v], and therefore dx = −(x − x̂) − [ϕ′∗ (x)ϕ′ (x) + αI]−1 [−ϕ′∗ (x)K(x, x̂ − x, x̂ − x) + αϕ′∗ (x)ϕ′ (x)v+ dt α(ϕ′∗ (x̂)ϕ′ (x̂) − ϕ′∗ (x)ϕ′ (x))v]. Since x̂ does not depend on t, it follows from (3.8) that dx d(x(t) − x̂) = = −(x − x̂) − [ϕ′∗ (x)ϕ′ (x) + αI]−1 [−ϕ′∗ (x)K(x, x̂ − x, x̂ − x) dt dt +αϕ′∗ (x)ϕ′ (x)v + α(ϕ′∗ (x̂) − ϕ′∗ (x))ϕ′ (x̂)v + αϕ′∗ (x)G(x, x̂ − x, v)]. Now let us derive an inequality for d dt ||x − x̂||2 . One has d ||x − x̂||2 = −2||x − x̂||2 + 2([ϕ′∗ (x)ϕ′ (x) + αI]−1 [ϕ′∗ (x)K(x, x̂ − x, x̂− x)], x − x̂)− dt 2α([ϕ′∗ (x)ϕ′ (x) + αI]−1 ϕ′∗ (x)ϕ′ (x)v, x − x̂) + 2α(ϕ′ (x̂)v, G(x, x̂ − x, [ϕ′∗ (x)ϕ′ (x) +αI]−1 (x − x̂))) + 2α([ϕ′∗ (x)ϕ′ (x) + αI]−1 ϕ′∗ (x)G(x, x̂ − x, v), x − x̂). Since the operator ϕ′∗ (x)ϕ′ (x) is selfadjoint and nonnegative we have the following spectral representation: Z ∞ λ ′∗ ′ −1 ′∗ ′ dEλ , [ϕ (x)ϕ (x) + αI] ϕ (x)ϕ (x) = λ + α 0 where Eλ is the resolution of the identity of the selfadjoint operator ϕ′∗ (x)ϕ′ (x). Since 0 ≤ λ/(λ + α) ≤ 1 for α > 0 and λ ≥ 0, it follows that ||[ϕ′∗ (x)ϕ′ (x) + α(t)I]−1 ϕ′∗ (x)ϕ′ (x)|| ≤ 1. (3.9) Also one has the following estimate: ||[ϕ′∗ (x)ϕ′ (x) + α(t)I]−1 || ≤ 1/α(t). (3.10) From (3.9),(3.10) and (3.6) one gets the following differential inequality for A(t) := ||x(t) − x̂|| : N1 N2 2 Ȧ ≤ −A + A + α||v|| + 2N1 N2 ||v||A. 2α In order to finish the proof of the Lemma 3.1, we derive from the last inequality the inequality for w(t) by taking into account that α̇(t)/α(t) is monotonically increasing function. ✷ The following lemma is a simple corollary of the more general results established in (6 ). Lemma 3.2. Let f (t, u) be a continuous function on [0, T ] × (−∞, +∞) such that the Cauchy problem u̇ = f (t, u(t)), u(0) = u0 (3.11) is uniquely solvable on [0, T ] and v(t) be a differentiable function defined on [0, T ] and satisfies the conditions v̇ ≤ f (t, v(t)), t ∈ [0, T ], v(0) = v0 . (3.12) If v0 ≤ u0 then v(t) ≤ u(t) for t ∈ [0, T ]. p It follows from inequality (2.5) that c = C22 − 4C1 C3 > 0 for constants C1 , C2 , C3 defined in (3.3). Let u1 and u2 be correspondingly the smaller and the larger roots of the equation C1 u2 − C2 u + C3 = 0. For u0 satisfying the inequality u2 < u0 < C2 /(2C1 ) the solution of the Cauchy problem u̇ = C1 u2 − C2 u + C3 , u(0) = u0 . (3.13) is given by the formula: u2 − u0 ct u − u2 = e . u − u1 u0 − u1 Let us show that u(t) is defined for all t ∈ [0, ∞). Indeed, u0 ∈ (u1 , u2 ), so for sufficiently small t one has u = u1 + u2 − u1 u2 −u0 ct u0 −u1 e + 1 . (3.14) Thus u1 < u(t) < u(0) < u2 . This means that u(t) does not leave the interval (u1 , u2 ) for all t ∈ [0, ∞) and u(t) is well defined for all t ∈ [0, ∞). From condition (i) of Theorem 2.3 one obtains the following estimate: w(0) = 1 − 2N1 N2 ||v|| + ||x0 − x̂|| < α(0) N1 N2 α̇(0) α(0) = C2 . 2C1 Therefore from Lemmas 3.1 and 3.2 it follows that ||x(t1 ) − x̂|| ≤ u1 + α(t1 ) Thus ||x(t1 ) − x̂|| < u2 − u1 u2 −u0 ct1 + u0 −u1 e 1 < u0 < C2 . 2C1 (3.15) C2 C2 α(t1 ) < α(0) ≤ R. 2C1 2C1 This contradicts the assumption ||x(t1 )− x̂|| = R. So the following lemma is proved. Lemma 3.3. If the assumptions of Theorem 2.3 hold and for an arbitrary positive T the solution of the problem (2.1) exists on the interval [0, T ], then the integral curve of the solution of (2.1) lies in the interior of the ball B(x̂, R) for all t from the interval [0, T ]. Now let us show that there exists the unique solution of (2.1) on [0, ∞)provided that x0 satisfies conditions ii) and iii) of Theorem 2.3. The Cauchy problem (2.1) is equivalent to the integral equation x(t) = x0 + Zt F (s, x(s))ds, (3.16) 0 where F (s, x(s)) := −[ϕ′∗ (x(s))ϕ′ (x(s)) + α(s)I]−1 [ϕ′∗ (x(s))ϕ(x(s)) + α(s)(x(s) − x0 )]. Let us fix an arbitrary large positive number T and use the successive approximation method to solve equation (3.16) on [0, T ]: xn+1 (t) = x0 + Zt F (s, xn (s))ds, x0 (t) = x0 , t ∈ [0, T ]. (3.17) 0 Since ϕ′ (x) and ϕ′′ (x) are assumed to be bounded in B(x̂, R), see (2.2), and α(t) is positive on [0, T ], for every t ∈ [0, T ] the function F (t, x) has bounded Fréchet derivative with respect to x in B(x̂, R). So one has: ||F (t, x1 ) − F (t, x2 )|| ≤ K(T )||x1 − x2 || for all t ∈ [0, T ] and x1 , x2 belong to B(x̂, R). Thus, one easily gets the estimate ||xn+1 (s) − xn (s)|| ≤ ||F (x0 )||K n (T ) Tn n! valid on the maximal subinterval [0, T1 ] = {t : t ∈ [0, T ] and x(t) ∈ B(x̂, R)}. Therefore iterative process (3.17) converges uniformly and determines the unique solution of equation (3.16) on [0, T1 ]. If T1 < T , it follows from the maximality of the subinterval [0, T1 ] that x(T1 ) is a boundary point of B(x̂, R). But this contradicts to Lemma 3.3. So the solution of the problem (2.1) exists and belongs to the interior of the ball B(x̂, R) on every interval [0, T ] and consequently on [0, ∞). To finish the proof of Theorem 2.3 it is sufficient to note that equation (3.15) implies estimate C2 α(t) ||x(t) − x̂|| ≤ 2C1 for all t ∈ [0, ∞). ✷ 4. Numerical Results To test numerically the method described above, we chose the inverse gravimetry problem (7 ). The goal of the numerical test is to illustrate the choice of the regularization function α(t) and to compare two methods of solving the Cauchy problem (2.1): the Euler method, which corresponds to the iterative scheme (1.4), and the Runge-Kutta method. Let the sources of a gravitational field with a constant density ρ be distributed in the domain D = {−l ≤ t ≤ l, −H ≤ z ≤ −H + x(t)}, where x(t) is an interface between two media, l and H are parameters of the domain. The potential V of such a field is given by the double integral: Z Z 1 1 ρ ln p dS = V (t, z) = 2π D (t − s)2 + (z − τ )2 ρ =− 4π Zl ds −H+x(s) Z ln[(t − s)2 + (z − τ )2 ]dτ. −H −l For the z - component of the gravitational field one has ∂V (t, z) ρ − =− ∂z 4π Zl ds −l ρ = 4π Zl ln −H+x(s) Z ∂ ln[t − s)2 + (x − τ )2 ]dτ = ∂τ −H (t − s)2 + (z + H)2 ds. (t − s)2 + (z + H − x(s))2 −l In particular, on the surface z = 0 we obtain the following nonlinear operator equation Zl ρ K(t, s, x(s))ds − y(t) = 0, (4.1) ϕ(x) ≡ 4π −l where K(t, s, x(s)) = ln (t − s)2 + H 2 . (t − s)2 + (H − x(s))2 (t,0) is given and the interface between two The gravity strength anomaly y(t) = − ∂V∂z media (with and without the sources of a gravimetry field) x(s) is to be determined. Let ϕ act between the pair of Hilbert spaces H1 and H2 . Assume that H1 = H 1 [−l, l] or L2 [−l, l] and H2 = L2 [−l, l]. The Fréchet derivative of this operator is the following one Zl 2(H − x(s))h(s) ′ ϕ (x)h = ds. (4.2) (t − s)2 + (H − x(s))2 −l For any fixed x ∈ {x ∈ L2 [−l, l], x ≤ H − ε, ε > 0} the kernel Kx′ (t, s, x(s)) ≡ 2(H − x(s)) (t − s)2 + (H − x(s))2 is a square integrable function on [−l, l]×[−l, l], therefore ϕ′ (x) in (4.2) is a compact linear operator in L2 [−l, l]. This means that the operators ϕ′ (x) and ϕ′∗ (x)ϕ′ (x) are not boundedly invertible. So, one can not use classical iterative schemes such as Newton or Gauss - Newton in the case of equation (4.1). We solve the Cauchy problem (2.1) for ϕ given by (4.1) with some regularization function α(t). The problem is numerically solved by means of two finite difference methods, namely, Euler’s method xk+1 = xk + τ F (tk , xk ), x0 = x(0) (4.3) and the Runge - Kutta method xk+ 12 = xk + τ F (tk , xk ), 2 xk+1 = xk + τ F (tk+ 21 , xk+ 12 ), (4.4) x0 = x(0). Here F (tk , xk ) ≡ −[ϕ′∗ (xk )ϕ′ (xk ) + α(tk )I]−1 (ϕ′∗ (xk )ϕ(xk ) + α(tk )(xk − x0 )) and an equal grid size τ > 0 defines the node points, tk = kτ, k = 0, 1, .. For the successful realization of Continuous Gauss-Newton method an appropriate regularization function should be chosen. At the beginning of the process values of α(t) should not be very small for the operator ϕ′∗ (x(t))ϕ′ (x(t)) + α(t)I to be stably invertible and at the same time α(t) should tend to zero sufficiently fast to ensure convergence of the function x(t) to the solution of problem (4.1). In the numerical experiments the functions α0 /(β + t)m , α0 e−βt and α0 2−βt were used. The experiments have shown that for all the considered functions the numerical solution is evaluated with appropriate accuracy for sufficiently large range of the parameters m, α0 and β. The following tables illustrate the dependence of the accuracy of the numerical results on parameters for the following data l = 1, H = 2, ρ = 1, x0 = 1. For the numerical tests the function y(t) in (4.1) was chosen as the solution of the direct problem for the model function xmod (t) = (1 − t2 )2 . The integral in (4.1) was calculated by Simpson’s formula with the number of node points equals 201 and with a step size equal to 0.01. In the tables below ∆E and ∆R are the absolute errors, σE and σR are the discrepancies ||ϕ(x(t))||, see (4.1), of the Euler and the Runge-Kutta methods respectively. The first table shows the dependence of the absolute errors and the discrepancies ϕ(x) on the regularization function. Table 1. α(t) α0 (1 + t)−2 α0 (1 + t)−4 α0 (1 + t)−6 α0 (1 + t)−8 α0 (1 + t)−10 α0 2−3.5t α0 e−3.5t N 38 25 43 61 79 127 85 α0 = 0.1, τ = 0.1 ∆E σE 0.26 3.75 · 10−2 0.25 0.25 0.12 0.12 5.30 · 10−2 2.87 · 10−3 2.25 · 10−2 4.90 · 10−4 1.07 · 10−2 1.14 · 10−5 1.08 · 10−2 2.84 · 10−4 ∆R 0.27 0.11 1.75 · 10−2 5.33 · 10−2 2.23 · 10−2 1.08 · 10−2 1.08 · 10−2 σR 4.09 · 10−2 0.11 2.03 · 10−2 3.58 · 10−3 6.22 · 10−4 3.03 · 10−5 3.60 · 10−4 Then it is assumed that the type of a function α(t) is chosen and the dependence of the absolute errors and the discrepancies on parameters β (Tables 2 and 3) and α0 (Table 4) is analyzed. Table 2. β 1 2 3 4 5 6 7 8 9 10 N 29 156 100 73 56 45 37 31 26 23 α(t) = α0 e−βt , α0 = 0.1, τ = 0.1 ∆E σE ∆R 0.26 7.60 · 10−2 0.26 1.09 · 10−2 4.58 · 10−6 1.10 · 10−2 1.06 · 10−2 5.78 · 10−5 1.07 · 10−2 1.15 · 10−2 7.54 · 10−4 1.14 · 10−2 1.52 · 10−2 4.50 · 10−3 1.50 · 10−2 2.08 · 10−2 1.40 · 10−2 2.17 · 10−2 3.90 · 10−2 3.06 · 10−2 3.89 · 10−2 0.11 5.82 · 10−2 8.59 · 10−2 0.18 0.10 0.15 0.24 0.14 0.20 σR 9.10 · 10−2 1.52 · 10−5 9.30 · 10−5 1.10 · 10−3 6.12 · 10−3 1.84 · 10−2 3.80 · 10−2 6.52 · 10−2 0.11 0.16 Table 3. β 1 2 3 4 5 6 7 8 9 10 α(t) = α0 e−βt , α0 = 0.1, τ = 0.6 ∆E σE ∆R 0.23 5.35 · 10−2 0.26 4.13 · 10−2 1.00 · 10−3 1.09 · 10−2 8.90 · 10−2 5.04 · 10−2 1.05 · 10−2 1.46 · 10−2 1.29 · 10−4 1.25 · 10−2 1.76 · 10−2 7.14 · 10−4 1.84 · 10−2 0.16 3.85 · 10−3 2.98 · 10−2 0.38 2.58 · 10−2 0.12 0.83 1.48 0.18 0.83 01.48 0.21 0.83 01.48 0.23 N 4 20 16 12 9 7 6 5 4 3 σR 0.12 6.11 · 10−5 1.73 · 10−4 1.79 · 10−3 1.08 · 10−2 3.22 · 10−2 5.61 · 10−2 8.60 · 10−2 0.15 0.19 Table 4. α0 10−3 10−2 10−1 N 71 79 85 α(t) = α0 e−βt , β = 3.5, τ = 0.1 ∆E σE ∆R 1.62 · 10−2 9.77 · 10−4 1.61 · 10−2 1.33 · 10−2 4.55 · 10−4 2.29 · 10−2 1.07 · 10−2 2.84 · 10−4 1.08 · 10−2 σR 1.33 · 10−2 2.90 · 10−4 3.60 · 10−4 Analyzing the results of the numerical experiments (a part of them is included in the Tables) one concludes the following: (i) The Runge-Kutta method is more stable with respect to changes of the regularization function α(t) and the step size τ , than the Euler method; (ii) for all the considered functions suitable parameters can be chosen, however in the case when α(t) = α0 e−βt the accuracy with which the solution x(t) is calculated is higher; (iii) an appropriate range of values of the parameter α0 is from 0.001 to 0.1, for larger values the accuracy is lower, and for smaller values the processes do not converge; (iv) the range of appropriate values of β is large enough: from 2 to 6 for τ = 0.6 and from 2 to 7 for τ = 0.1. Remark 4.1. The reason why result (i), formulated above, is emphasized can be understood if one remembers that problem (4.1) is ill-posed. If a problem is illposed then the usage of a higher-order accuracy difference scheme (or quadrature formula) may lead to less accurate results, as was observed in the literature (see e.g. (4 ), p. 155). On the other hand, if a problem is well-posed, then the usage of a higher-order accuracy scheme should lead to more accurate results. Acknowledgments The authors thank Dr. V.Protopopescu for useful remarks. References 1. R.G. Airapetyan and I.V. Puzynin, Newtonian iterative scheme with simultaneous iterations of inverse derivative, Comp. Phys. Comm. 102 (1997) 97-108. 2. A.B. Bakushinskii, Iterative methods for nonlinear operator equations without regularity. New approach, Dokl. Russian Acad. Sci. 330 (1993) 282-284. 3. M.K. Gavurin, Nonlinear functional equations and continuous analogies of iterative methods, Izv. Vuzov. Ser. Matematika. 5 (1958) 18–31. 4. V.A.Morozov and A.I.Grebennikov, Methods of the Solution of Ill-Posed Problems. Alghorithmic Approach. (Moscow University Press, Moscow, 1992). 5. J.M. Ortega and W.C. Rheinboldt, Iterative Solution of Nonlinear Equations in Several Variables (Academic Press, 1970). 6. J. Szarski, Differential Inequalities (PWN, Warszawa, 1967). 7. V.V. Vasin and A.L. Ageev, Ill-Posed Problems with a priori Information (Nauka. Ekaterinburg, 1993). 8. E.P. Zhidkov and I.V. Puzynin, The solution of the boundary problems for second order nonlinear differential equations by means of the stabilization method, Soviet Math. Dokl. 8 (1967) 614-616.