4.5.1 Spherical coordinates

4.5.1.1 [333] Chain reaction PDE

4.5.1.1 [333] Chain reaction PDE

problem number 333

Added May 7, 2019.

Assume \(\phi \) independence. Solve for \(u(r,\theta ,t) \)

\begin {align*} u_t = k ( \lambda u + \nabla ^2(u) ) \end {align*}

Where \(\nabla ^2(u)= \frac {1}{r^2} \frac {\partial }{\partial r}( r^2 u_r)+ \frac {1}{r^2 \sin \theta } \frac {\partial }{\partial \theta }(\sin \theta u_\theta )\) with \(k>0\).

Boundary conditions \(u(R,\theta ,t)=0\).

Mathematica

ClearAll["Global`*"]; 
U = u[r, theta, t]; 
pde = D[U, t] == k*(lambda*U + (1/r^2)*D[r^2*D[U, r], r] + (1/(r^2*Sin[theta]))*D[Sin[theta]*D[U, theta], theta]); 
bc  = u[R, theta, t] == 0; 
sol = AbsoluteTiming[TimeConstrained[DSolve[{pde, bc}, U, {r, theta, t}, Assumptions -> {k > 0, r < R}], 60*10]];
 

Failed

Maple

restart; 
U   := u(r,theta,t); 
pde := diff(U,t) = k*(lambda*U +1/r^2* diff(r^2*diff(U,r),r)+ 1/(r^2*sin(theta))*diff(sin(theta)*diff(U,theta),theta)); 
bc  := u(R,theta,t) = 0; 
cpu_time := timelimit(60*10,CodeTools[Usage](assign('sol',pdsolve([pde,bc], U) assuming k>0, r<R),output='realtime'));
 

\[u(r,theta,t) = 0\] trivial solution

Hand solution

Solve for \(u\left ( r,\theta ,t\right ) \) in spherical coordinates (assuming \(\phi \) independence) the chain reaction equation \(\frac {1}{k}u_{t}=\lambda u+\nabla ^{2}u\) with boundary conditions \(u\left ( R,\theta ,t\right ) =0\).\begin {align*} \frac {1}{k}u_{t} & =\lambda u+\nabla ^{2}u\\ & =\lambda u+\left ( \frac {1}{r^{2}}\frac {\partial }{\partial r}\left ( r^{2}u_{r}\right ) +\frac {1}{r^{2}\sin \theta }\frac {\partial }{\partial \theta }\left ( \sin \theta u_{\theta }\right ) \right ) \\ & =\lambda u+\frac {1}{r^{2}}\left ( 2ru_{r}+r^{2}u_{rr}\right ) +\frac {1}{r^{2}\sin \theta }\frac {\partial }{\partial \theta }\left ( \sin \theta u_{\theta }\right ) \end {align*}

Let \(u=R\left ( r\right ) \Theta \left ( \theta \right ) T\left ( t\right ) \). Substituting into the above gives\begin {align*} \frac {1}{k}T^{\prime }R\Theta & =\lambda TR\Theta +\frac {1}{r^{2}}\left ( 2rR^{\prime }\Theta T+r^{2}R^{\prime \prime }\Theta T\right ) +\frac {1}{r^{2}\sin \theta }\frac {\partial }{\partial \theta }\left ( \sin \theta \left ( \Theta ^{\prime }RT\right ) \right ) \\ \frac {1}{k}T^{\prime }R\Theta & =\lambda TR\Theta +\frac {2}{r}R^{\prime }\Theta T+R^{\prime \prime }\Theta T+\frac {RT}{r^{2}\sin \theta }\frac {\partial }{\partial \theta }\left ( \Theta ^{\prime }\sin \theta \right ) \end {align*}

Dividing by \(TR\Theta \) gives\[ \frac {1}{k}\frac {T^{\prime }}{T}=\lambda +\frac {2}{r}\frac {R^{\prime }}{R}+\frac {R^{\prime \prime }}{R}+\frac {1}{\Theta r^{2}\sin \theta }\frac {\partial }{\partial \theta }\left ( \Theta ^{\prime }\sin \theta \right ) \] The left side depends on \(t\) only and the right depends on \(r,\theta \) only. Let the separation variable be \(-n\). This gives the following 2 equations\begin {align} \frac {1}{k}\frac {T^{\prime }}{T} & =-n\tag {1}\\ \lambda +\frac {2}{r}\frac {R^{\prime }}{R}+\frac {R^{\prime \prime }}{R}+\frac {1}{\Theta r^{2}\sin \theta }\frac {\partial }{\partial \theta }\left ( \Theta ^{\prime }\sin \theta \right ) & =-n\tag {2} \end {align}

Now we consider (2). Multiplying both sides of (2) by \(r^{2}\) gives\begin {align*} \lambda r^{2}+2r\frac {R^{\prime }}{R}+r^{2}\frac {R^{\prime \prime }}{R}+\frac {1}{\Theta \sin \theta }\frac {\partial }{\partial \theta }\left ( \Theta ^{\prime }\sin \theta \right ) & =-nr^{2}\\ 2r\frac {R^{\prime }}{R}+r^{2}\frac {R^{\prime \prime }}{R}+r^{2}\left ( \lambda +n\right ) & =-\frac {1}{\Theta \sin \theta }\frac {\partial }{\partial \theta }\left ( \Theta ^{\prime }\sin \theta \right ) \end {align*}

The left side depends on \(r\) and the right side depends on \(\theta \). Let the separation variable be \(l\left ( l+1\right ) \) where \(l\) is integer. Hence we obtain the following two equations\begin {align} -\frac {1}{\Theta \sin \theta }\frac {\partial }{\partial \theta }\left ( \Theta ^{\prime }\sin \theta \right ) & =l\left ( l+1\right ) \tag {4}\\ 2r\frac {R^{\prime }}{R}+r^{2}\frac {R^{\prime \prime }}{R}+r^{2}\left ( \lambda +n\right ) & =l\left ( l+1\right ) \tag {5} \end {align}

Starting with (4)\begin {align*} \frac {\partial }{\partial \theta }\left ( \Theta ^{\prime }\sin \theta \right ) +l\left ( l+1\right ) \Theta \sin \theta & =0\\ \Theta ^{\prime \prime }\sin \theta +\Theta ^{\prime }\cos \theta +l\left ( l+1\right ) \Theta \sin \theta & =0 \end {align*}

Using the substitution \(z=\cos \theta \) the above becomes\[ \left ( 1-z^{2}\right ) \Theta ^{\prime \prime }-2z\Theta ^{\prime }+l\left ( l+1\right ) \Theta =0 \] This Legendre ODE. Solution is \(P_{l}\left ( \theta \right ) \). The other solution to the above ODE is ignored as not bounded. Now back to solving (5). Writing it as\begin {align} 2rR^{\prime }+r^{2}R^{\prime \prime }+r^{2}\left ( \lambda +n\right ) R & =l\left ( l+1\right ) R\nonumber \\ r^{2}R^{\prime \prime }+2rR^{\prime }+\left ( r^{2}\left ( \lambda +n\right ) -l\left ( l+1\right ) \right ) R & =0\tag {6} \end {align}

This can be converted to Bessel ODE using substitution. First let \(v=r\sqrt {\left ( \lambda +n\right ) }\). Then \(R^{\prime }\left ( r\right ) =\sqrt {\left ( \lambda +n\right ) }R^{\prime }\left ( v\right ) ,R^{\prime \prime }\left ( r\right ) =\left ( \lambda +n\right ) R^{\prime \prime }\left ( v\right ) \) and (6) becomes\begin {align} \left ( \lambda +n\right ) r^{2}R^{\prime \prime }\left ( v\right ) +2r\sqrt {\left ( \lambda +n\right ) }R^{\prime }\left ( v\right ) +\left ( r^{2}\left ( \lambda +n\right ) -l\left ( l+1\right ) \right ) R & =0\nonumber \\ v^{2}R^{\prime \prime }\left ( v\right ) +2vR^{\prime }\left ( v\right ) +\left ( v^{2}-l\left ( l+1\right ) \right ) R & =0\tag {7} \end {align}

Now, we apply second transformation \(R\left ( v\right ) =\frac {Z\left ( v\right ) }{\sqrt {v}}\) Then\begin {align*} R^{\prime }\left ( v\right ) & =\frac {Z^{\prime }\left ( v\right ) }{\sqrt {v}}-\frac {1}{2}Z\left ( v\right ) \frac {1}{v^{\frac {3}{2}}}\\ R^{\prime \prime }\left ( v\right ) & =\frac {Z^{\prime \prime }\left ( v\right ) }{\sqrt {v}}-\frac {1}{2}Z^{\prime }\left ( v\right ) \frac {1}{v^{\frac {3}{2}}}-\frac {1}{2}Z^{\prime }\left ( v\right ) \frac {1}{v^{\frac {3}{2}}}-\frac {1}{2}\left ( -\frac {3}{2}\right ) Z\left ( v\right ) \frac {1}{v^{\frac {5}{2}}}\\ & =\frac {Z^{\prime \prime }\left ( v\right ) }{\sqrt {v}}-Z^{\prime }\left ( v\right ) \frac {1}{r^{\frac {3}{2}}}+\frac {3}{4}Z\left ( v\right ) \frac {1}{v^{\frac {5}{2}}} \end {align*}

Hence (7) becomes\[ v^{2}\left ( \frac {Z^{\prime \prime }\left ( v\right ) }{\sqrt {v}}-Z^{\prime }\left ( v\right ) \frac {1}{r^{\frac {3}{2}}}+\frac {3}{4}Z\left ( v\right ) \frac {1}{v^{\frac {5}{2}}}\right ) +2v\left ( \frac {Z^{\prime }\left ( v\right ) }{\sqrt {v}}-\frac {1}{2}Z\left ( v\right ) \frac {1}{v^{\frac {3}{2}}}\right ) +\left ( v^{2}-l\left ( l+1\right ) \right ) \frac {Z\left ( v\right ) }{\sqrt {v}}=0 \] Multiplying by \(\sqrt {v}\) gives\begin {align*} v^{2}\left ( Z^{\prime \prime }\left ( v\right ) -Z^{\prime }\left ( v\right ) \frac {1}{v}+\frac {3}{4}Z\left ( v\right ) \frac {1}{v^{2}}\right ) +2v\left ( Z^{\prime }\left ( v\right ) -\frac {1}{2}Z\left ( v\right ) \frac {1}{v}\right ) +\left ( v^{2}-l\left ( l+1\right ) \right ) Z\left ( v\right ) & =0\\ \left ( v^{2}Z^{\prime \prime }\left ( v\right ) -vZ^{\prime }\left ( v\right ) +\frac {3}{4}Z\left ( v\right ) \right ) +\left ( 2vZ^{\prime }\left ( v\right ) -Z\left ( v\right ) \right ) +\left ( v^{2}-l\left ( l+1\right ) \right ) Z\left ( v\right ) & =0\\ v^{2}Z^{\prime \prime }\left ( v\right ) +vZ^{\prime }\left ( v\right ) +\frac {1}{4}Z\left ( v\right ) +\left ( v^{2}-l\left ( l+1\right ) \right ) Z\left ( v\right ) & =0\\ v^{2}Z^{\prime \prime }\left ( v\right ) +vZ^{\prime }\left ( v\right ) +\left ( v^{2}-l\left ( l+1\right ) +\frac {1}{4}\right ) Z\left ( v\right ) & =0\\ v^{2}Z^{\prime \prime }\left ( v\right ) +vZ^{\prime }\left ( v\right ) +\left ( v^{2}-l^{2}+l+\frac {1}{4}\right ) Z\left ( v\right ) & =0\\ v^{2}Z^{\prime \prime }\left ( v\right ) +vZ^{\prime }\left ( v\right ) +\left ( v^{2}-\left ( l+\frac {1}{2}\right ) ^{2}\right ) Z\left ( v\right ) & =0 \end {align*}

This is now in standard Bessel ODE form. Comparing it to \(v^{2}Z^{\prime \prime }\left ( v\right ) +vZ^{\prime }\left ( v\right ) +\left ( v^{2}-d^{2}\right ) Z\left ( v\right ) =0\) shows the order is \(d=l+\frac {1}{2}\). The solutions are \[ Z\left ( v\right ) =c_{1}J_{l+\frac {1}{2}}\left ( v\right ) +c_{2}Y_{l+\frac {1}{2}}\left ( v\right ) \] But \(R\left ( v\right ) =\frac {Z\left ( v\right ) }{\sqrt {v}}\), hence\[ R\left ( v\right ) =c_{1}\frac {J_{l+\frac {1}{2}}\left ( v\right ) }{\sqrt {v}}+c_{2}\frac {Y_{l+\frac {1}{2}}\left ( v\right ) }{\sqrt {v}}\] But \(v=r\sqrt {\lambda +n}\) then above becomes\[ R\left ( r\right ) =c_{1}\frac {J_{l+\frac {1}{2}}\left ( r\sqrt {\left ( \lambda +n\right ) }\right ) }{\sqrt {r\sqrt {\left ( \lambda +n\right ) }}}+c_{2}\frac {Y_{l+\frac {1}{2}}\left ( r\sqrt {\left ( \lambda +n\right ) }\right ) }{\sqrt {r\sqrt {\left ( \lambda +n\right ) }}}\] The solution assumed bounded at \(r=0\) hence \(c_{2}=0\) and the above becomes\[ R\left ( r\right ) =c_{1}\frac {J_{l+\frac {1}{2}}\left ( r\sqrt {\left ( \lambda +n\right ) }\right ) }{\sqrt {r\sqrt {\left ( \lambda +n\right ) }}}\] Let \(m^{2}=\left ( \lambda +n\right ) \), hence\begin {align*} R\left ( r\right ) & =c_{1}\frac {J_{l+\frac {1}{2}}\left ( mr\right ) }{\sqrt {mr}}\\ & =j_{l}\left ( mr\right ) \end {align*}

where \(j_{l}\left ( mr\right ) \) are the spherical Bessel functions. Boundary conditions at \(r=R\) gives \[ j_{l}\left ( mR\right ) =0 \] Hence \(mR\) or \(\sqrt {\lambda +n}R\) are the zeros of spherical Bessel functions \(j_{l}\left ( mR\right ) \). There are infinite zeros for each \(l\). Let the \(v^{th}\) zero of \(j_{l}\) be called \(Z_{l,v}\). Then \(mR_{l,v}=Z_{l,v}\) or \(\sqrt {\lambda +n}_{l,v}=\frac {Z_{l,v}}{R}\). or \[ n=\left ( \frac {Z_{l,v}}{R}\right ) ^{2}-\lambda \] The solution to the time ODE is therefore\begin {align*} T_{l,v} & =A_{l,v}e^{-nkt}\\ & =A_{l,v}e^{-\left ( \left ( \frac {Z_{l,v}}{R}\right ) ^{2}-\lambda \right ) kt} \end {align*}

Hence the complete solution is\begin {align*} u\left ( r,\theta ,t\right ) & =e^{-i\alpha kt}P_{l}\left ( \theta \right ) j_{l}\left ( mr\right ) \\ & =\sum _{l=1}^{\infty }\sum _{v=0}^{\infty }A_{l,v}e^{-nkt}P_{l}\left ( \theta \right ) j_{l}\left ( \frac {Z_{l,v}}{R}r\right ) \end {align*}

\(A_{l.v}\) constants still need to be found from initial conditions. For each \(l\), we have infinite sum over all \(v^{\prime }s\) zeros of \(j_{l}\).