Geodesics curve

https://faculty.etsu.edu/knisleyj/multicalc/Chap3/Chap3-8/part1.htm

Geodesic curve is a surface curve that coincides surface normal at each point. This means geodesic curve normal vector must be perpendicular to surface tangent plane. Any vector in surface tangent plane is perpendicular to geodesic curve normal vector. Therefore, \(\frac{dr}{du}\) and \(\frac{dr}{dv}\) are vectors in surface tangent plane, \(\rho''\) is surface curve normal vector, \(\frac{dr}{du} \cdot \rho''=0\) and \(\frac{dr}{dv}\cdot \rho''=0\).

To replicate this cylinder example, We have a cylinder \(\bar r(u,v)=\{\cos(u),\sin(u),v\}\). A surface curve \(\rho(t)=\bar r(u=t,v=2t+1)=\{\cos(t),\sin(t),2t+1\}\). Show \(\rho\) is a geodesic surface curve.

library(Deriv)
## Warning: package 'Deriv' was built under R version 4.0.5
#surface 
rbar=c("cos(u)","sin(u)","v")

#curve on surface 
u="t"
v="2*t+1"
rho=gsub("u", u, rbar,fixed=TRUE)
rho=gsub("v", v, rho,fixed=TRUE)

#degree 1 and 2 partial on rbar
drdu=sapply(rbar,function(m) Deriv(m,"u"))
drdv=sapply(rbar,function(m) Deriv(m,"v"))

#degree 1 and 2 partial on rho
rho1=sapply(rho, function(m) Deriv(m,"t"))
rho2=sapply(rho1,function(m) Deriv(m,"t"))
dot_prod<-function(a,b){
#paste0("(",c("1","2","3"),")")
 a=paste0("(",a,")")
 b=paste0("(",b,")")
 res=Simplify(paste(paste(a,b,sep="*"),collapse="+") )
 return(paste0("(",res,")"))
}

dot_prod(rho2,drdu) # u == t, => 0
## [1] "(cos(t) * sin(u) - cos(u) * sin(t))"
dot_prod(rho2,drdv) # => 0
## [1] "(0)"
#rho( t) = r( t, 2t+1) is a geodesic on the cylinder

Geodesic curvature

https://math.stackexchange.com/questions/1590579/geodesic-curvature-of-sphere-parallels?rq=1

The normal vector of a sphere is just the sphere surface point itself, and divided by the radius we can get unit normal vector. So we have the sphere with radius \(a\) and using unit-speed parameterization, the surface is \[\bar r = \{r \cos(\theta/r),r\sin(\theta/r), \sqrt{a^2 - r^2}\}\] And the surface normal is \[N = \{\frac{r}{a} \cos(\frac{\theta}{r}),\frac{r}{a}\sin(\frac{\theta}{r}), \frac{\sqrt{a^2 - r^2}}{a}\}\] But I went through the cross product just to find matching, it requires manual simplification. \[N=\frac{\frac{\partial \bar r}{\partial \theta}\times \frac{\partial\bar r}{\partial r}}{||\frac{\partial \bar r}{\partial \theta}\times \frac{\partial\bar r}{\partial r}||}\]

library(Deriv)

#rbar
x="r*cos(theta/r)"
y="r*sin(theta/r)"
z="sqrt(a^2-r^2)"
rbar=c(x,y,z)

#dr/du , dr/dv
drdu=sapply(rbar,function(m) Deriv(m,"theta"))
drdv=sapply(rbar,function(m) Deriv(m,"r"))

drdv=c("cos(theta/r) + theta * sin(theta/r)/r","sin(theta/r) - theta * cos(theta/r)/r","-(r/sqrt(a^2 - r^2))")

cross_prod<-function(x,y){
a1=paste0("(",x[2],")*(",y[3],")-(",y[2],")*(",x[3],")")
a2=paste0("(",x[3],")*(",y[1],")-(",y[3],")*(",x[1],")")
a3=paste0("(",x[1],")*(",y[2],")-(",y[1],")*(",x[2],")")

return(c(Simplify(a1),Simplify(a2),Simplify(a3)))
}

cross_prod(drdu,drdv)
## [1] "-(r * cos(theta/r)/sqrt(a^2 - r^2))"                                                                                    
## [2] "-(r * sin(theta/r)/sqrt(a^2 - r^2))"                                                                                    
## [3] "-(cos(theta/r) * (cos(theta/r) + theta * sin(theta/r)/r) + sin(theta/r) *      (sin(theta/r) - theta * cos(theta/r)/r))"
#it requires manual simplification.
j1="r * cos(theta/r)/sqrt(a^2 - r^2)"
j2="r * sin(theta/r)/sqrt(a^2 - r^2)"
j3="1"

#it requires manual simplification.
normal_length= paste0("(",j1,")^2+","(",j2,")^2+","(",j3,")^2")
normal_length="a/sqrt(a^2 - r^2)"


N1=paste0("(",j1,")/(",normal_length,")")
N2=paste0("(",j2,")/(",normal_length,")")
N3=paste0("(",j3,")/(",normal_length,")")

N=c(Simplify(N1),Simplify(N2),Simplify(N3))
N
## [1] "r * cos(theta/r)/a" "r * sin(theta/r)/a" "sqrt(a^2 - r^2)/a"

When we consider the parameter \(r\) to be fixed, we can plot a small circle. So the parameter \(r\) is the radius of the small circle \(\gamma(\theta)\). \[\gamma(\theta) = \left\{r \cos\frac{\theta}{r}, \;r \sin\frac{\theta}{r}, \; \sqrt{a^2 - r^2}\right\}\]

First order derivative is an unit tangent \[\gamma'(\theta) = \left\{- \sin\frac{\theta}{r}, \; \cos\frac{\theta}{r}, \; 0\right\}\]

gamma1=drdu

Second order derivative is curvature vector \(\bar k\) \[\gamma''(\theta) = \left\{-\frac{1}{r} \cos\frac{\theta}{r}, \;-\frac{1}{r} \sin\frac{\theta}{r}, \; 0\right\}\]

gamma2=sapply(drdu,function(m) Deriv(m,"theta"))

Now we compute that, use \(\cos^2(\theta/r) + \sin^2(\theta/r)=1\) \[k_g = \gamma'' \cdot (\mathbf{N}\times \gamma ') =\frac{h}{ra} = \frac{\sqrt{a^2 - r^2}}{ra}\]

dot_prod<-function(a,b){
#paste0("(",c("1","2","3"),")")
 a=paste0("(",a,")")
 b=paste0("(",b,")")
 res=Simplify(paste(paste(a,b,sep="*"),collapse="+") )
 return(paste0("(",res,")"))
}
dot_prod(gamma2,cross_prod(N,gamma1))
## [1] "((cos(theta/r)^2 + sin(theta/r)^2) * sqrt(a^2 - r^2)/(a * r))"
#((cos(theta/r)^2 + sin(theta/r)^2) == 1
kg="sqrt(a^2 - r^2)/(a * r))"

In the original stackexchange example we see an angle \(\psi\) between \(N\) and \(n\). We know this is called surface curvature decomposition. So we have \[k_g=k\cos(\psi)\] Please note the \(k_g\) and \(k\) also have vector forms \(\bar k_g\) and \(\bar k\). We had \(\bar k=\gamma''\).

Geodesic parallels

https://www.kristakingmath.com/blog/how-to-find-orthogonal-trajectories

We had this picture from our plane geometry studies. Blue lines are called family of lines. Black circles are called orthogonal trajectories, because the black circles cut the blue lines orthogonally. If we wrap this picture around a sphere, we can see that the blue lines become great circles (geodesic curves), and the black circles are called the geodesic parallels. orthogonal trajectory

So that on a curved surface, orthogonal trajectories are geodesic parallels. We use geodesic parallels to measure distance along geodesics.

Geodesic coordinates

https://math.stackexchange.com/questions/971482/about-geodesic-polar-coordinate https://tutorial.math.lamar.edu/Classes/CalcII/PolarCoordinates.aspx https://byjus.com/maths/polar-coordinates/ https://math.stackexchange.com/questions/2334570/how-can-you-calculate-distance-in-plane-polar-coordinates-using-the-metric-tenso

Let \(u\) and \(v\) be parameters of a curved surface, if \(u\) and \(v\) vary perpendicularly along the geodesic and geodesic parallels, we say \(u\) and \(v\) are geodesic parameters of of the surface. For example longitude and latitude are the geodesic parameters of the earth surface(polar coordinate system). And the earth radius \(R\) has a direction that is perpendicular to both longitude and latitude, so the earth radius \(R\) is another geodesic parameter. But we only consider 2 parameters surface in our discussion.

Geodesic torsion

https://encyclopediaofmath.org/wiki/Geodesic_torsion

Line of curvature is a curve on a surface whose tangents are always in the direction of principal curvature. https://mathworld.wolfram.com/LineofCurvature.html

https://math.stackexchange.com/questions/1076396/curvature-and-torsion-of-a-helix

Geodesic torsion means torsion of a geodesic. For a non-geodesic curve on surface we need a tangent direction \(T\) to choose a geodesic, to compute geodesic torsion, we need to use the geodesic normal vector which coincides to the surface normal vector. \[\tau_g=\frac{dN}{ds} \cdot(T \times N)=\frac{dN/dt}{ds/dt} \cdot(T \times N)\]

The \(\bar r(u,v)=\{x=r\cos(u),y=r\sin(u),z=v\}\) is a cylinder, and \(\gamma(t)=\bar r(at+b,ct+d)\) is a geodesic on cylinder surface, \(a,b,c,d\) are constants. It means let \[u=at+b\\ v=ct+d\] plug into \(\bar r\) we get a geodesic on cylinder surface.

So we let \(\gamma=\{r\cos(t),r\sin(t),pt\}\), we see \[u=t\\ v=pt\] On cylinder surface this is a circular helix of radius \(r\) and pitch \(2\pi p\). It has curvature \(\frac{|r|}{r^2+p^2}\) and torsion \(\frac{p}{r^2+p^2}.\)

Let’s compute torsion. Geodesic must coincide normal vector to the surface.

require(Deriv)

rbar=c("r*cos(u)","r*sin(u)","v")
drdu=sapply(rbar,function(m) Deriv(m,"u"))
drdv=sapply(rbar,function(m) Deriv(m,"v"))

cross_prod<-function(x,y){
a1=paste0("(",x[2],")*(",y[3],")-(",y[2],")*(",x[3],")")
a2=paste0("(",x[3],")*(",y[1],")-(",y[3],")*(",x[1],")")
a3=paste0("(",x[1],")*(",y[2],")-(",y[1],")*(",x[2],")")

return(c(Simplify(a1),Simplify(a2),Simplify(a3)))
}

cross_prod(drdu,drdv)
## [1] "r * cos(u)"     "--(r * sin(u))" "0"

Cross product we can get the surface unit normal. Parameter \(u\) is changed to \(t\) we can have the helix unit normal, because the normal vectors coincide.

N=c("cos(u)","sin(u)", "0")#cylinder unit normal 

N=c("cos(t)","sin(t)", "0")#plug in u=t, v=pt, helix unit normal

Compute \(\frac{dN}{dt}\)

dNdt=sapply(N,function(m) Deriv(m,"t"))

Compute unit tangent \(T\)

gamma=c("r*cos(t)","r*sin(t)","p*t")#gamma is the helix 

dgammadt=sapply(gamma,function(m) Deriv(m,"t"))

T=paste0(dgammadt,"/sqrt(r^2+p^2)")
dsdt="sqrt(r^2+p^2)"
dNds=paste0(dNdt,"/",dsdt)

The \(ds/dt\) we can get from the differential form of arc length. \[ds/dt=\sqrt{(dx/dt)^2+(dy/dt)^2+(dz/dt)^2}=\sqrt{r^2\sin^2(t)+r^2\cos^2(t)+p^2}=\sqrt{r^2+p^2}\]

Use formula \(\tau_g=\frac{dN/dt}{ds/dt} \cdot(T \times N)\). We can get \(\tau_g=p/(p^2 + r^2)\).

dot_prod<-function(a,b){
#paste0("(",c("1","2","3"),")")
 a=paste0("(",a,")")
 b=paste0("(",b,")")
 res=Simplify(paste(paste(a,b,sep="*"),collapse="+") )
 return(paste0("(",res,")"))
}

dot_prod(cross_prod(N,T),dNds)#cos(t)^2 + sin(t)^2==1
## [1] "(-(p * (cos(t)^2 + sin(t)^2)/(p^2 + r^2)))"
tau="p/(p^2 + r^2)"

Gauss-Bonnet Theorem

https://www.math.ucla.edu/~greene/Gauss-Bonnet.htm

https://mathworld.wolfram.com/Gauss-BonnetFormula.html

https://math.stackexchange.com/questions/4127316/gauss-bonnet-theorem-for-half-a-cone

http://www.rdrop.com/~half/Creations/Puzzles/cone.geodesics/index.html

https://math.stackexchange.com/questions/2656057/how-can-a-cone-have-non-zero-riemann-curvature-yet-can-be-made-out-of-a-piece-of

https://math.stackexchange.com/questions/2458251/find-the-equation-of-the-cone-z-sqrtx2y2-in-spherical-coordinates?msclkid=5ac880e3cf7b11ec8352fcd5d4d1a838

Compute Gauss Bonnet theorem \[\int_Ck_gds+\int \int_A K \, dA=2\pi-\sum\theta_i\] for half a cone \[x^2+y^2=z^2, y\ge0, 0\le z\le1\] We need to determine exterior angles \(\theta_1\), \(\theta_2\), \(\theta_3\). By looking at the picture below, we see \(\theta_1\) and \(\theta_2\) are right angles. And \(\theta_3\) can be anti-clockwise and clockwise, i.e. \(\theta_3=\pi \pm \pi/\sqrt{2}\). Gauss Bonnet

After unfold the cone, we see the Gauss curvature \(K=0\) everywhere. Therefore \(\int \int_A K \, dA=0\).

Next is to determine the curvature \(\int_Ck_gds\) on area boundaries. The boundaries vertex-A and vertex-B are straight lines and geodesics, their curvature is 0. The curve A-B has geodesic curvature \(k_g=1/\sqrt{2}\). The cone has coordinates \[\bar r=\{x=z\cos(\theta),y=z\sin(\theta),z\}\] Here is a image found online.

cone coordinates

N is the unit normal computed as a cross product of degree 1 partial derivatives. T is the tangent of the base circle. K is the derivative of T. Then we can use formula to get \[k_g=k\cdot(N\times T)=\frac{-1}{\sqrt{2}}\] But if we consider the base circle as part of circle with radius \(\sqrt{2}\), then it has curvature \(1/\sqrt{2}\). Anyway here is the full code.

library(Deriv)
## Warning: package 'Deriv' was built under R version 4.0.5
#rbar, cone 
x="z*cos(u)"
y="z*sin(u)"
z="z"

rbar=c(x,y,z)

# 1st order partial
drdu=sapply(rbar,function(m) Deriv(m,"u"))
drdz=sapply(rbar,function(m) Deriv(m,"z"))


cross_prod<-function(x,y){
a1=paste0("(",x[2],")*(",y[3],")-(",y[2],")*(",x[3],")")
a2=paste0("(",x[3],")*(",y[1],")-(",y[3],")*(",x[1],")")
a3=paste0("(",x[1],")*(",y[2],")-(",y[1],")*(",x[2],")")

return(c(Simplify(a1),Simplify(a2),Simplify(a3)))
}

dot_prod<-function(a,b){
#paste0("(",c("1","2","3"),")")
 a=paste0("(",a,")")
 b=paste0("(",b,")")
 res=Simplify(paste(paste(a,b,sep="*"),collapse="+") )
 return(paste0("(",res,")"))
}

cross_prod(drdu,drdz)
## [1] "z * cos(u)"                   "--(z * sin(u))"              
## [3] "-(z * (cos(u)^2 + sin(u)^2))"
#it requires manual simplification.
j1="z*cos(u)"
j2="z*sin(u)" 
j3="-z"

normal_length= Simplify(paste0("(",j1,")^2+","(",j2,")^2+","(",j3,")^2"))
normal_length="sqrt(2)*z"

N1=paste0("(",j1,")/(",normal_length,")")
N2=paste0("(",j2,")/(",normal_length,")")
N3=paste0("(",j3,")/(",normal_length,")")

N=c(Simplify(N1),Simplify(N2),Simplify(N3))
N
## [1] "cos(u)/1.4142135623731" "sin(u)/1.4142135623731" "-(1/1.4142135623731)"
#geodesic curvature on base
x="cos(u)"
y="sin(u)"
z="1"

base=c(x,y,z)

T=sapply(base,function(m) Deriv(m,"u"))
k=sapply(T,function(m) Deriv(m,"u"))

#kg=k.(NxT)
kg=dot_prod(k,cross_prod(N,T))
#==-0.707106781186545(cos(u)^2+sin(u)^2)==-1/sqrt(2)
kg=-1/sqrt(2)

Next, the half base circle has length \(\pi\), \(k_g=\pm \frac{1}{\sqrt{2}}\), so the integral \[\int_Ck_gds=\int^0_{\pi}\pm\frac{1}{\sqrt{2}}ds=\pm \frac{\pi}{\sqrt{2}}\]

Together \[\int_Ck_gds+\int \int_A K \, dA=2\pi-\sum\theta_i\\ \pm \frac{\pi}{\sqrt{2}}+0=2\pi-(2\pi\pm \frac{\pi}{\sqrt{2}})\]