# Wrong result for gradient in spherical coordinates [closed]

Hello,

I use the following code to calculate the divergence of a scalar field in spherical coordinates:

E.<x,y,z> = EuclideanSpace()
CA = E.cartesian_coordinates(); CA
BP.<r, θ, ϕ> = E.chart(r"r:[0,+oo)  θ:[0,pi]  ϕ:[0,2*pi]" )
BP_to_CA = BP.transition_map(CA, [r * sin(θ) * cos(ϕ), r * sin(θ) * sin(ϕ), r * cos(θ)])
g = E.metric()

E.set_default_chart(BP)
E.set_default_frame(BP.frame())
show( g.display() )
show( g.display_comp())

f = E.scalar_field({BP: function('F')(r, θ, ϕ)}, name='f')
f.display(BP)
f.laplacian().expr(BP).factor()

But the result is wrong:

It must be:

What am I doing wrong?

By the way: how can I define a vector field using BP ?

edit retag reopen merge delete

### Closed for the following reason the question is answered, right answer was accepted by pfeifhns close date 2023-11-26 19:02:24.431536

Sort by » oldest newest most voted

Thank you very much for your answer. You are right, I expected the answer in e_r, e_theta and e_phi. Is there a way to get the gradient in e_r, e_theta and e_phi?

The other problem I have is that I don't know how to define a vector space in the BP-chart. I also would like to calculate divergence and curl.

more

I have edited my answer to show how to a new vector frame and express the gradient in it.

( 2023-11-26 15:45:56 +0100 )edit

The answer provided by Sage for the gradient of f is correct:

grad(f) = d(F)/dr ∂/∂r + d(F)/dθ/r^2 ∂/∂θ + d(F)/dϕ/(r^2*sin(θ)^2) ∂/∂ϕ

Remember that (∂/∂r, ∂/∂θ, ∂/∂ϕ) is not an orthonormal basis. The standard orthonormal basis associated with spherical coordinates is $(e_r, e_\theta, e_\phi)$ with $e_r := \frac{\partial}{\partial r}$, $e_\theta := \frac{1}{r} \frac{\partial}{\partial \theta}$, $e_\phi := \frac{1}{r\sin\theta}\frac{\partial}{\partial \phi}$ and apparently you were expecting an answer in that basis.

You might be interested by taking a look at these tutorials:

https://nbviewer.org/github/sagemanif...

https://nbviewer.org/github/sagemanif...

https://nbviewer.org/github/sagemanif...

Update (26 nov. 2023): here is how to set up the orthonormal basis $(e_r, e_\theta, e_\phi)$ in your code and use it to express the gradient. You have to use the method vector_frame():

BPf = BP.frame()
e = E.vector_frame(('e_r' , 'e_θ', 'e_ϕ'), (BPf[1], 1/r*BPf[2], 1/(r*sin(θ))*BPf[3]),
symbol_dual=('e^r' , 'e^θ', 'e^ϕ'))
e

yields

Vector frame (E^3, (e_r,e_θ,e_ϕ))

You can check that $e$ is an orthornomal basis by asking for the components of the metric tensor $g$ with respecf to $e$:

g[e,:]

yields

[1 0 0]
[0 1 0]
[0 0 1]

Then if suffices to pass $e$ as the argument of display():

yields

grad(f) = d(F)/dr e_r + d(F)/dθ/r e_θ + d(F)/dϕ/(r*sin(θ)) e_ϕ

An alternative is to use the predefined orthonormal frame in EuclideanSpace, which is returned by the method spherical_frame(). For this you have to use the predefined spherical coordinates and not to construct them from scratch as you do. Accordingly, your code would become:

E.<x,y,z> = EuclideanSpace()
CA = E.cartesian_coordinates()
BP.<r, θ, ϕ> = E.spherical_coordinates()
E.set_default_chart(BP)
f = E.scalar_field({BP: function('F')(r, θ, ϕ)}, name='f')
e = E.spherical_frame()

Even shorter, if you don't need the Cartesian coordinate at all, you can introduce the Euclidean space with only spherical coordinates on it:

E.<r, θ, ϕ> = EuclideanSpace(coordinates='spherical')
BP = E.spherical_coordinates()
f = E.scalar_field({BP: function('F')(r, θ, ϕ)}, name='f')
e = E.spherical_frame()