laplaces method

The Laplace Method for ODE solving PhL 12.7.08 These are notes on a little 3-page web download that I have now printed and put title info on. From a book on the subject of "integral transforms" which would I guess includes Fourier transforms and Laplace transforms. 18.1 Laplace's Method We start with 18.1 which is an ODE of arbitrary order n, but each coefficient is at most linear in the coordinate x. The transform is shown in 18.2. There, if contour C were (0,∞), then y(-x) would be the usual Laplace transform and p is usually called s as in Schaum p 161. The author's comment on this seems wrong, but no matter. We are used to the regular Laplace Transform as solving ODE's with constant coefficients, but here this will be generalized to linear coefficients of the form a+bx. So here we go: (k sum runs 0 to n ) Σk (ak + bkx) ∂ky = 0 18.1 y(x) = ∫C S(p) epx dp 18.2 Σk (ak + bkx) ∂k { ∫C S(p) epx dp }= 0 Σk (ak + bkx) { ∫C S(p) ∂k epx dp }= 0 Σk (ak + bkx) { ∫C S(p) pk epx dp }= 0 And then we can write x epx = ∂p epx in the second term and we then have ∫C dp S(p) Σk ak pk epx + ∫C dp S(p) Σk bk pk ∂p epx = 0 and now all x reference is gone except in the exponent. Now define F and P as shown so we have = ∫C dp S(p) F(p) epx + ∫C dp S(p)G(p) ∂p epx = 0 Now in the second term, do parts integration to get = ∫C dp S(p) F(p) epx – ∫C dp ∂p [S(p)G(p)] epx + { S(p)G(p) epx }parts = 0 = ∫C dp { S(p) F(p) – ∂p [S(p)G(p)] } epx + { S(p)G(p) epx }parts = 0 18.3 And after a slow start, we arrive there at 18.3. Now assume that we select the contour so that the parts vanish, then we see that the equation is true if: S(p) F(p) – ∂p [S(p)G(p)] = 0 18.4 But we know how to solve this simple first order ODE: SF –∂p [SG] = 0 => SF – S∂pG – G∂pS = 0 => SF dp = SdG + GdS => Fdp = dG + G dS/S => F/Gdp – dG/G = dS/S => ln S = ∫dp (F/G) – ln G => ln(SG) = ∫dp (F/G) => SG = exp [∫dp (F/G)] => S(p) = [A/G(p)] exp [∫p dp' (F(p')/G(p')] 18.5 where we have added an arbitrary overall constant A since 18.4 is linear in S. So let's now review. Somebody hands you this ODE of order n Σk (ak + bkx) ∂ky(x) = 0 You apply the transform y(x) = ∫C S(p) epx dp and your new ODE equation in p-space is this S(p) F(p) – ∂p [S(p)G(p)] = 0 which you can solve to get S(p) = [A/G(p)] exp [∫p dp' (F(p')/G(p')] where F(p) = Σk ak pk G(p) = Σk bk pk so that F and G are polynomials of degree n, where n was the order of the original ODE. But this integral is of a ratio of polynomials and probably for small n it is a doable integral. So then we have to do this integral, and we have to find a contour C for our transform that makes parts = 0. 18.2 Application to the Hermites. We just follow the instructions, I did pencil notes, and we get a super simple result for S(p). This means that the solution to our original ODE must be this: y(x) = ∫C S(p) epx dp = – ½ ∫C {A [ exp(-p2/4)/pν+1 ]epx dp = – ½ ∫C A[ exp(-p2/4 + px)/pν+1 ]dp which is now an "integral representation" for a solution of the ODE. Pencil notes show F(p) = (2ν + p2) G(p) = -2p The constant is arbitrary anyway, so write w(x) = A' ∫C [ exp(-p2/4 + px)/pν+1 ]dp 18.7 Now let's write out those "parts" for this problem: { S(p)G(p) epx }parts = { exp(-p2/4 + px)/pν+1 * (-2p) * epx} parts ~ { exp(-p2/4 + 2px)/pν } parts Rescaling: Now comes a scaling thing which is not completely obvious to me. Go back to the original deal y(x) = ∫C S(p) epx dp ∫C dp S(p) F(p) epx – ∫C dp ∂p [S(p)G(p)] epx + { S(p)G(p) epx }parts = 0 F(p) = (2ν + p2) G(p) = -2p S(p) F(p) – ∂p [S(p)G(p)] = 0 Replace p everywhere with variable q = ½ p and p = 2q and we get y(x) = 2 ∫C S(2q) e2qx dq ∫C dq S(2q) F(2q) e2qx – ∫C dq ½ ∂q [S(2q)G(2q)] e2qx + 2{ S(2q)G(2q) e2qx }parts = 0 Now define new functions s(q) = S(2q) f(q) = F(2q) = (2ν+4q2) g(q) = G(2q) = -4q Then we have y(x) = 2 ∫C s(q) e2qx dq ∫C dq s(q) f(q) e2qx – ∫C dq ½ ∂q [s(q) g(q)] e2qx + 2{ s(q)g(q) e2qx }parts = 0 This then looks like our original equation but we have an extra ½ and "2" which came from making the ∂p term in this case. Our q-space ODE is then s(q) f(q) – ½ ∂q [s(q)g(q)] = 0 like replacing G with g/2 s(q) = [A/g(p)] exp [∫q dp' 2 (f(p')/g(p')] = A/ (-4q) * exp [∫q dp' ((2ν+4p'2)/(-2p')] = A/ (-4q) * exp [∫q dp [ -ν (1/p') - 2p' ] ] = A/ (-4q) * exp [ - ν ln q – q2 ] = A/ (-4qν+1) * exp [– q2 ] ~ A'/qν+1 exp(-q2) So review the problem. We start with the Hermite ODE in x. We do this transform: y(x) = ∫C s(q) e2qx dq The new ODE in q space is this: s(q) f(q) – ½ ∂q [s(q)g(q)] = 0 The solution to this equation is this: s(q) = A'/qν+1 exp(-q2) So we then get this integral representation for the Hermite solution w(x) = A∫C exp(-q2+2qx) /qν+1 dq 18.7 and then the parts we need to vanish is this: 2{ s(q)g(q) e2qx }parts ~ { A'/qν+1 exp(-q2) (-4q) e2qx }parts ~ { 1/qν exp(-q2+2qx) }parts but my power now disagrees with 18.8. I think the author is wrong here about the parts. Errata check? Brian Davies is retired at Canberra, gives no errata, nor do I find on the web with author, book title, and errata. So I would guess there is none, I might email him. The power does not matter much because of the strong exponential. Let's now look at his list of contours. (i) The real axis one seems reasonable, but what do you do at p=0? Is this a "tick" integral, or do you detour around it, and if so, which way? (ii) If origin is one end, cannot have ν be positive, need Re(ν) < 0, we agree and the error is confirmed (iii) circle around the origin? OK if there are no branch cuts from the power, otherwise all is single-valued. So need ν = integer and OK, we agree. (iv) a drawn contour p 306. Certainly the two end points have no parts since both at infinity. For the Hermite Polynomials, we have ν = n = integer, so the circle of (iii) works. This gives 18.9 with a certain normalization that he states: coefficient of xn is 2n, and I see this is how they are done in page 151 Schaum, so boom, we have an integral representation for the Hermites! How could I prove this normalization? Hn(x) = 2n xn + lower powers ∂n Hn(x) = 2n n! Apply these derivatives to the integral representation and use ∂n e2qx = (2q)n e2qx e2qx ∂n Hn(x) = ∂n A∫C dq exp(-q2+2qx) /qn+1 = = A∫C dq exp(-q2+2qx) (2q)n/qn+1 = A 2n ∫C dq exp(-q2+2qx)/q The residue is just 2πi I think, so I then get 2n n! = 2πi A 2n => A = n!/(2πi) 18.9 Very good. Now let's just do the whole thing as a residue, so that Hn(x) = A∫C dq exp(-q2+2qx) /qn+1 = 2πi (1/n!) A [ ∂qn+1 exp(-q2+2qx)]q=0 = [ ∂qn exp(-q2+2qx)]q=0 Test this for n = 1: H1(x) = [ ∂q exp(-q2+2qx)]q=0 = (-2q+2x)| = 2x = correct Plan A: try converting this result directly to Rodriquez' form (failed, go to Plan B) Now change from variable q to variable y where x is treated as a constant: -y2 = -q2+2qx q2 - 2xq – y2 = 0 q = variable -2ydy = -2qdq + 2xdq = 2(x-q)dq ydy = (q-x)dq dy/dq = (q-x)/y Now solve for q in terms of y: q = ( 2x ± )/2 = ( x ± ) = q(y) (q-x) = ± . But we shall be interested in the region near q = 0 for our parts evaluation, and we care about x > 0, so we need to choose the sign this way (q-x) = – dy/dq = (q-x)/y = – /y ∂q = d/dq = dy/dq * d/dy = – /y * ∂y Test this for n=1: H1(x) = [ ∂q exp(-q2+2qx)]q=0 = - [/y * ∂y exp(-y2) ]y=0 = - /y *(-2y) = +2 | = 2x // correct But sign is wrong, could fix by taking other square root sign above, but how justify that? Higher powers? ∂q2 = [ – /y * ∂y ] [ – /y * ∂y ] This just looks very messy. Write = r(y) so we have ∂yr = ½ 1/r 2y = (y/r) . Then ∂q2 = [ – r/y * ∂y ] [ – r/y * ∂y ] = + (r/y) { (∂yr) y ∂y - r/y2∂y + r/y ∂y2 } = (r/y) { (y/r) y ∂y - r/y2∂y + r/y ∂y2 } = (r/y) { [(y2/r) - (r/y2)]∂y + (r/y) ∂y2 } I guess I am missing something here, like we are expanding something around the wrong point. It probably works in the end, but lets instead do Plan B: Rodriquez' formula Start with Hn(x) = n!/(2πi) ∫C dq exp(-q2+2qx) /qn+1 Do the recommended change of variables to u = q-x and q = u+x . We get = n!/(2πi) ∫C du exp(-u2+x2) /(u+x)n+1 where now the contour is a circle around the point u = -x. But this is = n!/(2πi)exp(x2) ∫C du exp(-u2) /(u+x)n+1 // agrees with p 306 A = n!/(2πi)exp(x2) * 2πi /n! [∂un exp(-u2)]u=-x = exp(x2) [∂un exp(-u2)]u=-x Now take u → -u and the derivative picks up (-1)n and we have = (-1)n exp(x2) [∂un exp(-u2)]u=+x = (-1)n exp(x2) ∂xn exp(-x2) // agrees with p 306 B and this is the Rodriguez seen in Schaum p 151. So this method probably works for lots of ODE's Generating function. I have to do this on my own since the author blocked the next page! Let's go back to our result earlier Hn(x) = n!/(2πi) ∫C dq exp(-q2+2qx) /qn+1 // his 18.9 Suppose we were to expand exp(-q2+2qx) formally as a power series in q: exp(-q2+2qx) = Σk=0,∞ Ak qk Then we have Hn(x) = n!/(2πi) ∫C dq Σk Ak qk /qn+1 = n!/(2πi) ∫C dq Σk Ak /qn+1-k = n!/(2πi) Σk Ak ∫C dq /qn+1-k If k = n+1 or more, there is no pole so get 0. If k = n, we have a simple pole and pick up 2πi. If k < n, get 0 since derivative of 1 = 0. So this says Hn(x) = n!/(2πi) Σk Ak 2πi δn,k = n! An Therefore we have shown that exp(-q2+2qx) = Σn=0,∞ Hn(x)/n! qn and this is the generating formula shown Schiff page 13.10. Comment: This Laplace Method gives us a whole pile of results: (1) A set of integral representations for the solution to an ODE with up to linear coefficients, any order. Each choice of contour with vanishing parts gives a different integral representation with the same integrand. The parameters can in general be complex. (2) The "circle" choice for a parameter = integers then gives you both Rodriquez and the generating function formula for the polynomial solutions of the ODE. Try doing the Legendre ODE which has quadratic coefficients (does not succeed) The Legendre ODE on page 146 Schaum has (1-x2) as the y" coefficient, so the Laplace method as presented above does not work. That is probably why the Legendre generating function looks so different. I could try the Laplace method keeping the next term, but I will then end up with a second order ODE in p, but maybe it is one with a simple solution. Try next term idea: Σk (ak + bkx + ckx2) { ∫C S(p) pk epx dp }= 0 And then we can write x epx = ∂p epx in the second term and third terms and we then have ∫C dp S(p) Σk ak pk epx + ∫C dp S(p) Σk bk pk ∂p epx + ∫C dp S(p) Σk ck pk ∂p2 epx = 0 and now all x reference is gone except in the exponent. Now define F and G and H as shown so we have = ∫C dp S(p) F(p) epx + ∫C dp S(p)G(p) ∂p epx + ∫C dp S(p)H(p) ∂2 epx = 0 Now in the second term, do parts integration to get = ∫C dp S(p) F(p) epx – ∫C dp ∂p [S(p)G(p)] epx + { S(p)G(p) epx }parts = 0 = ∫C dp { S(p) F(p) – ∂p [S(p)G(p)] } epx + { S(p)G(p) epx }parts = 0 18.3 But now we have to do parts twice on the third term. Let's just ignore the parts for the moment and cut to the chase: = ∫C dp S(p) F(p) epx – ∫C dp ∂p [S(p)G(p)] epx + ∫C dp ∂2(S(p)H(p)) epx = 0 Then we get S(p) F(p) – ∂p [S(p)G(p)] + ∂2[S(p)H(p)] = 0 Legendre says (1-x2)y(2) - 2xy(1) + n(n+1)y(0) = 0 p 146 Schaum So we then have F(p) = Σk ak pk G(p) = Σk bk pk H(p) = Σk ck pk a0 = n(n+1) a2 = 1 b1= -2 c2 = -1 So we then have F(p) = n(n+1) + p2 G(p) = -2 p H(p) = -p2 so our p-space ODE is then S(p)[ n(n+1) + p2]+2∂p [S(p)p] – ∂2[S(p) p2] = 0 ∂2 [S(p) p2] - 2 ∂p [S(p)p] + [ n(n+1) + p2] S(p) = 0 But write ∂2 [S(p) p2] = ∂ [ S 2p + ∂S p2] = (∂S)2p + 2S + ∂2S p2 + 2p ∂S = p2∂2S + 4p ∂S + 2S ∂p [S(p)p] = (∂S)p + S so our ODE becomes p2∂2S + 4p ∂S + 2S - 2p (∂S)- 2S + [ n(n+1) + p2] S(p) = 0 p2∂2S + 2p ∂S + [ n(n+1) + p2] S(p) = 0 and as claimed, this thing is worse that what we started with. This just says that the Legendre situation is not going to involve those same expo forms that the Hermite does. Knowing the Legendre answer, we might try starting our Laplace Method with a different transform, but I don't think that is the right approach. So conclusion is that Laplace works in those cases specified, ie, ODE that have at worst linear coefficients. Confluent Hypergeometric ODE This baby is this: z F(2) + (b-z)F(1) - aF(0) = 0 Schiff page 139 and this is a Laplace candidate. We have a0 = -a a1 = b b1 = -1 b2 = 1 F(p) = Σk ak pk = -a + bp G(p) = Σk bk pk = -p + p2 = p(p-1) S(p) = [A/p(p-1)] exp [∫p dp' ((-a+bp')/p'(p'-1)] The integral in the exponent is ∫p dx (-a+bx)/x(x-1) I have Maple compute this and get a*ln(x)-ln(x-1)*a+ln(x-1)*b = a ln(p) + (b-a) ln(p-1) Then we find exp[a ln(p) + (b-a) ln(p-1)] = pa (p-1)b-a And we seem to find that S(p) = [A/p(p-1)] pa (p-1)b-a = A pa-1 (p-1)b-a-1 Then our solution to the ODE will have the form y(x) = ∫C S(p) epx dp = A ∫C dp pa-1 (p-1)b-a-1epx We can change (p-1) to (1-p) and absorb into the constant. Doing this and changing p to t we get y(x) = ∫C dt etx ta-1 (1-t)b-a-1 // agrees Messiah p 480 B.4 Since this agrees, we know we have the right scaling, so we can go look at the parts condition, { S(p)G(p) epx }parts = 0 {pa-1 (p-1)b-a-1 p(p-1) epx } = { pa (p-1)b-a epx}parts = 0 So let's take this to say { ta (1-t)b-a etx}parts = 0 Consider now the structure in the t-plane of our integrand for y(x) above [ b = integer! ] On the main sheet, I would say arg(t) = 0 and arg(1-t) = 0 so as Messiah claims, these phases are equal. You could draw the cut going from 0 to 1 if you wanted (see contours.doc! ). Now back to what we had y(z) = ∫C dt etz ta-1 (1-t)b-a-1 // agrees Messiah p 480 B.4 The claim is that this function with this contour is analytic in z. In the neighborhood of any z, nothing weird happens, so I buy this claim. Messiah calls it an "entire function". Since the usual 1F1 power series solution is also analytic in the plane including z=0, he claims these two functions must then be proportional. Notice the notation: F(α|β|z) [ Messiah p 480] = M(α,β,z) [ AS p 504] = F1(α; β; z) = Φ(α,β; z) [ both Bateman p 248, the second used by GR 1057] GR refer to the "degenerate" hypergeometric function instead of "confluent". So we have F(α|β|z) = K ∫C dt etz ta-1 (1-t)b-a-1 Now expand etz in a power series to get = K Σk zk/k! ∫C dt tk ta-1 (1-t)b-a-1 = K Σk zk/k! ∫C dt tk+a-1 (1-t)b-a-1 Looking at Messiah B.2, it would appear that K ∫C dt tk+a-1 (1-t)b-a-1 = Γ(a+k)Γ(b)/{ Γ(a) Γ(b+k) } and so on. See separate document on the Messiah appendix "Messiah W functions.doc" Comment: The Laplace Method says that a contour integral of a certain integrand gives you a solution of an ODE with linear coefficients. The "parts" at the two ends of the contour must add to zero, which is always true for a closed contour, and true for some unclosed contours. An important point to remember: as you change the contour C, you change the solution! That is, varying the contour gives you solutions to the ODE which are likely different from each other. Perhaps contour C1 gives solution u1 while contour C2 gives solution 3u1 + 5u2 .