Vector fields in spherical coordinates

Unit vectors and the divergence

Specification using the spherical unit vectors

In spherical coordinates, a vector field F can be specified as

F(ρ,φ,θ)=(Subscript[F, ρ])(ρ,φ,θ)(Subscript[e, ρ])+(Subscript[F, φ])(ρ,φ,θ)(Subscript[e, φ])+(Subscript[F, θ])(ρ,φ,θ)(Subscript[e, θ]),

where (Subscript[e, ρ]), (Subscript[e, φ]), and (Subscript[e, θ]) are the spherical unit vectors.  Each of the unit vectors is a vector of length 1 pointing in the direction induced by a small change in the corresponding variable.  Here’s an illustration:

pt[rho_,phi_,theta_] = rho{Sin[phi]Cos[theta],Sin[phi]Sin[theta],Cos[phi]};
erho[rho_,phi_,theta_] = {Sin[phi]Cos[theta],Sin[phi]Sin[theta],Cos[phi]};
etheta[rho_,phi_,theta_] = {-Sin[theta],Cos[theta],0};
ephi[rho_,phi_,theta_] = Simplify[Cross[etheta[rho,phi,theta],erho[rho,phi,theta]]];
With[{rho=2,phi=Pi/6,theta=23Pi/12},
point=pt[rho,phi,theta];
Show[{
Graphics3D[{{Opacity[0.05],Sphere[{0,0,0},rho]},
{PointSize[Large],Point[point]},
{Thin,Line[{{0,0,0},point}]},
Thick,
Arrow[{point,point+erho[rho,phi,theta]}],
Arrow[{point,point+etheta[rho,phi,theta]}],
Arrow[{point,point+ephi[rho,phi,theta]}]
}],
ParametricPlot3D[pt[rho,phi,thetaLoc],{thetaLoc,0,2Pi},
PlotStyle -> Directive[Opacity[0.2]]],
ParametricPlot3D[pt[rho,phiLoc,theta],{phiLoc,0,2Pi},
PlotStyle -> Directive[Opacity[0.2]]]
}]
]

Note that the sperical unit vectors can be expressed in terms of the Cartesian unit vectors

(Subscript[e, ρ])=sin(φ)cos(θ)i+sin(φ)sin(θ)j+cos(φ)k=((x i + y j + z k)/(Sqrt[(x^2)+(y^2)+(z^2)]))
(Subscript[e, θ])=-sin(θ)i+cos(θ)j = ((-y i+ x j)/(Sqrt[(x^2)+(y^2)]))
(Subscript[e, φ])=(Subscript[e, θ])×(Subscript[e, ρ])=cos(φ)cos(θ)i+cos(φ)sin(θ)j-sin(φ)k=((x z i+y z j-((x^2)+(y^2)+(z^2))k)/((Sqrt[(x^2)+(y^2)])(Sqrt[(x^2)+(y^2)+(z^2)])))

Going the other way, we have

i=sin(φ)cos(θ)(Subscript[e, ρ])+cos(φ)cos(θ)(Subscript[e, φ])-sin(θ)(Subscript[e, θ])
j=sin(φ)sin(θ)(Subscript[e, ρ])+cos(φ)sin(θ)(Subscript[e, φ])+cos(θ)(Subscript[e, θ])
k=cos(φ)(Subscript[e, ρ])-sin(φ)(Subscript[e, θ])

The divergence

The divergence is

div F = (1/(ρ^2))(∂/(∂ρ))((Subscript[F, ρ])(ρ^2))+(1/(ρ sin(φ)))(∂/(∂φ))((Subscript[F, φ])sin(φ))+(1/(ρ sin(φ)))((∂(Subscript[F, θ]))/(∂θ)).

Simple example

Suppose our field is given in Cartesian coordinates as

F(x,y,z)=⟨(x/(Sqrt[(x^2)+(y^2)+(z^2)])),(y/(Sqrt[(x^2)+(y^2)+(z^2)])),(z/(Sqrt[(x^2)+(y^2)+(z^2)]))⟩.

In spherical coordinates, this has the very simple representation F=(Subscript[e, ρ]). Thus, (Subscript[F, ρ])=1 while (Subscript[F, φ])=(Subscript[F, θ])=0.  The divergence computation is simply

div F = (1/(ρ^2))(∂/(∂ρ))((ρ^2))=((2ρ)/(ρ^2))=(2/ρ).

In Cartesian coordinates, the computation is

Simplify[D[(x/(Sqrt[(x^2)+(y^2)+(z^2)])),x]+D[(y/(Sqrt[(x^2)+(y^2)+(z^2)])),y]+D[(z/(Sqrt[(x^2)+(y^2)+(z^2)])),z]]

Derivation

To compute the divergence at a point given in spherical coordinates, we encase the point in a closed surface with six sides: two sides of constant ρ, two of constant φ, and two of constant θ.  We then integrate

(∫_S)F·dn

and allow the surrounding volume, which is (Subscript[ρ, 0]^2)sin((Subscript[φ, 0]))Δρ Δφ Δθ, to shrink to zero.  We should get three terms in the forumula - one for each variable.

pt[rho_,phi_,theta_] := rho{Sin[phi]Cos[theta],Sin[phi]Sin[theta],Cos[phi]};
rhoPatch[rho_,phi0_,theta0_,dphi_,dtheta_] := ParametricPlot3D[
pt[rho,phi,theta],
{phi,phi0-dphi/2,phi0+dphi/2},
{theta,theta0-dtheta/2,theta0+dtheta/2},
Mesh -> None, PlotStyle -> Directive[Opacity[0.4],ColorData[1,1]]];
phiPatch[phi_,rho0_,theta0_,drho_,dtheta_] := ParametricPlot3D[
pt[rho,phi,theta],
{rho,rho0-drho/2,rho0+drho/2},
{theta,theta0-dtheta/2,theta0+dtheta/2},
Mesh -> None, PlotStyle -> Directive[Opacity[0.4],ColorData[1,1]]];
thetaPatch[theta_,rho0_,phi0_,drho_,dphi_] := ParametricPlot3D[
pt[rho,phi,theta],
{rho,rho0-drho/2,rho0+drho/2},
{phi,phi0-dphi/2,phi0+dphi/2},
Mesh -> None, PlotStyle -> Directive[Opacity[0.4],ColorData[1,1]]];
{rho0,phi0,theta0} = {0.9,Pi/3,23Pi/12};
{drho,dphi,dtheta} = {0.2,0.2,0.3};
Show[{
rhoPatch[rho0+drho/2,phi0,theta0,dphi,dtheta],
rhoPatch[rho0-drho/2,phi0,theta0,dphi,dtheta],
phiPatch[phi0+dphi/2,rho0,theta0,drho,dtheta],
phiPatch[phi0-dphi/2,rho0,theta0,drho,dtheta],
thetaPatch[theta0+dtheta/2,rho0,phi0,drho,dphi],
thetaPatch[theta0-dtheta/2,rho0,phi0,drho,dphi],
Graphics3D[{
{Opacity[0.05],Sphere[],Thick,
Line[{{0,0,0},pt[rho0+drho/2,phi0+dphi/2,theta0+dtheta/2]}],
Line[{{0,0,0},pt[rho0+drho/2,phi0-dphi/2,theta0+dtheta/2]}],
Line[{{0,0,0},pt[rho0+drho/2,phi0+dphi/2,theta0-dtheta/2]}],
Line[{{0,0,0},pt[rho0+drho/2,phi0-dphi/2,theta0-dtheta/2]}],
PointSize[Large],Point[pt[rho0,phi0,theta0]]
}}]
}, PlotRange -> All]

The ρ term

Outward flux (away from the origin) through side with constant ρ=(Subscript[ρ, 0])+Δρ/2 is

(Subscript[F, ρ])((Subscript[ρ, 0])+((Δρ)/2),(Subscript[φ, 0]),(Subscript[θ, 0]))((((Subscript[ρ, 0])+((Δρ)/2)))^2)sin((Subscript[φ, 0]))Δφ Δθ

Inward flux (towards the origin) through side with constant ρ=(Subscript[ρ, 0])-Δρ/2 is

-(Subscript[F, ρ])((Subscript[ρ, 0])-((Δρ)/2),(Subscript[φ, 0]),(Subscript[θ, 0]))((((Subscript[ρ, 0])-((Δρ)/2)))^2)sin((Subscript[φ, 0]))Δφ Δθ

Adding these together, we get

(((Subscript[F, ρ])((Subscript[ρ, 0])+((Δρ)/2),(Subscript[φ, 0]),(Subscript[θ, 0]))((((Subscript[ρ, 0])+((Δρ)/2)))^2)-(Subscript[F, ρ])((Subscript[ρ, 0])-((Δρ)/2),(Subscript[φ, 0]),(Subscript[θ, 0]))((((Subscript[ρ, 0])-((Δρ)/2)))^2))/((Subscript[ρ, 0]^2)Δρ))(Subscript[ρ, 0]^2)sin((Subscript[φ, 0]))Δρ Δφ Δθ ∼ (1/(Subscript[ρ, 0]^2))(∂/(∂ρ))((Subscript[F, ρ])(ρ^2))Δ V

That leads to a term of

(1/(ρ^2))(∂/(∂ρ))((Subscript[F, ρ])(ρ^2))

in the spherical formula for divergence.

The φ term

The top is wedge with area Δ A=(Subscript[ρ, 0])sin((Subscript[φ, 0])+Δφ/2)Δρ Δθ and corresponds to constant φ=(Subscript[φ, 0])+Δφ/2.  The outward flux through this side is

(Subscript[F, φ])((Subscript[ρ, 0]),(Subscript[φ, 0])+Δφ/2,(Subscript[θ, 0]))(Subscript[ρ, 0])sin((Subscript[φ, 0])+Δφ/2)Δρ Δθ

Outward flux through the bottom, corresponding to constant φ=(Subscript[φ, 0])-Δφ/2 is

-(Subscript[F, φ])((Subscript[ρ, 0]),(Subscript[φ, 0])-Δφ/2,(Subscript[θ, 0]))(Subscript[ρ, 0])sin((Subscript[φ, 0])-Δφ/2)Δρ Δθ

Adding these together, we get

(((Subscript[F, φ])((Subscript[ρ, 0]),(Subscript[φ, 0])+Δφ/2,(Subscript[θ, 0]))(Subscript[ρ, 0])sin((Subscript[φ, 0])+Δφ/2)-(Subscript[F, φ])((Subscript[ρ, 0]),(Subscript[φ, 0])-Δφ/2,(Subscript[θ, 0]))(Subscript[ρ, 0])sin((Subscript[φ, 0])-Δφ/2))/((Subscript[ρ, 0])sin((Subscript[φ, 0]))Δφ))(Subscript[ρ, 0]^2)sin((Subscript[φ, 0]))Δρ Δφ Δθ ∼ (1/((Subscript[ρ, 0])sin((Subscript[φ, 0]))))(∂/(∂φ))((Subscript[F, φ])sin(φ))Δ V

That leads to a term of

(1/(ρ sin(φ)))(∂/(∂φ))((Subscript[F, φ])sin(φ))

in the spherical formula for divergence.

Activities

Activity 1

Consider the vector field F=(Subscript[e, φ])+(Subscript[e, θ]).

What does the field look like?

Geometrically, why should the flux of the field across any sphere centered at the origin be zero?

Verify the last part via the divergence.

Find the flux across the region inside the sphere of radius 2, outside the sphere of radius 1, and above the cone φ=π/4.

top = ParametricPlot3D[2{Cos[θ]Sin[φ],Sin[θ]Sin[φ],Cos[φ]},{φ,0,Pi/4},{θ,0,2Pi}];
bottom  = ParametricPlot3D[{Cos[θ]Sin[φ],Sin[θ]Sin[φ],Cos[φ]},{φ,0,Pi/4},{θ,0,2Pi}];
side = ParametricPlot3D[r{Cos[θ]Sin[Pi/4],Sin[θ]Sin[Pi/4],Cos[Pi/4]},{r,1,2},{θ,0,2Pi}];
Show[{top,bottom,side}, PlotRange->All]

Recall that, in spherical coordinates, the divergence takes the form

div F = (1/(ρ^2))(∂/(∂ρ))((Subscript[F, ρ])(ρ^2))+(1/(ρ sin(φ)))(∂/(∂φ))((Subscript[F, φ])sin(φ))+(1/(ρ sin(φ)))((∂(Subscript[F, θ]))/(∂θ)).

Activity 2

Consider the vector field F=sin(θ)(Subscript[e, θ]).

What does the field look like?

Geometrically, why should the flux of the field across any cylinder centered on the z axis be zero?

Verify the last part via the divergence.

Find the flux across the “cylindrical cube” with opposing (r,θ,z) vertices (1,0,0) and (2,π/2,1) as shown below.

bottom = ParametricPlot3D[{r*Cos[θ],r*Sin[θ],0},{r,1,2},{θ,0,Pi/2},
PlotPoints -> 40, BoundaryStyle->Thick];
top = ParametricPlot3D[{r*Cos[θ],r*Sin[θ],1},{r,1,2},{θ,0,Pi/2},
PlotPoints -> 40, BoundaryStyle->Thick];
inner = ParametricPlot3D[{Cos[θ],Sin[θ],z},{z,0,1},{θ,0,Pi/2},
PlotPoints -> 40, BoundaryStyle->Thick];
outer = ParametricPlot3D[{2Cos[θ],2Sin[θ],z},{z,0,1},{θ,0,Pi/2},
PlotPoints -> 40, BoundaryStyle->Thick];
xz = ParametricPlot3D[{r,0,z},{r,1,2},{z,0,1},
PlotPoints -> 40, BoundaryStyle->Thick];
yz = ParametricPlot3D[{0,r,z},{r,1,2},{z,0,1},
PlotPoints -> 40, BoundaryStyle->Thick];
Show[{top,bottom,inner,outer,xz,yz}, PlotRange->All]

Recall that, in cylindrical coordinates, the divergence takes the form

div F = (1/r)(∂/(∂r))(r (Subscript[F, r]))+(1/r)((∂(Subscript[F, θ]))/(∂θ))+((∂(Subscript[F, z]))/(∂z)).