stakgold chap 7 exercise 7_33

Exercise 7.33: Heat flow with spherical delta shell initial condition. (290) Exercise 7.33: Heat flow with spherical delta shell initial condition. (290) 1 Wrong Turn (ignore this section) 2 Back on the Track 4 Acrobat OCR Exercise 12 Redo a part of this problem after "properties of Bessel" written. 12 –––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––– Overview from Meta Notes Exercise 7.33 ( p 290) Problem Statement: Consider an infinite 3D uniform heat conducting medium with a spherical hole of radius r = a at the origin. The hole contains vacuum which we assume does not conduct heat. The "boundary" here then is the sphere at r = a and the great sphere at r = ∞ and we assume that u=0 on the boundary. We take the entire medium to be at u = 0 for time t< 0. At time t=0 we turn on a spherical shell initial condition at radius r0 > a lying inside the medium. This is not a continuous heat source, it is just an initial condition, like saying the left inch of a frozen rod is at u=1 at t = 0+ε. The heat energy in this initial condition region is going to spread out (with inwards and outwards), and that is our problem. This problem has both θ and φ symmetry, so is a 1D problem in variable r. Recall how a heat source initial condition turns into a Helmholtz equation driving source, and that happens here. To solve this problem, we first do a time Laplace to replace t by s, and we then get a 1D Helmholtz Green's problem in r for U(r,s): - (r2U')' -λr2U = δ(r-r0)/4π where λ = -s. We realize that this is just the special case n=0 of a Green's equation we encountered in Ex 7.31 (whose function was called an(r)) so at this point we know the general atomic form of the solution must be Z1/2(r)/. Also, we have certain BC's on U(r,s) at r=a and r=∞ (0 in both places). We then solve the Green's problem in r to get U(r,s) = (i/8) [ H1/2(1)(a)J1/2(r<) - J1/2(a) H1/2(1)(r<)] [H1/2(1)(r>)/ H1/2(1)(a)] But then since λ = -s this gets rewritten tat U(r,s) = (1/4π)[ K1/2(1)(a)I1/2(r<) - I1/2(a) K1/2(1)(r<)][K1/2(1)( r>)/ K1/2(1)( a)] But these simple K functions are in fact just elementary functions, allowing the above to be written as U(r,s) = - (1/[8πr r0)]){ exp[(2a-r-r0) ] - exp[-|r-r0|] } / I verify at this point that this meets both BC's. Each term above has a trivial Schaum inverse Laplace, and we then get this very simple result u(r,t) = (1/8πrr0){ exp(-|r-r0|2/4t)/ – exp(- (2a-r-r0)2/4t)/ } I then take limit t→0 to verify that we get the initial condition shell δ. I plot the solution in Maple. –––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––– Consider an infinite uniform 3D heat conducting medium with a spherical hole of radius r = a at the origin. The hole contains vacuum which we assume does not conduct heat. The "boundary" here is the sphere at r = a and the great sphere at r = ∞ and we assume that u=0 on the boundary. We take the entire medium to be at u = 0 for time t< 0. At time t=0 we turn on a spherical shell initial condition at radius r0 > a lying inside the medium. This is not a continuous heat source, it is just an initial condition, like saying the left inch of a frozen rod is at u=1 at t = 0+ε. The heat energy in this initial condition region is going to spread out (with inwards and outwards), and that is our problem. I suspect this will convert into a Helmholtz problem with a shell source on the RHS, but let's give a shot. The heat equation system is this (∂t - 2) u = 0 t>0 boundary condition: u(a,t) = 0 u(∞,t) = 0 initial condition: u(r,0) = f(r) = δ(r-r0)/4πr2 We know from symmetry that u(r,t) = u(r,t), so only one spatial coordinate to worry about. Let's check on the form of this source. If we integrate it over a thick shell containing the thin shell, we will get 4π from the dΩ and an r2 from the radial integral, and we will get 1 as the integral, so this is a "unit" initial condition in some sense. Notice that 2 = 2r only. The first step is to do our Laplace time transform U(r,s) = !Syntax Error, Idt e-stu(r,t) We just apply !Syntax Error, Idt e-st to both terms of our homo PDE and get !Syntax Error, Idt e-st ∂tu(r,t) - !Syntax Error, Idt e-st2u(r,t) = 0 !Syntax Error, Idt e-st ∂tu(r,t) - 2!Syntax Error, Idt e-stu(r,t) = 0 !Syntax Error, Idt e-st ∂tu(r,t) - 2 U(r,s) = 0 Now we do parts integration on the first term, !Syntax Error, Idt e-st ∂tu(r,t) = - !Syntax Error, Idt ∂te-st u(r,t) + [e-st u(r,t)]∞0 = + s !Syntax Error, Idt e-st u(r,t) + [ 0 - u(r,0)] = s U(r,s) - u(r,0) So our transformed equation is now s U(r,s) - u(r,0) - 2 U(r,s) = 0 or s U(r,s) - 2 U(r,s) = u(r,0) or (-2 - (-s) ) U(r,s) = u(r,0) = δ(r-r0)/4πr2 So this has the standard Helmholtz equation form with λ = -s and we see our source as expected. So we end up with this 1D Green's Function problem: -2f(r) = - r-2∂r(r2∂rf(r)) = - r-2 (r2f')' We think of s as a fixed parameter and U=U(r) and we then have - r-2(r2U')' +sU = δ(r-r0)/4πr2 - (r2U')' +sr2U = δ(r-r0)/4π // this is our 1D Green's Equation Now let's use λ = -s and stay in "standard form" so we write - (r2U')' -λr2U = δ(r-r0)/4π It is certainly useful to compare this to the 1D Green's equation we got in Ex 7.31, which was - (r2an'(r))' – [λr2– n(n+1)] an(r) = δ(r-r0) (2n+1)/4π Wrong Turn (ignore this section) Therefore, this Green's problem is the n=0 case of the problem we already solved! So we can just read off the solution W(r,λ) = i [ (2n+1)/[8] Jn+1/2(r<) H(1)n+1/2(r>) = U(r,s) Now if we had λ = -k2 we would write, again from Ex 7.31, that = +ik and Jν(ikr<) H(1)(ikr>) = (2/iπ) Iν(kr<)Kν(kr>) and in our present case we just set k = , so here is our solution U(r,s) = i [ (2n+1)/[8] (2/iπ) I n+1/2 (kr<)K n+1/2 (kr>) // but with n=0 and k= =( i /[8])(2/iπ) I 1/2 (r<)K 1/2 (r>) =( 1 /[4π]) I 1/2 (r<)K 1/2 (r>) So this is the s-space solution to our problem. How it happens that we know So we can write I 1/2 (r<)K 1/2 (r>) = (2/π)1/2 (r<)-1/2 sh(r<) (2/π)-1/2 (r>)-1/2 exp(-r>) = (r<)-1/2 sh(r<) (r>)-1/2 exp(-r>)/ // note that ()-1/2()-1/2 = 1/ = (rr0)-1/2 sh(r<) exp(-r>)/ Therefore our s-space solution may be written as U(r,s) =( 1 /[4π])(rr0)-1/2 sh(r<) exp(-r>)/ = (1/4πrr0) sh(r<) exp(-r>)/ In hopes of doing a lookup, we would do well to expand sh(r<) = (1/2) [exp(r<) - exp(-r<) ] and then we have U(r,s) =(1/8πrr0){ exp(r<) exp(-r>) - exp(-r<) exp(-r>)}/ =(1/8πrr0){ exp(-(r>- r<)) - exp(-(r>+ r<)) }/ Now (r>- r<) = a positive number = |r-r0| and r>+ r< = r+r0 so U(r,s) = (1/8πrr0) { exp(-|r-r0|)/ - exp(- (r+r0) )/ } Now from Schaum p 169 we have that exp(-b)/ ↔ exp(-b2/4t)/ So our solution should then be u(r,t) = (1/8πrr0) { exp(-|r-r0|2/4t)/ - exp(-(r+r0)2/4t)/ } But how does that respect the BC ? Oooooo.. I think I made a big mistake. I used the wrong left end boundary condition on the 1D Green's function! My solution might be right in the limit a→ 0. This is going to make things very much messier because in the low region we can have a mix of J and H1. Back on the Track We had our ODE for U = U(r,s) correctly I think obtained as - (r2U')' -λr2U = δ(r-r0)/4π And it is certainly useful to compare this to the 1D Green's equation we got in Ex 7.31, which was - (r2an'(r))' – [λr2– n(n+1)] an(r) = δ(r-r0) (2n+1)/4π and with n=0 this gives our current situation including the RHS stuff (but not the BC's). We know that our homo equation solutions will be of this form J1/2(r)/ and we know this because our homo is the Ex 7.31 homo with n=0. But, we have now to state our boundary conditions on U(s,t). That is where I went wrong above. We had: boundary condition: u(a,t) = 0 u(∞,t) = 0 and if we apply !Syntax Error, Idt e-st to both these we find that U(a,s) = 0 and U(∞,s) = 0 The section between the two lines is replaced at the end with a much smaller section! ____________________________________________________________________________ We know on the right side we still need H1/2(1)(r)/ . But on the left we can have fleft(r) = [ A J1/2(r) + B H1/2(1)(r)] / If we want this to vanish at a, we require that [ A J1/2(a) + B H1/2(1)(a)] / = 0 => A J1/2(a) + B H1/2(1)(a) => B/A = - J1/2(a)/ H1/2(1)(a) So we can frame things in terms only of constant A by writing fleft(r) = A [ J1/2(r) + B/A H1/2(1)(r)] / = A [ J1/2(r) - [J1/2(a)/ H1/2(1)(a)] H1/2(1)(r)] / But I know (from experience, thinking Smythian forms) that we should redefine this constant and write = C [ H1/2(1)(a)J1/2(r) - J1/2(a) H1/2(1)(r)] / and now it is obvious that this vanishes at r = a, so we then have U(r,s) = C (rr0)-1/2 [ H1/2(1)(a)J1/2(r<) - J1/2(a) H1/2(1)(r<)] H1/2(1)(r>) which is not so bad really. I won't bother to do the jump from scratch. we know jump = [ (2n+1)/4π ] /(-r02) = -1/(4πr02) // since n=0 for our current problem Our task now is to compute this jump from the U(r,s) form shown above, which will take some time and require a more compact notation perhaps. // Next day 8.3.11: // As before we write U(r+,s) = C (rr0)-1/2 [ H1/2(1)(a)J1/2(r0) - J1/2(a) H1/2(1)(r0)] H1/2(1)(r) U'(r+,s) = C (rr0)-1/2 [ H1/2(1)(a)J1/2(r0) - J1/2(a) H1/2(1)(r0)] H1/2(1)'(r) U'(r+,s)|r0 = C (r0)-1 [ H1/2(1)(a)J1/2(r0) - J1/2(a) H1/2(1)(r0)] H1/2(1)'(r0) U(r-,s) = C (rr0)-1/2 [ H1/2(1)(a)J1/2(r) - J1/2(a) H1/2(1)(r)] H1/2(1)(r0) U'(r-,s) = C (rr0)-1/2 [ H1/2(1)(a)J1/2'(r) - J1/2(a) H1/2(1)'(r)] H1/2(1)(r0) U'(r-,s) )|r0 = C (r0)-1 [ H1/2(1)(a)J1/2'(r0) - J1/2(a) H1/2(1)'(r0)] H1/2(1)(r0) jump = U'(r+,s)|r0 – U'(r-,s) )|r0 = C (r0)-1 * { [ H1/2(1)(a)J1/2(r0) - J1/2(a) H1/2(1)(r0)] H1/2(1)'(r0) – [ H1/2(1)(a)J1/2'(r0) - J1/2(a) H1/2(1)'(r0)] H1/2(1)(r0)} jump = U'(r+,s)|r0 – U'(r-,s) )|r0 = C (r0)-1 * { [ HaJ - Ja H] H' – [ HaJ' - Ja H'] H } // compact notation Now we have { [ HaJ - Ja H] H' - [ HaJ' - Ja H'] H } = HaJ H' - Ja H H' - HaJ' H + Ja H' H // cleared brackets = HaJ H' - HaJ' H // two terms cancelled (as usual) = Ha W[J,H] // Wronskian form = Ha2i/[π (r0)] // Wronskian taken from Ex 7.31 doc z = r0 So this tells us that jump = C (r0)-1 Ha2i/[π (r0)] = C (2i/πr02) H1/2(1)(a) We now set this to -1/(4πr02) as per above to get C (2i/πr02) H1/2(1)(a) = -1/(4πr02) C (2i) H1/2(1)(a) = -1/(4) C (8/i) H1/2(1)(a) = 1 C = (i/8) [H1/2(1)(a)]-1 Our Green's function solution is therefore U(r,s) = C (rr0)-1/2 [ H1/2(1)(a)J1/2(r<) - J1/2(a) H1/2(1)(r<)] H1/2(1)(r>) = (i/8) [ H1/2(1)(a)J1/2(r<) - J1/2(a) H1/2(1)(r<)] [H1/2(1)(r>)/ H1/2(1)(a)] This solves our Helmholtz equation where λ = -s, and meets the boundary conditions specified. Now we can use the idea that when λ = -k2 we chose = +ik so here we will use = i. So we could put that in everywhere above. Then we want to use our facts quoted ni Ex 7.31 which we write as Iν(z) = e-iπν/2 Jν(iz) z above the cut Kν(z) = (iπ/2) e+iπν/2Hν(1)(iz) z above the cut Jν(iz) = e+iπν/2 Iν(z) Hν(1)(iz) = (2/iπ) e-iπν/2 Kν(z) Hν(1)(iz) Jν(iz') = (2/iπ) Kν(z) Iν(z') Hν(1)(iz)/ Hν(1)(iz') = Kν(z)/ Kν(z') Applying these last two rules to our U(r,s) above we get (for example, iz = a = (+)ia = ia , so then z = a ) U(r,s) = (i/8) [ H1/2(1)(a)J1/2(r<) - J1/2(a) H1/2(1)(r<)] * [H1/2(1)(r>)/ H1/2(1)(a)] U(r,s) = (i/8) (2/iπ) [ K1/2(1)( a)I1/2(r<) - I1/2(a) K1/2(1)( r<)] * [K1/2(1)( r>)/ K1/2(1)( a)] = (1/4π)[ K1/2(1)( a)I1/2(r<) - I1/2(a) K1/2(1)( r<)] * [K1/2(1)( r>)/ K1/2(1)( a)] and we have everything now in terms of s. ____________________________________________________________________________ Next we note that these functions are just elementary functions so that that last factor in U becomes [K1/2(1)( r>)/ K1/2(1)( a)] = (r>)-1/2 exp(-r>)/ [(a)-1/2 exp(-a)] = (r>/a)-1/2 exp[-(r> - a)] and the first factor becomes [ K1/2(1)( a)I1/2(r<) - I1/2(a) K1/2(1)( r<)] = (ar<)-1/2 [ exp(-a) sh(r<) - exp(-r<) sh(a)] = (ar<)-1/2 * [ exp(-a) sh(r<) - exp(-r<) sh(a)] = (ar<)-1/2 (1/2)* [ exp(-a) {exp(r<)- exp(-r<)} - exp(-r<) {exp(a)- exp(-a)}] = (ar<)-1/2 (1/2)* [ exp(-a) {exp(r<) - exp(-r<)} - exp(-r<) {exp(a)- exp(-a)}] = (ar<)-1/2 (1/2)* [ exp(-a) exp(r<) - exp(-a)exp(-r<)} - exp(-r<)exp(a) + exp(-r<)exp(-a)}] = (ar<)-1/2 (1/2)* [ exp(-(a-r<) - exp(-(a+r<)} - exp(+(a-r<)+ exp(-(a+r<)} ] = (ar<)-1/2 (1/2)* [ exp(-(a-r<) - exp(+(a-r<) ] = - (ar<)-1/2 sh[(a-r<)] = - (ar<)-1/2 sh[(a-r<)] / To summarize, we have just evaluated our two factors: [ K1/2(1)( a)I1/2(r<) - I1/2(a) K1/2(1)( r<)] = - (ar<)-1/2 sh[(a-r<)] / [K1/2(1)( r>)/ K1/2(1)( a)] = (r>/a)-1/2 exp[-(r> - a)] If we install these now into our U expression we get U(r,s) = (i/8) (2/iπ) [ K1/2(1)( a)I1/2(r<) - I1/2(a) K1/2(1)( r<)] * [K1/2(1)( r>)/ K1/2(1)( a)] = (i/8) (2/iπ) {- (ar<)-1/2 sh[(a-r<)] / } {(r>/a)-1/2 exp[-(r> - a)]} = (i/8) (2/iπ) {- (r<)-1/2 sh[(a-r<)] / } {(r>)-1/2 exp[-(r> - a)]} = (i/8) (2/iπ) {- sh[(a-r<)] / } { exp[-(r> - a)]} = (1/[8r r0)] (2/π) {- sh[(a-r<)] / } { exp[-(r> - a)]} = - (1/[4πr r0)]) sh[(a-r<)] exp[-(r> - a)] / = - (1/[8πr r0)]){ exp[(a-r<)] - exp[-(a-r<)]} exp[-(r> - a)] / = - (1/[8πr r0)]){ exp[(a-r<)] exp[-(r> - a)] - exp[-(a-r<)] exp[-(r> - a)]} / = - (1/[8πr r0)]){ exp[(a-r<)] exp[- (r> - a)] - exp[-(a-r<)] exp[- (r> - a)]} / = - (1/[8πr r0)]){ exp[{(a-r<)-(r>-a)} ] - exp[-{(a-r<) + (r> - a)}] } / = - (1/[8πr r0)]){ exp[(2a-r<-r>) ] - exp[-(r> - r<)] } / = - (1/[8πr r0)]){ exp[(2a-r-r0) ] - exp[-|r-r0|] } / Comment: in the limit a→0, you would expect this result to match our "wrong turn" result above. If I take that limit on the above I get = - (1/[8πrr0)]){ exp[(-r-r0) ] - exp[-|r-r0|] } / And this agrees exactly with my "wrong turn" result. [ I ferreted out three errors by doing this check and they are all fixed above so the reader does not know they ever existed. ] [ Doing wrong turn paid off. ] Notice how small the difference is between the general a result and the a=0 result! Now here is our U(r,s) one more time with terms reversed U(r,s) = (1/8πrr0) { exp[-|r-r0|]/s – exp[(2a-r-r0) ]/ } Let's check our BC's here. First, U(a,s) = (1/8πrr0) { exp[-|a-r0|]/s – exp[(2a-a-r0) ]/ } = (1/8πrr0) { exp[-|a-r0|]/s – exp[(a-r0) ]/ } But we know that r0 where the initial condition was set is > a, the hole radius, so r0>a so r0-a > 0 so = (1/8πrr0) { exp[-(r0-a) ]/ – exp[-(r0-a)) ]/ } = 0 // correct Next we can see that as r→∞, both expos decay and again we get 0. So BC's are met! Now, once again, from Schaum p 169 we have that exp(-b)/ ↔ exp(-b2/4t)/ So our solution should then be U(r,s) = (1/8πrr0) { exp[-|r-r0|]/ – exp[(2a-r-r0) ]/ } u(r,t) =(1/8πrr0){ exp(-|r-r0|2/4t)/ – exp(- (2a-r-r0)2/4t)/ } If we take a→0, this replicates our "wrong turn" result. So this thing looks good. Now as t→∞, both terms approach e0/ = 0. This seems reasonable since our finite initial energy pulse in the shell gets absorbed by the infinite size medium. Now what about t→0 ? The second term → exp(-k∞)/ = 0 since k>0 for all r including r=r0. Note that we cannot have both r and r0 = a since r0 > a. So we have Limt→0 u(r,t) = (1/8πrr0) Limt→0 { exp(-|r-r0|2/4t)/ } As long as r ≠r0, we get 0, so it looks like a delta function. We then have to integrate and see what happens !Syntax Error, I dr limt→0 [u(r,t)] = (1/8πrr0)!Syntax Error, I dr Limt→0 { exp(-|r-r0|2/4t)/ } We keep ε and a small but finite value, then we can interchange order = (1/8πrr0) Limt→0 !Syntax Error, Idr { exp(-|r-r0|2/4t)/ } = (1/8πrr0) Limt→0 { (1/)!Syntax Error, Idr exp(-|r-r0|2/4t) Change variable to s = r-r0 to get = (1/8πrr0) Limt→0 { (1/)!Syntax Error, Ids exp(-|s|2/4t) } = (1/4πrr0) Limt→0 { (1/)!Syntax Error, Ids exp(-s2/4t) } = (1/4πrr0) Limt→0 { (1/) [erf(ε/2 )] } // according to Maple = (1/4πrr0) Limt→0 { erf(ε/2 )} = (1/4πrr0) erf(∞)} = (1/4πrr0) * 1 Therefore we have shown that Limt→0 u(r,t) = (1/4πrr0) δ(r-r0) = (1/4πr02) δ(r-r0) and this agrees with our initial condition. So I think my final answer is correct which I repeat here: U(r,s) = (1/8πrr0) { exp[-|r-r0|]/ – exp[(2a-r-r0) ]/ } u(r,t) =(1/8πrr0){ exp(-(r-r0)2/4t)/ – exp(- (2a-r-r0)2/4t)/ } **** Here is a little Maple treatment The left edge of the cube is r = 1. You see that the solution is glued down to 0 at r=1 and r=∞ as required. At t = 0 we start with a delta at r=2 and over time it spreads out. Heat solutions always look like this. Acrobat OCR Exercise By the way, I pasted the above Maple equation into Acrobat 8 and asked for OCR and got this u : = (1/(8*Pi *r*r0)) *(exp( -(r-r0)^2/( 4*t)) - exp( -(2* a-r-r0)^2/( 4*t))) ; As Marilyn said, the OCR is not very good. I highlight all the errors in red. Redo a part of this problem after "properties of Bessel" written. Above we arrive at this point - (r2U')' -λr2U = δ(r-r0)/4π λ = -s We apply the result of Example 1A of the properties doc with ν = 1/2 to conclude that U = (1/4π)[ I1/2(kr<) K1/2 (ka)) – K1/2 (kr<) I1/2 (ka)] [ K1/2 (kr>)/ K1/2 (ka)] where k=and this agrees with the result above just above the second line.