| 85 | sage: u, v = var('u, v') |
| 86 | sage: sphere = ParametrizedSurface3D([cos(u)*cos(v),sin(u)*cos(v),sin(v)],[u,v],'sphere') |
| 87 | sage: print latex(sphere) |
| 88 | \left[\cos\left(u\right) \cos\left(v\right), \sin\left(u\right) \cos\left(v\right), \sin\left(v\right)\right] |
| 89 | sage: print sphere._latex_() |
| 90 | \left[\cos\left(u\right) \cos\left(v\right), \sin\left(u\right) \cos\left(v\right), \sin\left(v\right)\right] |
| 91 | sage: print sphere |
| 92 | Parametrized surface ('sphere') with equation [cos(u)*cos(v), sin(u)*cos(v), sin(v)] |
| 93 | |
| 94 | We can plot the surface using the inner method plot:: |
| 95 | |
| 96 | sage: ellipsoid_parametric_equation_abc = ellipsoid_parametric_equation.substitute(a=2,b=1.5,c=1) |
| 97 | sage: ellipsoid_abc = ParametrizedSurface3D(ellipsoid_parametric_equation_abc,[u1,u2],'ellipsoid_abc') |
| 98 | sage: ellipsoid_abc_plot = ellipsoid_abc.plot((u1,0,2*pi),(u2,-pi/2,pi/2)) |
| 99 | sage: ellipsoid_abc_plot.show(aspect_ratio=(1,1,1)) |
| 100 | |
| 101 | Let us find the natural frame of the parametrization of the ellipsoid. We get the dictionary of vector fields of the natural frame of the ellipsoid:: |
| 102 | |
| 103 | sage: ellipsoid.natural_frame() |
| 104 | {1: (-a*sin(u1)*cos(u2), b*cos(u1)*cos(u2), 0), 2: (-a*sin(u2)*cos(u1), -b*sin(u1)*sin(u2), c*cos(u2))} |
| 105 | |
| 106 | Now find the normal vector field. |
| 107 | The normal vector field which is the vector product of the vectors of the natural frame is:: |
| 108 | |
| 109 | sage: ellipsoid.normal_vector() |
| 110 | (b*c*cos(u1)*cos(u2)^2, a*c*sin(u1)*cos(u2)^2, a*b*sin(u2)*cos(u2)) |
| 111 | |
| 112 | And the unit normal vector field of the elliptic paraboloid is:: |
| 113 | |
| 114 | sage: u, v = var('u,v') |
| 115 | sage: eparaboloid = ParametrizedSurface3D([u,v,u^2+v^2],[u,v],'elliptic paraboloid') |
| 116 | sage: eparaboloid.normal_vector(1) |
| 117 | (-2*u/sqrt(4*u^2 + 4*v^2 + 1), -2*v/sqrt(4*u^2 + 4*v^2 + 1), 1/sqrt(4*u^2 + 4*v^2 + 1)) |
| 118 | |
| 119 | Now let us find the coefficients of the first fundamental form of the torus:: |
| 120 | |
| 121 | u, v = var('u,v') |
| 122 | a, b = var('a,b') |
| 123 | sage: torus = ParametrizedSurface3D(((a + b*cos(u))*cos(v),(a + b*cos(u))*sin(v), b*sin(u)),[u,v],'torus') |
| 124 | sage: torus.first_fundamental_form_coefficients() |
| 125 | {(1, 2): 0, (1, 1): b^2, (2, 1): 0, (2, 2): b^2*cos(u)^2 + 2*a*b*cos(u) + a^2} |
| 126 | |
| 127 | We can find the length of a curve on the surface. |
| 128 | For example, let us find the length of the curve $u^1 = t$, $u^2 =t$, $t \in [0,2\pi]$, |
| 129 | on the ellipsoid with axes $a=1$, $b=1.5$ and $c=1$. So we take the curve:: |
| 130 | |
| 131 | sage: t = var('t') |
| 132 | sage: uu1 = t |
| 133 | sage: uu2 = t |
| 134 | |
| 135 | Then find the tangent vector:: |
| 136 | |
| 137 | sage: duu1 = diff(uu1,t) |
| 138 | sage: duu2 = diff(uu2,t) |
| 139 | sage: duu = vector([duu1,duu2]) |
| 140 | |
| 141 | Then calculate the symbolic expression for the length:: |
| 142 | |
| 143 | sage: L=sqrt(ellipsoid.first_fundamental_form(duu,duu).substitute(u1=uu1,u2=uu2)) |
| 144 | sage: integrate(L,t,0,2*pi) |
| 145 | integrate(sqrt(2*(a^2 - b^2)*sin(t)^2*cos(t)^2 + c^2*cos(t)^2 + (a^2*cos(t)^2 + b^2*sin(t)^2)*sin(t)^2 + (a^2*sin(t)^2 + b^2*cos(t)^2)*cos(t)^2), t, 0, 2*pi) |
| 146 | |
| 147 | And finally find the numerical value:: |
| 148 | |
| 149 | sage: print numerical_integral(L.substitute(a=2,b=1.5,c=1),0,1)[0] |
| 150 | 2.00127905972 |
| 151 | |
| 152 | Find the area of the sphere of radius $R$:: |
| 153 | |
| 154 | sage: R = var('R'); |
| 155 | sage: u, v = var('u,v'); |
| 156 | sage: assume(R>0) |
| 157 | sage: assume(cos(v)>0) |
| 158 | sage: sphere = ParametrizedSurface3D([R*cos(u)*cos(v),R*sin(u)*cos(v),R*sin(v)],[u,v],'sphere') |
| 159 | sage: integral(integral(sphere.area_form(),u,0,2*pi),v,-pi/2,pi/2) |
| 160 | 4*pi*R^2 |
| 161 | |
| 162 | We can find an orthonormal frame field $\{e_1, e_2\}$ of a surface and calculate it structure functions $c^k_{ij}$, |
| 163 | where $[e_i,e_j] = \sum_k c^k_{ij} e_k$. |
| 164 | Let us first find the orthonormal frame field on the elliptic paraboloid:: |
| 165 | |
| 166 | sage: u, v = var('u,v') |
| 167 | sage: eparaboloid = ParametrizedSurface3D([u,v,u^2+v^2],[u,v],'elliptic paraboloid') |
| 168 | sage: eparaboloid.orthonormal_frame() |
| 169 | {1: (1/sqrt(4*u^2 + 1), 0, 2*u/sqrt(4*u^2 + 1)), 2: (-4*u*v/(sqrt(4*u^2 + 1)*sqrt(4*u^2 + 4*v^2 + 1)), sqrt(4*u^2 + 1)/sqrt(4*u^2 + 4*v^2 + 1), 2*v/(sqrt(4*u^2 + 1)*sqrt(4*u^2 + 4*v^2 + 1)))} |
| 170 | |
| 171 | |
| 172 | Now find the inner (local) coordinates of the orthonormal frame field on the elliptic paraboloid:: |
| 173 | |
| 174 | sage: EE = eparaboloid.orthonormal_frame(coordinates='int') |
| 175 | sage: E1 = EE[1]; E2 = EE[2] |
| 176 | sage: CC = eparaboloid.frame_structure_functions(E1,E2) |
| 177 | sage: CC[1,2,1].simplify_full() |
| 178 | 4*sqrt(4*u^2 + 4*v^2 + 1)*v/(sqrt(4*u^2 + 1)*(16*u^4 + 4*(4*u^2 + 1)*v^2 + 8*u^2 + 1)) |
| 179 | |
| 180 | |
| 181 | We find the Gaussian and Mean Curvatures of the sphere:: |
| 182 | |
| 183 | sage: u, v = var('u,v') |
| 184 | sage: sphere = ParametrizedSurface3D((u,v,sqrt(1-u^2-v^2)),[u,v],'sphere') |
| 185 | sage: sphere.gauss_curvature() |
| 186 | 1 |
| 187 | sage: sphere.mean_curvature() |
| 188 | 1 |
| 189 | |
| 190 | Another example is to plot gaussian curvature of ellipsoid by color:: |
| 191 | |
| 192 | sage: u1, u2 = var('u1,u2'); |
| 193 | sage: u = [u1,u2] |
| 194 | sage: ellipsoid_equation(u1,u2) = [2*cos(u1)*cos(u2),1.5*cos(u1)*sin(u2),sin(u1)] |
| 195 | sage: ellipsoid = ParametrizedSurface3D(ellipsoid_equation(u1,u2),u,'ellipsoid') |
| 196 | sage: # set intervals for variables and the number of division points |
| 197 | sage: interval_1 = (-1.5,1.5) |
| 198 | sage: interval_2 = (0,6.28) |
| 199 | sage: plot_points = [10,20] |
| 200 | sage: # make the arguments array |
| 201 | sage: u1_array = [interval_1[0] + cnt*(interval_1[1]-interval_1[0])/plot_points[0] for cnt in range(0,plot_points[0])] |
| 202 | sage: u2_array = [interval_2[0] + cnt*(interval_2[1]-interval_2[0])/plot_points[1] for cnt in range(0,plot_points[1])] |
| 203 | sage: u_array = [ (uu1,uu2) for uu1 in u1_array for uu2 in u2_array] |
| 204 | sage: # Find the gaussian curvature |
| 205 | sage: K(u1,u2) = ellipsoid.gauss_curvature() |
| 206 | sage: # Make array of K values |
| 207 | sage: K_array = [K(uu[0],uu[1]) for uu in u_array] |
| 208 | sage: # Find minimum and max of the gauss curvature |
| 209 | sage: K_max = max(K_array) |
| 210 | sage: K_min = min(K_array) |
| 211 | sage: # Make the array of color coefficients |
| 212 | sage: cc_array = [ (ccc - K_min)/(K_max - K_min) for ccc in K_array ] |
| 213 | sage: points_array = [ellipsoid_equation(u_array[counter][0],u_array[counter][1]) for counter in range(0,len(u_array)) ] |
| 214 | sage: curvature_ellipsoid_plot = sum( point([xx for xx in points_array[counter]],color=hue(cc_array[counter]/2)) for counter in range(0,len(u_array)) ) |
| 215 | sage: curvature_ellipsoid_plot.show(aspect_ratio=1) |
| 216 | |
| 217 | We can find the principal curvatures and principal directions of the elliptic paraboloid:: |
| 218 | |
| 219 | sage: u, v = var('u,v') |
| 220 | sage: eparaboloid = ParametrizedSurface3D([u,v,u^2+v^2],[u,v],'elliptic paraboloid') |
| 221 | sage: pd = eparaboloid.principal_directions() |
| 222 | |
| 223 | We extract the principal curvatures:: |
| 224 | |
| 225 | sage: k1 = pd[0][0].simplify_full() |
| 226 | sage: k1 |
| 227 | 2*sqrt(4*u^2 + 4*v^2 + 1)/(16*u^4 + 8*u^2 + 16*v^4 + 8*(4*u^2 + 1)*v^2 + 1) |
| 228 | sage: k2 = pd[1][0].simplify_full() |
| 229 | sage: k2 |
| 230 | 2/sqrt(4*u^2 + 4*v^2 + 1) |
| 231 | |
| 232 | and let us check it:: |
| 233 | |
| 234 | sage: K = eparaboloid.gauss_curvature().simplify_full() |
| 235 | sage: K |
| 236 | 4/(16*u^4 + 8*u^2 + 16*v^4 + 8*(4*u^2 + 1)*v^2 + 1) |
| 237 | sage: H = eparaboloid.mean_curvature().simplify_full() |
| 238 | sage: H |
| 239 | 2*(2*u^2 + 2*v^2 + 1)/(4*u^2 + 4*v^2 + 1)^(3/2) |
| 240 | sage: (K - k1*k2).simplify_full() |
| 241 | 0 |
| 242 | sage: (2*H - k1 - k2).simplify_full() |
| 243 | 0 |
| 244 | |
| 245 | And we can find the intrinsic (local coordinates) of the principal directions:: |
| 246 | |
| 247 | sage: pd[0][1] |
| 248 | [(1, v/u)] |
| 249 | sage: pd[1][1] |
| 250 | [(1, -u/v)] |
| 251 | |
| 252 | Also with the ParametrizedSurface3D class on can find the coefficients of the second fundamental form, the shape operator, the rotation on the surface at a given angle, the connection coefficients, also one can calculate numerically the geodesics and the parallel translation along a curve (see the documentation below). |
| 253 | """ |
| 254 | |
957 | | @cached_method |
958 | | def principal_curvatures(self): |
959 | | """ |
960 | | Finds the principal curvatures of the surface |
961 | | |
962 | | OUTPUT: |
963 | | The two principal curvatures, given as a dictionary. |
964 | | |
965 | | EXAMPLES:: |
966 | | |
967 | | sage: R = var('R') |
968 | | sage: assume(R>0) |
969 | | sage: u, v = var('u,v') |
970 | | sage: assume(cos(v)>0) |
971 | | sage: sphere = ParametrizedSurface3D([R*cos(u)*cos(v),R*sin(u)*cos(v),R*sin(v)],[u,v],'sphere') |
972 | | sage: sphere.principal_curvatures() |
973 | | {1: -1/R, 2: -1/R} |
974 | | |
975 | | sage: u, v = var('u,v') |
976 | | sage: R, r = var('R, r') |
977 | | sage: assume(R>r,r>0) |
978 | | sage: torus = ParametrizedSurface3D([(R+r*cos(v))*cos(u),(R+r*cos(v))*sin(u),r*sin(v)],[u,v],'torus') |
979 | | sage: torus.principal_curvatures() |
980 | | {1: -cos(v)/(r*cos(v) + R), 2: -1/r} |
981 | | |
982 | | """ |
983 | | |
984 | | from sage.symbolic.assumptions import assume |
985 | | from sage.symbolic.relation import solve |
986 | | from sage.calculus.var import var |
987 | | |
988 | | KK = self.gauss_curvature() |
989 | | HH = self.mean_curvature() |
990 | | |
991 | | # jvkersch: when this assumption is uncommented, Sage raises an error stating that the assumption |
992 | | # is redundant... Can we safely omit this, based on some geometric reasoning? |
993 | | |
994 | | # assume(HH**2-KK>=0) |
995 | | # mikarm: This is a problem, I had a lot of trobles here. Of course, this inequality always hold true. |
996 | | # I insert this assumption because sage sometimes, in simplification, uses the complex numbers, though the roots |
997 | | # are, for sure, real. |
998 | | # This, in turn, causes problems when we substitute coordinates into the expression of principal curvatures. |
999 | | # Unfortunately, I did not manage to tell Sage that they are real (to declare the variables as real). |
1000 | | # Moreover, in a neighborhood of an umbilic point we even cannot "smoothly" order the set of principal curvatures. |
1001 | | # So, in general, at present this method is far from the final form. |
1002 | | |
1003 | | |
1004 | | |
1005 | | x = var('x') |
1006 | | sol = solve(x**2 -2*HH*x + KK == 0,x) |
1007 | | |
1008 | | #k1=var('k1') |
1009 | | #k2=var('k2') |
1010 | | |
1011 | | # jvkersch: when I tried to run the previous version of the code, I ran into the problem that if the equation for the principal curvatures had a double root (as in the case of the sphere example in the worksheet), solve returned only one root. Maybe this is a difference due to having different versions of sage. |
1012 | | |
1013 | | k1 = (x.substitute(sol[0])).simplify_full() |
1014 | | if len(sol) == 1: |
1015 | | k2 = k1 |
1016 | | else: |
1017 | | k2 = (x.substitute(sol[1])).simplify_full() |
1018 | | |
1019 | | return {1:k1,2: k2} |
1020 | | |
1103 | | """ |
1104 | | Finds the principal curvatures and principal directions of the surface |
1105 | | |
1106 | | OUTPUT: |
1107 | | |
1108 | | - The dictionary of lists [principal curvature, corresponding principal direction] |
1109 | | |
1110 | | If principal curvatures coincide, gives the warning that the surface is a sphere. |
1111 | | |
1112 | | EXAMPLES:: |
1113 | | |
1114 | | sage: u, v = var('u, v') |
1115 | | sage: R, r = var('R,r') |
1116 | | sage: assume(R>r,r>0) |
1117 | | sage: torus = ParametrizedSurface3D([(R+r*cos(v))*cos(u),(R+r*cos(v))*sin(u),r*sin(v)],[u,v],'torus') |
1118 | | sage: pdd = torus.principal_directions() |
1119 | | sage: pdd[1] |
1120 | | [-cos(v)/(r*cos(v) + R), (1, 0)] |
1121 | | sage: pdd[2] |
1122 | | [-1/r, (0, -(R*r*cos(v) + R^2)/r)] |
1123 | | |
1124 | | sage: R = var('R') |
1125 | | sage: assume(R>0) |
1126 | | sage: assume(cos(v)>0) |
1127 | | sage: sphere = ParametrizedSurface3D([R*cos(u)*cos(v),R*sin(u)*cos(v),R*sin(v)],[u,v],'sphere') |
1128 | | sage: sphere.principal_directions() |
1129 | | 'This is a sphere, so any direction is principal' |
1130 | | |
1131 | | sage: catenoid = ParametrizedSurface3D([R*cosh(v)*cos(u), R*cosh(v)*sin(u),v],[u,v],'catenoid') |
1132 | | sage: pd = catenoid.principal_directions() |
1133 | | sage: pd[1][1] |
1134 | | (1, 0) |
1135 | | sage: pd[2][1] |
1136 | | (0, 1/2*(sqrt(R^2*sinh(v)^2 + 1)*sqrt((4*sinh(v)^2*cosh(v)^2 + 1)*R^4 + 2*(2*cosh(v)^2 - 1)*R^2 + 1)*R*cosh(v) + sqrt(R^2*sinh(v)^2 + 1)*((2*sinh(v)^2*cosh(v) + cosh(v))*R^3 + R*cosh(v)))/(R^4*sinh(v)^4 + 2*R^2*sinh(v)^2 + 1)) |
1137 | | sage: pd[1][1]*pd[2][1] |
1138 | | 0 |
1139 | | """ |
1140 | | gg = self.first_fundamental_form_coefficients() |
1141 | | hh = self.second_fundamental_form_coefficients() |
1142 | | kk = self.principal_curvatures() |
1143 | | if kk[1]==kk[2]: |
1144 | | return "This is a sphere, so any direction is principal" |
1145 | | pd1 = simplify_vector([hh[(1,2)]-kk[1]*gg[(1,2)],-hh[(1,1)]+kk[1]*gg[(1,1)]]) |
1146 | | if pd1==vector([0,0]): |
1147 | | pd1 = vector([1,0]) |
1148 | | pd2 = simplify_vector([hh[(1,2)]-kk[2]*gg[(1,2)],-hh[(1,1)]+kk[2]*gg[(1,1)]]) |
1149 | | if pd2==vector([0,0]): |
1150 | | pd2 = vector([1,0]) |
1151 | | return {1:[kk[1],pd1],2:[kk[2],pd2]} |
1152 | | |
1153 | | |
1154 | | ### Jvkersch: doctests checked up to this point. |
1155 | | |
1156 | | # jvkersch: alternative definition of principal_directions. The behavior of this |
1157 | | # function is consistent with the eigenvector function in Sage (since it is just |
1158 | | # looking for eigenvectors of the shape operator). In particular, for a spherical |
1159 | | # surface it just returns an arbitrary pair of vectors. |
1160 | | def principal_directions_new(self): |