# Orthogonality of 3d diffusion equation solution base functions

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.

### Background

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.

### Question

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?

• 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. – reuns Jan 14 '16 at 12:54
• @user1952009 So you are saying that $\tan(\lambda_i R)=\lambda_i R$ for any boundary conditions imposed? – Robin Jan 14 '16 at 13:11
• 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. – 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.$$

• 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? – Robin Jan 15 '16 at 12:48
• @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! – Disintegrating By Parts Jan 15 '16 at 12:52