Dynamical systems method for solving nonlinear equations with monotone operators

A version of the Dynamical Systems Method (DSM) for solving ill-posed nonlinear equations with monotone operators in a Hilbert space is studied in this paper. An a posteriori stopping rule, based on a discrepancy-type principle is proposed and justified mathematically. The results of two numerical experiments are presented. They show that the proposed version of DSM is numerically efficient. The numerical experiments consist of solving nonlinear integral equations.


Introduction
In this paper we study a Dynamical Systems Method (DSM) for solving the equation where F is a nonlinear twice Fréchet differentiable monotone operator in a real Hilbert space H, and equation (1.1) is assumed solvable. Monotonicity means that Here, ·, · denotes the inner product in H. It is known (see, e.g., [8]), that the set N := {u : F (u) = f } is closed and convex if F is monotone and continuous. A closed and convex set in a Hilbert space has a unique minimal-norm element. This element in N we denote y, F (y) = f . We assume that where u 0 ∈ H is an element of H, R > 0 is arbitrary, and f = F (y) is not known; but f δ , the noisy data, are known and f δ − f ≤ δ. If F ′ (u) is not boundedly invertible, then solving for u given noisy data f δ is often (but not always) an illposed problem. When F is a linear bounded operator many methods for a stable solution of (1.1) were proposed (see [4]- [8] and the references therein). However, when F is nonlinear then the theory is less complete. The DSM for solving equation (1.1) was extensively studied in [8]- [15]. In [8] the following version of the DSM for solving equation (1.1) was studied: Here F is a monotone operator, and a(t) > 0 is a continuous function, defined for all t ≥ 0, strictly monotonically decaying, lim t→∞ a(t) = 0. These assumptions on a(t) hold throughout the paper and are not repeated. Additional assumptions on a(t) will appear later. Convergence of the above DSM was proved in [8] for any initial value u 0 with an a priori choice of stopping time t δ , provided that a(t) is suitably chosen. The theory of monotone operators is presented in many books, e.g., in [1], [7], [16]. Most of the results of the theory of monotone operators, used in this paper, can be found in [8]. In [6] methods for solving nonlinear equations in a finitedimensional space are discussed.
In this paper we propose and justify a stopping rule based on a discrepancy principle (DP) for the DSM (1.4). The main result of this paper is Theorem 3.1 in which a DP is formulated, the existence of the stopping time t δ is proved, and the convergence of the DSM with the proposed DP is justified under some natural assumptions apparently for the first time for a wide class of nonlinear equations with monotone operators.
These results are new from the theoretical point of view and very useful practically. The auxiliary results in our paper are also new and can be used in other problems of numerical analysis. These auxiliary results are formulated in Lemmas 2.2-2.4, 2.7, 2.10, 2.11, and in the remarks. In particular, in Remark 3.3 we emphasize that the trajectory of the solution stays in a ball of a fixed radius R for all t ≥ 0.
In Section 4 the results of two numerical experiments are presented. In the second experiment we demonstrate numerically that our method for solving equation (1.1) can be used even for a wider class of equations than the basic Theorem 3.1 guarantees.
Let a = a(t), 0 < a(t) ց 0, and assume a ∈ C 1 [0, ∞). Then the solution V δ (t) := V δ,a(t) of (2.1) is a function of t. From the triangle inequality one gets From Lemma 2.2 it follows that for large a(0) one has .
where ǫ > 0 is sufficiently small, for sufficiently large a(0) > 0. Below the words decreasing and increasing mean strictly decreasing and strictly increasing.
Proof. The uniqueness of t 1 follows from Lemma 2.3. We have F (y) = f , and Here the inequality V δ − y, F (V δ ) − F (y) ≥ 0 was used. Therefore, On the other hand, we have This implies (2.6) a V δ − y ≤ a y + δ.
Similarly, from the equation From (2.8) and (2.9), one gets the following estimate: Let us recall the following lemma, which is basic in our proofs.
Remark 2.8. In the proof of Lemma 2.7 a(0) and λ can be chosen so that a(0) λ is uniformly bounded as δ → 0 regardless of the rate of growth of the constant where we have assumed, without loss of generality, that 0 < δ < 1. With this choice of d and λ, the ratio a(0) λ is bounded uniformly with respect to δ ∈ (0, 1) and does not depend on R.
Indeed, with the above choice one has a(0) ≤c, wherẽ c > 0 is a constant independent of δ, and one can assume that λ ≥ 1 without loss of generality.
This remark is used in Remark 3.3, where we prove that the trajectory of u δ (t), defined by (3.1), stays in a ball B(u 0 , R) for all 0 ≤ t ≤ t δ , where the number t δ is defined by formula (3.3) (see below), and R > 0 is sufficiently large. An upper bound on R is given in Remark 3.3.
Indeed, if, for example, u 0 = 0, then by Lemmas 2.2 and 2.3 one gets Lemma 2.10 is proved.
Proof. Let p = 1 2 in Lemma 2.10. Then This implies Multiplying (2.34) by e s 2 V δ (s) , integrating from 0 to t, using inequality (2.33) and the fact that V δ (s) is nondecreasing, one gets This implies inequality (2.32). Lemma 2.11 is proved.

Main result
, A a := A + aI, where I is the identity operator, and u δ (t) solves the following Cauchy problem: where C 1 > 1 and ζ ∈ (0, 1] are some constants. We also assume, without loss of generality, that δ ∈ (0, 1).
Assume that equation F (u) = f has a solution, possibly nonunique, and y is the minimal norm solution to this equation. Let f be unknown, but f δ be given, where C 1 > 1 and 0 < ζ ≤ 1. If ζ ∈ (0, 1) and t δ satisfies (3.3), then Let We use Taylor's formula and get This t 0 exists and is unique since a(t) > 0 monotonically decays to 0 as t → ∞.

Let us prove (3.4).
From (3.15) with t = t δ , and from (2.10), one gets Thus, for sufficiently small δ, one gets whereC < C 1 is a constant. Therefore, We claim that Let us prove (3.20). Using (3.1), one obtains This and (2.1) imply Multiplying (3.21) by v, one obtains the following two inequalities:

Inequalities (3.23) and (3.26) imply
Inequality (3.28) implies From (3.29) and (3.26), one gets Therefore, (3.31) From the results in Section 2 (see Lemma 2.11), it follows that there exists an a(t) such that For example, one can choose where d, c, b > 0. Moreover, one can always choose u 0 such that If 2b < c, then (3.33) implies e − t 2 a(0) ≤ a(t). Therefore,

Thus, lim
Since V δ (t) increases (see Lemma 2.3), the above formula implies lim δ→0 a(t δ ) = 0. Since 0 < a(t) ց 0, it follows that lim δ→0 t δ = ∞, i.e., (3.20) holds. It is now easy to finish the proof of the Theorem 3.1. From the triangle inequality and inequalities (3.14) and (2.8) one obtains Here we have used the fact that t δ < t 0 (see the proof of Theorem 3.1). Since one can choose a(t) and λ so that a(0) λ is uniformly bounded as δ → 0 and regardless of the growth of M 1 (see Remark 2.8) one concludes that R can be chosen independent of δ and M 1 . Since the function u → arctan 3 u is increasing on R, one has Moreover, Thus, F is a monotone operator. Note that Therefore, the operator F , defined in (4.1), is injective and equation (1.1), with this F , has at most one solution.
The Fréchet derivative of F is If u(x) vanishes on a set of positive Lebesgue measure, then F ′ (u) is not boundedly invertible. If u ∈ C[0, 1] vanishes even at one point x 0 , then F ′ (u) is not boundedly invertible in H.
In numerical implementation of the DSM, one often discretizes the Cauchy problem (3.1) and gets a system of ordinary differential equations (ODEs). Then, one can use numerical methods for solving ODEs to solve the system of ordinary differential equations obtained from discretization. There are many numerical methods for solving ODEs (see, e.g., [2]).
In practice one does not have to compute u δ (t δ ) exactly but can use an approximation to u δ (t δ ) as a stable solution to equation (1.1). To calculate such an approximation, one can use, for example, the iterative scheme u n+1 = u n − (F ′ (u n ) + a n I) −1 (F (u n ) + a n u n − f δ ), (4.6) and stop iterations at n := n δ such that the following inequality holds: The existence of the stopping time n δ is proved in [3, p. 733] and the choice u 0 = 0 is also justified in this paper. Iterative scheme (4.6) and stopping rule (4.7) are used in the numerical experiments. We proved in [3, p. 733] that u n δ converges to u * , a solution of (1.1). Since F is injective as discussed above, we conclude that u n δ converges to the unique solution of equation (1.1) as δ tends to 0. The accuracy and stability are the key issues in solving the Cauchy problem. The iterative scheme (4.6) can be considered formally as the explicit Euler's method with the stepsize h = 1 (see, e.g., [2]). There might be other iterative schemes which are more efficient than scheme (4.6), but this scheme is simple and easy to implement. Integrals of the form 1 0 e −|x−y| h(y)dy in (4.1) and (4.5) are computed by using the trapezoidal rule. The noisy function used in the test is The noise level δ and the relative noise level are defined by the formulas In the test κ is computed in such a way that the relative noise level δ rel equals some desired value, i.e., We have used the relative noise level as an input parameter in the test. In all the figures the x-variable runs through the interval [0, 1], and the graphs represent the numerical solutions u DSM (x) and the exact solution u exact (x).
In the test we took h = 1, C = 1.01, and γ = 0.99. The exact solution in the test is otherwise, here x ∈ [0, 1], and the right-hand side is f = F (u e ). As mentioned above, F ′ (u) is not boundedly invertible in any neighborhood of u e .
It is proved in [3] that one can take a n = d 1+n , and d is sufficiently large. However, in practice, if we choose d too large, then the method will use too many iterations before reaching the stopping time n δ in (4.7). This means that the computation time will be large in this case. Since and we choose d = C 0 δ γ , C 0 > 0. In the experiments our method works well with C 0 ∈ [7,10]. In numerical experiments, we found out that the method diverged for smaller C 0 . In the test we chose a n by the formula a n := C 0 δ 0.99 n+1 . The number of nodal points, used in computing integrals in (4.1) and (4.5), was N = 100. The accuracy of the solutions obtained in the tests with N = 30 and N = 50 was slightly less accurate than the one for N = 100.
Numerical results for various values of δ rel are presented in Table 1. In this experiment, the noise function f noise is a vector with random entries normally distributed, with mean value 0 and variance 1. Table 1 shows that the iterative scheme yields good numerical results. Table 1. Results when C 0 = 7, N = 100 and u = u e .  Figure 2 presents the numerical results when N = 100 and C 0 = 7 with δ = 0.003 and δ = 0.001. In these cases, it took 58 and 59 iterations to get the numerical solutions for δ rel = 0.003 and δ rel = 0.001, respectively.
We also carried out numerical experiments with u(x) ≡ 1, x ∈ [0, 1], as the exact solution. Note that F ′ (u) is boundedly invertible at this exact solution. However, in any arbitrarily small (in L 2 norm) neighborhood of this solution, there are infinitely many elements u at which F ′ (u) is not boundedly invertible, because, as we have pointed out earlier, F ′ (u) is not boundedly invertible if u(x) is continuous and vanishes at some point x ∈ [0, 1]. In this case one cannot use the usual methods like Newton's method or the Newton-Kantorovich method. Numerical results for this experiment are presented in Table 2.   Therefore, the assumptions of Theorem 3.1 are not satisfied. Our goal is to show by this numerical example, that numerically our method may work for an even wider class of problems than that covered by Theorem 3.1.
The operator B is compact in H = L 2 [0, 1]. The operator u −→ u 3 is defined on a dense subset D of L 2 [0, 1], for example, on D := C[0, 1]. If u, v ∈ D, then (4.10) This and the inequality B(u − v), u − v ≥ 0, followed from equality (4.3), imply Note that the equal sign of inequality (4.10) happens iff u = v a.e. in Lebesgue measure. Thus, F is injective. Therefore, the element u n δ obtained from the iterative scheme (4.6) and the stopping rule The Fréchet derivative of F is If u(x) vanishes on a set of positive Lebesgue measure, then F ′ (u) is obviously not boundedly invertible. If u ∈ C[0, 1] vanishes even at one point x 0 , then F ′ (u) is not boundedly invertible in H.
We also use the iterative scheme (4.6) with the stopping rule (4.7). We use the same exact solution u e as in (4.8). The right-hand side f is computed by f = F (u e ). Note that F ′ is not boundedly invertible in any neighborhood of u e .
In experiments we found that our method works well with C 0 ∈ [1,4]. Indeed, in the test we chose a n by the formula a n := C 0 δ 0.9 n+6 . The number of node points used in computing integrals in (4.1) and (4.5) was N = 30. In the test, the accuracy of the solutions obtained when N = 30, N = 50 were slightly less accurate than the one when N = 100.
Numerical results for various values of δ rel are presented in Table 3. In this experiment, the noise function f noise is a vector with random entries normally distributed of mean 0 and variance 1. Table 3 shows that the iterative scheme yields good numerical results. Table 3. Results when C 0 = 2 and N = 100.  We have included the results of the numerical experiments with u(x) ≡ 1, x ∈ [0, 1], as the exact solution. The operator F ′ (u) is boundedly invertible in L 2 ([0, 1]) at this exact solution. However, in any arbitrarily small L 2 -neighborhood of this solution, there are infinitely many elements u at which F ′ (u) is not boundedly invertible as was mentioned above. Therefore, even in this case one cannot use the usual methods such as Newton's method or the Newton-Kantorovich method. Numerical results for this experiment are presented in Table 4. Table 4. Results when C 0 = 1, N = 30 and u(x) = 1, x ∈ [0, 1]. From the numerical experiments we can conclude that the method works well in this experiment. Note that the function F used in this experiment is not defined on the whole space H = L 2 [0, 1] but defined on a dense subset D = C[0, 1] of H.