Finite difference method (diffusion equation) for 3D spherical case

There is one 3D-diffusion IBVP for sphere (Dirichlet problem with zero conditions on the surface of sphere). The equation has the following view: $$\frac{\partial u}{\partial t}=\operatorname{div}\left(A(r,\varphi,\psi)\nabla u(r,\varphi,\psi,t)\right)+f(r,\varphi,\psi,t) ,$$ where the matrix $A$ is skew with variable coefficients: $$\left(\matrix{1 \quad a \quad c\\-a \quad 1 \quad b\\-c\quad -b \quad 1 }\right)$$ I've written a numerical scheme for this problem, but the program works strangely: if all non-diagonal coefficients are zeros, everything is cool, but when at least one coefficient is non-zero, my scheme begins to show great error in neighbourhood of points $\varphi=0$ and $\varphi=\pi$. I've checked program many times, moreover, I've shown it to expert --he said, that program code is correct. He also said that this problem with big errors is connected with singularities in coefficients of terms $\frac{\partial u}{\partial r},\frac{\partial u}{\partial \varphi},\frac{\partial u}{\partial \psi}$. Has anyone faced with similar problems? How can I avoid those singularities? Here is one example of error. The surfaces' equation is $u(r,\varphi,\psi,t)=\cos(t)(1-r^2)$

[upd] Here is my numerical scheme: $$\frac{u^{n+1}_{i,j,k}-u^{n}_{i,j,k}}{\Delta t}=AA_{i,j,k}\frac{u^{n+1}_{i+1,j,k}-u^{n+1}_{i-1,j,k}}{2\Delta r}+BB_{i,j,k}\frac{u^{n+1}_{i,j+1,k}-u^{n+1}_{i,j-1,k}}{2\Delta \varphi}+CC_{i,j,k}\frac{u^{n+1}_{i,j,k+1}-u^{n+1}_{i,k,k-1}}{2\Delta \psi}+\frac{u^{n+1}_{i+1,j,k}-2u^{n+1}_{i,j,k}+u^{n+1}_{i-1,j,k}}{\Delta r^2}+\frac{u^{n+1}_{i,j+1,k}-2u^{n+1}_{i,j,k}+u^{n+1}_{i,j-1,k}}{r_{i}^{2}\Delta \varphi^2}+\frac{u^{n+1}_{i,j,k+1}-2u^{n+1}_{i,j,k}+u^{n+1}_{i,j,k-1}}{r_{i}^{2}\sin^{2}\varphi_{j}\Delta \psi^2}+f_{i,j,k}$$ Coefficients $AA,BB,CC$ are the following: $$AA_{i,j,k}=\frac{2}{r_{i}}+a_{i,j,k}\frac{\cot \varphi_{j}}{r_{i}}+\frac{1}{r_{i}}\left(\frac{\partial a}{\partial \varphi}\right)_{i,j,k}-\frac{1}{r_{i}\sin\varphi_{j}}\left(\frac{\partial c}{\partial \psi}\right)_{i,j,k}$$ $$BB_{i,j,k}=-\frac{a_{i,j,k}}{r_{i}^2}+\frac{\cot\varphi_{j}}{r_{i}^2}-\frac{1}{r_{i}}\left(\frac{\partial a}{\partial r}\right)_{i,j,k}+\frac{1}{r_{i}^2\sin\varphi_{j}}\left(\frac{\partial b}{\partial \psi}\right)_{i,j,k}$$ $$CC_{i,j,k}=\frac{c_{i,j,k}}{r_{i}^2\sin\varphi_{j}}+\frac{1}{r_{i}\sin\varphi_{j}}\left(\frac{\partial c}{\partial r}\right)_{i,j,k}-\frac{1}{r_{i}^2\sin\varphi_{j}}\left(\frac{\partial b}{\partial \varphi}\right)_{i,j,k}$$

• How are you sampling the space in your scheme? Is it just evenly dividing the ranges of the 3 coordinates? – Phil H Nov 28 '13 at 18:19
• I use uniform grid for angles $\varphi,\psi$. For $r$ i used both logarithmic and uniform grids, but the results were the same – cool Nov 28 '13 at 19:21
• How does the error vary as you increase the number of samples along $r$ - does it get better to a point and then worse? – Phil H Nov 29 '13 at 10:45
• The error increases with increasing of the number of samples along $r$, but there is some upper bound, you can see it on the picture related to this topic, such picture is similar for any number of samples. – cool Nov 29 '13 at 16:55
• What form do you use for the divergence and the gradient in your expression? Since you don't show the program or talk about any of the details of your numerical scheme, there is little we can suggest here. In general, you will get into exactly the kind of trouble you show if you use finite differences, but things may work if you use an appropriate finite element formulation. – Wolfgang Bangerth Dec 2 '13 at 3:08