# Numerical gradient in spherical coordinates

Assume that we have a function $u$ defined in a ball in a discrete way: we know only the values of $u$ in the nodes $(i,j,k)$ of spherical grid, where $i$ is a radius coordinate, $j$ is a coordinate for angle $\varphi$, $k$ is a coordinate for angle $\psi$.

Consider a vector-function $$\nabla u_{i,j,k}=\left(\frac{\partial u}{\partial r}_{i,j,k},\frac{1}{r_i\sin\psi_k}\frac{\partial u}{\partial \varphi}_{i,j,k},\frac{1}{r_i}\frac{\partial u}{\partial \psi}_{i,j,k}\right)-$$ gradient of $u$.

I need to know the values of $\nabla u_{i,j,k}$ on z-axis in cartesian coordinates, which corresponds to $\psi=0$ -- axis in spherical coordinates, but we can not use the formula above, because in case $\psi=0$ the second term turns to infinity.

Actually, we can find the values of $\frac{\partial u}{\partial z}$ with a help of the formula of numerical derivative, but we have a problem with finding $\frac{\partial u}{\partial x}$,$\frac{\partial u}{\partial y}$, because the grid is not rectangular. Could you help me with this stuff and advise me what to do?

• Could you not calculate the numerical derivative in the $(r,\varphi,\psi)$ coordinate system? This should give you a vector which you can then project onto the $x$- $y$- and $z$-axes to give you the $x$- $y$- and $z$-components of the gradient Feb 18, 2014 at 16:48
• No, we can not calculate numerical derivative, because for $\psi=0$ $\varphi$-angle is not defined and we have singularity in formula.
– cool
Feb 18, 2014 at 17:09
• There is likely a solution using quaternions, but its kinda tricky to figure out what exactly that looks like wrt. your nomencalture. The effect you are trying to fight is commonly referred to as "gimbal-lock" Feb 19, 2014 at 6:20
• I found one variant to deal with this stuff correctly : here we can use least square method for gradient reconstruction, but I did not find exact explanation, how to use it
– cool
Feb 19, 2014 at 12:15
• Perhaps it's the case - for moderately well behaved $u$$- that \lim_{\psi\to 0} \frac{1}{\sin \psi} \frac{\partial u}{\partial \phi}=0. But I don't think the middle component of your vector has meaning, so if you're plugging this into a differential equation in polar coordinates the differential equation might also give no weight to the middle coordinate. – user121 Feb 22, 2014 at 22:36 ## 1 Answer There are 3 ways to avoid this situation, but before use one must check if this way is suitable due to computation error: 1) Green-Gauss cell method: here the definition of gradient is used:$$\nabla u_i \approx\frac{1}{V_{i}} \int\limits_{\partial V_i}u d\overline{S}\approx \sum\limits_{k=1}^{n}{u_{f_k} S_k \overline{n}_k} , $$where k - numbers of neighbours of cell V_{i} 2) Least squares method: the error$$\sum\limits_{k=1}^{n}{\frac{1}{d_{ik}}E_{i,k}^{2}}, E_{i,k}=\nabla u_i \cdot\Delta r_{i,k}+u_i-u_k$$must be minimized, hence we get the components of$\nabla u_i\$

3) Interpolation method. The value of gradient is interpolated from the values of gradient vector-function.