# Ticket #10132: trac_10132_diff_geom_sage_doctests.patch

File trac_10132_diff_geom_sage_doctests.patch, 55.9 KB (added by jvkersch, 10 years ago)

Latest version of the patch, with 100% coverage and doctests in place.

• ## doc/en/reference/index.rst

# HG changeset patch
# User Joris Vankerschaver <joris.vankerschaver@gmail.com>
# Date 1290619376 28800
# Node ID b834b692db4c5e7fd57c369e54f8763224d2c58a
# Parent  120c07be6358d93bcff503363d379c26b8342f2b
#10132: Differential Geometry via Sage

diff -r 120c07be6358 -r b834b692db4c doc/en/reference/index.rst
 a modabvar modmisc tensor riemann history_and_license
• ## new file doc/en/reference/riemann.rst

diff -r 120c07be6358 -r b834b692db4c doc/en/reference/riemann.rst
 - Differential Geometry of Curves and Surfaces ============================================ .. toctree:: :maxdepth: 2 sage/riemann/parametrized_surface3d
• ## sage/all.py

diff -r 120c07be6358 -r b834b692db4c sage/all.py
 a from sage.ext.fast_eval      import fast_float from sage.tensor.all     import * from sage.riemann.all    import * from copy import copy, deepcopy
• ## new file sage/riemann/__init__.py

diff -r 120c07be6358 -r b834b692db4c sage/riemann/__init__.py
 - # This comment is here so the file is non-empty (so Mercurial will check it in).
• ## new file sage/riemann/all.py

diff -r 120c07be6358 -r b834b692db4c sage/riemann/all.py
 - from parametrized_surface3d import ParametrizedSurface3D
• ## new file sage/riemann/parametrized_surface3d.py

diff -r 120c07be6358 -r b834b692db4c sage/riemann/parametrized_surface3d.py
 - """ Differential Geometry of Parametrized Surfaces jvkersch: we need to have some more examples of how to use this class here. AUTHORS: - Mikhail Malakhaltsev (2010-09-25): initial version """ #***************************************************************************** #       Copyright (C) 2010  Mikhail Malakhaltsev # #  Distributed under the terms of the GNU General Public License (GPL) #                  http://www.gnu.org/licenses/ #***************************************************************************** # The following is maybe too specific ... from vector_functions import * # mikarm: This is for the moment, at the final version will be replaced by import #attach('/home/prince/MY_SAGE_MODULES/vector_functions.py') from sage.structure.sage_object import SageObject from sage.calculus.functional import diff from sage.functions.other import sqrt from sage.misc.cachefunc import cached_method class GenericManifold(SageObject): """ Placeholder class representing a smooth manifold.  Currently this class serves as a base class for ParametrizedSurface3D, but does not have any other functionality. TODO: - Make this class aware of CoordinatePatch functionality. """ pass class ParametrizedSurface3D(GenericManifold): r""" This is a class for finding basic invariants of surfaces. The orientation of the surface is determined by the parametrization, that is, the natural frame with positive orientation is given by $\partial_1 \vec r$, $\partial_2 \vec r$. jvkersch: This description, as well as the doctests, should be expanded. INPUT: - surface_equation - list/vector specifying a parametric representation of the surface - variables - coordinates on the surface - name_of_surface - string with the name of the surface (optional). EXAMPLES:: sage: u, v = var('u,v') sage: paraboloid = ParametrizedSurface3D([u,v,u^2+v^2],[u,v],'paraboloid'); paraboloid Parametrized surface ('paraboloid') with equation [u, v, u^2 + v^2] """ def __init__(self,equation,variables,*name): """ See ParametrizedSurface3D for full documentation. EXAMPLES:: sage: u, v = var('u,v') sage: eq = [3*u + 3*u*v^2 - u^3, 3*v + 3*u^2*v - v^3, 3*(u^2-v^2)] sage: enneper = ParametrizedSurface3D(eq,[u,v],'enneper_surface') """ self.equation = equation self.variables_list = variables self.variables = {1:self.variables_list[0],2:self.variables_list[1]} self.name = name # Callable version of the underlying equation def eq_callable(u, v): u1, u2 = self.variables_list return [fun.subs({u1: u, u2: v}) for fun in self.equation] self.eq_callable = eq_callable # Various index tuples self.index_list   = [(i,j) for i in (1,2) for j in (1,2)] self.index_list_3 = [(i,j,k) for i in (1,2) for j in (1,2) for k in (1,2)] def _latex_(self): r""" Returns the LaTeX representation of this parametrized surface. EXAMPLES:: sage: u, v = var('u, v') sage: sphere = ParametrizedSurface3D([cos(u)*cos(v),sin(u)*cos(v),sin(v)],[u,v],'sphere') sage: latex(sphere) \left[\cos\left(u\right) \cos\left(v\right), \sin\left(u\right) \cos\left(v\right), \sin\left(v\right)\right] sage: sphere._latex_() \left[\cos\left(u\right) \cos\left(v\right), \sin\left(u\right) \cos\left(v\right), \sin\left(v\right)\right] """ from sage.misc.latex import latex return latex(self.equation) def _repr_(self): r""" Returns the string representation of this parametrized surface. EXAMPLES:: sage: u, v = var('u, v') sage: eq = [3*u + 3*u*v^2 - u^3, 3*v + 3*u^2*v - v^3, 3*(u^2-v^2)] sage: enneper = ParametrizedSurface3D(eq,[u,v],'enneper_surface') sage: print enneper Parametrized surface ('enneper_surface') with equation [-u^3 + 3*u*v^2 + 3*u, 3*u^2*v - v^3 + 3*v, 3*u^2 - 3*v^2] sage: enneper._repr_() "Parametrized surface ('enneper_surface') with equation [-u^3 + 3*u*v^2 + 3*u, 3*u^2*v - v^3 + 3*v, 3*u^2 - 3*v^2]" """ name = 'Parametrized surface' if self.name is not None: name += " ('%s')" % self.name return '%(designation)s with equation %(eq)s' % \ {'designation': name, 'eq': str(self.equation)} def plot3d(self, urange, vrange): """ Enable easy plotting directly from the surface class. INPUT: - urange - 3-tuple (u, u_min, u_max) - vrange - 3-tuple (v, v_min, v_max) EXAMPLES:: sage: u, v = var('u, v') sage: eq = [3*u + 3*u*v^2 - u^3, 3*v + 3*u^2*v - v^3, 3*(u^2-v^2)] sage: enneper = ParametrizedSurface3D(eq,[u,v],'enneper_surface') sage: enneper.plot3d((u, -5, 5), (v, -5, 5)) """ from sage.plot.plot3d.parametric_plot3d import parametric_plot3d P = parametric_plot3d(self.equation, urange, vrange) return P @cached_method def natural_frame(self): """ Returns the natural tangent frame on the parametrized surface.  The vectors of this frame are tangent to the coordinate lines on the surface. OUTPUT: - The natural frame as a dictionary. EXAMPLES:: sage: u, v = var('u, v') sage: eparaboloid = ParametrizedSurface3D([u,v,u^2+v^2],[u,v],'elliptic paraboloid') sage: eparaboloid.natural_frame() {1: (1, 0, 2*u), 2: (0, 1, 2*v)} """ dr1 = vector([diff(f,self.variables[1]).simplify_full() for f in self.equation]) dr2 = vector([diff(f,self.variables[2]).simplify_full() for f in self.equation]) return {1:dr1, 2:dr2} @cached_method def normal_vector(self, normalized=False): """ Returns the normal vector field of the parametrized surface. INPUT: - normalized - default False - specifies whether the normal vector should be normalized. OUTPUT: - Normal vector field. EXAMPLES:: sage: u, v = var('u, v') sage: eparaboloid = ParametrizedSurface3D([u,v,u^2+v^2],[u,v],'elliptic paraboloid') sage: eparaboloid.normal_vector(normalized=False) (-2*u, -2*v, 1) sage: eparaboloid.normal_vector(normalized=True) (-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)) """ dr = self.natural_frame() normal = dr[1].cross_product(dr[2]) if normalized: normal /= normal.norm() return simplify_vector(normal) @cached_method def _compute_first_fundamental_form_coefficient(self, index): """ Helper function to compute coefficients of the first fundamental form. Do not call this method directly; instead use first_fundamental_form_coefficient. This method expects its argument to be a list, and hence can be cached. EXAMPLES:: sage: u, v = var('u, v') sage: eparaboloid = ParametrizedSurface3D([u,v,u^2+v^2],[u,v]) sage: eparaboloid._compute_first_fundamental_form_coefficient((1,2)) 4*u*v """ dr = self.natural_frame() return (dr[index[0]]*dr[index[1]]).simplify_full() def first_fundamental_form_coefficient(self, index): r""" Compute a single component $g_{ij}$ of the first fundamental form.  If the parametric representation of the surface is given by the vector function $\vec r(u^i)$, where $u^i$, $i = 1, 2$ are curvilinear coordinates, then .. math:: g_{ij} = \frac{\partial \vec r}{\partial u^i} \cdot \frac{\partial \vec r}{\partial u^j}. INPUT: - index - length-2 index (i, j) of the component OUTPUT: - Component g_ij of the first fundamental form EXAMPLES:: sage: u, v = var('u, v') sage: eparaboloid = ParametrizedSurface3D([u,v,u^2+v^2],[u,v]) sage: eparaboloid.first_fundamental_form_coefficient((1,2)) 4*u*v When the index is invalid, an error is raised:: sage: u, v = var('u, v') sage: eparaboloid = ParametrizedSurface3D([u,v,u^2+v^2],[u,v]) sage: eparaboloid.first_fundamental_form_coefficient((1,5)) Traceback (most recent call last): ... ValueError: Index (1, 5) out of bounds. """ index = tuple(sorted(index)) if index not in self.index_list: raise ValueError, "Index %s out of bounds." % str(index) return self._compute_first_fundamental_form_coefficient(index) @cached_method def first_fundamental_form_coefficients(self): """ Returns the coefficients of the first fundamental form as a dictionary. The keys are tuples $(i, j)$, where $i$ and $j$ range over $1, 2$, while the values are the corresponding coefficients $g_{ij}$. OUTPUT: - Dictionary of first fundamental form coefficients. EXAMPLES:: sage: u, v = var('u,v') sage: sphere = ParametrizedSurface3D([cos(u)*cos(v),sin(u)*cos(v),sin(v)],[u,v],'sphere') sage: sphere.first_fundamental_form_coefficients() {(1, 2): 0, (1, 1): cos(v)^2, (2, 1): 0, (2, 2): 1} """ coefficients = {} for index in self.index_list: sorted_index = list(sorted(index)) coefficients[index] = \ self._compute_first_fundamental_form_coefficient(index) return coefficients def first_fundamental_form(self, vector1, vector2): r""" Evaluate the first fundamental form on two vectors expressed with respect to the natural coordinate frame on the surface. In other words, if the vectors are $v = (v^1, v^2)$ and $w = (w^1, w^2)$, calculate $g_{11} v^1 w^1 + g_{12}(v^1 w^2 + v^2 w^1) + g_{22} v^2 w^2$, with $g_{ij}$ the coefficients of the first fundamental form. INPUT: - vector1, vector2 - vectors on the surface. OUTPUT: - First fundamental form evaluated on the input vectors. EXAMPLES:: sage: u, v = var('u, v') sage: v1, v2, w1, w2 = var('v1, v2, w1, w2') sage: sphere = ParametrizedSurface3D([cos(u)*cos(v),sin(u)*cos(v),sin(v)],[u,v],'sphere') sage: sphere.first_fundamental_form(vector([v1,v2]),vector([w1,w2])) v1*w1*cos(v)^2 + v2*w2 sage: vv = vector([1,2]) sage: sphere.first_fundamental_form(vv,vv) cos(v)^2 + 4 sage: sphere.first_fundamental_form([1,1],[2,1]) 2*cos(v)^2 + 1 """ gg = self.first_fundamental_form_coefficients() return sum(gg[ind]*vector1[ind[0]-1]*vector2[ind[1]-1] for ind in self.index_list) @cached_method def area_form_squared(self): """ Returns the square of the coefficient of the area form on the surface. In terms of the coefficients $g_{ij}$ (where $i, j = 1, 2$) of the first fundamental form, this invariant is given by .. math:: A^2 = g_{11}g_{22} - g_{12}^2. See also area_form. OUTPUT: - Square of the area form EXAMPLES:: sage: u, v = var('u, v') sage: sphere = ParametrizedSurface3D([cos(u)*cos(v),sin(u)*cos(v),sin(v)],[u,v],'sphere') sage: sphere.area_form_squared() cos(v)^2 """ gg = self.first_fundamental_form_coefficients() return  (gg[(1,1)]*gg[(2,2)]-gg[(1,2)]**2).simplify_full() @cached_method def area_form(self): r""" Returns the coefficient of the area form on the surface.  In terms of the coefficients $g_{ij}$ (where $i, j = 1, 2$) of the first fundamental form, the coefficient of the area form is given by .. math:: A = \sqrt{g_{11}g_{22} - g_{12}^2}. See also area_form_squared. OUTPUT: - Coefficient of the area form EXAMPLES:: sage: u, v = var('u,v') sage: sphere = ParametrizedSurface3D([cos(u)*cos(v),sin(u)*cos(v),sin(v)],[u,v],'sphere') sage: sphere.area_form() abs(cos(v)) """ return sqrt(self.area_form_squared()).simplify_full() @cached_method def first_fundamental_form_inverse_coefficients(self): r""" Returns the coefficients $g^{ij}$ of the inverse of the fundamental form, as a dictionary.  The inverse coefficients are defined by .. math:: \sum_j g^{ij} g_{jk} = \delta^i_k with $\delta^i_k$ the Kronecker delta. OUTPUT: - Dictionary of the inverse coefficients. EXAMPLES:: sage: u, v = var('u, v') sage: sphere = ParametrizedSurface3D([cos(u)*cos(v),sin(u)*cos(v),sin(v)],[u,v],'sphere') sage: sphere.first_fundamental_form_inverse_coefficients() {(1, 2): 0, (1, 1): cos(v)^(-2), (2, 1): 0, (2, 2): 1} """ gg = self.first_fundamental_form_coefficients() DD = gg[(1,1)]*gg[(2,2)]-gg[(1,2)]**2 gi11 = (gg[(2,2)]/DD).simplify_full() gi12 = (-gg[(1,2)]/DD).simplify_full() gi21 = gi12 gi22 = (gg[(1,1)]/DD).simplify_full() return {(1,1): gi11, (1,2): gi12, (2,1): gi21, (2,2): gi22} def first_fundamental_form_inverse_coefficient(self, index): r""" Returns a specific component $g^{ij}$ of the inverse of the fundamental form. INPUT: - index - index (i, j) of the component g^{ij}. OUTPUT: - Component of the inverse of the fundamental form. EXAMPLES:: sage: u, v = var('u, v') sage: sphere = ParametrizedSurface3D([cos(u)*cos(v),sin(u)*cos(v),sin(v)],[u,v],'sphere') sage: sphere.first_fundamental_form_inverse_coefficient((1, 2)) 0 sage: sphere.first_fundamental_form_inverse_coefficient((1, 1)) cos(v)^(-2) """ index = tuple(sorted(index)) if index not in self.index_list: raise ValueError, "Index %s out of bounds." % str(index) return self.first_fundamental_form_inverse_coefficients()[index] @cached_method def rotation(self,theta): r""" Gives the matrix of the rotation operator over a given angle $\theta$ with respect to the natural frame. INPUT: - theta - rotation angle OUTPUT: - Rotation matrix with respect to the natural frame. ALGORITHM: The operator of rotation over $\pi/2$ is $J^i_j = g^{ik}\omega_{jk}$, where $\omega$ is the area form.  The operator of rotation over an angle $\theta$ is $\cos(\theta) I + sin(\theta) J$ EXAMPLES:: sage: u, v = var('u, v') sage: assume(cos(v)>0) sage: sphere = ParametrizedSurface3D([cos(u)*cos(v),sin(u)*cos(v),sin(v)],[u,v],'sphere') We first compute the matrix of rotation over pi/3:: sage: rotation = sphere.rotation(pi/3); rotation [                1/2 -1/2*sqrt(3)/cos(v)] [ 1/2*sqrt(3)*cos(v)                 1/2] We verify that three succesive rotations over pi/3 yield minus the identity:: sage: rotation^3 [-1  0] [ 0 -1] """ from sage.functions.trig import sin, cos gi = self.first_fundamental_form_inverse_coefficients() w12 = self.area_form() RR11 = (cos(theta) + sin(theta)*gi[1,2]*w12).simplify_full() RR12 = (- sin(theta)*gi[1,1]*w12).simplify_full() RR21 = (sin(theta)*gi[2,2]*w12).simplify_full() RR22 = (cos(theta) - sin(theta)*gi[2,1]*w12).simplify_full() rotation = matrix([[RR11,RR12],[RR21,RR22]]) return rotation @cached_method def orthonormal_frame(self, coordinates='ext'): r""" Returns the orthonormal frame field on the surface, expressed either in exterior coordinates (i.e. expressed as vector fields in the ambient space $\mathbb{R}^3$, the default) or interior coordinates (with respect to the natural frame) INPUT: - coordinates - either ext (default) or int. OUTPUT: - Orthogonal frame field as a dictionary. ALGORITHM: We normalize the first vector $\vec e_1$ of the natural frame and then get the second frame vector as $\vec e_2 = [\vec n, \vec e_1]$, where $\vec n$ is the unit normal to the surface. EXAMPLES:: sage: u, v = var('u,v') sage: assume(cos(v)>0) sage: sphere = ParametrizedSurface3D([cos(u)*cos(v), sin(u)*cos(v), sin(v)], [u, v],'sphere') sage: frame = sphere.orthonormal_frame(); frame {1: (-sin(u), cos(u), 0), 2: (-sin(v)*cos(u), -sin(u)*sin(v), cos(v))} sage: (frame[1]*frame[1]).simplify_full() 1 sage: (frame[1]*frame[2]).simplify_full() 0 sage: frame[1] == sphere.orthonormal_frame_vector(1) True We compute the orthonormal frame with respect to the natural frame on the surface:: sage: frame_int = sphere.orthonormal_frame(coordinates='int'); frame_int {1: (1/cos(v), 0), 2: (0, 1)} sage: sphere.first_fundamental_form(frame_int[1], frame_int[1]) 1 sage: sphere.first_fundamental_form(frame_int[1], frame_int[2]) 0 sage: sphere.first_fundamental_form(frame_int[2], frame_int[2]) 1 """ from sage.symbolic.constants import pi if coordinates not in ['ext', 'int']: raise ValueError, \ r"Coordinate system must be exterior ('ext') or interior ('int')." c  = self.first_fundamental_form_coefficient([1,1]) if coordinates == 'ext': f1 = self.natural_frame()[1] E1 = simplify_vector(f1/sqrt(c)) E2 = simplify_vector(self.normal_vector(normalized=True).cross_product(E1)) else: E1 =  vector([(1/sqrt(c)).simplify_full(), 0]) E2 = simplify_vector(self.rotation(pi/2)*E1) return  {1:E1, 2:E2} @cached_method def orthonormal_frame_vector(self, index, coordinates='ext'): r""" Returns a specific basis vector field of the orthonormal frame field on the surface, expressed in exterior or interior coordinates.  See orthogonal_frame for more details. INPUT: - index - index (i, j) of the basis vector; - coordinates - either ext (default) or int. OUTPUT: - Orthonormal frame vector field. EXAMPLES:: sage: u, v = var('u, v') sage: assume(cos(v)>0) sage: sphere = ParametrizedSurface3D([cos(u)*cos(v),sin(u)*cos(v),sin(v)],[u,v],'sphere') sage: V1 = sphere.orthonormal_frame_vector(1); V1 (-sin(u), cos(u), 0) sage: V2 = sphere.orthonormal_frame_vector(2); V2 (-sin(v)*cos(u), -sin(u)*sin(v), cos(v)) sage: (V1*V1).simplify_full() 1 sage: (V1*V2).simplify_full() 0 sage: from sage.riemann.vector_functions import simplify_vector sage: n = sphere.normal_vector(normalized=True) sage: simplify_vector(V1.cross_product(V2) - n) (0, 0, 0) """ return self.orthonormal_frame(coordinates)[index] def lie_bracket(self, v, w): """ Returns the Lie bracket of two vector fields that are tangent to the surface. The vector fields should be given in coordinates with respect to the natural frame. INPUT: - v and w - vector fields on the surface, expressed as tuples of functions. OUTPUT: - The Lie bracket [v, w] EXAMPLES:: sage: u, v = var('u, v') sage: assume(cos(v)>0) sage: sphere = ParametrizedSurface3D([cos(u)*cos(v),sin(u)*cos(v),sin(v)],[u,v],'sphere') sage: sphere.lie_bracket([u,v],[-v,u]) (0, 0) sage: EE_int = sphere.orthonormal_frame(coordinates='int') sage: sphere.lie_bracket(EE_int[1],EE_int[2]) (sin(v)/cos(v)^2, 0) """ vv = vector([xx for xx in v]) ww = vector([xx for xx in w]) uuu = self.variables_list Dvv = matrix([ [ diff(xxx,uu).simplify_full() for uu in uuu ] for xxx in vv]) Dww = matrix([ [ diff(xxx,uu).simplify_full() for uu in uuu ] for xxx in ww]) return simplify_vector(vector(Dvv*ww - Dww*vv)) def frame_structure_functions(self, e1, e2): r""" Returns the structure functions $c^k_{ij}$ for a frame field $e_1, e_2$, i.e. a pair of vector fields on the surface which are linearly independent at each point.  The structure functions are defined using the Lie bracket by .. math:: [e_i,e_j] = c^k_{ij}e_k. INPUT: - e1, e2 - vector fields on the surface, expressed as pairs of functions. OUTPUT: Dictionary of structure functions, where the key (i, j, k) refers to the function $c^k_{ij}$. EXAMPLES:: sage: u, v = var('u, v') sage: assume(cos(v) > 0) sage: sphere = ParametrizedSurface3D([cos(u)*cos(v),sin(u)*cos(v),sin(v)],[u,v],'sphere') sage: sphere.frame_structure_functions([u,v],[-v,u]) {(1, 2, 1): 0, (2, 1, 2): 0, (2, 2, 2): 0, (1, 2, 2): 0, (1, 1, 1): 0, (2, 1, 1): 0, (2, 2, 1): 0, (1, 1, 2): 0} We construct the structure functions of the orthonormal frame on the surface:: sage: EE_int = sphere.orthonormal_frame(coordinates='int') sage: CC = sphere.frame_structure_functions(EE_int[1],EE_int[2]); CC {(1, 2, 1): sin(v)/cos(v), (2, 1, 2): 0, (2, 2, 2): 0, (1, 2, 2): 0, (1, 1, 1): 0, (2, 1, 1): -sin(v)/cos(v), (2, 2, 1): 0, (1, 1, 2): 0} sage: sphere.lie_bracket(EE_int[1],EE_int[2]) - CC[(1,2,1)]*EE_int[1] - CC[(1,2,2)]*EE_int[2] (0, 0) """ ee1 = vector([xx for xx in e1]) ee2 = vector([xx for xx in e2]) lb = simplify_vector(self.lie_bracket(ee1,ee2)) trans_matrix = matrix([[ee1[0],ee2[0]],[ee1[1],ee2[1]]]) zz = simplify_vector(trans_matrix.inverse()*lb) return {(1,1,1):0, (1,1,2):0, (1,2,1):zz[0], (1,2,2):zz[1], (2,1,1):-zz[0], (2,1,2):-zz[1],(2,2,1):0, (2,2,2):0} @cached_method def _compute_second_order_frame_element(self, index): """ Compute an element of the second order frame of the surface.  See second_order_natural_frame for more details. This method expects its arguments in tuple form for caching.  As it does no input checking, it should not be called directly. EXAMPLES:: sage: u, v = var('u, v') sage: paraboloid = ParametrizedSurface3D([u, v, u^2 + v^2], [u,v], 'paraboloid') sage: paraboloid._compute_second_order_frame_element((1, 2)) (0, 0, 0) sage: paraboloid._compute_second_order_frame_element((2, 2)) (0, 0, 2) """ variables = [self.variables[i] for i in index] ddr_element = vector([diff(f, variables).simplify_full() for f in self.equation]) return ddr_element @cached_method def second_order_natural_frame(self): r""" Returns the second-order frame of the surface, i.e. computes the second-order derivatives (with respect to the parameters on the surface) of the parametric expression $\vec r = \vec r(u^1,u^2)$ of the surface. OUTPUT: - Dictionary where the keys are 2-tuples (i, j) and the values are the corresponding derivatives r_ij. EXAMPLES: We compute the second-order natural frame of the sphere:: sage: u, v = var('u, v') sage: sphere = ParametrizedSurface3D([cos(u)*cos(v),sin(u)*cos(v),sin(v)],[u,v],'sphere') sage: sphere.second_order_natural_frame() {(1, 2): (sin(u)*sin(v), -sin(v)*cos(u), 0), (1, 1): (-cos(u)*cos(v), -sin(u)*cos(v), 0), (2, 1): (sin(u)*sin(v), -sin(v)*cos(u), 0), (2, 2): (-cos(u)*cos(v), -sin(u)*cos(v), -sin(v))} """ vectors = {} for index in self.index_list: sorted_index = tuple(sorted(index)) vectors[index] = \ self._compute_second_order_frame_element(sorted_index) return vectors def second_order_natural_frame_element(self, index): r""" Returns a vector in the second-order frame of the surface, i.e. computes the second-order derivatives of the parametric expression $\vec{r}$ of the surface with respect to the parameters listed in the argument. INPUT: - index - a 2-tuple (i, j) specifying the element of the second-order frame. OUTPUT: - The second-order derivative r_ij. EXAMPLES:: sage: u, v = var('u, v') sage: sphere = ParametrizedSurface3D([cos(u)*cos(v),sin(u)*cos(v),sin(v)],[u,v],'sphere') sage: sphere.second_order_natural_frame_element((1, 2)) (sin(u)*sin(v), -sin(v)*cos(u), 0) """ index = tuple(sorted(index)) if index not in self.index_list: raise ValueError, "Index %s out of bounds." % str(index) return self._compute_second_order_frame_element(index) @cached_method def _compute_second_fundamental_form_coefficient(self, index): """ Compute a coefficient of the second fundamental form of the surface.  See second_fundamental_form_coefficient for more details. This method expects its arguments in tuple form for caching.  As it does no input checking, it should not be called directly. EXAMPLES:: sage: u, v = var('u,v') sage: paraboloid = ParametrizedSurface3D([u, v, u^2+v^2], [u, v], 'paraboloid') sage: paraboloid._compute_second_fundamental_form_coefficient((1,1)) 2/sqrt(4*u^2 + 4*v^2 + 1) """ NN = self.normal_vector(normalized=True) v  = self.second_order_natural_frame_element(index) return (v*NN).simplify_full() def second_fundamental_form_coefficient(self, index): r""" Returns the coefficient $h_{ij}$ of the second fundamental form corresponding to the index $(i, j)$.  If the equation of the surface is $\vec{r}(u^1, u^2)$, then .. math:: h_{ij} = \vec{r}_{u^i u^j} \cdot \vec{n}, where $\vec{n}$ is the unit normal. INPUT: - index - a 2-tuple (i, j) OUTPUT: - Component h_{ij} of the second fundamental form. EXAMPLES:: sage: u, v = var('u,v') sage: assume(cos(v)>0) sage: sphere = ParametrizedSurface3D([cos(u)*cos(v),sin(u)*cos(v),sin(v)],[u,v],'sphere') sage: sphere.second_fundamental_form_coefficient((1, 1)) -cos(v)^2 sage: sphere.second_fundamental_form_coefficient((2, 1)) 0 """ index = tuple(index) if index not in self.index_list: raise ValueError, "Index %s out of bounds." % str(index) return self._compute_second_fundamental_form_coefficient(index) @cached_method def second_fundamental_form_coefficients(self): """ Returns the coefficients $h_{ij}$ of the second fundamental form as a dictionary, where the keys are the indices $(i, j)$ and the values are the corresponding components $h_{ij}$. When only one component is needed, consider the function second_fundamental_form_coefficient. OUTPUT: Dictionary of second fundamental form coefficients. EXAMPLES:: sage: u, v = var('u, v') sage: assume(cos(v)>0) sage: sphere = ParametrizedSurface3D([cos(u)*cos(v),sin(u)*cos(v),sin(v)],[u,v],'sphere') sage: sphere.second_fundamental_form_coefficients() {(1, 2): 0, (1, 1): -cos(v)^2, (2, 1): 0, (2, 2): -1} """ coefficients = {} for index in self.index_list: coefficients[index] = \ self._compute_second_fundamental_form_coefficient(index) return coefficients def second_fundamental_form(self,vector1,vector2): r""" Evaluates the second fundamental form on two vectors on the surface. If the vectors are given by $v=(v^1,v^2)$ and $w=(w^1,w^2)$, the result of this function is $h_{11} v^1 w^1 + h_{12}(v^1 w^2 + v^2 w^1) + h_{22} v^2 w^2$. INPUT: - vector1, vector2 - 2-tuples representing the input vectors. OUTPUT: - Value of the second fundamental form evaluated on the given vectors. EXAMPLES: We evaluate the second fundamental form on two symbolic vectors:: sage: u, v = var('u, v') sage: v1, v2, w1, w2 = var('v1, v2, w1, w2') sage: assume(cos(v) > 0) sage: sphere = ParametrizedSurface3D([cos(u)*cos(v),sin(u)*cos(v),sin(v)],[u,v],'sphere') sage: sphere.second_fundamental_form(vector([v1, v2]), vector([w1, w2])) -v1*w1*cos(v)^2 - v2*w2 We evaluate the second fundamental form on vectors with numerical components:: sage: vect = vector([1,2]) sage: sphere.second_fundamental_form(vect, vect) -cos(v)^2 - 4 sage: sphere.second_fundamental_form([1,1], [2,1]) -2*cos(v)^2 - 1 """ hh = self.second_fundamental_form_coefficients() return sum(hh[ind]*vector1[ind[0]-1]*vector2[ind[1]-1] for ind in self.index_list) @cached_method def gauss_curvature(self): r""" Finds the gaussian curvature of the surface, given by $K = \frac{h_{11}h_{22} - h_{12}^2}{g_{11}g_{22} - g_{12}^2}$, where $g_{ij}$ and $h_{ij}$ are the coefficients of the first and second fundamental form, respectively. OUTPUT: - Gaussian curvature of the surface. EXAMPLES:: sage: R = var('R') sage: assume(R>0) sage: u, v = var('u,v') sage: assume(cos(v)>0) sage: sphere = ParametrizedSurface3D([R*cos(u)*cos(v),R*sin(u)*cos(v),R*sin(v)],[u,v],'sphere') sage: sphere.gauss_curvature() R^(-2) """ hh = self.second_fundamental_form_coefficients() return ((hh[(1,1)]*hh[(2,2)]-hh[(1,2)]**2)/self.area_form_squared()).simplify_full() @cached_method def mean_curvature(self): r""" Finds the mean curvature of the surface, given by $H = \frac{1}{2}\frac{g_{22}h_{11} - 2g_{12}h_{12} + g_{11}h_{22}}{g_{11}g_{22} - g_{12}^2}$, where $g_{ij}$ and $h_{ij}$ are the components of the first and second fundamental forms, respectively. OUTPUT: - Mean curvature of the surface EXAMPLES:: sage: R = var('R') sage: assume(R>0) sage: u, v = var('u,v') sage: assume(cos(v)>0) sage: sphere = ParametrizedSurface3D([R*cos(u)*cos(v),R*sin(u)*cos(v),R*sin(v)],[u,v],'sphere') sage: sphere.mean_curvature() -1/R """ gg = self.first_fundamental_form_coefficients() hh = self.second_fundamental_form_coefficients() denom = 2*self.area_form_squared() enum =(gg[(2,2)]*hh[(1,1)]-2*gg[(1,2)]*hh[(1,2)]+gg[(1,1)]*hh[(2,2)]).simplify_full() return (enum/denom).simplify_full() @cached_method def principal_curvatures(self): """ Finds the principal curvatures of the surface OUTPUT: The two principal curvatures, given as a dictionary. EXAMPLES:: sage: R = var('R') sage: assume(R>0) sage: u, v = var('u,v') sage: assume(cos(v)>0) sage: sphere = ParametrizedSurface3D([R*cos(u)*cos(v),R*sin(u)*cos(v),R*sin(v)],[u,v],'sphere') sage: sphere.principal_curvatures() {1: -1/R, 2: -1/R} sage: u, v = var('u,v') sage: R, r = var('R, r') sage: assume(R>r,r>0) sage: torus = ParametrizedSurface3D([(R+r*cos(v))*cos(u),(R+r*cos(v))*sin(u),r*sin(v)],[u,v],'torus') sage: torus.principal_curvatures() {1: -cos(v)/(r*cos(v) + R), 2: -1/r} """ from sage.symbolic.assumptions import assume from sage.symbolic.relation import solve from sage.calculus.var import var KK = self.gauss_curvature() HH = self.mean_curvature() # jvkersch: when this assumption is uncommented, Sage raises an error stating that the assumption # is redundant... Can we safely omit this, based on some geometric reasoning? # assume(HH**2-KK>=0) # mikarm: This is a problem, I had a lot of trobles here. Of course, this inequality always hold true. # I insert this assumption because sage sometimes, in simplification, uses the complex numbers, though the roots # are, for sure, real. # This, in turn, causes problems when we substitute coordinates into the expression of principal curvatures. # Unfortunately, I did not manage to tell Sage that they are real (to declare the variables as real). # Moreover, in a neighborhood of an umbilic point we even cannot "smoothly" order the set of principal curvatures. # So, in general, at present this method is far from the final form. x = var('x') sol = solve(x**2 -2*HH*x + KK == 0,x) #k1=var('k1') #k2=var('k2') # 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. k1 = (x.substitute(sol[0])).simplify_full() if len(sol) == 1: k2 = k1 else: k2 = (x.substitute(sol[1])).simplify_full() return {1:k1,2: k2} @cached_method def shape_operator_coefficients(self): r""" Returns the components of the shape operator of the surface as a dictionary. See shape_operator for more information. OUTPUT: - Dictionary where the keys are two-tuples (i, j), with values the corresponding component of the shape operator. EXAMPLES:: sage: R = var('R') sage: u, v = var('u,v') sage: assume(cos(v)>0) sage: sphere = ParametrizedSurface3D([R*cos(u)*cos(v),R*sin(u)*cos(v),R*sin(v)],[u,v],'sphere') sage: sphere.shape_operator_coefficients() {(1, 2): 0, (1, 1): -1/R, (2, 1): 0, (2, 2): -1/R} """ gi = self.first_fundamental_form_inverse_coefficients() hh = self.second_fundamental_form_coefficients() shop11 = (gi[(1,1)]*hh[(1,1)] + gi[(1,2)]*hh[(1,2)]).simplify_full() shop12 = (gi[(1,1)]*hh[(2,1)] + gi[(1,2)]*hh[(2,2)]).simplify_full() shop21 = (gi[(2,1)]*hh[(1,1)] + gi[(2,2)]*hh[(1,2)]).simplify_full() shop22 = (gi[(2,1)]*hh[(2,1)] + gi[(2,2)]*hh[(2,2)]).simplify_full() return {(1,1): shop11, (1,2): shop12, (2,1): shop21, (2,2): shop22} # jvkersch: I've changed the definition of the shape operator to return the matrix # of the shape operator, rather than a matrix times a vector.  This will make it easier # to compute eigenvalues and eigenvectors of the shape operator. def shape_operator(self): r""" Returns the shape operator of the surface as a matrix.  The shape operator is defined as the derivative of the Gauss map, and is computed here in terms of the first and second fundamental form by means of the Weingarten equations. OUTPUT: - Matrix of the shape operator EXAMPLES:: sage: R = var('R') sage: assume(R>0) sage: u, v = var('u,v') sage: assume(cos(v)>0) sage: sphere = ParametrizedSurface3D([R*cos(u)*cos(v),R*sin(u)*cos(v),R*sin(v)],[u,v],'sphere') sage: S = sphere.shape_operator(); S [-1/R    0] [   0 -1/R] The eigenvalues of the shape operator are the principal curvatures of the surface:: sage: u, v = var('u,v') sage: paraboloid = ParametrizedSurface3D([u, v, u^2+v^2], [u, v], 'paraboloid') sage: S = paraboloid.shape_operator(); S [2*(4*v^2 + 1)/(4*u^2 + 4*v^2 + 1)^(3/2)        -8*u*v/(4*u^2 + 4*v^2 + 1)^(3/2)] [       -8*u*v/(4*u^2 + 4*v^2 + 1)^(3/2) 2*(4*u^2 + 1)/(4*u^2 + 4*v^2 + 1)^(3/2)] sage: S.eigenvalues() [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), 2/sqrt(4*u^2 + 4*v^2 + 1)] sage: paraboloid.principal_curvatures() {1: 2/(4*u^2 + 4*v^2 + 1)^(3/2), 2: 2/sqrt(4*u^2 + 4*v^2 + 1)} """ #vv = vector([xx for xx in v]) shop = self.shape_operator_coefficients() shop_matrix=matrix([[shop[(1,1)],shop[(1,2)]],[shop[(2,1)],shop[(2,2)]]]) return simplify_matrix(shop_matrix) @cached_method def principal_directions(self): """ Finds the principal curvatures and principal directions of the surface OUTPUT: - The dictionary of lists [principal curvature, corresponding principal direction] If principal curvatures coincide, gives the warning that the surface is a sphere. EXAMPLES:: sage: u, v = var('u, v') sage: R, r = var('R,r') sage: assume(R>r,r>0) sage: torus = ParametrizedSurface3D([(R+r*cos(v))*cos(u),(R+r*cos(v))*sin(u),r*sin(v)],[u,v],'torus') sage: pdd = torus.principal_directions() sage: pdd[1] [-cos(v)/(r*cos(v) + R), (1, 0)] sage: pdd[2] [-1/r, (0, -(R*r*cos(v) + R^2)/r)] sage: R = var('R') sage: assume(R>0) sage: assume(cos(v)>0) sage: sphere = ParametrizedSurface3D([R*cos(u)*cos(v),R*sin(u)*cos(v),R*sin(v)],[u,v],'sphere') sage: sphere.principal_directions() 'This is a sphere, so any direction is principal' sage: catenoid = ParametrizedSurface3D([R*cosh(v)*cos(u), R*cosh(v)*sin(u),v],[u,v],'catenoid') sage: pd = catenoid.principal_directions() sage: pd[1][1] (1, 0) sage: pd[2][1] (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)) sage: pd[1][1]*pd[2][1] 0 """ gg = self.first_fundamental_form_coefficients() hh = self.second_fundamental_form_coefficients() kk = self.principal_curvatures() if kk[1]==kk[2]: return "This is a sphere, so any direction is principal" pd1 = simplify_vector([hh[(1,2)]-kk[1]*gg[(1,2)],-hh[(1,1)]+kk[1]*gg[(1,1)]]) if pd1==vector([0,0]): pd1 = vector([1,0]) pd2 = simplify_vector([hh[(1,2)]-kk[2]*gg[(1,2)],-hh[(1,1)]+kk[2]*gg[(1,1)]]) if pd2==vector([0,0]): pd2 = vector([1,0]) return {1:[kk[1],pd1],2:[kk[2],pd2]} ### Jvkersch: doctests checked up to this point. # jvkersch: alternative definition of principal_directions.  The behavior of this # function is consistent with the eigenvector  function in Sage (since it is just # looking for eigenvectors of the shape operator).  In particular, for a spherical # surface it just returns an arbitrary pair of vectors. def principal_directions_new(self): r""" Finds the principal curvatures and principal directions of the surface OUTPUT: For each principal curvature, returns a list of the form (rho, V, n), where rho is the principal curvature, V is the corresponding principal direction, and n is the multiplicity. EXAMPLES:: sage: u, v = var('u, v') sage: R, r = var('R,r') sage: assume(R>r,r>0) sage: torus = ParametrizedSurface3D([(R+r*cos(v))*cos(u),(R+r*cos(v))*sin(u),r*sin(v)],[u,v],'torus') sage: torus.principal_directions_new() [(-cos(v)/(r*cos(v) + R), [(1, 0)], 1), (-1/r, [(0, 1)], 1)] """ return self.shape_operator().eigenvectors_left() @cached_method def connection_coefficients(self): r""" Computes the connection coefficients or Christoffel symbols $\Gamma^k_{ij}$ of the surface. If the coefficients of the first fundamental form are given by $g_{ij}$ (where $i, j = 1, 2$), then .. math:: \Gamma^k_{ij} = \frac{1}{2} g^{kl} \left( \frac{\partial g_{li}}{\partial x^j} - \frac{\partial g_{ij}}{\partial x^l} + \frac{\partial g_{lj}}{\partial x^i} \right). Here, $(g^{kl})$ is the inverse of the matrix $(g_{ij})$, with $i, j = 1, 2$. OUTPUT: Dictionary of connection coefficients, where the keys are 3-tuples $(i,j,k)$ and the values are the corresponding coefficients $\Gamma^k_{ij}$. EXAMPLES:: sage: r = var('r') sage: assume(r > 0) sage: u, v = var('u,v') sage: assume(cos(v)>0) sage: sphere = ParametrizedSurface3D([r*cos(u)*cos(v),r*sin(u)*cos(v),r*sin(v)],[u,v],'sphere') sage: sphere.connection_coefficients() {(1, 2, 1): -sin(v)/cos(v), (2, 2, 2): 0, (1, 2, 2): 0, (2, 1, 1): -sin(v)/cos(v), (1, 1, 2): sin(v)*cos(v), (2, 2, 1): 0, (2, 1, 2): 0, (1, 1, 1): 0} """ x = self.variables gg = self.first_fundamental_form_coefficients() gi = self.first_fundamental_form_inverse_coefficients() dg = {} for i,j,k in self.index_list_3: dg[(i,j,k)]=gg[(j,k)].differentiate(x[i]).simplify_full() structfun={} for i,j,k in self.index_list_3: structfun[(i,j,k)]=sum(gi[(k,s)]*(dg[(i,j,s)]+dg[(j,i,s)]-dg[(s,i,j)])/2 for s in (1,2)).full_simplify() return structfun @cached_method def _create_geodesic_ode_system(self): r""" Helper method to create a fast floating-point version of the geodesic equations, used by geodesics_numerical. EXAMPLES:: sage: p, q = var('p,q') sage: sphere = ParametrizedSurface3D([cos(q)*cos(p),sin(q)*cos(p),sin(p)],[p,q],'sphere') sage: ode = sphere._create_geodesic_ode_system() sage: ode.function(0.0, (1.0, 0.0, 1.0, 1.0)) [1.00000000000000, 1.00000000000000, -0.45464871341284091, 3.1148154493098041] """ from sage.ext.fast_eval import fast_float from sage.calculus.var import var from sage.gsl.ode import ode_solver u1 = self.variables[1] u2 = self.variables[2] v1, v2 = var('v1, v2') C = self.connection_coefficients() dv1 = - C[(1,1,1)]*v1**2 - 2*C[(1,2,1)]*v1*v2 - C[(2,2,1)]*v2**2 dv2 = - C[(1,1,2)]*v1**2 - 2*C[(1,2,2)]*v1*v2 - C[(2,2,2)]*v2**2 fun1 = fast_float(dv1, str(u1), str(u2), str(v1), str(v2)) fun2 = fast_float(dv2, str(u1), str(u2), str(v1), str(v2)) geodesic_ode = ode_solver() geodesic_ode.function = \ lambda t, (u1, u2, v1, v2) : \ [v1, v2, fun1(u1, u2, v1, v2), fun2(u1, u2, v1, v2)] return geodesic_ode def geodesics_numerical(self, p0, v0, tinterval): r""" Numerical integration of the geodesic equations.  Explicitly, the geodesic equations are given by .. math:: \frac{d^2 u^i}{dt^2} + \Gamma^i_{jk} \frac{d u^j}{dt} \frac{d u^k}{dt} = 0. Solving these equations gives the coordinates $(u^1, u^2)$ on the surface of the geodesic.  The coordinates in space can then be found by substituting $(u^1, u^2)$ into the vector representating $\vec{r}(u^1, u^2)$ of the surface. ALGORITHM: The geodesic equations are integrated forward in time using the ode solvers from sage.gsl.ode.  See _create_geodesic_ode_system for more details. INPUT: - p0 - 2-tuple with coordinates of the initial point. - v0 - 2-tuple with components of the initial tangent vector to the geodesic. - tinterval - List [a,b,M], where (a,b) is the domain of the geodesic and M is the number of division points OUTPUT: List [t, [u1(t),u2(t)], [v1(t),v2(t)], [x1(t),x2(t),x3(t)]], where - t is a subdivision point; - [u1(t),u2(t)] is the list of coordinates of the geodesic point; - [v1(t),v2(t)] is the list of coordinates of the vector tanget to the geodesic; - [x1(t),x2(t),x3(t)] is the list of coordinates of the geodesic point in the three-dimensional space. EXAMPLES:: sage: p, q = var('p,q') sage: v = [p,q] sage: assume(cos(q)>0) sage: sphere = ParametrizedSurface3D([cos(q)*cos(p),sin(q)*cos(p),sin(p)],v,'sphere') sage: gg_array = sphere.geodesics_numerical([0.0,0.0],[1.0,1.0],[0,2*pi,5]) sage: gg_array[0] [0, [0.000000000000000, 0.000000000000000], [1.00000000000000, 1.00000000000000], [1, 0, 0]] sage: gg_array[1] [1.2566370614359172, [0.76440092815484362, 1.8586224292405702], [-0.28386842687533731, 1.9194187166389509], [-0.204895409519981, 0.692104714174527, 0.692104714458033]] """ u1 = self.variables[1] u2 = self.variables[2] solver = self._create_geodesic_ode_system() t_interval, n = tinterval[0:2], tinterval[2] solver.y_0 = [p0[0], p0[1], v0[0], v0[1]] solver.ode_solve(t_span=t_interval, num_points=n) parsed_solution = \ [[vec[0], vec[1][0:2], vec[1][2:], self.eq_callable(vec[1][0], vec[1][1])] \ for vec in solver.solution] return parsed_solution @cached_method def _create_pt_ode_system(self, curve, t): """ Helper method to create a fast floating-point version of the parallel transport equations, used by parallel_translation_numerical. INPUT: - curve - curve in intrinsic coordinates along which to do parallel transport. - t - curve parameter EXAMPLES:: sage: p, q = var('p,q') sage: sphere = ParametrizedSurface3D([cos(q)*cos(p),sin(q)*cos(p),sin(p)],[p,q],'sphere') sage: s = var('s') sage: ode = sphere._create_pt_ode_system((s, s), s) sage: ode.function(0.0, (1.0, 1.0)) [-0.0, 0.0] """ from sage.ext.fast_eval import fast_float from sage.calculus.var import var from sage.gsl.ode import ode_solver u1 = self.variables[1] u2 = self.variables[2] v1, v2 = var('v1, v2') du1 = diff(curve[0], t) du2 = diff(curve[1], t) C = self.connection_coefficients() for coef in C: C[coef] = C[coef].subs({u1: curve[0], u2: curve[1]}) dv1 = - C[(1,1,1)]*v1*du1 - C[(1,2,1)]*(du1*v2 + du2*v1) - C[(2,2,1)]*du2*v2 dv2 = - C[(1,1,2)]*v1*du1 - C[(1,2,2)]*(du1*v2 + du2*v1) - C[(2,2,2)]*du2*v2 fun1 = fast_float(dv1, str(t), str(v1), str(v2)) fun2 = fast_float(dv2, str(t), str(v1), str(v2)) pt_ode = ode_solver() pt_ode.function = lambda t, (v1, v2): [fun1(t, v1, v2), fun2(t, v1, v2)] return pt_ode def parallel_translation_numerical(self,curve,t,v0,tinterval): r""" Numerically solves the equations for parallel translation of a vector along a curve on the surface.  Explicitly, the equations for parallel translation are given by .. math:: \frac{d u^i}{dt} + u^j \frac{d c^k}{dt} \Gamma^i_{jk} = 0, where $\Gamma^i_{jk}$ are the connection coefficients of the surface, the vector to be transported has components $u^j$ and the curve along which to transport has components $c^k$. ALGORITHM: The parallel transport equations are integrated forward in time using the ode solvers from sage.gsl.ode.  See _create_pt_ode_system for more details. INPUT: - curve - 2-tuple of functions which determine the curve with respect to the local coordinate; - t - symbolic variable denoting the curve parameter; - v0 - 2-tuple representing the initial vector; - tinterval - list [a,b,N], where (a,b) is the domain of the curve and N is the number of subdivision points. OUTPUT: The list consisting of lists [t, [v1(t),v2(t)]], where - t is a subdivision point; - [v1(t),v2(t)] is the list of coordinates of the vector parallel translated along the curve. EXAMPLES:: sage: p, q = var('p,q') sage: v = [p,q] sage: assume(cos(q)>0) sage: sphere = ParametrizedSurface3D([cos(q)*cos(p),sin(q)*cos(p),sin(p)],v,'sphere') sage: ss = var('ss') sage: vv_array = sphere.parallel_translation_numerical([ss,ss],ss,[1.0,1.0],[0.0,pi/4,20]) sage: vv_array[0] (0.000000000000000, [1.00000000000000, 1.00000000000000]) sage: vv_array[5] (0.19634954084936207, [0.98060186565184926, 1.0389928973994338]) """ u1 = self.variables[1] u2 = self.variables[2] solver = self._create_pt_ode_system(tuple(curve), t) t_interval, n = tinterval[0:2], tinterval[2] solver.y_0 = v0 solver.ode_solve(t_span=t_interval, num_points=n) return solver.solution
• ## new file sage/riemann/vector_functions.py

diff -r 120c07be6358 -r b834b692db4c sage/riemann/vector_functions.py
 - ### TODO: Some of these functions have Sage equivalents -- need to use these instead. from sage.modules.free_module_element import vector from sage.matrix.constructor import matrix # vector functions #def scalar_product(a,b): #    return(a[0]*b[0]+a[1]*b[1]+a[2]*b[2]) #def norm(a): #    return(sqrt(scalar_product(a,a))) def vector_product(a,b): return(vector([a[1]*b[2]-a[2]*b[1],a[2]*b[0]- a[0]*b[2],a[0]*b[1]-a[1]*b[0]])) def mixed_product(a,b,c): AA=matrix([[a[0],a[1],a[2]],[b[0],b[1],b[2]],[c[0],c[1],c[2]]]) return AA.determinant() # standard vector functions def vector_function_e(u): return vector([cos(u),sin(u),0]) def vector_function_g(u): return vector([-sin(u),cos(u),0]) vector_k=vector([0,0,1]) # simplify list of functions def simplify_list(a): return([i.simplify_full() for i in a]) def simplify_vector(a): return(vector([i.simplify_full() for i in a])) def simplify_matrix(A): return matrix([[ xx.simplify_full() for xx in row ] for row in A])
• ## setup.py

diff -r 120c07be6358 -r b834b692db4c setup.py