202010_surfaces.sagews

%md 
# Surfaces. First and second fundamental form, curvatures.

Surfaces. First and second fundamental form, curvatures.

%md
## Paraboloid

Paraboloid

u1,u2 = var('u1, u2')
paraboloid_eq(u1, u2) = (u1*cos(u2), u1*sin(u2), u1^2)
paraboloid = ParametrizedSurface3D(paraboloid_eq(u1,u2), [u1,u2], 'paraboloid')
var('v1, v2, w1, w2')
v = vector([v1, v2]) 
w = vector([w1, w2]) 
print paraboloid.first_fundamental_form(v, w)
print paraboloid.first_fundamental_form(v, w).substitute(u1=sqrt(2), u2=pi/4)

(v1, v2, w1, w2) u1^2*v2*w2 + (4*u1^2 + 1)*v1*w1 9*v1*w1 + 2*v2*w2

print 'K= ', paraboloid.gauss_curvature()
print 'K(sqrt(2) = ', paraboloid.gauss_curvature().substitute(u1=sqrt(2))

K= 4/(16*u1^4 + 8*u1^2 + 1) K(sqrt(2) = 4/81

%md 
## Cone

Cone

u1,u2 = var('u1, u2')
cone_eq(u1, u2) = (u1*cos(u2), u1*sin(u2), u1)
cone = ParametrizedSurface3D(cone_eq(u1,u2), [u1,u2], 'cone')
print cone.first_fundamental_form_coefficients()
print cone.second_fundamental_form_coefficients()
print cone.mean_curvature()
print cone.principal_directions()

{(1, 2): 0, (1, 1): 2, (2, 1): 0, (2, 2): u1^2} {(1, 2): 0, (1, 1): 0, (2, 1): 0, (2, 2): 1/2*sqrt(2)*u1} 1/4*sqrt(2)/u1 [(1/2*sqrt(2)/u1, [(0, 1)], 1), (0, [(1, 0)], 1)]

%md
## Helicoid

Helicoid

u1,u2 = var('u1, u2')
helicoid_eq(u1, u2) = (u1*cos(u2), u1*sin(u2), u2)
helicoid = ParametrizedSurface3D(helicoid_eq(u1,u2), [u1,u2], 'helicoide')
print helicoid.first_fundamental_form_coefficients()
print helicoid.second_fundamental_form_coefficients()
print helicoid.mean_curvature()
print helicoid.principal_directions()
print helicoid.shape_operator()
print helicoid.gauss_curvature()
print helicoid.connection_coefficients()

{(1, 2, 1): 0, (2, 2, 2): 0, (1, 2, 2): 1/2/u1, (2, 1, 1): 0, (1, 1, 2): 0, (2, 2, 1): -2*u1/(4*u1 + 1), (2, 1, 2): 1/2/u1, (1, 1, 1): -1/2/(4*u1^2 + u1)}

%md 
## Paraboloid
 

Paraboloid

u1,u2 = var('u1, u2')
paraboloid_eq(u1, u2) = (sqrt(u1)*cos(u2), sqrt(u1)*sin(u2), u1)
paraboloid = ParametrizedSurface3D(paraboloid_eq(u1,u2), [u1,u2], 'paraboloid')
paraboloid.connection_coefficients()

{(1, 2, 1): 0, (2, 2, 2): 0, (1, 2, 2): 1/2/u1, (2, 1, 1): 0, (1, 1, 2): 0, (2, 2, 1): -2*u1/(4*u1 + 1), (2, 1, 2): 1/2/u1, (1, 1, 1): -1/2/(4*u1^2 + u1)}

%md
# Surfaces of revolution

Surfaces of revolution

u1, u2 = var('u1, u2')
f = function('f')(u1)
surf_rev_eq(u1, u2) = (f(u1)*cos(u2), f(u1)*sin(u2), u1)
surf_rev = ParametrizedSurface3D(surf_rev_eq(u1, u2), [u1, u2], 'surface of revolution')
print surf_rev.first_fundamental_form_coefficients()
print surf_rev.second_fundamental_form_coefficients()

{(1, 2): 0, (1, 1): diff(f(u1), u1)^2 + 1, (2, 1): 0, (2, 2): f(u1)^2} {(1, 2): 0, (1, 1): -diff(f(u1), u1, u1)/sqrt(diff(f(u1), u1)^2 + 1), (2, 1): 0, (2, 2): f(u1)/sqrt(diff(f(u1), u1)^2 + 1)}

{(1, 2): 0, (1, 1): diff(f(u1), u1)^2 + 1, (2, 1): 0, (2, 2): f(u1)^2}

%md 
# Geodesic curves

Geodesic curves

%md
## Helicoid geodesics

Helicoid geodesics

u1,u2 = var('u1, u2')
helicoid_eq(u1, u2) = (u1*cos(u2), u1*sin(u2), u2)
helicoid = ParametrizedSurface3D(helicoid_eq(u1,u2), [u1,u2], 'helicoide')
psi = pi/4 #starting angle
initial_point = helicoid_eq(2, 0)
ort_frame = [helicoid.orthonormal_frame()[indx].substitute(u1 = 2, u2 = 0) for indx in [1, 2]]
v_int = (sin(psi), cos(psi)/3)
initial_data_graph = point(initial_point,size=20) + line([initial_point, initial_point + ort_frame[0]],color='black') + line([initial_point, initial_point + ort_frame[1]],color='black')
pg_list = helicoid.geodesics_numerical([2, 0],v_int,[0, 10, 1000])
graph = points([ xxx[3] for xxx in pg_list ], color = 'red') + parametric_plot3d(helicoid_eq(u1,u2, 1, 2), (u1, 1, 20), (u2, -pi, pi), opacity=0.3)
show(graph)

%md 
## Torus geodesics

Torus geodesics

var('u1, u2')
var('a, b')
torus_eq(u1, u2, a, b) = [ (a*cos(u1) + b)*cos(u2), (a*cos(u1) + b)*sin(u2), a*sin(u1) ]
torus = ParametrizedSurface3D(torus_eq(u1,u2,a,b) , [u1, u2])
torus_1 = ParametrizedSurface3D(torus_eq(u1,u2,1,2) , [u1, u2])
psi = pi/6 #starting angle
initial_point = torus_eq(0, 0, 1, 2)
ort_frame = [torus_1.orthonormal_frame()[indx].substitute(u1 = 0, u2 = 0) for indx in [1, 2]]
v_int = (sin(psi), cos(psi)/3)
initial_data_graph = point(initial_point,size=20) + line([initial_point, initial_point + ort_frame[0]],color='black') + line([initial_point, initial_point + ort_frame[1]],color='black')
pg_list = torus_1.geodesics_numerical([0,0],v_int,[0, 30, 1000])
graph = points([ xxx[3] for xxx in pg_list ], color = 'red') + parametric_plot3d(torus_eq(u1,u2, 1, 2), (u1, 0, 2*pi), (u2, 0, 2*pi), opacity=0.3)
u1_lim = arccos(3*cos(psi) - 2).n()
print u1_lim
circle_limit = parametric_plot3d(torus_eq(u1_lim, u2, 1, 2), (u2, 0, 2*pi), color='green', thickness=2)
show(graph + initial_data_graph + circle_limit, aspect_ratio=(1, 1, 1))

(u1, u2) (a, b) 0.929697791487669

3D rendering not yet implemented

%md 
## Catenoid geodesics

Catenoid geodesics

u1,u2 = var('u1, u2')
catenoid_eq(u1, u2) = (cosh(u1)*cos(u2), cosh(u1)*sin(u2), u1)
catenoid = ParametrizedSurface3D(catenoid_eq(u1,u2), [u1,u2], 'catenoid')
psi = 0 #starting angle
initial_point = catenoid_eq(0.2, 0)
ort_frame = [catenoid.orthonormal_frame()[indx].substitute(u1 = 2, u2 = 0) for indx in [1, 2]]
v_int = (sin(psi), cos(psi))
initial_data_graph = point(initial_point,size=20) + line([initial_point, initial_point + ort_frame[0]],color='black') + line([initial_point, initial_point + ort_frame[1]],color='black')
pg_list = catenoid.geodesics_numerical([0.2, 0],v_int,[0, 4, 1000])
graph = points([ xxx[3] for xxx in pg_list ], color = 'red') + parametric_plot3d(catenoid_eq(u1,u2, 1, 2), (u1, -2, 2), (u2, 0, 2*pi), opacity=0.3)
show(graph)

3D rendering not yet implemented

%md 
# Minimal surfaces

Minimal surfaces

%md 
## Catenoid

u1,u2 = var('u1, u2')
catenoid_eq(u1, u2) = (cosh(u1)*cos(u2), cosh(u1)*sin(u2), u1)
catenoid = ParametrizedSurface3D(catenoid_eq(u1,u2), [u1,u2], 'catenoid')
print "mean curvature = ", catenoid.mean_curvature()
print catenoid.area_form()
print integrate(catenoid.area_form(), (u1, -1, 1)).expand()

mean curvature = 0 cosh(u1)^2 1/4*e^2 - 1/4*e^(-2) + 1

%md 
## Enneper surface

Enneper surface

var('u, v')
enneper_eq = [u - u^3/3 + u*v^2, v - v^3/3 + u^2*v, u^2 - v^2]
enneper = ParametrizedSurface3D(enneper_eq, [u, v], "enneper surface")
print  "mean curvature =", enneper.mean_curvature()
parametric_plot3d(enneper_eq, (u, -2, 2), (v, -2, 2))

(u, v) mean curvature = 0

3D rendering not yet implemented

%md 
## Scherk surface

Scherk surface

var('u, v')
scherk_eq = [u, v, log(cos(u)/cos(v))]
scherk = ParametrizedSurface3D(scherk_eq, [u, v], "scherk surface")
print "mean curvature=", scherk.mean_curvature()
parametric_plot3d(scherk_eq, (u, -5, 5), (v, -5, 5))

(u, v) mean curvature= 0

3D rendering not yet implemented

%md
# Sufaces of constant positive curvature

Sufaces of constant positive curvature

%md
## Sievert surface

Sievert surface

u, v = var('u, v')
C = var('C')
phi(u, C) = -u/sqrt(C + 1) + arctan(sqrt(C + 1)*tan(u))
a(u, v, C) = 2/(C + 1 - C*sin(v)^2*cos(u)^2)
r(u, v, C) = a(u, v, C)*sqrt( (C + 1)*(1 + C*sin(u)^2) )*sin(v)/sqrt(C)
x(u, v, C) = r(u, v, C)*cos(phi(u, C))
y(u, v, C) = r(u, v, C)*sin(phi(u, C))
z(u, v, C) = (log(tan(v/2)) + a(u, v)*(C + 1)*cos(v))/sqrt(C)
sievert_eq(u, v, C) = [x(u, v, C), y(u, v, C), z(u, v, C)] 
sievert = ParametrizedSurface3D(sievert_eq(u, v, 1), [u, v], "Sievert surface")

parametric_plot3d(sievert_eq(u, v, 1), (u, -pi/2, pi/2), (v, 0, pi/2))

# print sievert.gauss_curvature()

3D rendering not yet implemented

%md
## Remb surface

Remb surface

var('u, v')
var('c')
fU(u,c) = cosh(u*sqrt(c))/sqrt(c)
fV(v,c) = cos(v*sqrt(c+1))/sqrt(c+1)
gU(u,c) = diff(fU(u,c), u)
gV(v,c) = diff(fV(v,c), v)
a(u, v, c) = 2*fV(v)/( (c+1)*(fU(u)^2 - fV(v)^2) )
x(u, v, c) = a(u, v, c)*(fU(u, c)*cos(u) - gU(u, c)*sin(u))
y(u, v, c) = -a(u, v, c)*(fU(u, c)*sin(u) + gU(u, c)*cos(u))
z(u, v, c) = v - a(u, v, c)*gV(v, c)
remb_eq(u, v, c) = [x(u, v, c), y(u, v, c), z(u, v, c)] 
remb = ParametrizedSurface3D(remb_eq(u, v, 1), [u, v], "Rhemb's surface")

parametric_plot3d(remb_eq(u, v, 1), (u, -1, 1), (v, -1, 1))
print remb.gauss_curvature()

(u, v) c

3D rendering not yet implemented

1

print remb.mean_curvature()
plot3d(remb.mean_curvature(), (u, -1, 1), (v, -1, 1)) 

1/8*(sqrt(2)*cos(sqrt(2)*v)^4 + 12*sqrt(2)*cos(sqrt(2)*v)^2*cosh(u)^2 + 4*sqrt(2)*cosh(u)^4)/(cos(sqrt(2)*v)^3*cosh(u) + 2*cos(sqrt(2)*v)*cosh(u)^3)

3D rendering not yet implemented

1.76207903005194

%md
# Geodesics of hyperbolic plane

Geodesics of hyperbolic plane

s = var('s')
a, R = var('a, R')
x(s, a, R) = a + R*(1 - exp(2*s))/(1 + exp(2*s))
y(s, a, R) = 2*R*exp(s)/(1 + exp(2*s))
parametric_plot([x(s, 1, 1), y(s, 1, 1)], (s, -10, 10))

xx = a + R*(1 - exp(2*s))/(1 + exp(2*s))
yy = 2*R*exp(s)/(1 + exp(2*s))
show((diff(xx, s, 2) - 2*diff(xx, s)*diff(yy, s)/yy).simplify_full())
show((diff(yy, s, 2) + (diff(xx, s)^2 - diff(yy, s)^2)/yy).simplify_full())

$\displaystyle 0$

$\displaystyle 0$

︠52ef27ac-cfdc-4a14-acdd-55ee2d607b20︠