# 3d-diffusion equation in spherical coordinates (numerical), boundary problem

There is one boundary problem

$$\frac{\partial u}{\partial t}= \operatorname{div}\left(a^2 E \nabla u\left(r,\varphi,\psi \right) \right)$$ in a ball $$B_{1}(0)=\left\{x \in \mathbb{R^3}: \left\| x \right\| \le 1 \right\}$$

(boundary and initial values are now not principal). It's obvious that this equation can be expanded to: \begin{gather} \frac{\partial u}{\partial t}=a^2\left(\frac{1}{r^{2}} \frac{\partial}{\partial r}\left(r^2 \frac{\partial u}{\partial r}\right)+\frac{1}{r\sin\varphi}\frac{\partial }{\partial \varphi}\left(\frac{\sin\varphi}{r} \frac{\partial u}{\partial \varphi} \right)+\frac{1}{r\sin\varphi}\frac{\partial}{\partial \psi}\left(\frac{1}{r\sin\varphi}\frac{\partial u}{\partial \psi} \right)\right).\end{gather} This boundary problem is solved numerically. So, we write an implicit finite-difference scheme

$$\frac{u^{n+1}_{i,j,k}-u^{n}_{i,j,k}}{\Delta t}=a^2 \left( \frac{2}{r_{1}+(i-1)\Delta r}\frac{u^{n+1}_{i+1,j,k}-u^{n+1}_{i-1,j,k}}{2\Delta r}+\frac{u^{n+1}_{i+1,j,k}-2u^{n+1}_{i,j,k}+u^{n+1}_{i-1,j,k}}{\Delta r^2}\right) +$$ $$+a^2\left(\frac{\cot \left(\varphi_{1}+(j-1)\Delta \varphi \right) ) }{\left(r_{1}+(i-1)\Delta r\right)^2}\frac{u^{n+1}_{i,j+1,k}-u^{n+1}_{i,j-1,k}}{2\Delta \varphi}+\frac{1}{\left(r_{1}+(i-1)\Delta r\right)^2}\frac{u^{n+1}_{i,j+1,k}-2u^{n+1}_{i,j,k}+u^{n+1}_{i,j-1,k}}{\Delta \varphi^{2}} \right) +$$ $$+a^2\left(\frac{1}{\left(r_{1}+(i-1)\Delta r\right)^2 \sin^2\left(\varphi_{1}+(j-1)\Delta \varphi \right)}\frac{u^{n+1}_{i,j,k+1}-2u^{n+1}_{i,j,k}+u^{n+1}_{i,j,k-1}}{\Delta \psi^2} \right).$$ The main problem appears in realization of this scheme. The values $r_{1},\varphi_{1}$ -- can be equal to zero in some moments of evaluation ,that's why we face with singularities during evaluation. Could you help me, what should I do with this scheme (perhaps, change it in the moments with singularities) in order to avoid singularities during evaluation?

• Depending on the boundary conditions, you may be able to identify $i = 0$ with $i = N-1$. What are the boundary conditions? Also, can you explain to me where the $i-1$ terms in the denominator come from in the discretization? – Eric Kightley Jul 23 '13 at 10:28
• The boundary conditions are Dirichlet conditions: $u(\partial D)=0$. They are defined on $r_{n}$'s layer -- on the surface of $B_{1}(0)$. $i−1$ comes from the next fact: there is a summand 1r...∂∂.. in the equation, in the process of discretization we replace r by the value $r_1+(i−1)\Delta r$, where i -- index of the position in a ball (coordinate). When $i=1$, $r=r_{1}=0$ - here we are in the center of the ball. In other cases the same situation occurs with $φ_{1}=0$ – cool Jul 29 '13 at 18:57