# Solving the diffusion equation in a ball with Neumann BC

I was hoping to ask for finding $\lambda_n$ in the following problem, I need to jog my memory as to how to solve these using fourier series. I am only interested in $\lambda_n$ as these terms will represent my decay rates.

Suppose we have a ball, radius $r_0$, with the diffusion equation $$\frac{\partial c(r,t)}{\partial t} = \frac{D}{r^2} \frac{\partial }{\partial r} \left( r^2 \frac{\partial }{\partial r} c(r,t) \right), \quad r\in[0,r_0],t\geq 0$$

with neumann boundary conditions $$D\frac{\partial c}{\partial r} = -pc ,\quad \text{on}\ \ r=r_0$$

with initial conditions $$c(r,0)= \begin{cases} c_0 & r\in[0,R] \\ 0 & r\in(R,r_0] \end{cases}$$

where $D,p>0$.

If we let $$c = R(r)T(t)$$

then we find $$\frac{1}{D}\frac{T'}{T} = \frac{\partial^2 R}{\partial r^2} + \frac{2}{r}\frac{\partial^2 R}{\partial r^2} = -\lambda^2$$

Therefore $T(t)=\exp(-D \lambda^2 t)$. I'm interested in describing the half-life of this system by examining the decay rates of the series solution to it.

Your last displayed equation is wrong. You get $$r^2R''+2rR'+\lambda^2r^2R=0\ .$$ Introduce a new unknown function $Q(r):=r\>R(r)$. You then obtain $Q'=R+rR'$, $Q''=2R'+rR''$ and therefore $$r(Q''+\lambda^2 Q)=r^2R''+2rR'+\lambda^2 r^2R=0\ .$$ It follows that $$Q(r)=a\cos(\lambda r)+b\sin(\lambda r)\ .$$ Since $R(r)={1\over r}Q(r)$ one necessarily has $a=0$, so that we are left with $R(r)={1\over r}\sin(\lambda r)$. We now have to determine the $\lambda$s for which the boundary condition $$D\>R'(r_0)=-p\>R(r_0)\tag{1}$$ is fulfilled. To this end we compute $$R'(r)=-{1\over r^2}\sin(\lambda r)+{\lambda\over r}\cos(\lambda r)\ .$$ Plugging this and $R(r)$ into $(1)$ we obtain the transcendental equation (typical for this kind of problem) $$\tan\alpha={D\over D-r_0p}\>\alpha$$for the dimensionless unknown $\alpha:=\lambda r_0$. Now determine numerically the smallest positive solution of this equation, and you are done.
• Thanks for the update! This is really helpful. So the smallest +ve solution for $\lambda$ will give the power of the dominant exponential, and hence the half-life? Commented Apr 3, 2017 at 13:37
• Replacing a $\lambda$ by its negative just multiplies the corresponding "basis function" $R_\lambda(r)$ by $-1$, hence gives nothing new. Commented Apr 3, 2017 at 14:04