PDE, inhomogeneous diffusion equation

[fluidistic]

Homework Statement

Mathews and Walker problem 8-2 (page 253):
Assume that the neutron density n inside [itex]U_{235}[/itex] obeys the differential equation [itex]\nabla ^2 n+\lambda n =\frac{1}{\kappa } \frac{\partial n }{\partial t}[/itex] (n=0 on surface).
a)Find the critical radius [itex]R_0[/itex] such that the neutron density insdie a [itex]U_{235}[/itex] sphere of radius [itex]R_0[/itex] or greater is unstable and increases exponentially with time.

b)Suppose two hemispheres, each just barely stable, are brought together to form a sphere. This sphere is unstable, worth n ~ [itex]e^t/\tau[/itex]. Find the "time-constant" tau of the resulting explosion.

Homework Equations

Already given.

The Attempt at a Solution

I tried to solve this PDE via a Laplace transform, but now I realize I do not think it is well appropriated due to the nabla operator.
Applying Laplace transform gave me [itex]\nabla ^2 F(n)+\lambda F(n)=\frac{1}{\kappa } [n(t=0)+sF(n)][/itex]. Also I do not know the initial condition n(t=0).

I guess I'll have to try another method and that it might involve taking the Laplacian in spherical coordinates, I'm open to use any method to solve this PDE. I'd like to know whether I'm right in thinking that it's not solvable using a Laplace transform and what other method could I use.
Thank you very much in advance.

 

[nlsherrill]

Sorry, I really wish I could help but I don't think I know. This sounds like such a cool problem. What book is this question in?

I am assuming the neutron is to be considered a sphere? I think spherical coordinates would be your best bet, but I have not had experience with PDE's in spherical... Isn't the laplace transform usually used for solving PDE's on semi infinite boundaries?

 

[fluidistic]

Sorry, I really wish I could help but I don't think I know. This sounds like such a cool problem. What book is this question in?

I am assuming the neutron is to be considered a sphere? I think spherical coordinates would be your best bet, but I have not had experience with PDE's in spherical... Isn't the laplace transform usually used for solving PDE's on semi infinite boundaries?

The book is from Mathews and Walker -Mathematical methods of physics-. I think we do not consider a single neutron, but a sphere (nucleous of uranium) with a "nicely behaved" neutron density function.
I don't really know about the Laplace transform.
Yeah it looks like I'll have to take the Laplacian in spherical coordinates but... I'm unsure what method to use. A kind of separation of variables maybe.

 

[Mindscrape]

Hmm, yeah, in my experience, every PDE is really it's own beast that you just have to try different methods on. I think separation of variables should work out fine if the density is constant, is it?

 

[nlsherrill]

The book is from Mathews and Walker -Mathematical methods of physics-. I think we do not consider a single neutron, but a sphere (nucleous of uranium) with a "nicely behaved" neutron density function.
I don't really know about the Laplace transform.
Yeah it looks like I'll have to take the Laplacian in spherical coordinates but... I'm unsure what method to use. A kind of separation of variables maybe.

I think maybe try the separation of variables for spherical polar coordinates and see what you get...it will surely involve the tough to deal with bessel functions though

 

[fluidistic]

Ok guys I'll try separation of variables. But the neutron density isn't constant. It vanishes at the boundaries of the sphere (its surface). And inside the sphere I have no idea how it is, probably very complicated. It's the function that satisfies the PDE.

 

[Mindscrape]

My bad, I meant lambda whatever it is... I'm in EM mode. x-(

 

[fluidistic]

My bad, I meant lambda whatever it is... I'm in EM mode. x-(

No problem.
Indeed this method looks promising. I'm not math ninja though and I think I have to do some tricks.
So I indeed took the Laplacian in spherical coordinates. After simplifications I reached [tex]\frac{1}{R}\left (\frac{2}{r} R'+R'' \right )+\frac{1}{\Theta r^2 } \left ( \cot (\theta ) \Theta ' + \Theta '' \right )+\frac{\Psi ''}{\Psi r^2 \sin ^2 \theta }+\lambda =\frac{T'}{\kappa T} [/tex].
Where [itex]T(t)[/itex], [itex]R(r)[/itex], [itex]\Psi (\phi )[/itex] and [itex]\Theta (\theta )[/itex]. I assumed [itex]n (r, \theta, \phi , t)[/itex] to be the product of the 4 functions.

After a few considerations (like a function depending only on positon is equal to a function depending only on time for all positions+time must be that these functions are constants. Also, if I multiply the eq. by [itex]r^2 \sin \theta[/itex] or by [itex]r^2[/itex] I get other conditions leading me to what follows), I get that
(1)[itex]\frac{1}{\Psi } \Psi ''=C_1[/itex]. I could solve this equation but there will be 3 cases, depending upon the sign of [itex]C_1[/itex] and if it's worth 0.
(2)[itex]T'=\kappa T C_2 \Rightarrow T(t)= \tilde C_2 e ^{\kappa t}[/itex] and [itex]C_2=\frac{T'}{\kappa t }[/itex].
(3)[itex]C_3=R''+\frac{2R'}{r}+R(\lambda - C_2 )[/itex].
(4)[itex]C_4= \Theta ''+\Theta '+\frac{C_1 \Theta }{\sin \theta }[/itex].
Now it's "only" a matter of solving all these equations and then multiply each solutions together to form the general solution to the PDE.

 

[Mindscrape]

Yep, I think you're good so far. You just have to look up the solutions to the ODEs. They should be pretty similar to simply solving the helmholtz equation in spherical coordinates.

 

[sgd37]

what you have there is a modified schrodinger equation i.e. it doesn't give oscillating solutions but you can solve it just the same way

 

[fluidistic]

Ok I've followed a book for Hemlholtz equation (Mathews and Walker page 228) but it's not that easy to me.
I realize I shouldn't have calculated the derivatives in the Laplacian. So I've redone all the algebra keeping terms like in the book. But there's still a difference from the Helmholtz equation [itex]\nabla ^2 X +k^2 X=0[/itex]. In my case the equation isn't homogeneous but equal to a constant that I called [itex]C_1[/itex].
[tex]\frac{1}{Rr^2}\frac{1}{dr}(r^2R')+\frac{1}{r^2 \sin \theta }\frac{d}{d\theta } (\sin (\theta )\Theta )+\frac{1}{r^2 \sin ^2 \theta }\frac{\Psi'}{\Psi }+\lambda=\frac{1}{\kappa} \frac{T'}{T}=C_1[/tex].
So that [itex]T(t)=Ae^{c_1\kappa t }[/itex]. (I made a small mistake in my last post I think).
I also get that [itex]\frac{\Psi '}{\Psi}=\text{constant}=-m^2 \Rightarrow \Psi (\phi )=Be^{im\phi }+Ce^{-im\phi}[/itex]. This equation is the same for Helmholtz equation. I do not really know why the book chose the constant as "[itex]-m^2[/itex]", probably for convenience in the algebra. Nor do I know why it considered only a negative constant instead of getting 3 cases (positive, negative and equal to 0).
The next equation is [itex]\frac{1}{\sin \theta} \frac{1}{d\theta }(\sin (\theta ) \Theta ' )-\frac{m^2}{\sin ^2 \theta }\Theta = \text{constant}\cdot \Theta=-l(l+1)\Theta[/itex]. Which is the same as in Helmholtz equation.
Lastly the radial equation is [itex]\frac{1}{r^2}\frac{d}{dr} (r^2R')+\left [ \lambda -C_1-\frac{l(l+1)}{r^2} \right ]R=0[/itex]. Which differs slightly from Helmholtz equation. Indeed, I have [itex]\lambda-C_1[/itex] instead of [itex]k^2[/itex].
The book then gives the solution to the theta equation, [itex]P_l^m (x)[/itex], [itex]Q_l^m(x)[/itex]. Where [itex]x=\cos \theta[/itex].
While for the radial equation, adjusting the given result but substituting their k by [itex]\sqrt {\lambda - C_1}[/itex], I get [itex]R=\frac{J_{l+1/2}(\sqrt {\lambda -C_1}r)}{\sqrt r}[/itex] and [itex]\frac{Y_{l+1/2}(\sqrt {\lambda -C_1}r)}{\sqrt r}[/itex]. I do not know however if I have to take [itex]\pm \sqrt {\lambda -C_1}[/itex] or just the positive square root, etc.
So to answer to part a) I think I must NOT only look into the radial equation solution. For a particular value of [itex]r=R_0[/itex], when I multiply all functions together, I should get that the solution diverges when t tends to infinity.
Unfortunately multiplying all these functions together will create a too complicated function to analyze...

 

[sgd37]

m varies from +l to -l so you have all three considerations in your solution unless l is half integer of course

 

[fluidistic]

I'm back at this problem. I found some extra help, basically telling me to use the symmetry of the problem.
So I tried to redo all the exercise from start.
I'll use separation of variables and due to the symmetry of the problem I'll assume a solution of the form [itex]n(r,t)=R(r)T(t)[/itex].
My first doubt comes as soon as I take the Laplacian.
If I take the Laplacian in spherical coordinates, the PDE simplifies to [itex]\frac{2R'}{R}+r\frac{R''}{R}+\lambda = \frac{1}{\kappa } \frac{T'}{T}[/itex]. This gives me 2 DE's in which one isn't that simple.
However if I use the logic that the Laplacian only acts over the space coordinates (in this case r, and not on t), then the PDE simplifies to [itex]\frac{R''}{R}+\lambda = \frac{1}{\kappa } \frac{T'}{T}[/itex]. This gives rise to 2 "simple" DE's.
Why on Earth I get 2 different DE's to solve whether I consider the Laplacian in spherical coordinates or if I consider it only acts over r?
Which way is the correct one and why does the other way fail?
So I've kept following the easiest way (in other words the way of [itex]\nabla ^2 (RT)=T R''[/itex].)
First DE to solve: [itex]\frac{1}{\kappa} \frac{T'}{T}=m[/itex], with solution [itex]T(t)=Ae^{m \kappa t}[/itex].
The second DE to solve: [itex]\frac{R''}{R}+\lambda =m[/itex]. Note that m is the constant of separation of the method of separation of variables used in the problem.
Here I used logic and I assumed that only the case [itex]m>\lambda[/itex] holds true. Because if [itex]m<\lambda[/itex] I'd get a periodic solution in r, which I don't think to be true and if [itex]m=0[/itex] I'd get the constant solution equal to 0 which is also to be disgarted.
Hence [itex]R(r)=Be^{r\sqrt{m-\lambda }} + Ce^{-r\sqrt {m-\lambda}}[/itex]. Here I guess it's safe to assume [itex]C=0[/itex] so that [itex]n(r,t)[/itex] remains finite when [itex]r=0[/itex].
Thus [itex]R(r)=Be^{r\sqrt {m-\lambda}}[/itex].
Now I apply the boundary condition that [itex]n(\tilde r , t )=0[/itex] where [itex]\tilde r[/itex] is the radius of the sphere. But this implies that [itex]B=0[/itex] and thus [itex]n(r,t)=0[/itex] for all r and all t's. Which is obviously wrong.
Can someone point me out my error(s) and/or help me with this problem?
Thanks in advance.

 

[Mute]

I'm back at this problem. I found some extra help, basically telling me to use the symmetry of the problem.
So I tried to redo all the exercise from start.
I'll use separation of variables and due to the symmetry of the problem I'll assume a solution of the form [itex]n(r,t)=R(r)T(t)[/itex].
My first doubt comes as soon as I take the Laplacian.
If I take the Laplacian in spherical coordinates, the PDE simplifies to [itex]\frac{2R'}{R}+r\frac{R''}{R}+\lambda = \frac{1}{\kappa } \frac{T'}{T}[/itex]. This gives me 2 DE's in which one isn't that simple.
However if I use the logic that the Laplacian only acts over the space coordinates (in this case r, and not on t), then the PDE simplifies to [itex]\frac{R''}{R}+\lambda = \frac{1}{\kappa } \frac{T'}{T}[/itex]. This gives rise to 2 "simple" DE's.
Why on Earth I get 2 different DE's to solve whether I consider the Laplacian in spherical coordinates or if I consider it only acts over r?
Which way is the correct one and why does the other way fail?
So I've kept following the easiest way (in other words the way of [itex]\nabla ^2 (RT)=T R''[/itex].)
First DE to solve: [itex]\frac{1}{\kappa} \frac{T'}{T}=m[/itex], with solution [itex]T(t)=Ae^{m \kappa t}[/itex].
The second DE to solve: [itex]\frac{R''}{R}+\lambda =m[/itex]. Note that m is the constant of separation of the method of separation of variables used in the problem.
Here I used logic and I assumed that only the case [itex]m>\lambda[/itex] holds true. Because if [itex]m<\lambda[/itex] I'd get a periodic solution in r, which I don't think to be true and if [itex]m=0[/itex] I'd get the constant solution equal to 0 which is also to be disgarted.
Hence [itex]R(r)=Be^{r\sqrt{m-\lambda }} + Ce^{-r\sqrt {m-\lambda}}[/itex]. Here I guess it's safe to assume [itex]C=0[/itex] so that [itex]n(r,t)[/itex] remains finite when [itex]r=0[/itex].
Thus [itex]R(r)=Be^{r\sqrt {m-\lambda}}[/itex].
Now I apply the boundary condition that [itex]n(\tilde r , t )=0[/itex] where [itex]\tilde r[/itex] is the radius of the sphere. But this implies that [itex]B=0[/itex] and thus [itex]n(r,t)=0[/itex] for all r and all t's. Which is obviously wrong.
Can someone point me out my error(s) and/or help me with this problem?
Thanks in advance.

You have to use the 3d Laplacian. Your "simple" case amounted to using the 1d Laplacian, which is incorrect. You really have a 3d problem, you are just assuming that the solution you are seeking is the spherically symmetric solution, such that the angular derivatives in the Laplacian vanish, which leaves you with only the radial term from the 3d Laplacian. This symmetry does not allow you to treat the problem as if it were just a 1d problem. Suppose you were a masochist and wanted to try to solve this problem in cartesian coordinates: you would have to keep all three derivatives. What you did when you only used [itex]\partial^2/\partial r^2[/itex] would amount to only keeping the [itex]\partial^2/\partial x^2[/itex] term in the cartesian Laplacian, for instance.

 

[EulersFormula]

This problem is very similar to the Shrodinger equation. I recommend taking a look at it. I also wouldn't throw away solutions that suggest the neutron density to vary harmonically with distance. This is not uncommon.

 

[fluidistic]

Ok thank you very much guys, now it gets clearer to me.
So basically the 2 DE's to solve are:
(1): [itex]\frac{T'}{\kappa T }=m[/itex].
(2): [itex]\frac{2R'}{rR}+\frac{R''}{R}+\lambda =m[/itex].
The solution to (1) is I said already, [itex]T(t)=T_0e^{\kappa m t}[/itex].
I'm now trying to figure out the solution to (2). I tried the method of Frobenius, I proposed a solution of the form [itex]R(r)=\sum _{n=0} ^{\infty } a_n r^{n+c}[/itex]. Deriving with respect to r twice and plugging in the DE (2) I reached the secular equation [itex]c(c+1)=0[/itex] so that [itex]c=0[/itex] or [itex]c=-1[/itex]. Since the 2 roots differ by an integer, I can "only" reach 1 type of solution of the form I assumed, not the 2. I'd have to work out the second form of the solution by variation of parameters once I'm done with [itex]R_1 (r)[/itex].
So that I took [itex]c=-1[/itex] (no choice), this gave me the recurrence relation [itex]a_n = - \frac{a_{n-2} (\lambda - m )}{n ^2 -n}[/itex]. Where [itex]a_1[/itex] is arbitrary and [itex]a_0[/itex] too (except that it can't be 0).
For convenience I chose [itex]a_0=1[/itex] but even though I calculated the first terms ([itex]a_2[/itex], [itex]a_4[/itex] and [itex]a_6[/itex]) I don't really see a pattern to express the general form of [itex]a_n[/itex] independently of the previous [itex]a_{n-p}[/itex]'s. Also, should I choose [itex]a_1=0[/itex] for convenience so that all odd terms would be worth 0?
By the way this might be similar to the hydrogen atom in a way yes, but the method to solve the DE isn't the same I guess because one usually "factors out" the behavior the DE when r tends to infinity (in the hydrogen atom case) while here I'm only interested in solving the PDE inside the region of a sphere of radius a or [itex]\tilde r[/itex]; so that it's not really the same path to take. I can be wrong of course.
Any idea on how to solve the DE (2)?

 

[Chestermiller]

Ok thank you very much guys, now it gets clearer to me.
So basically the 2 DE's to solve are:
(1): [itex]\frac{T'}{\kappa T }=m[/itex].
(2): [itex]\frac{2R'}{rR}+\frac{R''}{R}+\lambda =m[/itex].
The solution to (1) is I said already, [itex]T(t)=T_0e^{\kappa m t}[/itex].
I'm now trying to figure out the solution to (2). I tried the method of Frobenius, I proposed a solution of the form [itex]R(r)=\sum _{n=0} ^{\infty } a_n r^{n+c}[/itex]. Deriving with respect to r twice and plugging in the DE (2) I reached the secular equation [itex]c(c+1)=0[/itex] so that [itex]c=0[/itex] or [itex]c=-1[/itex]. Since the 2 roots differ by an integer, I can "only" reach 1 type of solution of the form I assumed, not the 2. I'd have to work out the second form of the solution by variation of parameters once I'm done with [itex]R_1 (r)[/itex].
So that I took [itex]c=-1[/itex] (no choice), this gave me the recurrence relation [itex]a_n = - \frac{a_{n-2} (\lambda - m )}{n ^2 -n}[/itex]. Where [itex]a_1[/itex] is arbitrary and [itex]a_0[/itex] too (except that it can't be 0).
For convenience I chose [itex]a_0=1[/itex] but even though I calculated the first terms ([itex]a_2[/itex], [itex]a_4[/itex] and [itex]a_6[/itex]) I don't really see a pattern to express the general form of [itex]a_n[/itex] independently of the previous [itex]a_{n-p}[/itex]'s. Also, should I choose [itex]a_1=0[/itex] for convenience so that all odd terms would be worth 0?
By the way this might be similar to the hydrogen atom in a way yes, but the method to solve the DE isn't the same I guess because one usually "factors out" the behavior the DE when r tends to infinity (in the hydrogen atom case) while here I'm only interested in solving the PDE inside the region of a sphere of radius a or [itex]\tilde r[/itex]; so that it's not really the same path to take. I can be wrong of course.
Any idea on how to solve the DE (2)?

The solution is expressible as Bessel functions. See Eqn. 9.1.52 in Abramowitz and Stegan.

 

[fluidistic]

The solution is expressible as Bessel functions. See Eqn. 9.1.52 in Abramowitz and Stegan.

Thanks a lot for the reference (it's a free ebook).
So it states that the solution to the DE [itex]w''-\frac{2\nu -1 }{z}w'+\lambda ^2 w =0[/itex] is [itex]w =z^{\nu } C_{\nu} (\lambda z )[/itex]. Where C could be any linear combination of the Bessel functions of any kind.
My problem is... if I compare the DE to solve (2) or a slightly simplified but equivalent to DE: (3): [itex]R''+ \frac{2}{r}R'+R(\lambda - m )=0[/itex], one finds that [itex]\nu =-\frac{1}{2}[/itex]. So clearly nu isn't the smaller roots of the indicial equation as I thought.
What is very strange is that Wolfram alpha gives a totally different answer which seems to me unrelated to the Bessel function (http://www.wolframalpha.com/input/?i=y''+(2/x)y'+y*a=0) and it seems to confirm that [itex]c=-1[/itex] is the smaller root of the indicial equation.
According to the solution in the book of Abramowitz and Stegun the solution would be of the form of [itex]\frac{C_{-\frac{1}{2}}(r \sqrt {\lambda -m})}{\sqrt r}[/itex].
I'm totally confused.

 

[Mute]

Thanks a lot for the reference (it's a free ebook).
So it states that the solution to the DE [itex]w''-\frac{2\nu -1 }{z}w'+\lambda ^2 w =0[/itex] is [itex]w =z^{\nu } C_{\nu} (\lambda z )[/itex]. Where C could be any linear combination of the Bessel functions of any kind.
My problem is... if I compare the DE to solve (2) or a slightly simplified but equivalent to DE: (3): [itex]R''+ \frac{2}{r}R'+R(\lambda - m )=0[/itex], one finds that [itex]\nu =-\frac{1}{2}[/itex]. So clearly nu isn't the smaller roots of the indicial equation as I thought.
What is very strange is that Wolfram alpha gives a totally different answer which seems to me unrelated to the Bessel function (http://www.wolframalpha.com/input/?i=y''+(2/x)y'+y*a=0) and it seems to confirm that [itex]c=-1[/itex] is the smaller root of the indicial equation.
According to the solution in the book of Abramowitz and Stegun the solution would be of the form of [itex]\frac{C_{-\frac{1}{2}}(r \sqrt {\lambda -m})}{\sqrt r}[/itex].
I'm totally confused.

For half integer orders the Bessel functions are basically just combinations of signs and cosines divided by some power of x.

See http://en.wikipedia.org/wiki/Spherical_Bessel_function#Spherical_Bessel_functions:_jn.2C_yn

You should be able to mash the solution from Abramowitz and Stegun into looking like the one wolfram alpha gave you.