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)
