Bessel.Frobeinus astronomy

Notes on Bessel’s Equation and the Gamma Function Charles Byrne (Charles [email protected]) Department of Mathematical Sciences University of Massachusetts at Lowell Lowell, MA 01854, USA April 8, 2009 1 Bessel’s Equations For each non-negative constant p, the associated Bessel Equation is x2d2y dx2+xdy dx+ (x2−p2)y= 0, (1.1) which can also be written in the form y/prime/prime+P(x)y/prime+Q(x)y= 0, (1.2) with P(x) =1 xandQ(x) = 1−p2 x2. Solutions of Equation (1.1) are Bessel functions . These functions first arose in Daniel Bernoulli’s study of the oscillations of a hanging chain, and now play important roles in many areas of applied mathematics [1]. We begin this note with Bernoulli’s problem, to see how Bessel’s Equation be- comes involved. We then consider Frobenius-series solutions to second-order linear differential equations with regular singular points; Bessel’s Equation is one of these. Once we obtain the Frobenius-series solution of Equation (1.1), we discover that it involves terms of the form p!, for (possibly) non-integer p. This leads to the Gamma Function , which extends the factorial function to such non-integer arguments. The Gamma Function, defined for x >0 by the integral Γ(x) =/integraldisplay∞ 0e−ttx−1dt, (1.3) is a higher transcendental function that cannot be evaluated by purely algebraic means, and can only be approximated by numerical techniques. With clever changes 1 of variable, a large number of challenging integration problems can be rewritten and solved in terms of the gamma function. We prepare for our discussion of Bernoulli’s hanging chain problem by recalling some important points in the derivation of the one-dimensional wave equation for the vibrating string problem. 2 The Vibrating String Problem In the vibrating string problem, the string is fixed at end-points (0 ,0) and (1 ,0). The position of the string at time tis given by y(x, t), where xis the horizontal spatial variable. It is assumed that the string has a constant mass density, m. Consider the small piece of the string corresponding to the interval [ x, x+ ∆x]. Its mass is m∆x, and so, from Newton’s equating of force with mass times acceleration, we have that the force fon the small piece of string is related to acceleration by f=m(∆x)∂2y ∂t2. (2.1) In this problem, the force is not gravitational, but comes from the tension applied to the string; we denote by T(x) the tension in the string at x. This tensile force acts along the tangent to the string at every point. Therefore, the force acting on the left end-point of the small piece is directed to the left and is given by −T(x) sin( θ(x)); at the right end-point it is T(x+ ∆x) sin( θ(x+ ∆x)), where θ(x) is the angle the tangent line at xmakes with the horizontal. For small-amplitude oscillations of the string, the angles are near zero and the sine can be replaced by the tangent. Since tan(θ(x)) =∂y ∂x(x), we can write the net force on the small piece of string as f=T(x+ ∆x)∂y ∂x(x+ ∆x)−T(x)∂y ∂x(x). (2.2) Equating the two expressions for fin Equations (2.1) and (2.2) and dividing by ∆ x, we obtain T(x+ ∆x)∂y ∂x(x+ ∆x)−T(x)∂y ∂x(x) ∆x=m∂2y ∂t2. (2.3) Taking limits, as ∆ x→0, we arrive at the Wave Equation ∂ ∂x/parenleftBig T(x)∂y ∂x(x)/parenrightBig =m∂2y ∂t2. (2.4) For the vibrating string problem, we also assume that the tension function is constant, that is, T(x) =T, for all x. Then we can write Equation (2.4) as the more familiar T∂2y ∂x2=m∂2y ∂t2. (2.5) 2 We could have introduced the assumption of constant tension earlier in this dis- cussion, but we shall need the wave equation for variable tension Equation (2.4) when we consider the hanging chain problem. 3 The Hanging Chain Problem Imagine a flexible chain hanging vertically. Assume that the chain has a constant mass density m. Let the origin (0 ,0) be the bottom of the chain, with the positive x- axis running vertically, up through the chain. The positive y-axis extends horizontally to the left, from the bottom of the chain. As before, the function y(x, t) denotes the position of each point on the chain at time t. We are interested in the oscillation of the hanging chain. This is the vibrating string problem turned on its side, except that now the tension is not constant. 3.1 The Wave Equation for the Hanging Chain The tension at the point xalong the chain is due to the weight of the portion of the chain below the point x, which is then T(x) =mgx. Applying Equation (2.4), we have ∂ ∂x/parenleftBig mgx∂y ∂x(x)/parenrightBig =m∂2y ∂t2. (3.1) As we normally do at this stage, we separate the variables, to find potential solutions. 3.2 Separating the Variables We consider possible solutions having the form y(x, t) =u(x)v(t). (3.2) Inserting this y(x, t) into Equation (3.1), and doing a bit of algebra, we arrive at gxu/prime/prime(x) +gu/prime(x) +λu(x) = 0 , (3.3) and v/prime/prime(t) +λv(t) = 0 , (3.4) where λis the separation constant. It is Equation (3.3), which can also be written as d dx(gxu/prime(x)) +λu(x) = 0 , (3.5) that interests us here. 3 3.3 Obtaining Bessel’s Equation With a bit more work, using the change of variable z= 2/radicalBig λ g√xand the Chain Rule (no pun intended!), we find that we can rewrite Equation (3.3) as z2d2u dz2+zdu dz+ (z2−02)u= 0, (3.6) which is Bessel’s Equation (1.1), with the parameter value p= 0. 4 Solving Bessel’s Equations Second-order linear differential equations with the form y/prime/prime(x) +P(x)y/prime(x) +Q(x)y(x) = 0 , (4.1) with neither P(x) nor Q(x) analytic at x=x0, but with both ( x−x0)P(x) and (x−x0)2Q(x) analytic, are said to be equations with regular singular points . Writing Equation (1.1) as y/prime/prime(x) +1 xy/prime(x) + (1−p2 x2)y(x) = 0 , (4.2) we see that Bessel’s Equation is such a regular singular point equation, with the singular point x0= 0. Solutions to such equations can be found using the technique of Frobenius series. 4.1 Frobenius-series solutions A Frobenius series associated with the singular point x0= 0 has the form y(x) =xm(a0+a1x+a2x2+...), (4.3) where mis to be determined, and a0/negationslash= 0. Since xP(x) and x2Q(x) are analytic, we can write xP(x) =p0+p1x+p2x2+..., (4.4) and x2Q(x) =q0+q1x+q2x2+..., (4.5) with convergence for |x|< R. Inserting these expressions into the differential equa- tion, and performing a bit of algebra, we arrive at ∞/summationdisplay n=0/braceleftBigg an[(m+n)(m+n−1) + ( m+n)p0+q0] +n−1/summationdisplay k=0ak[(m+k)pn−k+qn−k]/bracerightBigg xn= 0. 4 (4.6) Setting each coefficient to zero, we obtain a recursive algorithm for finding the an. To start with, we have a0[m(m−1) +mp0+q0] = 0. (4.7) Since a0/negationslash= 0, we must have m(m−1) +mp0+q0= 0; (4.8) this is called the Indicial Equation . We solve the quadratic Equation (4.8) for m=m1 andm=m2. 4.2 Bessel Functions Applying these results to Bessel’s Equation (1.1), we see that P(x) =1 x,Q(x) = 1−p2 x2, and so p0= 1 and q0=−p2. The Indicial Equation (4.8) is now m2−p2= 0, (4.9) with solutions m1=p, and m2=−p. The recursive algorithm for finding the anis an=−an−2/n(2p+n). (4.10) Since a0/negationslash= 0 and a−1= 0, it follows that the solution for m=pis y=a0xp/bracketleftBigg 1−x2 22(p+ 1)+x4 242!(p+ 1)( p+ 2)−.../bracketrightBigg . (4.11) Setting a0= 1/2pp!, we get the pth Bessel function, Jp(x) =∞/summationdisplay n=0(−1)n/parenleftBigx 2/parenrightBig2n+p/n!(p+n)!. (4.12) The most important Bessel functions are J0(x) and J1(x). We have a Problem! So far, we have allowed pto be any real number. What, then, do we mean by p! and ( n+p)!? To answer this question, we need to investigate the gamma function. 5 The Gamma Function We want to define p! for pnot a non-negative integer. The Gamma Function is the way to do this. 5 5.1 Extending the Factorial Function As we said earlier, the Gamma Function is defined for x >0 by Γ(x) =/integraldisplay∞ 0e−ttx−1dt. (5.1) Using integration by parts, it is easy to show that Γ(x+ 1) = xΓ(x). (5.2) Using Equation (5.2) and the fact that Γ(1) =/integraldisplay∞ 0e−tdt= 1, (5.3) we obtain Γ(n+ 1) = n!, (5.4) forn= 0,1,2, .... 5.2 Extending Γ(x)to negative x We can use Γ(x) =Γ(x+ 1) x(5.5) to extend Γ( x) to any x <0, with the exception of the non-negative integers, at which Γ(x) is unbounded. 5.3 An Example We have Γ(1 2) =/integraldisplay∞ 0e−tt−1/2dt. (5.6) Therefore, using t=u2, we have Γ(1 2) = 2/integraldisplay∞ 0e−u2du. (5.7) Squaring, we get Γ(1 2)2= 4/integraldisplay∞ 0/integraldisplay∞ 0e−u2e−v2dudv. (5.8) 6 In polar coordinates, this becomes Γ(1 2)2= 4/integraldisplayπ 2 0/integraldisplay∞ 0e−r2rdrdθ = 2/integraldisplayπ 2 01dθ=π. (5.9) Consequently, we have Γ(1 2) =√π. (5.10) 6 An Application of the Bessel Functions in As- tronomy In remote sensing applications, it is often the case that what we measure is the Fourier transform of what we really want. This is the case in medical imaging, for example, in both x-ray tomography and magnetic-resonance imaging. It is also often the case in astronomy. Consider the problem of determining the size of a distant star. We model the star as a distance disk of uniform brightness. Viewed as a function of two variables, it is the function that, in polar coordinates, can be written as f(r, θ) = g(r), that is, it is a radial function that is a function of ronly, and independent of θ. The function g(r) is, say, one for 0 ≤r≤R, where Ris the radius of the star, and zero, otherwise. From the theory of Fourier transform pairs in two-dimensions, we know that the two-dimensional Fourier transform of fis also a radial function; it is the function H(ρ) = 2 π/integraldisplayR 0rJ0(rρ)dr, where J0is the zero-th order Bessel function. From the theory of Bessel functions, we learn that d dx[xJ1(x)] =xJ0(x), so that H(ρ) =2π ρRJ1(Rρ). When the star is viewed through a telescope, the image is blurred by the atmosphere. It is commonly assumed that the atmosphere performs a convolution filtering on the light from the star, and that this filter is random and varies somewhat from one observation to another. Therefore, at each observation, it is not H(ρ), but H(ρ)G(ρ) that is measured, where G(ρ) is the filter transfer function operating at that particular time. 7 Suppose we observe the star Ntimes, for each n= 1,2, ..., N measuring values of the function H(ρ)Gn(ρ). If we then average over the various measurements, we can safely say that the first zero we observe in our measurements is the first zero of H(ρ), that is, the first zero of J1(Rρ). The first zero of J1(x) is known to be about 3.8317, so knowing this, we can determine R. Actually, it is not truly Rthat we are measuring, since we also need to involve the distance Dto the star, known by other means. What we are measuring is the perceived radius, in other words, half the subtended angle. Combining this with our knowledge of D, we get R. 7 Orthogonality of Bessel Functions As we have seen previously, the orthogonality of trigonometric functions plays an important role in Fourier series. A similar notion of orthogonality holds for Bessel functions. We begin with the following theorem. Theorem 7.1 Letu(x)be a non-trivial solution of u/prime/prime(x) +q(x)u(x) = 0 . If /integraldisplay∞ 1q(x)dx=∞, then u(x)has infinitely many zeros on the positive x-axis. Bessel’s Equation x2y/prime/prime(x) +xy/prime(x) + (x2−p2)y(x) = 0 , (7.1) can be written in normal form as y/prime/prime(x) +/parenleftBigg 1 +1−4p2 4x2/parenrightBigg y(x) = 0 , (7.2) and, as x→ ∞ , q(x) = 1 +1−4p2 4x2→1, so, according to the theorem, every non-trivial solution of Bessel’s Equation has infinitely many positive zeros. Now consider the following theorem, which is a consequence of the Sturm Com- parison Theorem to be discussed later. Theorem 7.2 Letyp(x)be a non-trivial solution of Bessel’s Equation x2y/prime/prime(x) +xy/prime(x) + (x2−p2)y(x) = 0 , forx >0. If0≤p <1 2, then every interval of length πcontains at least one zero of yp(x); ifp=1 2, then the distance between successive zeros of yp(x)is precisely π; and ifp >1 2, then every interval of length πcontains at most one zero of yp(x). 8 It follows from these two theorems that, for each fixed p, the function yp(x) has an infinite number of positive zeros, say λ1< λ 2< ..., with λn→ ∞ . For fixed p, letyn(x) =yp(λnx).We have the following orthogonality theorem. Theorem 7.3 Form/negationslash=n,/integraltext1 0xym(x)yn(x)dx= 0. Proof: Letu(x) =ym(x) and v(x) =yn(x). Then we have u/prime/prime+1 xu/prime+ (λ2 m−p2 x2)u= 0, and v/prime/prime+1 xv/prime+ (λ2 n−p2 x2)v= 0. Multiplying on both sides by xand subtracting one equation from the other, we get x(uv/prime/prime−vu/prime/prime) + (uv/prime−vu/prime) = (λ2 m−λ2 n)xuv. Since d dx/parenleftBig x(uv/prime−vu/prime)/parenrightBig =x(uv/prime/prime−vu/prime/prime) + (uv/prime−vu/prime), it follows, by integrating both sides over the interval [0 ,1], that x(uv/prime−vu/prime)|1 0= (λ2 m−λ2 n)/integraldisplay1 0xu(x)v(x)dx. But x(uv/prime−vu/prime)|1 0=u(1)v/prime(1)−v(1)u/prime(1) = 0 . References [1] Simmons, G. (1972) Differential Equations, with Applications and Historical Notes . New York: McGraw-Hill. 9