Imagine I have the following situation: inside a sphere containing mostly water, solute A is diffusing from the center to the end of the sphere:

enter image description here

At r=0, mass of A is continuously generated through the equation:

$\frac{\partial C(r=0,t)}{\partial t} = \dot{r}_{g}$

where $\dot{r}_{g}$ is a constant giving the mass of A formed per volume, per time.


$\frac{\partial C(r=0,t)}{\partial r} = 0$

due to symmetry reasons. In my opinion, at the end, this translates into the boundary condition $C(r=0,t) = \dot{r}_{g}\cdot t$

Moreover, at $t = 0$, $C(r,t=0) = 0$, because no A was still generated inside the sphere. So this is my initial condition.

From $r = 0$ to $r = R$, the diffusion equation for A applies:

$\frac{\partial C(r,t)}{\partial t} = D[\frac{\partial^2 C(r,t)}{\partial r^2} + \frac{2}{r}\frac{\partial C(r,t)}{\partial r}]$

What I would like to know is if there is a solution to this partial differential equation and how can we obtain it.

I tried to use separation of variables:

$C(r,t) = R^{*}(r)\cdot T(t)$

converting the original PDE into:

$\frac{1}{T}\frac{\partial T}{\partial t} = D[\frac{1}{R^{*}}\frac{\partial^2 R^{*}}{\partial r^2} + \frac{2}{r\cdot R^{*}}\frac{\partial R^{*}}{\partial r}]$

Although the following equation can be solved:

$\frac{1}{T}\frac{\partial T}{\partial t} = \lambda$

where $\lambda$ is the separation constant, I don't know how to proceed with the other equation:

$r\frac{\partial^2 R^{*}}{\partial r^2} + 2\frac{\partial R^{*}}{\partial r} - \frac{\lambda}{D}rR^{*} = 0$

in order to make the boundary and initial conditions that I specified valid.

  • $\begingroup$ Might Mathematics be better suited for this math question? $\endgroup$
    – Kyle Kanos
    Jan 10, 2017 at 18:39

1 Answer 1


Drop the unnecessary asterisk for simplicity:


After separation, drop the partials (separation separates the PDE into two ODEs), move $D$ to the lefthand side and use separation constant $-\lambda^2$, so you get two ODEs:

$$\frac{1}{DT}\frac{T'}{T}=-\lambda^2\tag{1}$$ $$rR''+2R'+\lambda^2rR=0\tag{2}$$

(Note: a sign error has been corrected here)

The latter is a Sturm-Liouville equation and it solves to:

$$R(r)=c_1\frac{\sin \lambda r}{r}+c_2\frac{\cos \lambda r}{r}\tag{3}$$

This offers the possibility of a fairly simple solution, at least for certain types of boundary conditions.

For example we can make it a bound problem by specifying: $$C(0,t)=C_0$$

Simplify a lot by a simple substitution:


So that the first boundary condition becomes homogeneous:


As second boundary condition we can choose:

$$C_r(R,t)=0\implies u_r(R,t)=0\implies R'(R)=0,$$

where $R$ is the radius of the sphere.

Using the first BC in $(3)$, we can see that necessarily:


or else $R(0)\to +\infty$.

So we have:

$$R(r)=c_1\frac{\sin \lambda r}{r}$$

Applying the second BC gives us:

$$c_1\Big(-\frac{\sin \lambda R}{R^2}+\lambda\frac{\cos \lambda R}{R}\Big)=0$$

Which, assuming $c_1\neq 0$ and $R\neq0$, reworks to:

$$\boxed{\lambda=\frac {\tan \lambda R}{R}}$$

This transcendental equation has no analytical solutions but does have an infinite number of numerical solutions: $\lambda_1$, $\lambda_2$, $\lambda_3$... $\lambda_n$. These are the eigenvalues of the originl PDE.

A particular solution of $R$ is thus:

$$R_n(r)=c_1\frac{\sin \lambda_n r}{r}$$

Equation $(1)$ simply solves to:

$$T_n(t)=c_3e^{-D\lambda_n^2 t}$$

The particular solution to $u(r,t)$ is thus:

$$u_n(r,t)=A_n\frac{\sin \lambda_n r}{r}\times e^{-D\lambda_n^2 t}$$

Due to the Superposition Principle the overal solution becomes:

$$u(r,t)=\displaystyle\sum_{n=1}^{+\infty}\Big(A_n\frac{\sin \lambda_n r}{r}\times e^{-D\lambda_n^2 t}\Big)$$

The coefficients $A_n$ are found from the initial condition for $r>0$ and using the substitution higher up:

$$C(r,0)=0\implies u(r,0)=-C_0$$

$$\implies -C_0=\displaystyle\sum_{n=1}^{+\infty}A_n\frac{\sin \lambda_n r}{r}$$

Using Fourier we then get the coefficients $A_n$ from:

$$A_n=\frac2R\int_0^R\Big(-C_0\frac{\sin \lambda_n r}{r}\Big)dr$$

Finally, back-substituting we get:


Note: a different and probably simpler solution caould be obtained by using the second BC:


Physically this corresponds to a sphere constantly being washed with fresh ($C=0$) solvent.

in order to make the boundary and initial conditions that I specified valid.

That's putting the horse before the cart: you need to choose the right boundary conditions first. Your ODE will then solve correctly.

And I'm not sure you've chosen the right BCs:

due to symmetry reasons. In my opinion, at the end, this translates into the boundary condition

$C(r=0,t) = \dot{r}_{g}\cdot t$

Not necessarily. You seem to have simply integrated:

$\frac{\partial C(0,t)}{\partial t} = \dot{r}_{g}$

But have omitted an integration constant.

I think you also lack a BC for $r=R$. Often in diffusion problems we set:

$$\Big(\frac{\partial C(r,t)}{\partial r}\Big)_{r=R}=0$$

Physically this means the solute cannot diffuse beyond $r=R$.

  • $\begingroup$ Thank you Gert for the answer. What I am thinking now, is that probably I will never be able to get a finite concentration at r = 0. Since this is a sphere and at the middle point the area is zero, concentration has to be infinite in order to have some flux: J = A.dC/dr $\endgroup$
    – Marcelo
    Jan 6, 2017 at 18:40
  • $\begingroup$ One way to avoid that is to specify a small central area ($r \leq r_0$) where $C=C_0$ (constant). Then change to a new variable: $u=C-C_0$, this way the boundary condition is homogeneous. $\endgroup$
    – Gert
    Jan 6, 2017 at 19:37
  • $\begingroup$ Your initial condition becomes $u(r,0)=-C_0$. When all is done and dusted your solution will be: $C(r,t)=C_0+u(r,t)$. $\endgroup$
    – Gert
    Jan 6, 2017 at 20:02
  • $\begingroup$ I was just looking around on the internet and what I found closest to the solution I was looking for is the Green's function that satisfies: $\frac{\partial C}{\partial t} - k(\frac{\partial^{2} C}{\partial x^{2}}+\frac{\partial^{2} C}{\partial y^{2}}+\frac{\partial^{2} C}{\partial z^{2}}) = \delta(r-r_{0})$. The solution to this equation involves the gaussian function (en.wikipedia.org/wiki/Green's_function). However, that solution is for a non-bounded domain and for a point source that happens at t = 0 $\endgroup$
    – Marcelo
    Jan 10, 2017 at 10:19
  • $\begingroup$ What I don't know is how to specify that $\frac{\partial C}{\partial r} = 0$ for $r = R$ (because diffusion does not proceed further that $r=R$) and for $r = 0$ (if function is C1 and for symmetry reasons). I believe this is done with the method of images. And also how to specify that point source is continuous. $\endgroup$
    – Marcelo
    Jan 10, 2017 at 10:26

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.