Ask Your Question
0

Wrong result for gradient in spherical coordinates [closed]

asked 2023-11-25 13:44:12 +0100

pfeifhns gravatar image

updated 2023-11-25 14:15:34 +0100

FrédéricC gravatar image

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()

show(f.gradient().display())

But the result is wrong:

grad(f)=∂F∂r∂∂r+∂F∂θr^2∂∂θ+∂F∂ϕr^2sin⁡(θ)^22∂∂ϕ

It must be:

grad(f)=∂F∂r∂∂r+∂F∂θr∂∂θ+∂F∂ϕrsin⁡(θ)∂∂ϕ

What am I doing wrong?

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

edit retag flag offensive 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

2 Answers

Sort by » oldest newest most voted
0

answered 2023-11-25 18:44:02 +0100

pfeifhns gravatar image

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.

edit flag offensive delete link more

Comments

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

eric_g gravatar imageeric_g ( 2023-11-26 15:45:56 +0100 )edit
0

answered 2023-11-25 17:59:34 +0100

eric_g gravatar image

updated 2023-11-26 15:44:38 +0100

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():

f.gradient().display(e)

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()
f.gradient().display(e)

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()
f.gradient().display(e)
edit flag offensive delete link more

Question Tools

Stats

Asked: 2023-11-25 13:44:12 +0100

Seen: 53 times

Last updated: Nov 26 '23