On the Solution of a Class of Cauchy Integral Equations

arXiv:math/0512029v1 [math.CA] 1 Dec 2005 ON THE SOLUTION OF A CLASS OF CAUCHY INTEGRAL EQUATIONS E. DE MICHELI AND G. A. VIANO Abstract. In a previous paper we have presented a new method for solving a class of Cauchy integral equations. In this work we discuss in detail how to manage this method numerically, when only a finite and noisy data set is available: particular attention is focused on the question of the numerical stability. 1. Introduction In a previous paper [3] we have proved the following theorems. Theorem 1. Let us consider the following series ∞ 1 X an z n 2π n=0 (z = x + iy; x, y ∈ R), (1) p and suppose that the set of numbers {fn }∞ 0 , fn = n an (p ≥ 0, n = 0, 1, 2, . . .) satisfies the following Hausdorff conditions n  2+ǫ X n 2+ǫ (1+ǫ) ∆i f(n−i) <M (n = 0, 1, 2, . . . ; ǫ > 0), (2) (n + 1) i i=0 where M is a positive constant, and
∆ is the difference operator defined by: ∆fn = Pk k fn+1 − fn , ∆k fn = m=0 (−1)m m fn+k−m . Then: (1) series (1) converges uniformly to a function f (z) analytic in the unit disk D0 = {z | |z| < 1}; (2) f (z) admits a holomorphic extension to the “cut-plane” {z ∈ C \ (1, +∞)}; (3) the jump function F (x) = −i[f+ (x) − f− (x)], (f± (x) = limη→0 f (x ± iη)), η>0 belongs to L2 (1, +∞), and, moreover, if p ≥ 1, it is a function of class C (p−1) . The ǫ in (2) (missing in Ref. [3]) is needed to guarantee the continuity of ã(λ) (Carlsonian interpolation of the an ’s) at λ = −1/2 + iν, ν ∈ R (see next formula (15)). Let us remind the reader that these results in the case p = 0 (i.e. fn = an ) are due to Stein and Wainger [9]; however, these authors do not formulate the requirements on the Taylor coefficients an by the use of the Hausdorff condition (2). Theorem 2. If in series (1) the coefficients an satisfy the Hausdorff condition (2) (i.e., fn = an , n = 0, 1, 2, . . .), then the jump function F (x) can be represented by 1 2 E. DE MICHELI AND G. A. VIANO the following expansion, that converges in the sense of L2 -norm: F (x) = ∞ X cm φm (x) m=0 (x ∈ [1, +∞)), (3) where the coefficients cm are given by: cm =    ∞ √ X 1 (−1)n , an Pm −i n + 2 n! 2 n=0 (4) Pm being the Pollaczek polynomials. The functions φm (x) form a basis in L2 (0, +∞), and are expressed by: e−1/x φm (x) = Bm (x) , (5) x √ where Bm (x) = im 2Lm (2/x), Lm being the Laguerre polynomials. From these results it derives that we can formally solve the following integral equation of Cauchy type: Z +∞ F (x) f (z) = dx, (6) x−z 1 if the infinite set of Taylor coefficients an = f (n) (0)/n! (n = 0, 1, 2, . . .) is known, and, in addition, if they are supposed to satisfy the Hausdorff condition (2) with fn = an . But, at this point, a very serious problem is met: How can this procedure be managed numerically? In fact, it must be noted that, in practice, only a finite number of Taylor coefficients can be known, and, moreover, they are usually affected by numerical errors. Furthermore, if these coefficients are the results of experimental measurements the situation is even worse because also the error proper of any measurement must be considered. Finally, let us remark that the determination of the jump function F (x) starting from an approximate knowledge of the function f (z) is a typical example of improperly posed problem in the sense of Hadamard: the solution does not depend continuously on the data. In order to make this point more clear notice that the problem of solving the Cauchy integral equation (6) is strictly connected to the problem of the analytic continuation up to the cut. We can conformally map the cut z-plane onto the unit disk in the ζ-plane (i.e. the domain |ζ| < 1); in this map the upper (lower) lip of the cut is mapped in the upper (lower) half of the unit circle. Therefore the continuation up to the cut corresponds to the continuation up to the unit circle |ζ| = 1. As is well known, the uniqueness of the analytic continuation does not guarantee its continuity in the L2 or in the uniform topology. Furthermore, performing the analytic continuation up to the boundary of the analyticity domain (i.e. up to the unit circle in the ζ-plane geometry) is a severely ill-posed problem [5]. In the standard methods of regularization, whose credit is generally due to Tikhonov [4, 11], normally some appropriate global bound on the solution are imposed, by assuming some prior knowledge on the solution itself. In this way a subspace of the solution space is determined, and then a solution belonging to this subspace is looked for. If this subspace is compact, then continuity follows from compactness. In the numerical analytic continuation in the unit disk (ζ-plane geometry), if one avoids going up to the boundary, this subspace is obtained by imposing a global bound at the boundary. This point can be easily understood by CAUCHY INTEGRAL EQUATIONS 3 observing that from a bound on the functions at |ζ| = 1 a bound on their derivatives can be derived (in any point inside the unit circle), via the Cauchy integral representation. Conversely, if the analytic continuation is required up to the cut (i.e. up to the boundary), then a global bound only on the functions is not sufficient, and a uniform bound on the first derivative becomes necessary. Then, in this case, the regularization by the use of a-priori bounds is rather cumbersome. Furthermore, it must be observed that in several problems of physical interest the prior knowledge which allows for imposing bounds on the solution is missing, or can be very poor. See on this point Refs. [7, 12], where the difference between synthesis and/or inverse problems is widely discussed: while in the first class of problems the prior knowledge is intrinsic to the problem itself, being a part of its formulation, for the second class (i.e., inverse problems) this is not true. We are thus motivated to go beyond the Tikhonov and Tikhonov–like methods, and, accordingly, to look for regularization procedures that do not make use of prior knowledge. The method that we present in this paper is very far from the Tikhonov–like regularization procedures, and it is actually based on the following observation: when in expansion (3) we use the actual data set, which is usually composed by only a finite number of data affected by noise, then, although the formal series (3) diverges, nevertheless the effect of the errors (and of the fact that the number of data is finite) remains quite small in the beginning of the expansion, and there will exist a point (i.e. a certain value of m) where the divergence sets in. Thus, the basic idea is to stop the expansion at the point where it turns to diverge. This rough and qualitative description can be put in rigorous form by proving that even though series (3) diverges, nevertheless it converges (in the sense of the L2 -norm) as the number of data tends to infinity and the noise vanishes. This is the main point of the paper and it will be proved in Section 2. A very delicate point of the method consists in determining a procedure for the determination of the truncation point in the approximation which will be proved to converge to the unknown jump function F (x). This procedure, that will be illustrated in Section 3, which is devoted to the numerical analysis, is based on the asymptotic behavior of the Pollaczek polynomial Pm [−i(n + 1/2)] for large values of m (at fixed n), whose proof will be given in the Appendix. 2. Solution of the Cauchy integral equation when the data are a finite number of noisy Taylor coefficients (ǫ) Let us assume that only a finite number of noisy Taylor coefficients an are (ǫ) (ǫ,N ) known, and that |an − an | ≤ ǫ (n = 0, 1, 2, . . . , N ; ǫ > 0). Let cm be the following finite sums    N √ X (−1)n (ǫ) 1 (ǫ,N ) . (7) an Pm −i n + cm = 2 n! 2 n=0 (0,∞) Accordingly, we have cm lemma. = cm (see formula (4)). We can then prove the following Lemma 1. The following statements hold true: (i) ∞ X 2 c(0,∞) = kF k2L2 [1,+∞) = C; m m=0 (8) 4 E. DE MICHELI AND G. A. VIANO (ii) ∞ X (ǫ,N ) cm 2 = +∞; (9) m=0 (iii) (ǫ,N ) lim cm = c(0,∞) = cm m (m = 0, 1, 2, . . .); N →∞ ǫ→0 (10) (iv) if k0 (ǫ, N ) is defined as k0 (ǫ, N ) = max{k ∈ N : then k X (ǫ,N ) cm m=0 2 ≤ C}, lim k0 (ǫ, N ) = +∞; (11) (12) N →∞ ǫ→0 (v) the sum (ǫ,N ) Mk = k X (ǫ,N ) cm 2 m=0 (k ∈ N) (13) satisfies the following properties: (a) it increases for increasing values of k; (b) the following relationships hold true (ǫ,N ) Mk (ǫ,N ) ≥ ck 2 ∼ k→+∞ 1 (2k)2N (N !)2 (N fixed). (14) Proof. (i) In Ref. [3] we proved that the coefficients an in expansion (1), which are supposed to satisfy condition (2), admit a unique interpolating function ã(λ) (λ = σ + iν) which belongs to the Hardy space H 2 (C−1/2 ) with norm kãk2 = R +∞ supσ>−1/2 ( −∞ |ã(σ + iν)|2 dν)1/2 (C−1/2 = {λ ∈ C | Re λ > −1/2}). Furthermore, we proved that ã(−1/2 + iν), i.e., the restriction of ã(σ + iν) to the line Re λ ≡ σ = −1/2, belongs to L2 (−∞, +∞), and can be represented by the following expansion which converges in the sense of the L2 –norm:   X ∞ 1 ã − + iν = dm ψm (ν), (15) 2 m=0 √ where dm = 2πcm (see formula (4)), and the functions ψm (ν) are the Pollaczek functions (see [3])   1 1 ψm (ν) = √ Γ + iν Pm (ν), (16) 2 π (see Refs. [2, 10] for what concerns the Pollaczek polynomials). In view of the Parseval equality we have: 2 Z +∞  ∞ X 1 dν. (17) |dm |2 = ã − + iν 2 −∞ m=0 Next, we recall that ã(−1/2 + iν) is the Fourier transform of F (v) exp(v/2), where F (v) represents the jump function associated with the cut located at θ = iv (v ∈ (0, +∞)) in the complex θ-plane geometry (see Fig. 1A in Ref. [3] which refers to CAUCHY INTEGRAL EQUATIONS 5 P∞ the series (1/2π) n=0 an e−inθ ). Therefore, in view of the Plancherel theorem and equality (17) we get (see formula (46) of Ref. [3]): 2 Z +∞ Z +∞  ∞ X 2 1 v/2 dν = 2π e F (v) dv = |dm |2 . (18) ã − + iν 2 0 −∞ m=0 Next we pass from the θ-plane geometry of Theorem 2 of Ref. [3] to the z-plane geometry of Theorem 2’ of Ref. [3] through the transformation z = exp(−iθ). Then, to the cut at θ = iv there corresponds the cut in the z-plane geometry located on the real x-axis (x ≡ Re z) from x = 1 up to x = +∞. Accordingly, the jump function associated with the cut will be F (ln x). Adopting the convention used in [3] we still denote this latter jump function by F (x) in order to avoid a useless proliferation of symbols. Therefore, by the substitution: ev = x, from equality (18) we obtain: Z +∞ ∞ 1 X |F (x)|2 dx = C. (19) |dm |2 = 2π m=0 1 √ (0,∞) Since cm = dm / 2π, from (19) equality (8) follows. (ǫ,N ) (ii) Let us rewrite the sums cm as follows:    N X 1 (ǫ,N ) (ǫ) cm = bn Pm −i n + , (20) 2 n=0 (ǫ) √ (ǫ) 2(−1)n an /n!. Now, we can write the following inequality:    N X 1 bn(ǫ) Pm −i n + = 2 n=0
  PN −1 (ǫ)    bn Pm −i n + 21 n=0 1 (ǫ) . ≥ bN Pm −i N + · 1−
  (ǫ) 2 b P −i N + 1 where bn = (ǫ,N ) cm m N (21) 2 In the lemma proved in the Appendix we show that the asymptotic behavior of the Pollaczek polynomials Pm [−i(n + 1/2)] for large values of m (at fixed n) is given by:    1 (−1)m im ∼ Pm −i n + (2m)n . (22) m→∞ 2 n! Therefore, we have: 
 PN −1 (ǫ) 1 n=0 bn Pm −i n + 2
  (ǫ) bN Pm −i N + 21 (23) 
 PN −1 (ǫ) 1 N −1 (ǫ) −i n + b P X n m n=0 2 bn N ! ∼ ≤ (2m)n−N −−−−→ 0. m→∞
  (ǫ) n! (ǫ) m→∞ 1 b Pm −i N + n=0 b N N 2 From (21), (22) and (23) it follows that for m sufficiently large: (ǫ) (ǫ,N ) ∼ cm m→∞ bN N! (2m)N . (24) 6 E. DE MICHELI AND G. A. VIANO (ǫ,N ) Therefore, limm→∞ |cm | = +∞, and statement (ii) follows. (0,∞) (ǫ,N ) (iii) We can write the difference [cm − cm ] as (N    X (−1)n √ 1 (0,∞) (ǫ,N ) (ǫ) cm − cm = 2 (an − an )Pm −i n + n! 2 n=0 )    ∞ X (−1)n 1 + . an Pm −i n + n! 2 (25) n=N +1 √ P∞ n (0,∞) 1 , it In view of the fact that 2 n=0 (−1) n! an Pm [−i(n + 2 )] converges to cm follows that the second term in bracket (25) tends to zero as N → ∞. Concerning the first term we may write the inequality:       N N X X 1 1 (−1)n 1 (ǫ) ≤ǫ , (26) (an − an )Pm −i n + Pm −i n + n! 2 n! 2 n=0 n=0 (ǫ) where the inequalities |an − an | ≤ ǫ (n = 0, 1, 2, . . . , N ) have been used. Next, by rewriting the Pollaczek polynomials Pm [−i(n + 1/2)] as     X j m 1 1 (m) n+ pj = , (27) Pm −i n + 2 2 j=0 and substituting this expression in inequality (26) we obtain    j m N X X 1  1  (m) n+ pj . ǫ n! 2 n=0 j=0 (28) Pm (m) Next, we perform the limit for N → ∞. In view of the fact that j=0 pj (n+1/2)j j P∞ converges, we can exchange the order of the is finite, and the series n=0 (n+1/2) n! sums and write:  j ∞ m X X 1 1 (m) pj n+ ǫ . (29) n! 2 n=0 j=0 Finally, performing the limit for ǫ → 0, and recalling equality (25), statement (iii) is obtained. (iv) From definition (11) it follows that for k1 = k0 + 1, we have: k1 X (ǫ,N ) cm 2 > C. (30) m=0 Statement (iv) (formula (12)) is proved if we can show that limN →∞ k1 (ǫ, N ) = +∞. ǫ→0 Let us suppose that limN →∞ k1 (ǫ, N ) is finite. Then there should exist a finite ǫ→0 number K (independent of ǫ and N ), such that limN →∞ k1 (ǫ, N ) ≤ K. Then from ǫ→0 inequality (30) we have: k1 (ǫ,N ) C< X m=0 (ǫ,N ) cm 2 ≤ K X m=0 2 (ǫ,N ) cm . (31) CAUCHY INTEGRAL EQUATIONS 7 But as N → ∞, ǫ → 0 we have (recalling also statement (iii) formula (10)) K X C< c(0,∞) m 2 < ∞ X 2 c(0,∞) m = C, (32) m=0 m=0 which leads to a contradiction. Then statement (iv) follows. (ǫ,N ) (v) Concerning statement (a), it follows obviously from definition (13) of Mk . Finally, the first one of relationships (14) is obvious; the second one follows from the asymptotic behavior of Pm [−i(n + 1/2)] at large m (for fixed n), proved in the lemma of the Appendix.  Remark. From statement (v) and formula (12) it follows that, for large values of (ǫ,N ) N and small values of ǫ, the sum Mk presents a plateau in a neighborhood of k = k0 , if the function F (x) being reconstructed is sufficiently regular (see Section 3 and figs. 1B, 2B, 3A, 4A). Now we may introduce the following approximation of the jump function F (x) (see expansion (3)): k0 (ǫ,N ) F (ǫ,N ) (x) = X (ǫ,N ) ψm (x) cm m=0 (x ∈ (1, +∞)). (33) Approximation F (ǫ,N ) (x) is defined through the truncation number k0 (ǫ, N ); the (ǫ,N ) latter can be numerically determined by plotting the sum Mk versus k, and exploiting properties (a) and (b), proved in statement (v) of the previous lemma, and the property stated in the remark above. For a more detailed analysis of this point see Section 3 where some numerical examples are given. Now, we want to prove that the approximation F (ǫ,N ) (x) converges asymptotically to F (x) in the sense of the L2 -norm, as N → ∞ and ǫ → 0. We can prove the following theorem. Theorem 3. The equality lim N →∞ ǫ→0 F − F (ǫ,N ) L2 [1,+∞) =0 (34) holds true. Proof. From the Parseval equality it follows that: ( ∞ ) k0 X X 2 2 2 (ǫ,N ) (0,∞) (ǫ,N ) (0,∞) F −F = cm + . cm − cm L2 [1,+∞) m=k0 +1 2 (0,∞) = C and limN →∞ m=0 cm ǫ→0 P∞ (ǫ,N ) 2 limN →∞ m=k0 +1 |cm | = 0. ǫ→0 Since P∞ (35) m=0 k0 (ǫ, N ) = +∞, it follows that It is then convenient to rewrite the second term of the r.h.s. of formula (35) as follows. Let us denote by ( (0,∞) cm if m is even; (0,∞) gm = (36) (0,∞) if m is odd; −icm 8 E. DE MICHELI AND G. A. VIANO (ǫ,N ) gm = Notice that (ǫ,N ) gm (ǫ,N ) Pk0 cm m=0 ( (0,∞) − cm (ǫ,N ) cm if m is even; (37) (ǫ,N ) −icm if m is odd. 2 2 Pk0  (ǫ,N ) (0,∞) (0,∞) = g − g , and gm and m m m=0 are real. Next, we introduce the following functions: ∞ X (0,∞) I[m,m+1[ (x), G(0,∞) (x) = gm G(ǫ,N ) (x) m=0 ∞ X = (ǫ,N ) gm I[m,m+1[ (x), (38a) (38b) m=0 where IE is the characteristic function of the set E. From statements (i), (ii), and (iii) of the previous lemma (formulae (8), (9) and (10)) we obtain Z +∞  ∞  2 2 X (0,∞) = C, (39) G(0,∞) (x) dx = gm 0 Z m=0 +∞ 0 ∞  2  2 X (ǫ,N ) gm = +∞, G(ǫ,N ) (x) dx = (40) m=0 G(ǫ,N ) (x) −−−−→ G(0,∞) (x) N →∞ ǫ→0 (x ∈ [0, +∞)). (41) Hereafter we assume, for the sake of simplicity and without loss of generality, that (ǫ,N ) every term gm is different from zero. Next, let X(ǫ, N ) be the unique root of the RX
2 RX
2 equation 0 G(ǫ,N ) (x) dx = C. Let us indeed observe that 0 G(ǫ,N ) (x) dx is a continuous non–decreasing function which is zero for X = 0, and +∞ for X → +∞. Furthermore, from statement (iv) of the previous lemma (formula (12)) we have limN →∞ X(ǫ, N ) = +∞. ǫ→0 Then we can write Z X(ǫ,N ) h i2 G(ǫ,N ) (x) − G(0,∞) (x) dx 0 = Z +∞ X(ǫ,N )  (0,∞) G Z 2 (x) dx − 2 X(ǫ,N ) 0 h i G(0,∞) (x) G(ǫ,N ) (x) − G(0,∞) (x) dx. (42) Next, we perform the limit for N → ∞ and ǫ → 0. Concerning the first term at the r.h.s. of formula (42) we have Z +∞  2 G(0,∞) (x) dx = 0. (43) lim N →∞ ǫ→0 X(ǫ,N ) For what concerns the second term we introduce the following function: ( G(ǫ,N ) (x) − G(0,∞) (x) if 0 ≤ x ≤ X(ǫ, N ), (ǫ,N ) H (x) = 0 if x > X(ǫ, N ). Then, we have by the use of the Schwarz inequality Z +∞ 2 H (ǫ,N ) (x) dx ≤ 4C (N < ∞, ǫ > 0). 0 (44) (45) CAUCHY INTEGRAL EQUATIONS 9 Moreover, from (41) we have H (ǫ,N ) (x) −−−−→ 0 N →∞ ǫ→0 x ∈ [0, +∞). (46) The family of functions {H (ǫ,N ) } is bounded in L2 [0, +∞), therefore it has a subsequence which is weakly convergent in L2 [0, +∞). The limit of this subsequence is zero. In fact, let us observe that |H (ǫ,N ) (x)| ≤ 2C; then we consider the function H (ǫ,N ) (x)Φ(x), where Φ(x) is an arbitrary element of the class of functions Cc∞ (R+ ). Then, we have |H (ǫ,N ) (x)Φ(x)| ≤ 2C|Φ(x)|, and this inequality does not depend on N and ǫ. In view of the Lebesgue dominated convergence theorem we can then write (see also the limit (46)) Z +∞ lim sup H (ǫ,N ) (x)Φ(x) dx = 0. (47) N →∞ ǫ→0 0 Since the set of functions Cc∞ (R+ ) is everywhere dense in L2 [0, +∞), given an arbitrary function ψ ∈ L2 [0, +∞) and an arbitrary number η > 0, there exists a function Φk ∈ Cc∞ (R+ ) such that kψ − Φk kL2 [0,+∞) < η. Furthermore, through the Schwarz inequality we have Z +∞ H (ǫ,N ) (x) [ψ(x) − Φk (x)] dx 0 ≤ Z +∞ H (ǫ,N ) (x) 2 1/2 Z +∞ 0 0 |ψ(x) − Φk (x)| 2 1/2 √ ≤ 2 Cη. From equality (47) and inequalities (48) we thus conclude that Z +∞ lim sup H (ǫ,N ) (x)ψ(x) dx = 0, N →∞ ǫ→0 (48) (49) 0 for any ψ ∈ L2 [0, +∞). Next, by using the same type of arguments, we can state that if there is an arbitrary subsequence belonging to the family {H (ǫ,N )}, that converges weakly in L2 [0, +∞), then the weak limit of this subsequence is necessarily zero. Finally, from the uniqueness of the (weak) limit point it follows that the whole family {H (ǫ,N ) } converges weakly to zero in L2 [0, +∞). We can thus write Z +∞ lim G(0,∞) (x)H (ǫ,N ) (x) dx = 0, (50) N →∞ ǫ→0 0 and from inequality (42) it follows Z X(ǫ,N ) h i2 G(ǫ,N ) (x) − G(0,∞) (x) dx = 0. lim N →∞ ǫ→0 Since Pk0 m=0 (ǫ,N ) cm (51) 0 (0,∞) − cm 2 lim ≤ R X(ǫ,N )  0 k0 X N →∞ ǫ→0 m=0 2 G(ǫ,N ) (x) − G(0,∞) (x) dx, we have (ǫ,N ) cm − c(0,∞) m 2 and, in view of equality (35), the theorem is proved. = 0, (52)  10 E. DE MICHELI AND G. A. VIANO 3. Numerical analysis In this section various numerical aspects of the issues discussed in the previous section will be illustrated. In particular, the effectiveness and the accuracy of the reconstruction of the jump function F (x) across the cut will be tested when only a finite number of coefficients an of the Taylor series (1) is known. Moreover, tests performed with coefficients an corrupted by random noise will be illustrated. The results hereafter presented summarize a large number of numerical tests performed on a variety of sample functions, including smooth functions, regular oscillating functions and discontinuous functions. We shall illustrate only the somehow extreme cases, that is, smooth (see Figs. 1 and 2) and discontinuous functions (see Fig. 3), since the behavior of the oscillating ones is very similar to the case of the smooth functions. In Fig. 1 the basic steps for the reconstruction of jump function in the presence of only round-off noise are shown (see the legend for the numerical details). In this 2 example the test function is F1 (x) = x4 +x(x−1) 3 +x2 +x+1 (x ∈ [1, +∞)). This function is characterized to be smooth so that we expect the corresponding reconstruction to be quite satisfactory. Fig. 1A illustrates a sample of the coefficients an , which are computed as Z +∞ an = x−n−1 F1 (x) dx. (53) 1 Then, the coefficients an do satisfy condition (2) of Theorem 1. In Fig. 1B the plot (0,N ) of Mk versus k (see Eq. (13)) is shown. According to the properties summarized in Lemma 1 (see also the remark after the proof of Lemma 1), a quite extended (0,N ) plateau is clearly present before Mk starts increasing rapidly (see Eq. (14)). This fact allows us to determine the truncation point k0 (0, N ) which is needed (0,N ) for constructing the regularized approximation F1 (see Eq. (33)). A simple algorithm for the automatic determination of the index k0 has been implemented; it exploits both the properties provided by Lemma 1 that k0 lies approximately on (0,N ) a plateau, and, moreover, that such a plateau must be located before Mk starts growing as a power of 2N (see formula (14)). (ǫ,N ) In the general case with ǫ 6= 0, the knowledge of the asymptotic behavior of Mk allows for restricting the range of k0 by defining an upper limit ka (k0 < ka ), that represents the value of k where approximately the asymptotic behavior sets in; in (ǫ,N ) practice, ka is set as the value of k where Mk and its asymptotic behavior starts being close enough. Then, the candidate plateaux are located by selecting the extended intervals of k < ka where the modulus of the first numerical derivative of (ǫ,N ) Mk is sufficiently small. Finally, k0 is chosen as the largest value of k belonging to the interval which is closest, but inferior, to ka . Figure 1B shows a quite simple situation, where only one plateau is present in the interval from about k = 15 to k = 65, and whose value, as expected from Eq. (11), is approximately the squared norm of the function F1 (x) (see the horizontal dashed line). Panels C and D show two reconstructions of the function F1 (x); in Fig. 1C we have k0 = 61, which is well within the plateau of Fig. 1B, and the comparison between the actual function F1 (x) and its approximation F (0,30) (x) (dotted line) shows the good quality of the reconstruction. Similar results (not displayed), obtained for different values of k0 (15 < k0 < 65), indicate that the actual choice CAUCHY INTEGRAL EQUATIONS an A M k(0,N) 11 B 0.009 0.04 0.008 0.007 0.03 0.006 0.02 0.005 || F1|| 0.01 2 0.003 0 0 2 4 6 8 10 12 14 16 0 10 20 30 40 n F(x) 1 50 60 70 80 k C F(x) 1 D 0.04 0.03 0.035 0.025 0.03 0.02 0.025 0.02 0.015 0.015 0.01 0.01 0.005 0 0.005 0 2 4 6 8 10 12 x 14 16 18 20 0 0 2 4 6 8 10 12 14 16 18 20 x Figure 1: Reconstruction of a smooth jump function. F1 (x) = (x−1)2 /(x4 +x3 +x2 +x+1), x ∈ [1, +∞); the number of coefficients an used for the computations is N = 30. (A) (0,30) Noiseless coefficients an computed according to Eq. (53). (B) Mk versus k. The dashed horizontal line indicates the squared norm of the true function. The algorithm for recovering the truncation index gave k0 (0, 30) = 61. (C) Reconstruction of the jump function. The solid line represents the actual function F1 (x). The dots are samples of the (0,30) computed by using k0 (0, 30) = 61. (D) Comparison reconstructed jump function F1 (0,30) between F1 (x) and F1 computed by using k0 (0, 30) = 75. of k0 is not critical, provided that it is located in the interval corresponding to the plateau. Figure 1D shows that at k0 = 75, that is just after the correct plateau, the reconstruction deteriorates evidently. Figure 2 illustrates another example of a smooth sample function that gives rise to a good reconstruction; in this case F2 (x) = (x − 1)2 exp(−x/2) (x ∈ [1, +∞)). PN In Fig. 2D we compare the approximation of f (z) given by fN (z) = n=0 an z n (Re z ∈ [−1/2, 1/2], Im z = 0) with the approximation of f (z) obtained by the 12 E. DE MICHELI AND G. A. VIANO an A M k(0,N) B 45 1.6 40 1.4 35 1.2 30 1 25 0.8 20 0.6 15 0.4 || F2|| 0.2 2 5 0 0 2 4 6 8 10 12 14 16 0 0 20 40 60 n F(x) 2 80 100 k C f(z) D 2 1.2 1.95 1 1.9 1.85 0.8 1.8 0.6 1.75 1.7 0.4 1.65 0.2 1.6 0 2 4 6 8 10 12 x 14 16 18 20 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 Re(z) Figure 2: Reconstruction of a smooth jump function. F2 (x) = (x − 1)2 exp(−x/2), (0,30) x ∈ [1, +∞); N = 30. (A) Noiseless coefficients an . (B) Mk versus k. The dashed 2 horizontal line indicates kF2 k . The algorithm for recovering the truncation index gave (0,30) k0 (0, 30) = 82. (C) Comparison between the actual solution F2 (x) (solid line) and F2 (dashed dotted line) computed by using k0 (0, 30) = 82. (D) Comparison between the P30 truncated Taylor series a z n (solid line) and the approximation computed through n=0 n (0,30) the integral of Cauchy type (see Eq. (6)) by using the function F2 (dots), for Re z ∈ [−1/2, 1/2] and Im z = 0. Cauchy integral (see Eq. (6)) representing the jump function F2 (x) with its regu(0,30) larized approximation F2 (x). For the sake of completeness, it should be mentioned that the erratic behavior of the noise can rarely produce very short plateaux located between the true value of (ǫ,N ) k0 and before Mk starts following its asymptotic behavior (i.e. for k ≃ ka ). In this case our procedure could fail to recover the correct value of k0 , and this minor drawback has been solved heuristically by simply rejecting the plateaux shorter than a given threshold L = 5. CAUCHY INTEGRAL EQUATIONS M k(0,N) A F(x) 3 13 B 1.4 30 1.2 25 1 0.8 20 0.6 15 0.4 0.2 || F3|| 2 5 0 0 10 20 30 40 50 60 70 80 0 5 10 k F(x) 3 15 20 x C f(z) D 1.4 2.8 1.2 1 2.6 0.8 2.4 0.6 0.4 2.2 0.2 2 0 0 5 10 x 15 20 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 Re(z) Figure 3: Reconstruction of a discontinuous jump function. F3 (x) = 1 if x ∈ [1, 10] and (0,30) F3 (x) = 0 elsewhere; N = 30. (A) Mk versus k. The dashed horizontal line indicates 2 kF3 k . The algorithm for recovering the truncation index gave k0 (0, 30) = 23. (B) Com(0,30) parison between the actual solution F3 (x) (solid line) and F3 (dashed line) computed by using k0 (0, 30) = 23. (C) Comparison between the actual solution F3 (x) (solid line) (0,30) and F3 (dashed line) computed by using k0 (0, 30) = 44. (D) Comparison between the P30 truncated Taylor series a z n (solid line) and the approximation computed through n=0 n (0,30) the integral of Cauchy type (see Eq. (6)) by using the function F3 with k0 (0, 30) = 23 (dashed line), for Re z ∈ [−1/2, 1/2] and Im z = 0. However, the situation can become more difficult, and an example is given in Fig. 3. Here the jump function F3 (x) is a discontinuous one, being constant in the range [1, 10] and null elsewhere (see Fig. 3B). In this case the lack of smoothness (ǫ,N ) of F3 causes the plateau of Mk to be hardly visible, and, consequently, the definition of the truncation point k0 becomes uncertain (see Fig. 3A). In this case, our algorithm gave k0 = 23, which belongs to the only detectable plateau, and the corresponding reconstruction is shown in Fig. 3B. Fig. 3C displays the 14 E. DE MICHELI AND G. A. VIANO reconstruction of the jump function obtained with k0 = 44, that represents the (ǫ,N ) index k such that Mk ≃ kF3 k2 (see formula (11)). In spite of the fact that, according to definition (11), k0 = 44 is the true truncation index, the corresponding approximation is evidently worse than the one obtained with k0 = 23 (see Fig. 3B). This fact can be easily explained by recalling that our reconstruction strategy rests on a finite truncation of a series expansion (k0 < ka < ∞ for N < ∞), and that the convergence to the true function is only asymptotic. Then, we pay entirely the penalty of the Gibbs-like phenomenon, so that our procedure can lead to inaccurate reconstructions in the case of discontinuous jump functions. An example of the dependence of the reconstruction on the noise is sketched in Fig. 4. Here the test function is again F2 (x), and the data illustrated in this figure are representative of an extensive numerical analysis performed over several different test functions that produced similar results. The coefficients an have been noised by adding white noise, simulated by computer generated random numbers (ǫ,N ) uniformly distributed in the interval [−ǫ, ǫ]. Fig. 4A groups the functions Mk computed with different values of ǫ in the range from ǫ = 10−10 through ǫ = 10−2 . It is clear that the plateaux become shorter, and, correspondingly, k0 decreases as ǫ increases. This fact indicates how, when the noise increases, the number (ǫ,N ) of coefficients cm that carry reliable information for reconstructing the jump function decreases. Examples of such a reconstructions are plotted in Fig. 4B. A crude measure of the dependence of the reconstruction stability on the noise ǫ is shown in Fig. 4C, where the truncation index k0 (ǫ, N ), chosen in this case by hand according to its definition (11), is plotted versus the noise bound ǫ. The plot indicates that the stability estimate is only logarithmic, as expected for the problem of the analytic continuation up to the boundary (see [5, 8]). However, in spite of this extremely poor stability, nevertheless it is still possible to obtain acceptable reconstructions of the jump function, as shown in Fig. 4B. Of course, this is not always the case (see for instance the example illustrated in Fig. 3), since the quality of the reconstruction depends strongly on the regularity of the jump function, but, at least for some classes of functions, approximation (33) yields useful results even in the presence of quite noisy data. This numerical examples make milder some radical statements that can be found in the literature (see [5, 8]) discouraging any attempt of numerical analytic continuation up to the boundary in physical situations, and indicate that the numerical procedure here proposed can be successfully applied whenever the function being reconstructed is sufficiently regular. To this purpose it is worth recalling statement (3) of Theorem 1. In Fig. 4D the plot of the reconstruction error, defined as the mean square error of F (ǫ,N ) (x) with respect to F (x) against the global signal–to–noise ratio (SNR) is shown. It can be seen that the reconstruction error is still acceptable for SNR up to about 50 dB, whereas it becomes quite large for smaller SNR. Appendix In this Appendix we prove the following lemma on the asymptotic behavior of the Pollaczek polynomials (see also [6]). CAUCHY INTEGRAL EQUATIONS M k(ε,N) A F(x) 2 15 B 1.2 35 -2 -3 -4 -5 -6 -7 -8 -9 -10 1 30 0.8 -6 -3 0.6 25 0.4 20 0.2 15 0 -0.2 || F2 || 2 -0.4 5 0 20 40 60 80 100 -0.6 2 4 6 k k 0(ε,N) 8 10 12 14 16 18 20 x C Eε 80 D 0.3 70 0.25 60 0.2 50 40 0.15 30 0.1 20 0.05 10 0 0 2 4 6 8 log( ε -1 ) 10 12 14 0 0 20 40 60 80 100 120 140 160 180 200 SNR [dB] (ǫ) Figure 4: Reconstruction of a smooth jump function from noisy coefficients an . The sample function is F2 (x); N = 30. The coefficients an have been noised by adding white (ǫ,30) noise uniformly distributed in the interval [−ǫ, ǫ]. (A) Mk versus k computed for ǫ −2 −10 −1 ranging from ǫ = 10 through ǫ = 10 with step 10 . The rightmost solid line indi(0,30) cates the noiseless Mk . (B) Comparison between the actual solution F2 (x) (solid line) (10−6 ,30) (10−3 ,30) and the reconstructions F2 and F2 (dashed lines). (C) Plot of k0 (ǫ, 30), set by hand according to its definition (11), against log(ǫ−1 ). Each k0 (ǫ, 30) was obtained (ǫ,30) as the average over 300 realizations of the corresponding Mk . (D) Plot of the reconstruction error Eǫ versus the signal–to–noise ratio SNR. Eǫ is defined as the mean square (ǫ,N) error of the reconstruction F2 (x) with respect to the actual function F2 (x), computed in the interval x ∈ [1, 20]. The signal–to–noise ratio is defined as the ratio of the mean power of the sequence of noiseless coefficients {an } to the noise variance. Lemma 2. The Pollaczek polynomials Pm [−i(n + 1/2)] satisfy the following asymptotic behavior for large values of m (at fixed n),    1 (−1)m im Pm −i n + ∼ (2m)n . (54) 2 n! 16 E. DE MICHELI AND G. A. VIANO Im z -5/2 -3/2 -1/2 0 Re z Figure 5: Contour used for evaluating integral (59). (1) (2) Proof. Let us introduce the following functions Qm (y) and Qm (y), that we call the associate Pollaczek functions: Z i +∞ Pm (x)w(x) (1) dx (Im y < 0), (55) Qm (y) = 2 −∞ x−y and Q(2) m (y) i = 2 Z +∞ −∞ Pm (x)w(x) dx x−y (Im y > 0), (56) where w(x) = π1 Γ(1/2 + ix)Γ(1/2 − ix). Let us now consider the following function H(t, y) defined by P∞ Z i +∞ { m=0 tm Pm (x)} w(x) dx (Im y < 0; |t| < 1). (57) H(t, y) = 2 −∞ x−y P∞ The series m=0 tm Pm (x) converges to the generating function of the Pollaczek polynomials for |t| < 1 [2, 10]. The generating function is given by: G(t, x) = (1 − it)ix−1/2 (1 + it)−ix−1/2 . Then, we can rewrite integral (57) as  Z +∞ (1 − it)ix−1/2 (1 + it)−ix−1/2 Γ( 21 + ix)Γ( 12 − ix) i H(t, y) = dx 2π −∞ x−y (Im y < 0; |t| < 1). (58) We can then integrate the r.h.s. of (58) by the contour integration method. Let us, indeed, consider the following integral for Re t = 0, 0 < Im t < 1 and Im y < 0, z      I  1 (1 − it) [(1 − it)(1 + it)]−1/2 1 1 Γ +z Γ −z dz, I(t, y) = 2π γ1 (1 + it) (−iz − y) 2 2 (59) where γ1 is the path shown in Fig. 5. CAUCHY INTEGRAL EQUATIONS 17 (1−it) ), it can be easily checked Rewriting the term [(1 − it)/(1 + it)]z as exp(z ln (1+it) that the contributions along the quarters of circle belonging to the path γ1 vanishes as Re z → −∞. Therefore, through the Cauchy theorem applied to integral (59), we obtain, for Re t = 0, 0 < Im t < 1, Im y < 0:  k ∞ X it + 1 1 1 , (60) H(t, y) = (1 − it) (k + iy + 1/2) it − 1 k=0 in view of the fact that the only singularities which play a role are the simple poles of Γ(z + 1/2) at z = −k − 1/2. Let us now recall that 2 F1 (1, a, a + 1; z) = ∞ X k=0 a zk (a + k) (|z| < 1, a ∈ C), (61) where 2 F1 denotes the Gauss hypergeometric function; we can thus rewrite the expansion (60) as   1 3 it + 1 1 , (62) H(t, y) = 2 F1 1, iy + , iy + ; (1 − it)(iy + 1/2) 2 2 it − 1 which shows that H(t, y) is analytic in the half–plane Im t > 0. By the use of the following relationship, which holds true for the Gauss hypergeometric function [1]:     z z −a |z| < 1, < 1 , (63) 2 F1 (a, b, c; z) = (1−z) 2 F1 a, c − b, c; z−1 z−1 we can rewrite the r.h.s. of formula (62) as   1 3 1 + it H(t, y) = 1, 1, iy + F . ; 2 1 2 2 2(iy + 12 ) (64) The hypergeometric series at the r.h.s. of formula (64) converges inside a circle of radius 2 and center t = i, belonging to the complex t-plane. We can then calcu1 (dm H(t, y)/dtm )t=0 . By the use of the standard late the following derivatives: m! formula for the derivative of the hypergeometric function we obtain [1]   1 dm H(t, y) im Γ(m + 1)Γ(iy + 21 ) = m! dtm 2m+1 Γ(m + iy + 32 ) t=0 (65)   3 1 . × 2 F1 m + 1, m + 1, m + iy + ; 2 2 Next, returning to the integral representation of H(t, y) given by formula (57), we exchange the integral with the sum. This is legitimate for Re t = 0, |t| < 1, since 1−it ix ) | = | exp[ix ln 1−it |( 1+it 1+it ]| = 1, (x ∈ R); then we write, for Re t = 0, |t| < 1 and Im y < 0:  X  Z +∞ ∞ ∞ X Pm (x)w(x) i m dx = tm Q(1) (66) H(t, y) = t m (y), 2 x − y −∞ m=0 m=0 and therefore  dm H 1 (t, y) im Γ(m + 1)Γ(iy + 21 ) = dtm 2m+1 Γ(m + iy + 23 ) t=0   3 1 . × 2 F1 m + 1, m + 1, m + iy + ; 2 2 Q(1) m (y) = 1 m!  (67) 18 E. DE MICHELI AND G. A. VIANO Let us now observe that the function 2 F1 (a, b, c; z)/Γ(c) is an entire function of (1) the variable c; therefore Qm (y) can be analytically continued in the complex yplane, through formula (67) and the only singularities that it presents in this (open) complex plane are the simple poles of the gamma function Γ(iy + 1/2): i.e. y = i(N + 1/2), (N = 0, 1, 2, . . .). Now we can employ the following relationship which holds true for the hypergeometric function: 2 F1 (a, b, c; z) and we arrive at = (1 − z)c−a−b 2 F1 (c − a, c − b, c; z) (|z| < 1), (68)   Γ(m + 1)Γ(iy + 12 ) 1 1 3 1 iy + F . , iy + , m + iy + ; 2 1 2 2 2 2 Γ(m + iy + 32 ) (69) Finally, recalling the asymptotic behavior of 2 F1 (a, b, c; z) for large |c|, at fixed a, b, z (see [1, formula 10, p. 76]), and the asymptotic behavior of the gamma function, we obtain   1 m (2m)−iy−1/2 . (70) Q(1) (y) ∼ i Γ iy + m 2 m −iy−1/2 Q(1) m (y) = i 2 (2) Let us now focus our attention on the function Qm (y), defined for Im y > 0. If we change x into −x, by observing that Pm (−x) = (−1)m Pm (x), (71a) w(−x) = w(x), (71b) we obtain m+1 Q(2) m (y) = (−1) i 2 Z +∞ −∞ Pm (x)w(x) dx = (−1)m+1 Q(1) m (−y). x+y (72) This latter relationship, proved for Im y > 0, can be extended by analytic continuation. Next, by the use of the Cauchy formula we obtain I Pm (x)w(x) i (1) dx = −πw(y)Pm (y), (73) Q(2) (y) − Q (y) = m m 2 γ2 x−y where γ2 is a path which encircles the singularity located at x = y in a counterclockwise sense. From (72) and (73) we finally obtain m (1) πw(y)Pm (y) = Q(1) m (y) + (−1) Qm (−y). From the asymptotic behavior (70) and by the use of (74) we obtain:    im 1 Pm (y) ∼ + iy (2m)−iy−1/2 Γ 2 Γ( 21 + iy)Γ( 12 − iy)    1 + (−1)m Γ − iy (2m)iy−1/2 . 2 Finally, by putting y = −i(n + 1/2), we get:      1 1 1 m −n m n ∼i (2m) + (−1) (2m) Pm −i n + 2 Γ(−n) Γ(n + 1) im (−1)m = (2m)n , n! which is exactly the result we want to prove. (74) (75) (76)  CAUCHY INTEGRAL EQUATIONS 19 References [1] Bateman Manuscript Project, “Higher Trascendental Functions” (A. Erdelyi, Director), Vol. 1, McGraw–Hill, New York, 1953. [2] Bateman Manuscript Project, “Higher Trascendental Functions” (A. Erdelyi, Director), Vol. 2, McGraw–Hill, New York, 1953. [3] E. De Micheli and G. A. Viano, Hausdorff moments, Hardy spaces and power series, J. Math. Anal. Appl. 234 (1999), 265–286. [4] C. W. Groetsch, “The Theory of Tikhonov Regularization for Fredholm Equations of the First Kind,” Pitman, Boston, 1984. [5] F. John, Continuous dependence on data for solutions of partial differential equations with a prescribed bound, Comm. Pure Appl. Math. 13 (1960), 551–585. [6] M. E. H. Ismail, Asymptotics of Pollaczek polynomials and their zeros, SIAM J. Math. Anal. 25 (1994), 462–473. [7] N. Magnoli and G. A. Viano, The source identification problem in electromagnetic theory, J. Math. Phys. 38 (1997), 2366–2388. [8] K. Miller and G. A. Viano, On the necessity of nearly–best–possible methods for analytic continuation of scattering data, J. Math. Phys. 14 (1973), 1037–1048. [9] E. M. Stein and S. Wainger, Analytic properties of expansions, and some variants of the Parseval–Plancherel formulas, Ark. Mat. 37 (1965), 553–567. [10] G. Szegö, “Orthogonal Polynomials,” Amer. Math. Soc., Providence, 1959. [11] A. Tikhonov and V. Arsenine, “Mèthodes de Résolution de Probléms Mal Posès,” Mir, Moscow, 1976. [12] G. A. Viano, On the regularization of the antenna synthesis problem, in “Partial Differential Equations and Applications” (P. Marcellini, G. T. Talenti, and E. Visentini, Eds.), Dekker, New York, 1996. (E. De Micheli) IBF – Consiglio Nazionale delle Ricerche, Via De Marini, 6 - 16149 Genova, Italy E-mail address, E. De Micheli: demicheli@ge.cnr.it (G. A. Viano) Dipartimento di Fisica - Università di Genova, Istituto Nazionale di Fisica Nucleare - sez. di Genova, Via Dodecaneso, 33 - 16146 Genova, Italy E-mail address, G.A. Viano: viano@ge.infn.it