I have a question regarding the orthogonality of the base functions for the solution of the 3d diffusion equation. My motivation is that I want to calculate the diffusion process inside a 3d hollow sphere with rotational invariant initial conditions and no-flux conditions at the surface. Since I have not found the solution process anywhere in one place, I'm posting a sketch of the process here as introduction.


The diffusion equation in 3d for rotational invariant geometries is $$ \frac{\partial c}{\partial t} = D \left( \frac{\partial^2 c}{\partial r^2} + \frac 2 r \frac{\partial c}{\partial r} \right) $$ with concentration $c(x,t)$, diffusion constant $D>0$, time $t \ge 0$ and radial coordinate $r \ge 0$. Using $u=r \, c$ this can be transformed to the 1d diffusion equation $$ \frac{\partial u}{\partial t} = D \frac{\partial^2 u}{\partial r^2} \quad . $$ This can be solved (amongst other) by separation of variables with the solution being $$ u(r,t)=\rho_0 + \rho_1 r + \sum\limits_{k=0}^\infty \exp\left(-D \lambda_k^2 t \right) \big( \kappa_k \cos\left( \lambda_k r \right) + \sigma_k \sin\left( \lambda_k r \right) \big) $$ with $\rho_0$, $\rho_1$, $\lambda_k$, $\kappa_k$ and $\sigma_k$ calculated to match the boundary and initial conditions (with $\rho_0$ being ambiguous to $\kappa_i$ for $\lambda_i=0$).

Since $u=r \, c$, $u=0$ for $r=0$ at all times, leading to $\rho_0=0$ and $\kappa_k=0$. The zero flux condition means $\partial c/\partial r = 0$ at $r=R$, and using the derivative of $c=u/r$, $$ \frac{\partial u}{\partial r} - \frac{u}{R} = 0 \quad (r=R) $$ is the second boundary condition. Evaluating this for the solution yields $$ \tan\left(\lambda_k R\right) = \lambda_k R $$ which I can solve numerically for the $\lambda_k$ ($\lambda_0 = 0$, $\lambda_1 \approx 4.5/R$, $\lambda_2 \approx 7.7/R$ and so on).

Since $\lambda_0=0$, the position-related base functions can be defined as $$ X_i(r):=\begin{cases} r/R & i = 0 \\ \sin(\lambda_i r) & i \ge 1 \end{cases} $$ The $X_i(r)$ are orthorgonal, i.e. $$ \frac 2 R \int\limits_0^R X_i(r) X_j(r) = \begin{cases} \frac 2 3 & i=j=0 \\ 1 - \frac{1}{2 \lambda_i R} \sin(2 \lambda_i R) & i=j\ge1 \\ 0 & i \neq j \end{cases} \quad , $$ so the $\sigma_i$ can be easily computed from initial conditions $c(r,t=0)=c_0(r)$ and the problem is solved.


From my (engineering) point of view, the orthogonality of the solution's base functions heavily depends on the boundary conditions, since $$ \int\limits_0^R r \, \sin(\lambda_i r) \, \mathrm{d}r = \frac{1}{\lambda_i^2} \cos(\lambda_i R) \big(\tan(\lambda_i R) - \lambda_i R \big) $$ and $$ \int\limits_0^R \sin(\lambda_i r) \sin(\lambda_j r) \, \mathrm{d}r = \frac{ R \lambda_i \lambda_j } {\lambda_i^2-\lambda_j^2} \cos(\lambda_i R) \cos(\lambda_j R) \left( \frac{ \tan(\lambda_i R) }{\lambda_i R} - \frac{\tan(\lambda_j R)}{\lambda_j R} \right) \quad . $$ So in order to become orthorgonal, either $\lambda_i R = \pi k$ with $k \in \mathbb{Z}$, or $\tan(\lambda_i R) = \lambda_i R$ (or, when ignoring $X_0$, $\tan(\lambda_i R) / \lambda_i R = \mathrm{const}$).

So my question is: What is the pattern behind this? Or is it just mere luck that my boundary conditions evaluate in a way that the base functions are orthorgonal?

  • $\begingroup$ no, the $\lambda_i$ depend on $R$ so the basis is still orthogonal for any chosen $R$. when the space is unbounded, the possible $\lambda$ are continuous and we can still find a subset of them being orthogonal (the real ones). any homogeneous linear differential equation has a simple eigen-mode orthogonal basis subset of the complex exponentials. $\endgroup$ – reuns Jan 14 '16 at 12:54
  • $\begingroup$ @user1952009 So you are saying that $\tan(\lambda_i R)=\lambda_i R$ for any boundary conditions imposed? $\endgroup$ – Robin Jan 14 '16 at 13:11
  • $\begingroup$ no but with a different $R$ it would be still orthogonal. if the shape of a boundary $u(f(x),0) = 0$ we won't get single frequency as solution, nor in the in-homogeneous case neither. $\endgroup$ – reuns Jan 14 '16 at 13:24

There is a natural weight factor that should be incorporated into the analysis, which is $r^2$ coming from the Jacobian of the spherical coordinate system. The corresponding inner product is $$ (f,g)_{r^2} = \int_{0}^{R}f(r)\overline{g(r)}r^2dr. $$ Using this inner product, the operator $$ Lf = \frac{d^2f}{d r^2}+\frac{2}{r}\frac{df}{dr} $$ is selfadjoint, once standard endpoint conditions are imposed. This is because $$ Lf = \frac{1}{r^2}\frac{d}{dr}\left(r^2\frac{df}{dr}\right), $$ which gives \begin{align} \int (Lf)g r^{2}dr & = \int \frac{d}{dr}\left(r^2\frac{df}{dr}\right)gdr \\ & = r^{2}\frac{df}{dr}g-\int r^{2}\frac{df}{dr}\frac{dg}{dr}dr \\ & = r^{2}\frac{df}{dr}g-r^{2}f\frac{dg}{dr}+\int f\frac{d}{dr}\left(r^2\frac{dg}{dr}\right)dr \\ & = r^{2}\{ f'g-fg' \}+\int f(Lg)r^2dr \end{align} Therefore, $$ (Lf,g)_{r^2} - (f,Lg)_{r^2} = r^2\{f'\overline{g}-f\overline{g}'\}|_{0}^{R} $$ Once proper endpoint conditions at the endpoints for the functions $f,g$ in the domain, $$ (Lf,g)_{r^2}=(f,Lg)_{r^2}. $$ When dealing with a weighted space like this, $rf \in L^2(0,R)$ with unit weight iff $f\in L^2_{r^2}(0,R)$, which motivates the substitution of the original equation that leads to an unweighted space. The resulting operator is then selfadjoint in $L^2(0,R)$ with unit weight. The transformation $f\in L^2_{r^2}(0,R)\mapsto rf \in L^2(0,R)$ is a unitary substitution between the Hilbert spaces $L^2_{r^2}(0,R)$ and $L^2(0,R)$. The resulting operator is also selfadjoint in $L^2(0,R)$ once the endpoint conditions are translated. The new ODE in $L^2(0,R)$ is sandwiched between multiplication by $1/r$ and multiplication by $r$. The new operator $M : \mathcal{D}(M)\subset L^2(0,R)\rightarrow L^2(0,R)$ is $$ M = rL(\frac{1}{r}f),\;\;\; f \in L^2(0,R) $$ Explicitly, \begin{align} Mf & = r\left[\frac{1}{r^2}\frac{d}{dr}\left(r^2\frac{d(\frac{1}{r}f)}{dr}\right)\right] \\ & = r\left[\frac{1}{r^2}\frac{d}{dr}\left(r\frac{df}{dr}-f\right)\right] \\ & = r\left[\frac{1}{r^2}\left(r\frac{d^2f}{dr^2}+\frac{df}{dr}-\frac{df}{dr}\right)\right] \\ & = \frac{d^2f}{dr^2}. \end{align} This also exposes the possible endpoint conditions at $r=0$. In the transformed coordinates, the left endpoint condition functionals are $\Phi(f)=f(0+)$ and $\Psi(f)=f'(0+)$. That is, $\{\Phi,\Psi\}$ is a basis of left endpoint functionals. These functions may be applied to functions on the original space to obtain the transformed basis of endpoint functionals $\{\Phi_{r^2},\Psi_{r^2}\}$ defined by \begin{align} \Phi_{r^2}(g) & = \lim_{r\rightarrow 0}(rg(r)) \\ \Psi_{r^2}(g) & = \lim_{r\rightarrow 0}(\frac{d}{dr}(rg)) =\lim_{r\rightarrow 0}(rg'(r)+g(r)). \end{align} To obtain a selfadjoint problem for $L^2(0,R)$ with separated endpoint conditions, you can impose one endpoint homogeneous condition at each endpoint. This is standard. The only natural condition at $r=0$ is $f(0+)=0$ because $g,Lg \in L^2_{r^2}$ implies that $\Phi_{r^2}(g)$ exists and $$ g(r) \sim \Phi_{r^2}(g)\frac{1}{r}. $$ Because $L$ is non-singular at $r=R$, the evaluation functionals for $g$ and $g'$ at $r=R$ is a basis of endpoint functionals, and these are equivalent to the the evaluations of $f$ and $f'$ at $r=R$. Let $L_{\alpha}$ denote the restriction of $L$ to the domain with $f,Lf\in L^2_{r^2}(0,R)$ satisfying endpoint conditions $$ \lim_{r\rightarrow 0}(rf(r))=0,\;\;\; \cos\alpha f(R)+\sin\alpha f'(R)=0. $$ Then $(Lf,g)_{r^2} = (f,Lg)_{r^2}$ for all $f,g \in \mathcal{D}(L_{\alpha})$. It is such symmetry that guarantees the orthogonality of eigenfunctions with different eigenvalues, just as it does for Hermitian matrices; more precisely, if $Lf=\lambda f$, $Lg=\mu g$ with $f,g\in\mathcal{D}(L_{\alpha})$ and $\lambda\ne \mu$, then $\mu,\lambda$ are real because $$ (\lambda-\overline{\lambda})\|f\|_{r^2}=(Lf,f)_{r^2}-(f,Lf)_{r^{2}}=0,\\ (\mu-\overline{\mu})\|g\|_{r^2}=(Lg,g)_{r^2}-(g,Lg)_{r^2}=0, $$ and $(f,g)_{r^2}=0$ because $$ (\lambda-\mu)(f,g)_{r^2}=(Lf,g)_{r^2}-(f,Lg)_{r^2}=0. $$

  • $\begingroup$ Thanks! Do I understand your answer correctly, that if the endpoint conditions ensure orthogonality in the original $L^2_{r^2}$ space (where $c$ lives), then they also ensure orthogonality in the transformed $L^2$ space (where $u$ lives)? Is that the core of your answer? $\endgroup$ – Robin Jan 15 '16 at 12:48
  • $\begingroup$ @Robin : Yes, unitary operators preserve inner products. That is $U : H_1 \rightarrow H_2$ is unitary iff $(Uh,Uk)_{2}=(h,k)_{1}$ for all $h,k\in H_1$. And the technique I've given you is a general technique for transforming to an unweighted Hilbert space. Please note that Fourier came up with this substitution method for the cooling problem of the sphere in his original ~1810 manuscript on Heat Conduction. He found the eigenvalues, too, as a transcendtal equation! $\endgroup$ – Disintegrating By Parts Jan 15 '16 at 12:52

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.