Документ взят из кэша поисковой машины. Адрес оригинального документа : http://imaging.cs.msu.ru/pub/2011.ITSF.Krylov_Mizotin.Laguerre.en.pdf Дата изменения: Mon Jul 4 19:20:44 2011 Дата индексирования: Mon Oct 1 19:49:05 2012 Кодировка: IBM-866

Pro jection Method for the Hankel Transform
Andrey S. Krylov Maxim M. Mizotin 11 2011 .
Numerical pro jection method for the zero-order Hankel transform inversion for the case of the data given on a finite interval has been developed. The justification of the convergence of the pro jection method has been done for the general case that also includes the sine-Fourier transform inversion. An inequality on the norms of the zeroth order Laguerre functions has been proved. Efficiency of the method was illustrated with the test data and for the problem of cylindrical distribution function calculation using melt surface layer diffraction data.

1

Introduction

The Hankel transform is widely used in optics [1, 7], hydrodynamics [2, 9], acoustics [3], image processing [16], etc. The Hankel transform is often used to solve axially symmetric problems in Fourier analysis. In this case, the problem dimensionality can be lowered and the zero-order Hankel transform is used: H [z (t)] = z (x)J0 (xt)x dx.
0

Usually the inversion of the Hankel transform is considered for the case H [z ] = u, H : L2 [0, ) L2 [0, ).

Nevertheless, in practical problems the right part u of this equation (the experimental data) is known only on a finite interval. Thus, instead of the original well-posed inverse problem on a whole half-line we solve an ill-posed inverse problem for the equation H [z ] = u, H : L2 [0, ) L2 [0, a]. (1) Usually the finite Hankel transform is used to solve this problem [15]. Here we have two difficulties. First, the inversion of the Hankel transform on a finite interval is an ill-posed problem. Second, we must assume that the right part and the solution are finite in this case. It is fulfilled only for zero function if we still consider the problem as the Hankel transform inversion problem. Moreover, due to the ill-posedness of the problem, known analytical formulae of inversion do not provide satisfactory accuracy [5]. An alternative way is to expand the given function into a series of orthonormal eigenfunctions of the Hankel transform. These functions are in fact the Laguerre 1


functions. Once the function expansion is known, one can easily evaluate its Hankel transform using the obtained expansion coefficients. However, numerical study of the method revealed some details which should be taken in account [18]. The Laguerre functions appreciably differ from zero only on a finite interval (we call it as "numerically concentrated"). The experiments showed that the expansion convergence slows down significantly when this interval goes beyond the region where the data is given. The authors propose an estimate for the "optimal" number of terms in expansion depending on the scale of the Laguerre functions and on the length of the given data interval. In [14] the numerical concentration property was used in conjunction with the a priori knowledge of numerical concentration interval of the solution. The scale for the Laguerre functions was chosen taking into account the concentration domains of both the data and the solution. This allowed to solve the problem without additional proposition on the zero continuation of the right part of the equation (1). In this paper we propose computational pro jection method for the case of noisy data given on a finite interval. The convergence of the method is theoretically proved and the optimal number of functions is determined, depending on the interval length and the noise level.

2

The general problem statement and the related works

Consider the general problem of inversion of the integral transform A on a halfline Az = u, A : L2 [0, ) L2 [0, ). (2) Usually, the exact right part u L2 [0, ) is unknown and only the function u and the error norm are given. Thus the inverse problem is defined by the a equations Az = u , A : L2 [0, ) L2 [0, a], (3) a u L2 [0, a] : u - u a a
L2 [0, a]

.

(4)

The problem (3) is ill-posed because its solution is not unique. The Sine-Fourier transform operator in the case when data is given on a finite interval is the most frequently considered problem of this type [4, 8, 19]: z (y ) sin(xy ) dy = u (x), S : L2 [0, ) L2 [0, a]. (5) Sz = a
0

There are two ma jor groups of methods used to solve this problems. The first one uses the zero continuation of u to the half-line (we call it "the standard a method"). The second one reduces the original problem to the Fredholm equation of the first kind by changing the integration limit to a finite value and applies regularization techniques. The ma jor drawback of the first method is so called "tear effect", while the second group does not use properties of the transform on the whole half-line, thus losing important information about the initial wellposed problem. In [11, 13] the computational pro jection method was suggested, which combines the regularization technique and the use of the integral transform properties. It 2


also takes into account the length of the interval where the data is given. This method was successfully applied for solution of the problem (5). In this work the computational pro jection method theory is further developed. The new general convergence theorem is proved. It generalizes the result of [11] for the case of noise addition (4) to the right part of (3). The pro jection method is applied to the numerical solution of the problem (1).

3

Computational pro jection metho d of the integral transform inversion

We consider the inverse problem (2). Let the operator A has a full, orthonormal in L2 [0, ) system of eigenfunctions {i }, and the eigenvalues {i } are finite and isolated from zero: > |i | > > 0. (6) Then the exact right part can be expanded in a Fourier series using the system of eigenfunctions u(x) = ci i (x), ci = u(x)i (x) dx,
i=0 0

and the exact solution z is given by: z (y ) =
i=0

di i (y ),

di =

ci . i

To solve numerically the problem (3), (4) we construct the pro jection method using the system of eigenfunctions of the well-posed problem (2). Denote
za (y ) = N -1 i=0 a

d

a, i

i (y ), i = 0, N - 1,

(7) (8)

d

a, i

1 ca, =i = i i


0

u (x)i (x) dx, a

where N is the number of terms in the Fourier series expansion. One of the key points of the pro jection method is the choice of the number N . In most cases, when no a priori information about function za (x) is available, the number N is taken as the minimal value, that gives an approximation error of u (x) lower than . But this choice of N has a drawback because it does not take a into account the value of a. When is small we need a large N to approximate u (x), and if a is finite, the coefficients ci and ca, can differ significantly for a i large i. Thus we find N as the maximal value of n that satisfies the condition i L
2

[a, )

K1 , n

i = 0, n - 1,

and K1 > 0 is a constant.

(9)

It selects functions i (x), that are numerically concentrated on the interval [0, a]. In practice it is important to take into account simultaneously the measurement error value and the interval length a to choose N . It can be done using the following theorem: 3


Theorem 3.1. Let eigenvalues i of the ful l, orthonormal in L2 [0, ) system of eigenfunctions {i } of the operator A (3) satisfy the condition (6) and the eigenfunctions are ordered so that: i L Let
2

[a, )



i+j L2 [a, )

,

i = 0, , j = k , , a > 0, k N. } ,

(10)

{ Na = max n N : n-1 (x)L2 [a, { N = min n N : u (x) - u (x) ~a a

)

K1 n

L2 [0,a]

K1 (0, 1], } K2 , K2 > 0,

(11) (12)

where u = Aza . If ~a

N = min {N , Na } , then
za

(13) (14) ) for i < Na . (15)

-z (

L2 [0,)

- 0,
a [a,)
0

and

da, i

1 - di |i |

K1 + uL Na

2

older inequality . The estimate (15) can be obtained using Hи for the right part of da, i 1 - di = |i | a u (x)i (x) dx -
a 0 0

u(x)i (x) dx

a 1 u(x)i (x) dx . (u (x) - u(x))i (x) dx + a |i |
0 a

Then, da, - di i 1( u(x) - u (x) a |i | + u(x)L2 [a,
L2 [0, a] )

i (x)L2

[0, a]

+ ) .

i (x)L2

( 1 + u(x)L |i |
[a, )

)

2

[a, )

K 1 Na

Let

{ u (x), x [0, a], a u (x) = u(x), x > a,

c i =
0



u (x)i dx

and z is the corresponding solution of the equation (2). Consider the approximation error: za - z L2 [0, ) za - z L2 [0, ) + z - z L2 [0, ) . The second term converges to zero when 0. It is necessary to show that za - z L2 [0, ) also converges to zero when 0 and a .

4


If N = N , then
za - z 2 L2 [0, ) 2 L2 [



+ u - u ~a The value
za

1 1( 2 2 u - u L2 [0, ) = 2 u - u L2 [0, a] + ~a ~a a 2 ) 1( 2 2 2 (K2 )2 + 2 u L [a, ) + 2 uL2 [ ~a a, ) 2 + uL
2
2

)
a, )

.

u ~

- z L2 [0, ) 0 when 0 and a . Consider the case N = Na . Let us estimate the first term:
za - z L2 [0, )

2 a L2 [a, )

[a, )

0 when 0 and a . Thus,

1

N -1 1 a, (ci - c )i - c i i i=0 i=N N -1 i=0

i L2 [0, )





(ca, - c ) i i

i L2 [0, )

1 + ci i i=N

. (16)
L2 [0, )

{ } It follows from (10) that max n-1 L2 [a, ) , n L2 [a, ) is ascendant function of n as (K1 / n) 0 while n . This guarantees the existence of the finite Na in (11) for the sufficiently large a. Now we show that the sequence Na does not have the upper bound when a . Suppose that this is not the case. Thus N : Because i L2
[a, )

a a > a :

Na < N .

0 when a and i is fixed, then for ai : a > ai i L2
[a, )

i =

K

1

N +1

K1 . N +1

For any a max{a0 , a1 , . . . , aN } we have i L2
[a, )

K1 N +1

i = 0, 1, . . . , N .

Consequently Na > N , and this contradicts to our assumption. Thus N when 0 and a , and the second term in (16) tends to zero. Consider the first term: N -1 N -1 1 a, 1 (c - c )i u ( )i ( ) d i = = i i=0 i i=0 2 1/2 [N -1 N -1 1 ( 1 u ( )i ( ) d uL2 = i=0 i=0
a



L2 [0, )

a

L2 [0, )

[a, )

i L2

[a, )

)

]
2

1/2



uL

2

[a, )

[( )2 ]1/2 1 c1 1 N = c1 uL N

2

[a, )

a

- 0.

Hence, the theorem is proved. 5


The theorem 3.1 enables us to construct the following algorithm for numerical solution of the problem (3)-(4): 1. Choose the number of the Fourier series terms N according to the formulae (11)(13). 2. Compute the Fourier coefficients (8) ca, , i = 0, N for the given function i u . a 3. Compute the solution as a partial sum of the Fourier series with the coefficients da, (7). i

4

Structure analysis applications of the pro jection metho d

Diffraction methods to study the structure of matter with a monochromatic coherent rays (X-rays, electron beams, etc.) measure the diffraction scattering on the ob ject's atoms. At present time various beam types and experimental setups allow to study the structure of both crystalline and non-crystalline solids and liquids. It is also possible to extract the structure of the surface layers. Under the kinematic approximation of the scattering theory the diffraction pattern acquired by the scattering of coherent monochromatic beam on the potential p(r) is given by [20]: f (s) = K p(r)ei(sr) dr, where coefficient K depends on the beam type. The measured function f (s) is often called as structure factor. In the inner regions (that are apart of the surface) of liquid and amorphous systems the scattering potential is radially symmetric, while near the surface the potential is cylindrically symmetric. The corresponding scattering potentials are called radial and cylindrical atomic distribution function respectively. Due to the symmetry of atom distribution functions in the case of the inner regions we obtain the problem of inversion of the sine-Fourier transform and for the surface layers we have the problem of the Hankel transform inversion. In both cases there are experimental constraints, thus the structure factor can be only measured on a finite interval and it is measured with an experimental error. So we come to the problem statements (5) and (1) respectively.

5

The Sine-Fourier transform

The pro jection method for the inversion of the sine-Fourier transform (5) with the data given on a finite interval was developed in [11]. The eigenfunctions of the sine-Fourier transform are the odd Hermite functions n , where n = e
-x2 /2

h2n+1 (x),

hn (x) = (-1)n e

-x2

d n e- x are the Hermite polynomials. dxn
2

6


Using the equivalence for the odd Hermite function norms n L2
2 [0, R]

= 2n

2 n-1 L2 [0, R] 2

- n (R)

n-1 2

(R),

it was shown in [12] that norms n L2 [0, R] and n+1 L2 [0, R] as functions of R intersects even for R > 1 (see fig. 1). Nevertheless, the theorem 3.1 can be used as the condition (10) is satisfied for k = 2.
1 0.9 0.8
2 n L2 [0,R]

0

2

1

2

0.7 0.6 0.5 0.4 0.3 0.2 0.1 0 0 1 2 3 R 4 5 6



(a)
0.302 0.301
2 n L2 [0,R]

0.3 0.299 0.298 0.297 2.26 2.265 2.27 R 2.275
2 5 L2 [0,R] 2 6 L2 [0,R]

2.255

2.28

2.285

(b)

. 1: (a) Norms of the first 15 odd Hermite functions on interval [0, (b) Intersection of 5-th and 6-th odd Hermite function L2 [0, R] norms

R];

6

Zero-order Hankel transform

The justification of the pro jection method for the zero-order Hankel transform inversion was not done before. To solve the problem (1) by the pro jection method we need to prove that eigenfunctions of the Hankel transform meet the requirements of the theorem 3.1. The zeroth order Laguerre functions are defined as n (x) = e where Ln = 1 x dn ( n e xe n! dxn
-x -x/2

Ln (x),

)

are the zeroth order Laguerre polynomials [6].

7


~ The Laguerre functions n (x) are the eigenfunctions of the transform H ~ H [z ] ( t) = z (x)J0 (xt) xt dx,
0

which is associated with the Hankel transform ] 1 ~[ H [z ] (t) = H z (x) x . t The eigenfunctions of the Hankel transform are [17]: n (x) = 2x n (x2 ), H [n ] (t) = (-1)n n (t). The examples of the eigenfunctions n (x) for n = 0, 5, and 10 are given in fig. 2. This system of functions is full and orthonormal in L2 [0, ) and the condition (6) is fulfilled. In this article we prove that n (x) can be ordered according to the condition (10).
1 0.75 0.5 n ( x) 0.25 0 -0.25 -0.5 0 1 2 3 4 x 5 6 7 8 0 ( x) 5 ( x) 10 ( x)

(17)

. 2: 0-th (solid line), 5-th (dashed line), 10-th (dash-dotted line) eigenfunctions of the Hankel transform

6.1

Inequality for the norms of the Hankel transform eigenfunctions

Theorem 6.1. Norms of the eigenfunctions (17) of the zero-order Hankel transform satisfy the inequality n L
2

[R, )

n

+1 L2 [R, )

,

n N, R > 0.

(18)

. The norms of n (x) and the norms of n (x) are connected as n (x)
2 L2 [0, R2 ] R 0
2


2 2xn

R 0

=

(x ) dx =

2

2 n (y ) dy = n (x)L

2

2

[0, R]

,

8


1 0.9 0.8
2 L2 [0,R]

0.7 0.6 0.5 0.4 0.3 0.2 0.1 0 0 2 4 6 8 10 R 12 14 16 18 20

. 3: Norms of the first 15 Laguerre functions (x) on the interval [0, R]

and (18) can be reduced to n (x)
2 L2 [0, R]

( x)

-

n+1

(x)

2 L2 [0, R]

R 2 n (x) - 0 2 n+1

=

(x) dx 0.

(19)

Consider the derivative of the expression d[ e dx
-x

(Ln (x) - Ln

+1

(x))

2

]

= 2e

-x

(

) )( Ln (x) - Ln+1 (x) L (x) - L +1 (x) - n n -e
-x

(Ln (x) - Ln+1 (x))2 . (20)

Using the relation [6]: L n
+1

(x) - L (x) + Ln (x) = 0, n

the expression (20) can be rewritten as : )2 ] ( d [ -x ( e Ln (x) - Ln+1 (x) = e-x 2L2 (x) - 2Ln (x)Ln+1 (x) - L2 (x) + n n dx ) 2 -x 2 2 2 2 + 2Ln (x)Ln+1 (x) - Ln+1 = e (Ln (x) - Ln+1 (x)) = n (x) - n+1 (x). Thus, we have:
0 R 2 n (x) - 2 n+1

(x) dx = e

-x

(

Ln (x)-Ln

+1

) (x)

2

R

=e
0

-R

(

)2 Ln (R)-Ln+1 (R) .

It proves the inequality (19) and the theorem. Moreover, the equality is possible ^ ^ ^ only at the points R, where Ln (R) = Ln+1 (R). The behavior of the norms of the zeroth order Laguerre functions is illustrated in fig. 3

9


2 1.5 1 z (t) 0.5 0 -0.5 Exact solution Standard method Projection method

0

5 t

10

15

. 4: Reconstruction of model function from u0,5 10

0.6 0.4 (r) 0.2 0 -0.2 -0.4 0 2 4 6 r Projection method Standard method 8 10 12 14

. 5: Reconstruction of model function z (t) =

1 -t2 /4a2 2a2 e

from u0,01 10

7
7.1

Numerical exp eriments
Mo del function

The pro jection method gives better results comparing with the "standard method" of the zero-order Hankel transform inversion (zero continuation of the data to the half-line). It can be illustrated with the model function { 1, t 2.5; z (t) = H [u] = 0, t > 2.5. Its zero-order Hankel transform is u(x) = 2.5 J1 (2.5x) . x To model the real data situation, we added uniformly distributed noise with = 0.5 to the function u(x), given on a uniform grid with a step 0.02. The 10


result for the cutoff value a = 14 is shown in fig. 4. The obtained discrepancies for the right part of the equation (3) in this case were 0.89 for the standard method and 0.32 for pro jectional method. The number of the functions in the pro jection method was chosen according to the interval length and the error value. It gave N = 51 for u0.5 . It can be 14 seen that the use of the "standard method" gave a more oscillating solution comparing with the pro jection method. 2 2 1 The function z (t) = 22 e-t /4 with the exact zero-order Hankel transform 22 u(x) = e- x is an illustration of the case of continuous function restoration. In this case we did not have the Gibbs phenomenon so the advantage of the pro jection method was more clear (see fig. 5). Here we took = 2, = 0.01, grid step 0.02 and the cutoff value a = 10.

7.2

Computation of cylindrical atomic distribution function

The advantage of the pro jection method comparing with "standard method" was also shown for the practical problem of computation of the cylindrical atomic distribution function of liquid germanium surface from measured electron diffraction data [21]. The source data was given on the interval [0, 9.8]. The problem H [ ] (s) = a (s), (21) was solved using the Hankel transform inversion formula with the zero continuation of the experimental data ("standard method") (r) and using the pro jection inv method (7)-(8) oj (r). In this experiment the noise level was small enough so pr that N > Na = 34. The solutions are shown in the fig. 6. To estimate the quality of the acquired solutions, the experimental structure factor a (s) was recalculated from (r) and oj (r) using (21), and the pr inv obtained residual norms in L2 were [ ] a (s) - H [ ] (s) = 3.49, a (s) - H oj (s) = 0.40. inv pr

8

Conclusion

It was illustrated with the zero-order Hankel transform that the pro jection method is effective for numerical inversion of a wide class of the Fourier integral transforms with the data given on a finite interval. The pro jection method enables one to take into consideration the length of the given data interval and the value of the data error. It looks promising to construct fast implementations of pro jection methods by analogy with the fast Hermite pro jection method [10].

Acknowledgements
The work was supported by the federal target program "Scientific and scientificpedagogical personnel of innovative Russia" in 2009-2013 and RFBR grant 1001-00535.

11


1 Projection method Standard method 0.5 (r)

0

-0.5

-1

0

5 r

10

15

(a)
Projection method Standard method

0.06 0.04 0.02 (r) 0 -0.02 -0.04 -0.06 -0.08 8 9 10 r

11

12

(b)

. 6: (a) Cylindrical atomic distribution function (r) of liquid germanium, calculated by the "standard method" and pro jection method; (b) Zoomed fragment

12



[1] Li Yu, M. Huang, M. Chen, W. Huang, and Z. Zhu. Quasi-discrete Hankel transform. Optics Letters, 23(6):409н411, 1998. [2] M. El-Shahed and M. Shawkey. Generalized finite Hankel transform. Integral Transforms and Special Functions, 17:39н44, 2006. [3] D. R. Mook. An algorithm for numerical evaluation of Hankel and Abel transform. IEEE Transactions on Acoustics, Speech and Signal Processing, 31(4):979н985, 1983. [4] V. F. Uhov, N. A. Vatolin, B. R. Gel'chinskiy, V.P. Beskachko, and O. A. Esin. Interparticle interaction in melted metals. Nauka, Moscow, 1979. (In Russian). [5] P. C. Shah and R. K. M. Thambynayagam. Application of the finite Hankel transform to a diffusion problem without azimuthal symmetry. Transport in Porous Media, 14(3):247н264, 2003. [6] G. Szogи Orthogonal Polynomials. American Mathematical Society, New o. York, 1939. [7] V. Magni, G. Cerullo, and De Silvestri. High-accuracy fast Hankel transform for optical beam propagation. J. Opt. Soc. Am. A, 9(11):2031н 2033, 1992. [8] V. V. Vasin and A. N. Ageev. Il l-posed problems with a priori information. Nauka, Ekaterinburg, 1993. (In Russian). [9] N. T. Eldabe, M. El-Shahed, and M. Shawkey. An extension of the finite Hankel transform. Applied Mathematics and Computation, 151:713нн717, 2004. [10] A. S. Krylov and D. N Korchagin. Fast Hermite Pro jection Method. Lecture Notes in Computer Science, 4141:329н338, 2006. [11] A. S. Krylov and A. V. Liakishev. Numerical Pro jection Method For Inverse Fourier Transform and its Application. Numerical Functional Analysis and optimization, 21:205н216, 2000. [12] A. S. Krylov and A. V. Liakishev. Inequality for Hermite function norms on finite interval. Vestnik Moskovskogo Universiteta, 1:17н19, 1999. (In Russian). [13] A.S. Krylov and A.V. Vvedenskii. Software Package for Radial Distribution Function Calculation. Journal of Non-Crystal line Solids, 192:683н687, 1995. [14] M. Mizotin, A. Krylov, and M. Spiridonov. Application of pro jection method for analysis of structural dependencies in surface layers of melts. Rasplavy, 5:18н26, 2009. [15] I. N. Sneddon. On finite Hankel transforms. Philosophical Magazine, 7:17н 25, 1946. 13


[16] W. E. Higgins and D. C. Munson. A Hankel transform approach to tomographic image reconstruction. IEEE Transactions on Medical Imaging, 7(1):59н72, 1988. [17] Arthur Erdelyi, editor. Tables of Integral Transforms, volume 2. McGrawHill, New York, 1954. [18] E. C. Cavanagh and B. D. Cook. Numerical evaluation of Hankel Transform via Gaussian-Laguerre polynomial expressions. IEEE Transactions on Acoustics, Speech and Signal Processing, 27:361н366, 1979. [19] G. M. Zhidomirov, editor. X-ray spectral method of studying amorphous objects structure: EXAFS spectroscopy. Nauka, Novosibirsk, 1988. (In Russian). [20] B. K. Veinshtein. Structural electronography. Academy of sciences, Moscow, 1956. (In Russian). [21] M. A. Spiridonov, A.V. Lavrov, and S. I. Popel. Atomic ordering in surface layers of cuprum-germanium melts. Metals, 2(3):49н53, 1990. (In Russian).

14