Ticket #10132: trac_10132_differential_geometry_sage.patch

File trac_10132_differential_geometry_sage.patch, 40.6 KB (added by jvkersch, 10 years ago)

Code of Mikhail as a Sage patch

  • sage/all.py

    # HG changeset patch
    # User Joris Vankerschaver <joris.vankerschaver@gmail.com>
    # Date 1287377049 25200
    # Node ID 930770c04a6d7628e9bcff5b0bd48a114c678903
    # Parent  7ae7b5175304fea5b08c9ee313c372a3936baa1e
    #10132: Differential geometry via Sage
    
    diff -r 7ae7b5175304 -r 930770c04a6d sage/all.py
    a b  
    139139from sage.ext.fast_eval      import fast_float
    140140
    141141from sage.tensor.all     import *
     142from sage.riemann.all    import *
    142143
    143144from copy import copy, deepcopy
    144145
  • new file sage/riemann/__init__.py

    diff -r 7ae7b5175304 -r 930770c04a6d sage/riemann/__init__.py
    - +  
     1# This comment is here so the file is non-empty (so Mercurial will check it in).
  • new file sage/riemann/all.py

    diff -r 7ae7b5175304 -r 930770c04a6d sage/riemann/all.py
    - +  
     1from parametrized_surface3d import ParametrizedSurface3D
  • new file sage/riemann/parametrized_surface3d.py

    diff -r 7ae7b5175304 -r 930770c04a6d sage/riemann/parametrized_surface3d.py
    - +  
     1"""
     2Basic invariants of parametrized surfaces.
     3
     4
     5AUTHORS:
     6        - Mikhail Malakhaltsev (2010-09-25): initial version
     7
     8"""   
     9#*****************************************************************************
     10#       Copyright (C) 2010  Mikhail Malakhaltsev <mikarm@gmail.com>
     11#
     12#  Distributed under the terms of the GNU General Public License (GPL)
     13#                  http://www.gnu.org/licenses/
     14#*****************************************************************************
     15
     16
     17# The following is maybe too specific ...
     18from vector_functions import *
     19
     20from sage.structure.sage_object import SageObject
     21from sage.calculus.functional import diff
     22from sage.functions.other import sqrt
     23
     24
     25# jvkersch: the following class can be expanded as needed (and should probably live somewhere else than in this file), but for now this is sufficient.
     26
     27class GenericManifold(SageObject):
     28    pass
     29
     30
     31class ParametrizedSurface3D(GenericManifold):
     32    """
     33    This is a class for finding basic invariants of surfaces.
     34
     35    USAGE:
     36    parametrized_surface3d(surface_equation,variables,name_of_surface)
     37    where
     38      surface_equation is a vector or list of three components
     39      variables is a list of two variables
     40      name_of_surface is a string (optional)
     41
     42    Warning: The orientation on the surface is given by the parametrization, that is the positive frame is
     43    $\partial_1 \\vec r$, $\partial_2 \\vec r$.
     44       
     45    This class includes the following methods:
     46      natural_frame
     47      normal_vector
     48      first_fundamental_form_coefficients
     49      first_fundamental_form
     50      first_fundamental_form_inverse_coefficients
     51      first_fundamental_form_inverse_coefficients
     52      area_form_squared
     53      area_form
     54      rotation
     55      orthonormal_frame
     56      lie_bracket
     57      frame_structure_functions
     58      second_order_natural_frame
     59      second_fundamental_form_coefficients
     60      second_fundamental_form
     61      gauss_curvature
     62      mean_curvature
     63      principal curvatures
     64      principal directions
     65      connection_coefficients
     66      geodesics_numerical
     67      parallel_translation_numerical
     68     
     69
     70    EXAMPLES::
     71
     72      sage: var('u,v')
     73      sage: paraboloid = parametrized_surface3d([u,v,u^2+v^2],[u,v],'paraboloid')
     74      sage: paraboloid = parametrized_surface3d([u,v,u^2+v^2],[u,v])
     75           
     76    """
     77
     78    def __init__(self,equation,variables,*name):
     79       
     80        # put here
     81       
     82        # attach('/Users/Joris/tmp/ticket_files/vector_functions.sage')
     83       
     84
     85        self.equation = equation
     86        self.variables_list = variables
     87        self.variables = {1:self.variables_list[0],2:self.variables_list[1]}
     88        self.name = name
     89       
     90       
     91        self.index_list = [(i,j) for i in (1,2) for j in (1,2)]
     92        self.index_list_3=[(i,j,k) for i in (1,2) for j in (1,2) for k in (1,2)]
     93
     94               
     95    def natural_frame(self,vector_number=0):
     96        """
     97        This functions find the natural frame of the parametrized surface
     98        INPUT:
     99        the argument can be empty, or equal to 1, or to 2
     100
     101        OUTPUT:
     102        if argument is equal to 1 or 2, the output is the corresponding vector of the natural frame
     103        if argument is empty, the output is the dictionary of the vector fields of the natural frame
     104
     105        EXAMPLES:: 
     106
     107            sage:  eparaboloid = parametrized_surface3d([u,v,u^2+v^2],[u,v],'elliptic paraboloid')
     108            sage:  eparaboloid.natural_frame()
     109            {1: (1, 0, 2*u), 2: (0, 1, 2*v)}
     110            sage:  eparaboloid.natural_frame(1)
     111            (1, 0, 2*u)
     112            sage:  eparaboloid.natural_frame(2)
     113            (0, 1, 2*v)
     114        """
     115
     116        if vector_number == 1:
     117           return vector([diff(f,self.variables[1]).simplify_full() for f in self.equation])
     118        elif vector_number == 2:
     119           return vector([diff(f,self.variables[2]).simplify_full() for f in self.equation])
     120        elif vector_number == 0:
     121           dr1 = vector([diff(f,self.variables[1]).simplify_full() for f in self.equation])
     122           dr2 = vector([diff(f,self.variables[2]).simplify_full() for f in self.equation])
     123           return  {1:dr1,2:dr2}
     124        else:
     125            print "You set an inappropriate argument. Read the doc on this function"
     126
     127    def normal_vector(self,length=0):
     128        """
     129        This functions find the normal vector of the parametrized surface
     130       
     131        INPUT:
     132        empty argument, of a nonzero real number l
     133       
     134        OUTPUT:
     135        without arguments gives the normal vector which is the vector product of the
     136        natural frame vectors;
     137
     138        with the argument l gives the normal vector of length l (the orientation determined
     139        by the sign of l)
     140       
     141        EXAMPLES::
     142
     143          sage: eparaboloid = parametrized_surface3d([u,v,u^2+v^2],[u,v],'elliptic paraboloid')
     144          sage: eparaboloid.normal_vector()
     145          (-2*u,-2*v,1)
     146          sage: eparaboloid.normal_vector(1)
     147          (-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))
     148
     149        """     
     150        dr = self.natural_frame()
     151        if length == 0:
     152            return simplify_vector(vector_product(dr[1],dr[2]))
     153        else:
     154            return simplify_vector(length*vector_product(dr[1],dr[2])/self.area_form())
     155
     156 
     157    def first_fundamental_form_coefficients(self,index1=0,index2=0):
     158        """
     159        This function gives the coefficients of the first fundamental form
     160       
     161        INPUT:
     162        empty argument, or one of the the following pair of indices: 1,1; 1,2; 2,1; 2,2.
     163
     164        OUTPUT:
     165        with empty argument the output is the dictionary of coefficients of the first fundamental form
     166
     167        with given indices the output is the corresponding coefficient of the first fundamental form,
     168        for example index1 = 1, index2 = 2 gives the coefficient $g_{12}$
     169
     170        EXAMPLES::
     171
     172            sage: var('u,v')
     173            sage: sphere = parametrized_surface3d([cos(u)*cos(v),sin(u)*cos(v),sin(v)],[u,v],'sphere')
     174            sage: sphere.first_fundamental_form_coefficients()
     175            {(1, 2): 0, (1, 1): cos(v)^2, (2, 1): 0, (2, 2): 1}
     176
     177            sage: sphere.first_fundamental_form_coefficients(1,1)
     178            cos(v)^2
     179
     180            sage: sphere.first_fundamental_form_coefficients(2,1)
     181            0
     182
     183            sage: sphere.first_fundamental_form_coefficients(3,2)
     184            'The argument is not appropriate. Read doc'
     185        """
     186        dr = self.natural_frame()
     187        indices=[index1,index2]
     188        if indices == [0,0]:
     189            g11 = (dr[1]*dr[1]).simplify_full()
     190            g12 = (dr[1]*dr[2]).simplify_full()
     191            g21 = g12
     192            g22 = (dr[2]*dr[2]).simplify_full()
     193            return {(1,1):g11,(1,2):g12,(2,1):g21,(2,2):g22}
     194        elif indices == [1,1]:
     195            return (dr[1]*dr[1]).simplify_full()
     196        elif indices == [1,2]:
     197            return (dr[1]*dr[2]).simplify_full()
     198        elif indices == [2,1]:
     199            return (dr[1]*dr[2]).simplify_full()
     200        elif indices == [2,2]:
     201            return (dr[2]*dr[2]).simplify_full()
     202        else:
     203            return 'The argument is not appropriate. Read doc'
     204
     205     
     206
     207    def first_fundamental_form(self,vector1,vector2):
     208        """
     209        Finds the value of first fundamental form on two vectors
     210
     211        INPUT:
     212        Two vectors $v=(v^1,v^2)$ and $w=(w^1,w^2)$
     213
     214        OUTPUT:
     215        $g_{11} v^1 w^1 + g_{12}(v^1 w^2 + v^2 w^1) + g_{22} v^2 w^2$
     216
     217        #EXAMPLES::
     218         
     219            sage: var('u,v')
     220            sage: var ('v1,v2,w1,w2')
     221            sage: sphere = parametrized_surface3d([cos(u)*cos(v),sin(u)*cos(v),sin(v)],[u,v],'sphere')
     222
     223            sage: sphere.first_fundamental_form(vector([v1,v2]),vector([w1,w2]))
     224            v1*w1*cos(v)^2 + v2*w2
     225
     226            sage: vv = vector([1,2])
     227            sage: sphere.first_fundamental_form(vv,vv)
     228            cos(v)^2 + 4
     229
     230            sage: sphere.first_fundamental_form([1,1],[2,1])
     231            2*cos(v)^2 + 1
     232        """
     233        gg = self.first_fundamental_form_coefficients()
     234        return sum(gg[ind]*vector1[ind[0]-1]*vector2[ind[1]-1] for ind in self.index_list)
     235
     236 
     237    def area_form_squared(self):
     238        """
     239        Gives $g_{11}g_{22} - g_{12}^2$, which is the square of the coefficient of the area form
     240       
     241        INPUT:
     242        No arguments
     243
     244        OUTPUT: 
     245        $g_{11}g_{22} - g_{12}^2$
     246
     247        EXAMPLES::
     248
     249            sage: var('u,v')
     250            sage: sphere = parametrized_surface3d([cos(u)*cos(v),sin(u)*cos(v),sin(v)],[u,v],'sphere')
     251            sage: sphere.area_form_squared()
     252            cos(v)^2
     253
     254        """
     255        gg = self.first_fundamental_form_coefficients()
     256        return  (gg[(1,1)]*gg[(2,2)]-gg[(1,2)]**2).simplify_full()
     257   
     258
     259    def area_form(self):
     260        """
     261        Gives $\sqrt{g_{11}g_{22} - g_{12}^2}$, which is the coefficient of the area form
     262       
     263        INPUT:
     264        No arguments
     265
     266        OUTPUT: 
     267        $\sqrt{g_{11}g_{22} - g_{12}^2}$
     268
     269        EXAMPLES::
     270
     271            sage: var('u,v')
     272            sage: sphere = parametrized_surface3d([cos(u)*cos(v),sin(u)*cos(v),sin(v)],[u,v],'sphere')
     273            sage: sphere.area_form_squared()
     274            abs(cos(v))
     275
     276        """
     277        return sqrt(self.area_form_squared()).simplify_full()
     278
     279
     280    def first_fundamental_form_inverse_coefficients(self):
     281        """
     282        Gives $g^{ij}$
     283       
     284        INPUT:
     285        No arguments
     286
     287        OUTPUT: 
     288        dictionary {(1,1):$g^{11}$, (1,2):$g^{12}$, (2,1):$g^{21}$, (2,2):$g^{22}$
     289
     290        EXAMPLES::
     291
     292            sage: var('u,v')
     293            sage: sphere = parametrized_surface3d([cos(u)*cos(v),sin(u)*cos(v),sin(v)],[u,v],'sphere')
     294            sage: sphere.first_fundamental_form_inverse_coefficients()
     295            {(1, 2): 0, (1, 1): cos(v)^(-2), (2, 1): 0, (2, 2): 1}
     296        """
     297
     298        gg = self.first_fundamental_form_coefficients()
     299        DD = gg[(1,1)]*gg[(2,2)]-gg[(1,2)]**2
     300        gi11 = (gg[(2,2)]/DD).simplify_full()
     301        gi12 = (-gg[(1,2)]/DD).simplify_full()
     302        gi21 = gi12
     303        gi22 = (gg[(1,1)]/DD).simplify_full()
     304        return {(1,1):gi11,(1,2):gi12,(2,1):gi21,(2,2):gi22}
     305
     306
     307    def rotation(self,theta):
     308        """
     309        Gives the matrix of the operator of rotation on the given angle $\\theta$ with respect to the natural frame
     310       
     311        INPUT:
     312        given angle $\\theta$
     313
     314        OUTPUT: 
     315        matrix of the operator of rotation on $\\theta$ with respect to the natural frame
     316
     317        ALGORITHM:
     318        The operator of rotation on $\pi/2$ is $J^i_j = g^{ik}\omega_{jk}$, where $\omega$ is the area form
     319        The operator of rotation on angle $\\theta$ is $\cos(\\theta) I + sin(\\theta) J$
     320
     321        EXAMPLES::
     322
     323            sage: var('u,v')
     324            sage: assume(cos(v)>0)
     325            sage: sphere = parametrized_surface3d([cos(u)*cos(v),sin(u)*cos(v),sin(v)],[u,v],'sphere')
     326            sage: rotation = sphere.rotation(pi/3)
     327            sage: rotation^3
     328            [-1  0]
     329            [ 0 -1]
     330            # it is true because the rotation at $\pi/3$ applied three times gives rotation a $\pi$
     331        """
     332
     333        from sage.functions.trig import sin, cos
     334
     335        gi = self.first_fundamental_form_inverse_coefficients()
     336        w12 = self.area_form()
     337        RR11 = (cos(theta) + sin(theta)*gi[1,2]*w12).simplify_full()
     338        RR12 = (- sin(theta)*gi[1,1]*w12).simplify_full()
     339        RR21 = (sin(theta)*gi[2,2]*w12).simplify_full()
     340        RR22 = (cos(theta) - sin(theta)*gi[2,1]*w12).simplify_full()       
     341        rotation = matrix([[RR11,RR12],[RR21,RR22]])
     342        return rotation
     343
     344
     345    def orthonormal_frame(self,vector_number=0,coordinates='ext'):
     346        """
     347        Gives the orthonormal frame field of the surface 
     348       
     349        INPUT:
     350        (vector_number,coordinates), where
     351        vector number is 0 (default), 1 or 2
     352        coordinates is 'ext'(default) or 'int'
     353               
     354
     355        OUTPUT:
     356        If vector_number is equal to 1 or 2, the output is the corresponding vector of an orthonormal frame.
     357        If argument vector_number is empty, the output is the dictionary of the vector fields of an orthonormal frame.
     358       
     359        If coordinates='ext', output is coordinates in $\mathbb{R}^3$
     360        If coordinates='int', output is coordinates with respect to the natural frame
     361
     362        ALGORITHM:
     363        We normalize the first vector $\\vec e_1$ of the natural frame and then get the second frame vector
     364        $\\vec e_2 = [\\vec n, \\vec e_1]$.
     365
     366        EXAMPLES::
     367
     368            sage: var('u,v')
     369            sage: assume(cos(v)>0)
     370            sage: sphere = parametrized_surface3d([cos(u)*cos(v),sin(u)*cos(v),sin(v)],[u,v],'sphere')
     371            sage: EE = sphere.orthonormal_frame()
     372            sage: (EE[1]*EE[1]).simplify_full()
     373            1
     374            sage: (EE[1]*EE[2]).simplify_full()
     375            0
     376            sage: simplify_vector(vector_product(EE[1],EE[2])-sphere.normal_vector(1))
     377            [0,0,0]
     378
     379            sage: EE_int = sphere.orthonormal_frame(coordinates='int')
     380            sage: EE_int
     381            {1: (1/cos(v), 0), 2: (0, 1)}
     382            sage: sphere.first_fundamental_form(EE_int[1],EE_int[1])
     383            1
     384            sage: sphere.first_fundamental_form(EE_int[1],EE_int[2])
     385            0
     386            sage: sphere.first_fundamental_form(EE_int[2],EE_int[2])
     387            1
     388           
     389        """
     390
     391        from sage.symbolic.constants import pi
     392
     393       
     394        if (vector_number,coordinates) not in [(0,'ext'),(0,'int'),(1,'ext'),(1,'int'),(2,'ext'),(2,'int')]:
     395           return "You set inappropriate arguments. Read the doc on this function"
     396
     397        if coordinates == 'ext':
     398            E1 = simplify_vector(self.natural_frame(1)/sqrt(self.first_fundamental_form_coefficients(1,1)))
     399            if vector_number == 0:
     400                E2 = simplify_vector(vector_product(self.normal_vector(1),E1))
     401                return  {1:E1, 2:E2}
     402            elif vector_number == 1:
     403                return E1
     404            elif vector_number == 2:
     405                return simplify_vector(vector_product(self.normal_vector(1),E1))
     406        else:
     407            E1 =  vector([(1/sqrt(self.first_fundamental_form_coefficients(1,1))).simplify_full(),0])
     408            if vector_number == 0:
     409               E2 = simplify_vector(self.rotation(pi/2)*E1)
     410               return  {1:E1, 2:E2}
     411            elif vector_number == 1:
     412               return E1
     413            elif vector_number == 2:
     414               E2 = simplify_vector(self.rotation(pi/2)*E1)
     415               return E2
     416       
     417    def lie_bracket(self,v,w):
     418        """
     419        Gives the Lie bracket of two vector fields which are given by coordinates with respect to the natural frame
     420         
     421       
     422        INPUT:
     423        v and w are vectors
     424               
     425
     426        OUTPUT:
     427        vector [v,w]
     428
     429
     430        EXAMPLES::
     431
     432            sage: var('u,v')
     433            sage: assume(cos(v)>0)
     434            sage: sphere = parametrized_surface3d([cos(u)*cos(v),sin(u)*cos(v),sin(v)],[u,v],'sphere')
     435            sage: sphere.lie_bracket([u,v],[-v,u])
     436            (0, 0)
     437
     438            sage: EE_int = sphere.orthonormal_frame(coordinates='int')
     439            sage: sphere.lie_bracket(EE_int[1],EE_int[2])
     440            (sin(v)/cos(v)^2, 0)
     441        """
     442        vv = vector([xx for xx in v])
     443        ww = vector([xx for xx in w])
     444        uuu = [self.variables[1],self.variables[2]]
     445        Dvv = matrix([ [ diff(xxx,uu).simplify_full() for uu in uuu ] for xxx in vv])   
     446        Dww = matrix([ [ diff(xxx,uu).simplify_full() for uu in uuu ] for xxx in ww])   
     447        return simplify_vector(vector(Dvv*ww - Dww*vv))
     448
     449
     450    def frame_structure_functions(self,e1,e2):
     451        """
     452        Gives the structure functions $c^k_{ij}$ for a frame field $e_1$, $e_2$, where
     453        $[e_i,e_j] = c^k_{ij}e_k$
     454         
     455       
     456        INPUT:
     457        Frame vectors $e_1$, $e_2$
     458               
     459
     460        OUTPUT:
     461        The dictionary of $c^k_{ij}$.
     462        Warning: the tuple (i,j,k) corresponds to $c^k_{ij}$
     463
     464
     465        EXAMPLES::
     466
     467            sage: var('u,v')
     468            sage: assume(cos(v)>0)
     469            sage: sphere = parametrized_surface3d([cos(u)*cos(v),sin(u)*cos(v),sin(v)],[u,v],'sphere')
     470            sage: sphere.frame_structure_functions([u,v],[-v,u])
     471            {(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}
     472            sage: structfun = sphere.frame_structure_functions(EE_int[1],EE_int[2])
     473            sage: structfun
     474            {(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}
     475            sage: sphere.lie_bracket(EE_int[1],EE_int[2]) - CC[(1,2,1)]*EE_int[1] - CC[(1,2,2)]*EE_int[2]
     476            (0, 0)
     477            """
     478        ee1 = vector([xx for xx in e1])
     479        ee2 = vector([xx for xx in e2])
     480       
     481        lb = simplify_vector(self.lie_bracket(ee1,ee2))
     482        trans_matrix = matrix([[ee1[0],ee2[0]],[ee1[1],ee2[1]]])
     483        zz = simplify_vector(trans_matrix.inverse()*lb)
     484        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}
     485
     486    def second_order_natural_frame(self,index1=0,index2=0):
     487        """
     488        Gives the second derivatives of the equation $\\vec r = \\vec r(u^1,u^2)$ of parametrized surface
     489         
     490       
     491        INPUT:
     492        empty argument, or one of the the following pair of indices: 1,1; 1,2; 2,1; 2,2.
     493               
     494
     495        OUTPUT:
     496        With empty argument the output is the dictionary of $\partial_{ij}\\vec r(u^1,u^2)$.
     497
     498        With given indices the output is the corresponding second derivative of the surface parametric equation,
     499        For example, index1 = 1, index2 = 2 gives $\partial_{12}\\vec r(u^1,u^2)$.
     500
     501
     502        EXAMPLES::
     503
     504            sage: var('u,v')
     505            sage: sphere = parametrized_surface3d([cos(u)*cos(v),sin(u)*cos(v),sin(v)],[u,v],'sphere')
     506            sage: sphere.second_order_natural_frame()
     507            {(1, 2): (sin(u)*sin(v), -sin(v)*cos(u), 0), (1, 1): (-cos(u)*cos(v),
     508            -sin(u)*cos(v), 0), (2, 1): (sin(u)*sin(v), -sin(v)*cos(u), 0), (2, 2):
     509            (-cos(u)*cos(v), -sin(u)*cos(v), -sin(v))}
     510           
     511            sage: sphere.second_order_natural_frame(1,1)
     512            (-cos(u)*cos(v), -sin(u)*cos(v), 0)
     513            sage: sphere.second_order_natural_frame(1,2)
     514            (sin(u)*sin(v), -sin(v)*cos(u), 0)
     515            sage: sphere.second_order_natural_frame(2,2)
     516            (-cos(u)*cos(v), -sin(u)*cos(v), -sin(v)
     517            """
     518        u = self.variables[1]
     519        v = self.variables[2]
     520        indices=[index1,index2]
     521       
     522        if indices == [0,0]:
     523            dr11 = vector([diff(f,u,2).simplify_full() for f in self.equation])
     524            dr12 = vector([diff(f,u,v).simplify_full() for f in self.equation])
     525            dr22 = vector([diff(f,v,2).simplify_full() for f in self.equation])
     526            return {(1,1):dr11,(1,2):dr12,(2,1):dr12,(2,2):dr22}
     527
     528        elif indices == [1,1]:
     529            return vector([diff(f,u,2).simplify_full() for f in self.equation])
     530        elif indices == [1,2]:
     531            return vector([diff(f,u,v).simplify_full() for f in self.equation])
     532        elif indices == [2,1]:
     533            return vector([diff(f,u,v).simplify_full() for f in self.equation])
     534        elif indices == [2,2]:
     535            return vector([diff(f,v,2).simplify_full() for f in self.equation])
     536        else:
     537            return 'The argument is not appropriate. Read doc'
     538       
     539       
     540    def second_fundamental_form_coefficients(self,index1=0,index2=0):
     541        """
     542        This function gives the coefficients $h_{ij}$ of the second fundamental form.
     543       
     544        INPUT:
     545        empty argument, or one of the the following pair of indices: 1,1; 1,2; 2,1; 2,2.
     546
     547        OUTPUT:
     548        with empty argument the output is the dictionary of coefficients of the second fundamental form
     549
     550        with given indices the output is the corresponding coefficient of the second fundamental form,
     551        for example index1 = 1, index2 = 2 gives the coefficient $h_{12}$
     552
     553        EXAMPLES::
     554
     555            sage: var('u,v')
     556            sage: assume(cos(v)>0)
     557            sage: sphere = parametrized_surface3d([cos(u)*cos(v),sin(u)*cos(v),sin(v)],[u,v],'sphere')
     558
     559            sage: sphere.second_fundamental_form_coefficients()
     560            {(1, 2): 0, (1, 1): -cos(v)^2, (2, 1): 0, (2, 2): -1}
     561            sage: sphere.second_fundamental_form_coefficients(1,1)
     562            -cos(v)^2
     563            sage: sphere.second_fundamental_form_coefficients(2,1)
     564            0
     565
     566            sage: sphere.second_fundamental_form_coefficients(3,2)
     567            'The argument is not appropriate. Read doc'
     568        """
     569        NN = self.normal_vector(1)
     570        indices=[index1,index2]
     571
     572        if indices == [0,0]:
     573            ddr = self.second_order_natural_frame()
     574            h11 = (ddr[1,1]*NN).simplify_full()
     575            h12 = (ddr[1,2]*NN).simplify_full()
     576            h21 = h12
     577            h22 = (ddr[2,2]*NN).simplify_full()
     578            return {(1,1):h11,(1,2):h12,(2,1):h21,(2,2):h22}
     579        elif indices == [1,1]:
     580            return (self.second_order_natural_frame(1,1)*NN).simplify_full()
     581        elif indices == [1,2]:
     582            return (self.second_order_natural_frame(1,2)*NN).simplify_full()
     583        elif indices == [2,1]:
     584            return (self.second_order_natural_frame(2,1)*NN).simplify_full()
     585        elif indices == [2,2]:
     586            return (self.second_order_natural_frame(2,2)*NN).simplify_full()
     587        else:
     588            return 'The argument is not appropriate. Read doc'
     589
     590    def second_fundamental_form(self,vector1,vector2):
     591        """
     592        Finds the value of second fundamental form on two vectors
     593
     594        INPUT:
     595        Two vectors $v=(v^1,v^2)$ and $w=(w^1,w^2)$
     596
     597        OUTPUT:
     598        $h_{11} v^1 w^1 + h_{12}(v^1 w^2 + v^2 w^1) + h_{22} v^2 w^2$
     599
     600        EXAMPLES::
     601
     602           sage: var('u,v')
     603           sage: var ('v1,v2,w1,w2')
     604           sage: assume(cos(v) > 0)
     605           sage: sphere = parametrized_surface3d([cos(u)*cos(v),sin(u)*cos(v),sin(v)],[u,v],'sphere')
     606           sage: sphere.second_fundamental_form(vector([v1,v2]),vector([w1,w2]))
     607           -v1*w1*cos(v)^2 - v2*w2
     608           sage: vv = vector([1,2])
     609           sage: sphere.second_fundamental_form(vv,vv)
     610           -cos(v)^2 - 4
     611           sage: sphere.second_fundamental_form([1,1],[2,1])
     612           -2*cos(v)^2 - 1
     613                   
     614        """
     615        hh = self.second_fundamental_form_coefficients()
     616        return sum(hh[ind]*vector1[ind[0]-1]*vector2[ind[1]-1] for ind in self.index_list)
     617
     618    def gauss_curvature(self):
     619        """
     620        Finds the gaussian curvature $K = \\frac{h_{11}h_{22} - h_{12}^2}{g_{11}g_{22} - g_{12}^2}$.
     621
     622        INPUT:
     623        No arguments
     624
     625        OUTPUT:
     626        $K = \\frac{h_{11}h_{22} - h_{12}^2}{g_{11}g_{22} - g_{12}^2}$
     627
     628        EXAMPLES::
     629
     630           sage: var('R')
     631           sage: assume(R>0)
     632           sage: var('u,v')
     633           sage: assume(cos(v)>0)
     634           sage: sphere = parametrized_surface3d([R*cos(u)*cos(v),R*sin(u)*cos(v),R*sin(v)],[u,v],'sphere')
     635           sage: sphere.gauss_curvature()
     636           R^(-2)
     637                   
     638        """
     639        hh = self.second_fundamental_form_coefficients()
     640        return ((hh[(1,1)]*hh[(2,2)]-hh[(1,2)]**2)/self.area_form_squared()).simplify_full()
     641
     642    def mean_curvature(self):
     643        """
     644        Finds the mean curvature $H = \\frac{1}{2}\\frac{g_{22}h_{11} - 2g_{12}h_{12} + g_{11}h_{22}}{g_{11}g_{22} - g_{12}^2}$.
     645
     646        INPUT:
     647        No arguments
     648
     649        OUTPUT:
     650        $H = \\frac{1}{2}\\frac{g_{22}h_{11} - 2g_{12}h_{12} + g_{11}h_{22}}{g_{11}g_{22} - g_{12}^2}$.
     651
     652        EXAMPLES::
     653           
     654           sage: var('R')
     655           sage: assume(R>0)
     656           sage: var('u,v')
     657           sage: assume(cos(v)>0)
     658           sage: sphere = parametrized_surface3d([R*cos(u)*cos(v),R*sin(u)*cos(v),R*sin(v)],[u,v],'sphere')
     659           sage: sphere.mean_curvature()
     660           -1/R
     661
     662        """
     663        gg = self.first_fundamental_form_coefficients()
     664        hh = self.second_fundamental_form_coefficients()
     665        denom = 2*self.area_form_squared()
     666        enum =(gg[(2,2)]*hh[(1,1)]-2*gg[(1,2)]*hh[(1,2)]+gg[(1,1)]*hh[(2,2)]).simplify_full()
     667        return (enum/denom).simplify_full()
     668       
     669
     670    def principal_curvatures(self):
     671        """
     672        Finds the principal curvatures of the surface
     673
     674        INPUT:
     675        No arguments
     676
     677        OUTPUT:
     678        The dictionary of principal curvatures         
     679
     680        EXAMPLES::
     681                     
     682           sage: var('R')
     683           sage: assume(R>0)
     684           sage: var('u,v')
     685           sage: assume(cos(v)>0)
     686           sage: sphere = parametrized_surface3d([R*cos(u)*cos(v),R*sin(u)*cos(v),R*sin(v)],[u,v],'sphere')
     687           sage: sphere.principal_curvatures()
     688           {1: -1/R, 2: -1/R}
     689
     690           sage: var('u,v')
     691           sage: var('R,r')
     692           sage: assume(R>r,r>0)
     693           sage: torus = parametrized_surface3d([(R+r*cos(v))*cos(u),(R+r*cos(v))*sin(u),r*sin(v)],[u,v],'torus')
     694           sage: torus.principal_curvatures()
     695           {1: -cos(v)/(r*cos(v) + R), 2: -1/r}
     696
     697        """
     698
     699        from sage.symbolic.assumptions import assume
     700        from sage.symbolic.relation import solve
     701        from sage.calculus.var import var
     702
     703        KK = self.gauss_curvature()
     704        HH = self.mean_curvature()
     705
     706        # jvkersch: when this assumption is uncommented, Sage raises an error stating that the assumption
     707        # is redundant... Can we safely omit this, based on some geometric reasoning?
     708       
     709        # assume(HH**2-KK>=0)
     710
     711        x = var('x')
     712        sol = solve(x**2 -2*HH*x + KK == 0,x)
     713       
     714        #k1=var('k1')
     715        #k2=var('k2')
     716
     717        # 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.  Hence the modification below...
     718
     719        k1 = (x.substitute(sol[0])).simplify_full()
     720        if len(sol) == 1:
     721            k2 = k1
     722        else:
     723            k2 = (x.substitute(sol[1])).simplify_full()
     724       
     725        #k1 = (x.substitute(sol[0])).simplify_full()
     726        #k2 = (x.substitute(sol[1])).simplify_full()
     727
     728        # Why return a dictionary like this?  Is there a tensor whose components are the principal curvatures?  Maybe it is better to return a list [k1, k2]
     729
     730        return {1:k1,2: k2}
     731   
     732    def principal_directions(self):
     733        """
     734        Finds the principal curvatures and principal directions of the surface
     735
     736        INPUT:
     737        No arguments
     738
     739        OUTPUT:
     740        The dictionary of lists [a principal curvature, the corresponding principal direction]         
     741
     742        If principal curvatures coincide, gives the warning that the surface is a sphere.
     743
     744        EXAMPLES::
     745
     746           sage: var('u,v')
     747           sage: var('R,r')
     748           sage: assume(R>r,r>0)
     749           sage: torus = parametrized_surface3d([(R+r*cos(v))*cos(u),(R+r*cos(v))*sin(u),r*sin(v)],[u,v],'torus')
     750           sage: pdd = torus.principal_directions()
     751           sage: pdd[1]
     752           [-cos(v)/(r*cos(v) + R), (1, 0)]
     753           sage: pdd[2]
     754           [-1/r, (0, -(R*r*cos(v) + R^2)/r)]
     755
     756           sage: var('RR')
     757           sage: assume(RR>0)
     758           sage: var('u,v')
     759           sage: assume(cos(v)>0)
     760           sage: sphere = parametrized_surface3d([RR*cos(u)*cos(v),RR*sin(u)*cos(v),RR*sin(v)],[u,v],'sphere')
     761           sage: sphere.principal_directions()
     762           'This is a sphere, so any direction is principal'
     763
     764           sage: var('aa')
     765           sage: assume(aa>0)
     766           sage: catenoid = parametrized_surface3d([aa*cosh(v)*cos(u),aa*cosh(v)*sin(u),v],[u,v],'catenoid')
     767           sage: pd = catenoid.principal_directions()
     768           sage: pd[1][1]
     769           (0, 1/2*(2*aa^3*sinh(v)^2*cosh(v) + sqrt(4*aa^4*sinh(v)^2*cosh(v)^2 + aa^4 + 4*aa^2*cosh(v)^2 - 2*aa^2 + 1)*aa*cosh(v) + (aa^3 + aa)*cosh(v))/(aa^2*sinh(v)^2 + 1)^(3/2))
     770           sage: pd[2][1]
     771           (1, 0)
     772           sage: pd[1][1]*pd[2][1]
     773           0
     774        """
     775        gg = self.first_fundamental_form_coefficients()
     776        hh = self.second_fundamental_form_coefficients()
     777        kk = self.principal_curvatures()
     778        if kk[1]==kk[2]:
     779            return "This is a sphere, so any direction is principal"
     780        pd1 = simplify_vector([hh[(1,2)]-kk[1]*gg[(1,2)],-hh[(1,1)]+kk[1]*gg[(1,1)]])
     781        if pd1==vector([0,0]):
     782           pd1 = vector([1,0])
     783        pd2 = simplify_vector([hh[(1,2)]-kk[2]*gg[(1,2)],-hh[(1,1)]+kk[2]*gg[(1,1)]])
     784        if pd2==vector([0,0]):
     785           pd2 = vector([1,0])
     786        return {1:[kk[1],pd1],2:[kk[2],pd2]}
     787       
     788    def connection_coefficients(self):
     789        """
     790        Finds the connection coefficients of the surface
     791
     792        INPUT:
     793        No arguments
     794
     795        OUTPUT:
     796        The dictionary of connection coefficients.
     797
     798        Warning: the triple $(i,j,k)$ corresponds to $\Gamma^k_{ij}$. 
     799
     800        EXAMPLES::
     801                                 
     802           sage: var('r')
     803           sage: assume(r > 0)
     804           sage: var('u,v')
     805           sage: assume(cos(v)>0)
     806           sage: sphere = parametrized_surface3d([r*cos(u)*cos(v),r*sin(u)*cos(v),r*sin(v)],[u,v],'sphere')
     807           sage: sphere.connection_coefficients()
     808           {(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}
     809           
     810        """
     811        x = self.variables
     812        gg = self.first_fundamental_form_coefficients()
     813        gi = self.first_fundamental_form_inverse_coefficients()
     814
     815        dg = {}
     816        for kkk in self.index_list_3:
     817            dg[kkk]=gg[(kkk[1],kkk[2])].differentiate(x[kkk[0]]).simplify_full()
     818        structfun={}
     819
     820        for kkk in self.index_list_3:
     821            structfun[kkk]=sum(gi[(kkk[2],s)]*(dg[(kkk[0],kkk[1],s)]+dg[(kkk[1],kkk[0],s)]-dg[(s,kkk[0],kkk[1])])/2 for s in (1,2)).full_simplify()
     822        return structfun
     823
     824
     825    def geodesics_numerical(self,p0,v0,tinterval):
     826        """
     827        This method gives the numerical solution for the geodesic equations 
     828
     829        INPUT:
     830        p0 is the list of the coordinates of the initial point
     831       
     832        v0 is the  list of the coordinates of the initial vector
     833
     834        tinterval is the list [a,b,M], where (a,b) is the domain of the geodesic, M is the number of division points
     835
     836        OUTPUT:
     837        The list consisting of lists [t, [u1(t),u2(t)], [v1(t),v2(t)], [x1(t),x2(t),x3(t)]], where
     838
     839        t is a subdivision point
     840
     841        [u1(t),u2(t)] is the list of coordinates of the geodesic point
     842         
     843        [v1(t),v2(t)] is the list of coordinates of the vector tanget to the geodesic 
     844
     845        [x1(t),x2(t),x3(t)] is the list of coordinates of the geodesic point in the three-dimensional space
     846
     847
     848        EXAMPLES::
     849
     850           sage: var('p,q')
     851           sage: v = [p,q]
     852           sage: assume(cos(q)>0)
     853           sage: sphere = parametrized_surface3d([cos(q)*cos(p),sin(q)*cos(p),sin(p)],v,'sphere')
     854           sage: gg_array = sphere.geodesics_numerical([0,0],[1,1],[0,2*pi,5])
     855           sage: gg_array[0]
     856           [0.0, [0.0, 0.0], [1.0, 1.0], (1, 0, 0)]
     857           sage: gg_array[1]
     858           [1.2566370614359172, [0.76440104189216407, 1.8586223516062499], [-0.2838683571264714, 1.9194187087799912], (-0.204895333443, 0.692104654602, 0.692104796553)]
     859                                           
     860         
     861        """
     862
     863        # jvkersch: rewrote a few things in this routine to make it pure python, i.e. without the sage extensions.
     864
     865        #from sage.calculus.var import var
     866
     867        u1 = self.variables[1]
     868        u2 = self.variables[2]
     869        #uu1, uu2 = var('uu1,uu2')
     870
     871        # jvkersch: hack for now
     872        def f(ux, vx):
     873            return [self.equation[0].subs({u1: ux, u2: vx}),
     874                    self.equation[1].subs({u1: ux, u2: vx}),
     875                    self.equation[2].subs({u1: ux, u2: vx})]
     876
     877        C111 = self.connection_coefficients()[(1,1,1)]
     878        C121 = self.connection_coefficients()[(1,2,1)]
     879        C221 = self.connection_coefficients()[(2,2,1)]
     880        C112 = self.connection_coefficients()[(1,1,2)]
     881        C122 = self.connection_coefficients()[(1,2,2)]
     882        C222 = self.connection_coefficients()[(2,2,2)]
     883
     884        def ode_system(unknown_functions,t):
     885            uu1,uu2,v1,v2 = unknown_functions
     886
     887            # jvkersch: this is just a temporary hack to evaluate the connection coefficients at the relevant point.  I don't think this is the best way to do things, or the fastest.
     888            c111 = C111.subs({u1: uu1, u2: uu2})
     889            c121 = C121.subs({u1: uu1, u2: uu2})
     890            c221 = C221.subs({u1: uu1, u2: uu2})
     891            c112 = C112.subs({u1: uu1, u2: uu2})
     892            c122 = C122.subs({u1: uu1, u2: uu2})
     893            c222 = C222.subs({u1: uu1, u2: uu2})
     894
     895            #print type(C221)
     896            #print u1, u2, uu1, uu2
     897            #print c111, c121, c221, c112, c122, c222
     898
     899            du1 = v1
     900            du2 = v2
     901            dv1 = - c111*v1**2 - 2*c121*v1*v2 - c221*v2**2
     902            dv2 = - c112*v1**2 - 2*c122*v1*v2 - c222*v2**2
     903            return float(du1), float(du2), float(dv1), float(dv2)
     904
     905        step = float((tinterval[1]-tinterval[0])/tinterval[2])
     906        ttt =  float(tinterval[0])
     907        tarray = [ ttt ]
     908        for counter in range(1,tinterval[2]+1):
     909           ttt = ttt + step
     910           tarray = tarray + [ ttt ]
     911
     912        # import ode solving routine
     913        import scipy.integrate
     914
     915        initial_data = (p0[0],p0[1],v0[0],v0[1])
     916        solution = scipy.integrate.odeint(ode_system,initial_data,tarray)
     917        return [[ tarray[counter],[solution[counter,0], solution[counter,1]], [solution[counter,2], solution[counter,3]], f(solution[counter,0],solution[counter,1])]   for counter in range(0,len(tarray))]
     918
     919
     920    def parallel_translation_numerical(self,curve,t,v0,tinterval):
     921        """
     922        This method gives the numerical solution to the equation of parallel translation of a vector
     923
     924        INPUT:
     925        curve equation = list of functions which determine the curve wrt the local coordinate
     926
     927        t - curve parameter
     928
     929        v0 - initial vector
     930
     931        tinterval = [a,b,N], (a,b) is the domain of the curve, N is the number of subdivision points
     932
     933        OUTPUT:
     934        The list consisting of lists [t, [v1(t),v2(t)]], where
     935
     936        t is a subdivision point
     937
     938        [v1(t),v2(t)] is the list of coordinates of the vector translated parallely along the curve
     939
     940        EXAMPLES::
     941
     942           sage: var('p,q')
     943           sage: v = [p,q]
     944           sage: assume(cos(q)>0)
     945           sage: sphere = parametrized_surface3d([cos(q)*cos(p),sin(q)*cos(p),sin(p)],v,'sphere')
     946           sage: var('ss')
     947           sage: vv_array = sphere.parallel_translation_numerical([ss,ss],ss,[1,1],[0,pi/4,20])
     948           sage: vv_array[0]
     949           [0.0, [1.0, 1.0]]
     950           sage: vv_array[5]
     951           [0.19634954084936207, [0.98060182522638917, 1.0389930474425113]]
     952
     953        """
     954        u1 = self.variables[1]
     955        u2 = self.variables[2]
     956        #var('uu1,uu2')
     957       
     958        #f(u1,u2) = self.equation
     959
     960        def f(ux, vx):
     961            return [self.equation[0].subs({u1: ux, u2: vx}),
     962                    self.equation[1].subs({u1: ux, u2: vx}),
     963                    self.equation[2].subs({u1: ux, u2: vx})]
     964
     965        #from sage.calculus.var import var
     966        #t = var('t')
     967        tt = t
     968
     969        C111 = self.connection_coefficients()[(1,1,1)]
     970        C121 = self.connection_coefficients()[(1,2,1)]
     971        C221 = self.connection_coefficients()[(2,2,1)]
     972        C112 = self.connection_coefficients()[(1,1,2)]
     973        C122 = self.connection_coefficients()[(1,2,2)]
     974        C222 = self.connection_coefficients()[(2,2,2)]
     975
     976        du1=diff(curve[0],tt)
     977        du2=diff(curve[1],tt)
     978        uu1=curve[0]
     979        uu2=curve[1]
     980
     981
     982        def ode_system(unknown_functions,t1):
     983            du1p = du1.subs(tt==t1)
     984            du2p = du2.subs(tt==t1)
     985            uu1p = uu1.subs(tt==t1)
     986            uu2p = uu2.subs(tt==t1)
     987
     988            c111 = C111.subs({u1: uu1p, u2: uu2p})
     989            c121 = C121.subs({u1: uu1p, u2: uu2p})
     990            c221 = C221.subs({u1: uu1p, u2: uu2p})
     991            c112 = C112.subs({u1: uu1p, u2: uu2p})
     992            c122 = C122.subs({u1: uu1p, u2: uu2p})
     993            c222 = C222.subs({u1: uu1p, u2: uu2p})
     994
     995            #print c111, c121, c221, c112, c122, c222
     996           
     997            v1,v2 = unknown_functions
     998            f1 = - c111*du1p*v1 - c121*(du1p*v2 + du2p*v1) - c221*du2p*v2
     999            f2 = - c112*du1p*v1 - c122*(du1p*v2 + du2p*v1) - c222*du2p*v2
     1000            dv1 = f1.subs(tt==t1)
     1001            dv2 = f2.subs(tt==t1)
     1002
     1003            #print dv1, dv2
     1004           
     1005            return float(dv1), float(dv2)
     1006
     1007        step = float((tinterval[1]-tinterval[0])/tinterval[2])
     1008        ttt =  float(tinterval[0])
     1009        tarray = [ ttt ]
     1010        for counter in range(1,tinterval[2]+1):
     1011           ttt = ttt + step
     1012           tarray = tarray + [ ttt ]
     1013
     1014
     1015        # import ode solving routine
     1016        import scipy.integrate
     1017
     1018        initial_data = (v0[0],v0[1])
     1019        solution = scipy.integrate.odeint(ode_system,initial_data,tarray)
     1020        return [[ tarray[counter],[solution[counter,0], solution[counter,1]]]   for counter in range(0,len(tarray))]
     1021
  • new file sage/riemann/vector_functions.py

    diff -r 7ae7b5175304 -r 930770c04a6d sage/riemann/vector_functions.py
    - +  
     1### TODO: Some of these functions have Sage equivalents -- need to use these instead.
     2
     3from sage.modules.free_module_element import vector
     4from sage.matrix.constructor import matrix
     5
     6# vector functions
     7def scalar_product(a,b):
     8    return(a[0]*b[0]+a[1]*b[1]+a[2]*b[2])
     9def norm(a):
     10    return(sqrt(scalar_product(a,a)))
     11def vector_product(a,b):
     12  return(vector([a[1]*b[2]-a[2]*b[1],a[2]*b[0]-
     13  a[0]*b[2],a[0]*b[1]-a[1]*b[0]]))
     14def mixed_product(a,b,c):
     15    AA=matrix([[a[0],a[1],a[2]],[b[0],b[1],b[2]],[c[0],c[1],c[2]]])
     16    return AA.determinant()
     17
     18# standard vector functions
     19def vector_function_e(u):
     20    return vector([cos(u),sin(u),0])
     21def vector_function_g(u):
     22    return vector([-sin(u),cos(u),0])
     23vector_k=vector([0,0,1])
     24
     25   
     26# simplify list of functions
     27def simplify_list(a):
     28    return([i.simplify_full() for i in a])
     29   
     30def simplify_vector(a):
     31    return(vector([i.simplify_full() for i in a]))
     32
     33def simplify_matrix(A):
     34    return matrix([[ xx.simplify_full() for xx in row ] for row in A])
     35
     36
  • setup.py

    diff -r 7ae7b5175304 -r 930770c04a6d setup.py
    a b  
    927927                     'sage.quadratic_forms',
    928928                     'sage.quadratic_forms.genera',
    929929
     930                     'sage.riemann',
     931
    930932                     'sage.rings',
    931933                     'sage.rings.finite_rings',
    932934                     'sage.rings.number_field',