Ticket #10132: trac_10132_differential_geometry_sage_v2.patch

File trac_10132_differential_geometry_sage_v2.patch, 51.9 KB (added by jvkersch, 10 years ago)

New version of the patch

  • sage/all.py

    # HG changeset patch
    # User Joris Vankerschaver <joris.vankerschaver@gmail.com>
    # Date 1288567720 25200
    # Node ID 38ffdf19521972372901d91e130fbe3735a5ef0b
    # Parent  f667e86a25fbe7138d4933ecf9b972ba8f40087d
    #10132: Differential Geometry via Sage
    
    diff -r f667e86a25fb -r 38ffdf195219 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 f667e86a25fb -r 38ffdf195219 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 f667e86a25fb -r 38ffdf195219 sage/riemann/all.py
    - +  
     1from parametrized_surface3d import ParametrizedSurface3D
  • new file sage/riemann/parametrized_surface3d.py

    diff -r f667e86a25fb -r 38ffdf195219 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
     24from sage.misc.cachefunc import cached_method
     25
     26
     27
     28# 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.
     29
     30class GenericManifold(SageObject):
     31    pass
     32
     33
     34class ParametrizedSurface3D(GenericManifold):
     35    """
     36    This is a class for finding basic invariants of surfaces.
     37
     38    USAGE:
     39    parametrized_surface3d(surface_equation,variables,name_of_surface)
     40    where
     41      surface_equation is a vector or list of three components
     42      variables is a list of two variables
     43      name_of_surface is a string (optional)
     44
     45    Warning: The orientation on the surface is given by the parametrization, that is the positive frame is
     46    $\partial_1 \\vec r$, $\partial_2 \\vec r$.
     47       
     48    This class includes the following methods:
     49      natural_frame
     50      normal_vector
     51      first_fundamental_form_coefficients
     52      first_fundamental_form
     53      first_fundamental_form_inverse_coefficients
     54      first_fundamental_form_inverse_coefficients
     55      area_form_squared
     56      area_form
     57      rotation
     58      orthonormal_frame
     59      lie_bracket
     60      frame_structure_functions
     61      second_order_natural_frame
     62      second_fundamental_form_coefficients
     63      second_fundamental_form
     64      gauss_curvature
     65      mean_curvature
     66      principal curvatures
     67      principal directions
     68      connection_coefficients
     69      geodesics_numerical
     70      parallel_translation_numerical
     71     
     72
     73    EXAMPLES::
     74
     75      sage: var('u,v')
     76      sage: paraboloid = parametrized_surface3d([u,v,u^2+v^2],[u,v],'paraboloid')
     77      sage: paraboloid = parametrized_surface3d([u,v,u^2+v^2],[u,v])
     78           
     79    """
     80
     81    def __init__(self,equation,variables,*name):
     82
     83        self.equation = equation
     84        self.variables_list = variables
     85        self.variables = {1:self.variables_list[0],2:self.variables_list[1]}
     86        self.name = name
     87
     88        # jvkersch:
     89
     90        # Since we don't have the Sage extensions at this point,
     91        # self.equation is just a vector of symbolic expressions
     92        # and is not automatically callable.  Therefore, we also
     93        # define a callable version of self.equation.
     94
     95        def eq_callable(u, v):
     96            u1, u2 = self.variables_list
     97            return [fun.subs({u1: u, u2: v}) for fun in self.equation]
     98           
     99        self.eq_callable = eq_callable
     100
     101        # Various index tuples       
     102        self.index_list = [(i,j) for i in (1,2) for j in (1,2)]
     103        self.index_list_3=[(i,j,k) for i in (1,2) for j in (1,2) for k in (1,2)]
     104
     105
     106    @cached_method           
     107    def natural_frame(self,vector_number=0):
     108        """
     109        This functions find the natural frame of the parametrized surface
     110        INPUT:
     111        the argument can be empty, or equal to 1, or to 2
     112
     113        OUTPUT:
     114        if argument is equal to 1 or 2, the output is the corresponding vector of the natural frame
     115        if argument is empty, the output is the dictionary of the vector fields of the natural frame
     116
     117        EXAMPLES:: 
     118
     119            sage:  eparaboloid = parametrized_surface3d([u,v,u^2+v^2],[u,v],'elliptic paraboloid')
     120            sage:  eparaboloid.natural_frame()
     121            {1: (1, 0, 2*u), 2: (0, 1, 2*v)}
     122            sage:  eparaboloid.natural_frame(1)
     123            (1, 0, 2*u)
     124            sage:  eparaboloid.natural_frame(2)
     125            (0, 1, 2*v)
     126        """
     127
     128        # jvkersch: I would suggest just returning a list with the basis vectors of the natural
     129        # frame, rather than distinguishing between the cases where an individual vector should
     130        # be returned.  The reasons are as follows:
     131        #
     132        #  1.  If the member function is called `natural_frame()' it makes sense for it to return
     133        #      a frame, and not a vector
     134        #  2.  By returning the frame, we don't have to check for the index being out of bounds
     135        #  3.  This behavior is more in line with the `basis()' member function of Sage vector
     136        #      spaces, which also returns a basis as a list of vectors.
     137        #
     138
     139        dr1 = vector([diff(f,self.variables[1]).simplify_full() for f in self.equation])
     140        dr2 = vector([diff(f,self.variables[2]).simplify_full() for f in self.equation])
     141
     142        # jvkersch: it is necessary to return a dictionary here?  Couldn't we just return
     143        # a list [dr1, dr2]?
     144       
     145        return {1:dr1, 2:dr2}
     146
     147
     148        # Old code:
     149
     150        #if vector_number == 1:
     151        #   return vector([diff(f,self.variables[1]).simplify_full() for f in self.equation])
     152        #elif vector_number == 2:
     153        #   return vector([diff(f,self.variables[2]).simplify_full() for f in self.equation])
     154        #elif vector_number == 0:
     155        #   dr1 = vector([diff(f,self.variables[1]).simplify_full() for f in self.equation])
     156        #   dr2 = vector([diff(f,self.variables[2]).simplify_full() for f in self.equation])
     157        #   return  {1:dr1,2:dr2}
     158        #else:
     159        #    print "You set an inappropriate argument. Read the doc on this function"
     160
     161    @cached_method
     162    def normal_vector(self, normalized=False):
     163        """
     164        This functions find the normal vector of the parametrized surface
     165       
     166        INPUT:
     167        empty argument, of a nonzero real number l
     168       
     169        OUTPUT:
     170        without arguments gives the normal vector which is the vector product of the
     171        natural frame vectors;
     172
     173        with the argument l gives the normal vector of length l (the orientation determined
     174        by the sign of l)
     175       
     176        EXAMPLES::
     177
     178          sage: eparaboloid = parametrized_surface3d([u,v,u^2+v^2],[u,v],'elliptic paraboloid')
     179          sage: eparaboloid.normal_vector()
     180          (-2*u,-2*v,1)
     181          sage: eparaboloid.normal_vector(1)
     182          (-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))
     183
     184        """
     185
     186        # jvkersch: I have changed the default argument in this function, so that this function
     187        # takes a boolean normalized=True/False,  since I guess that the
     188        # `length' argument, which was previously there, will mostly be used to normalize the
     189        # vector (the docstrings also reflect this).  Please revert if this is a bad assumption.
     190
     191        # I've also rewritten the code below somewhat to reflect this change and to use Sage
     192        # standard functions (cross_product, norm, etc)
     193       
     194        dr = self.natural_frame()
     195        normal = dr[1].cross_product(dr[2])
     196
     197        if normalized:
     198            normal /= normal.norm()
     199        return simplify_vector(normal)
     200       
     201        #if length == 0:
     202        #    return simplify_vector(vector_product(dr[1],dr[2]))
     203        #else:
     204        #    return simplify_vector(length*vector_product(dr[1],dr[2])/self.area_form())
     205
     206
     207
     208    # jvkersch: this is the pattern that I have used for most of the functions that
     209    # previously return either a single component or a list of components.
     210    #
     211    # I've defined a cached (private) method which does not check its arguments and
     212    # computes one component at a time.  Note that caching is just a matter of
     213    # prepending "@cached_method" to the method signature.  For this to work, the
     214    # arguments of the method must be hashable (e.g tuple, string, ...) and the method
     215    # must not change the internal state of the object.  For all of the methods in this
     216    # class, this is the case.
     217    #
     218    # I then added a method that does simple argument checking, converting the arguments to
     219    # tuples, and invokes the private method to do the work.  Note that this method should
     220    # not be cached, as the user might specify the index of the component as a list, which is
     221    # not hashable.
     222    #
     223    # A second method returns a dictionary of all the components, again invoking the hidden
     224    # helper method.
     225    #
     226    # I found that on first invocation, this methods are about as fast as the old implementation.
     227    # The speedup arises when calling the method several times, as then the results are
     228    # cached.
     229    #
     230    # I've left the old method in the class, and added "_new" to the name of the new methods, so
     231    # that they can easily be compared.
     232   
     233    @cached_method
     234    def _compute_first_fundamental_form_coefficient(self, index):
     235        dr = self.natural_frame()
     236        return (dr[index[0]]*dr[index[1]]).simplify_full()       
     237   
     238    def first_fundamental_form_coefficient_new(self, index):
     239        index = tuple(sorted(index))
     240        if index not in self.index_list:
     241            raise ValueError, "Index %s out of bounds." % str(index)
     242        return self._compute_first_fundamental_form_coefficient(index)
     243
     244    @cached_method
     245    def first_fundamental_form_coefficients_new(self):
     246        coefficients = {}
     247        for index in self.index_list:
     248            sorted_index = list(sorted(index))
     249            coefficients[index] = \
     250                self._compute_first_fundamental_form_coefficient(index)
     251        return coefficients
     252
     253
     254    def first_fundamental_form_coefficients(self,index1=0,index2=0):
     255        """
     256        This function gives the coefficients of the first fundamental form
     257       
     258        INPUT:
     259        empty argument, or one of the the following pair of indices: 1,1; 1,2; 2,1; 2,2.
     260
     261        OUTPUT:
     262        with empty argument the output is the dictionary of coefficients of the first fundamental form
     263
     264        with given indices the output is the corresponding coefficient of the first fundamental form,
     265        for example index1 = 1, index2 = 2 gives the coefficient $g_{12}$
     266
     267        EXAMPLES::
     268
     269            sage: var('u,v')
     270            sage: sphere = parametrized_surface3d([cos(u)*cos(v),sin(u)*cos(v),sin(v)],[u,v],'sphere')
     271            sage: sphere.first_fundamental_form_coefficients()
     272            {(1, 2): 0, (1, 1): cos(v)^2, (2, 1): 0, (2, 2): 1}
     273
     274            sage: sphere.first_fundamental_form_coefficients(1,1)
     275            cos(v)^2
     276
     277            sage: sphere.first_fundamental_form_coefficients(2,1)
     278            0
     279
     280            sage: sphere.first_fundamental_form_coefficients(3,2)
     281            'The argument is not appropriate. Read doc'
     282        """
     283
     284        dr = self.natural_frame()
     285       
     286        indices=[index1,index2]
     287        if indices == [0,0]:
     288            g11 = (dr[1]*dr[1]).simplify_full()
     289            g12 = (dr[1]*dr[2]).simplify_full()
     290            g21 = g12
     291            g22 = (dr[2]*dr[2]).simplify_full()
     292            return {(1,1):g11,(1,2):g12,(2,1):g21,(2,2):g22}
     293        elif indices == [1,1]:
     294            return (dr[1]*dr[1]).simplify_full()
     295        elif indices == [1,2]:
     296            return (dr[1]*dr[2]).simplify_full()
     297        elif indices == [2,1]:
     298            return (dr[1]*dr[2]).simplify_full()
     299        elif indices == [2,2]:
     300            return (dr[2]*dr[2]).simplify_full()
     301        else:
     302            return 'The argument is not appropriate. Read doc'
     303
     304     
     305
     306    def first_fundamental_form(self,vector1,vector2):
     307        """
     308        Finds the value of first fundamental form on two vectors
     309
     310        INPUT:
     311        Two vectors $v=(v^1,v^2)$ and $w=(w^1,w^2)$
     312
     313        OUTPUT:
     314        $g_{11} v^1 w^1 + g_{12}(v^1 w^2 + v^2 w^1) + g_{22} v^2 w^2$
     315
     316        #EXAMPLES::
     317         
     318            sage: var('u,v')
     319            sage: var ('v1,v2,w1,w2')
     320            sage: sphere = parametrized_surface3d([cos(u)*cos(v),sin(u)*cos(v),sin(v)],[u,v],'sphere')
     321
     322            sage: sphere.first_fundamental_form(vector([v1,v2]),vector([w1,w2]))
     323            v1*w1*cos(v)^2 + v2*w2
     324
     325            sage: vv = vector([1,2])
     326            sage: sphere.first_fundamental_form(vv,vv)
     327            cos(v)^2 + 4
     328
     329            sage: sphere.first_fundamental_form([1,1],[2,1])
     330            2*cos(v)^2 + 1
     331        """
     332        gg = self.first_fundamental_form_coefficients()
     333        return sum(gg[ind]*vector1[ind[0]-1]*vector2[ind[1]-1] for ind in self.index_list)
     334
     335
     336    # jvkersch: cached this method
     337
     338    @cached_method
     339    def area_form_squared(self):
     340        """
     341        Gives $g_{11}g_{22} - g_{12}^2$, which is the square of the coefficient of the area form
     342       
     343        INPUT:
     344        No arguments
     345
     346        OUTPUT: 
     347        $g_{11}g_{22} - g_{12}^2$
     348
     349        EXAMPLES::
     350
     351            sage: var('u,v')
     352            sage: sphere = parametrized_surface3d([cos(u)*cos(v),sin(u)*cos(v),sin(v)],[u,v],'sphere')
     353            sage: sphere.area_form_squared()
     354            cos(v)^2
     355
     356        """
     357        gg = self.first_fundamental_form_coefficients()
     358        return  (gg[(1,1)]*gg[(2,2)]-gg[(1,2)]**2).simplify_full()
     359
     360
     361    # jvkersch: cached this method
     362   
     363    @cached_method
     364    def area_form(self):
     365        """
     366        Gives $\sqrt{g_{11}g_{22} - g_{12}^2}$, which is the coefficient of the area form
     367       
     368        INPUT:
     369        No arguments
     370
     371        OUTPUT: 
     372        $\sqrt{g_{11}g_{22} - g_{12}^2}$
     373
     374        EXAMPLES::
     375
     376            sage: var('u,v')
     377            sage: sphere = parametrized_surface3d([cos(u)*cos(v),sin(u)*cos(v),sin(v)],[u,v],'sphere')
     378            sage: sphere.area_form_squared()
     379            abs(cos(v))
     380
     381        """
     382        return sqrt(self.area_form_squared()).simplify_full()
     383
     384
     385    # jvkersch: The following two methods compute the components of the
     386    # inverse metric in a cached way.  Same way as before.
     387
     388    @cached_method
     389    def first_fundamental_form_inverse_coefficients_new(self):
     390
     391        gg = self.first_fundamental_form_coefficients()
     392        DD = gg[(1,1)]*gg[(2,2)]-gg[(1,2)]**2
     393       
     394        gi11 = (gg[(2,2)]/DD).simplify_full()
     395        gi12 = (-gg[(1,2)]/DD).simplify_full()
     396        gi21 = gi12
     397        gi22 = (gg[(1,1)]/DD).simplify_full()
     398
     399        return {(1,1):gi11,(1,2):gi12,(2,1):gi21,(2,2):gi22}
     400
     401    def first_fundamental_form_inverse_coefficient_new(self, index):
     402        index = tuple(sorted(index))
     403        if index not in self.index_list:
     404            raise ValueError, "Index %s out of bounds." % str(index)
     405        return self._compute_first_fundamental_form_inverse_coefficients()[index]
     406
     407
     408    # jvkersch: This is the original uncached method.
     409   
     410    def first_fundamental_form_inverse_coefficients(self):
     411        """
     412        Gives $g^{ij}$
     413       
     414        INPUT:
     415        No arguments
     416
     417        OUTPUT: 
     418        dictionary {(1,1):$g^{11}$, (1,2):$g^{12}$, (2,1):$g^{21}$, (2,2):$g^{22}$
     419
     420        EXAMPLES::
     421
     422            sage: var('u,v')
     423            sage: sphere = parametrized_surface3d([cos(u)*cos(v),sin(u)*cos(v),sin(v)],[u,v],'sphere')
     424            sage: sphere.first_fundamental_form_inverse_coefficients()
     425            {(1, 2): 0, (1, 1): cos(v)^(-2), (2, 1): 0, (2, 2): 1}
     426        """
     427
     428        gg = self.first_fundamental_form_coefficients()
     429        DD = gg[(1,1)]*gg[(2,2)]-gg[(1,2)]**2
     430        gi11 = (gg[(2,2)]/DD).simplify_full()
     431        gi12 = (-gg[(1,2)]/DD).simplify_full()
     432        gi21 = gi12
     433        gi22 = (gg[(1,1)]/DD).simplify_full()
     434        return {(1,1):gi11,(1,2):gi12,(2,1):gi21,(2,2):gi22}
     435
     436
     437    # jvkersch: cached this method.  No differences otherwise
     438   
     439    @cached_method
     440    def rotation(self,theta):
     441        """
     442        Gives the matrix of the operator of rotation on the given angle $\\theta$ with respect to the natural frame
     443       
     444        INPUT:
     445        given angle $\\theta$
     446
     447        OUTPUT: 
     448        matrix of the operator of rotation on $\\theta$ with respect to the natural frame
     449
     450        ALGORITHM:
     451        The operator of rotation on $\pi/2$ is $J^i_j = g^{ik}\omega_{jk}$, where $\omega$ is the area form
     452        The operator of rotation on angle $\\theta$ is $\cos(\\theta) I + sin(\\theta) J$
     453
     454        EXAMPLES::
     455
     456            sage: var('u,v')
     457            sage: assume(cos(v)>0)
     458            sage: sphere = parametrized_surface3d([cos(u)*cos(v),sin(u)*cos(v),sin(v)],[u,v],'sphere')
     459            sage: rotation = sphere.rotation(pi/3)
     460            sage: rotation^3
     461            [-1  0]
     462            [ 0 -1]
     463            # it is true because the rotation at $\pi/3$ applied three times gives rotation a $\pi$
     464        """
     465
     466        from sage.functions.trig import sin, cos
     467
     468        gi = self.first_fundamental_form_inverse_coefficients()
     469        w12 = self.area_form()
     470        RR11 = (cos(theta) + sin(theta)*gi[1,2]*w12).simplify_full()
     471        RR12 = (- sin(theta)*gi[1,1]*w12).simplify_full()
     472        RR21 = (sin(theta)*gi[2,2]*w12).simplify_full()
     473        RR22 = (cos(theta) - sin(theta)*gi[2,1]*w12).simplify_full()       
     474        rotation = matrix([[RR11,RR12],[RR21,RR22]])
     475        return rotation
     476
     477
     478    # jvkersch: original method -- cached version below.
     479
     480    def orthonormal_frame(self,vector_number=0,coordinates='ext'):
     481        """
     482        Gives the orthonormal frame field of the surface 
     483       
     484        INPUT:
     485        (vector_number,coordinates), where
     486        vector number is 0 (default), 1 or 2
     487        coordinates is 'ext'(default) or 'int'
     488               
     489
     490        OUTPUT:
     491        If vector_number is equal to 1 or 2, the output is the corresponding vector of an orthonormal frame.
     492        If argument vector_number is empty, the output is the dictionary of the vector fields of an orthonormal frame.
     493       
     494        If coordinates='ext', output is coordinates in $\mathbb{R}^3$
     495        If coordinates='int', output is coordinates with respect to the natural frame
     496
     497        ALGORITHM:
     498        We normalize the first vector $\\vec e_1$ of the natural frame and then get the second frame vector
     499        $\\vec e_2 = [\\vec n, \\vec e_1]$.
     500
     501        EXAMPLES::
     502
     503            sage: var('u,v')
     504            sage: assume(cos(v)>0)
     505            sage: sphere = parametrized_surface3d([cos(u)*cos(v),sin(u)*cos(v),sin(v)],[u,v],'sphere')
     506            sage: EE = sphere.orthonormal_frame()
     507            sage: (EE[1]*EE[1]).simplify_full()
     508            1
     509            sage: (EE[1]*EE[2]).simplify_full()
     510            0
     511            sage: simplify_vector(vector_product(EE[1],EE[2])-sphere.normal_vector(1))
     512            [0,0,0]
     513
     514            sage: EE_int = sphere.orthonormal_frame(coordinates='int')
     515            sage: EE_int
     516            {1: (1/cos(v), 0), 2: (0, 1)}
     517            sage: sphere.first_fundamental_form(EE_int[1],EE_int[1])
     518            1
     519            sage: sphere.first_fundamental_form(EE_int[1],EE_int[2])
     520            0
     521            sage: sphere.first_fundamental_form(EE_int[2],EE_int[2])
     522            1
     523           
     524        """
     525
     526        from sage.symbolic.constants import pi
     527
     528       
     529        if (vector_number,coordinates) not in [(0,'ext'),(0,'int'),(1,'ext'),(1,'int'),(2,'ext'),(2,'int')]:
     530           return "You set inappropriate arguments. Read the doc on this function"
     531
     532       # jvkersch: changed the calling convention of natural_frame and normal_vector
     533        if coordinates == 'ext':
     534            E1 = simplify_vector(self.natural_frame()[1]/sqrt(self.first_fundamental_form_coefficients(1,1)))
     535            if vector_number == 0:
     536                E2 = simplify_vector(vector_product(self.normal_vector(normalized=True),E1))
     537                return  {1:E1, 2:E2}
     538            elif vector_number == 1:
     539                return E1
     540            elif vector_number == 2:
     541                return simplify_vector(vector_product(self.normal_vector(normalized=True),E1))
     542        else:
     543            E1 =  vector([(1/sqrt(self.first_fundamental_form_coefficients(1,1))).simplify_full(),0])
     544            if vector_number == 0:
     545               E2 = simplify_vector(self.rotation(pi/2)*E1)
     546               return  {1:E1, 2:E2}
     547            elif vector_number == 1:
     548               return E1
     549            elif vector_number == 2:
     550               E2 = simplify_vector(self.rotation(pi/2)*E1)
     551               return E2
     552
     553
     554    # jvkersch: cached version of the code to compute the orthonormal frame
     555
     556    @cached_method
     557    def orthonormal_frame_vector(self, vector_number=0, coordinates='ext'):
     558       
     559        from sage.symbolic.constants import pi
     560        if (vector_number,coordinates) not in [(0,'ext'),(0,'int'),(1,'ext'),(1,'int'),(2,'ext'),(2,'int')]:
     561           return "You set inappropriate arguments. Read the doc on this function"
     562
     563
     564       # jvkersch: I made some changes to the definition of natural_frame and normal_vector,
     565       # so the following had to be modified slightly to use this new calling convention.
     566       
     567        if coordinates == 'ext':
     568            E1 = simplify_vector(self.natural_frame()[1]/sqrt(self.first_fundamental_form_coefficients(1,1)))
     569            if vector_number == 0:
     570                E2 = simplify_vector(vector_product(self.normal_vector(),E1))
     571                return  {1:E1, 2:E2}
     572            elif vector_number == 1:
     573                return E1
     574            elif vector_number == 2:
     575                return simplify_vector(vector_product(self.normal_vector(),E1))
     576        else:
     577            E1 =  vector([(1/sqrt(self.first_fundamental_form_coefficients(1,1))).simplify_full(),0])
     578            if vector_number == 0:
     579               E2 = simplify_vector(self.rotation(pi/2)*E1)
     580               return  {1:E1, 2:E2}
     581            elif vector_number == 1:
     582               return E1
     583            elif vector_number == 2:
     584               E2 = simplify_vector(self.rotation(pi/2)*E1)
     585               return E2
     586
     587    @cached_method
     588    def orthonormal_frame_new(self, coordinates='ext'):
     589        E1 = self.orthonormal_frame_vector(1, coordinates)
     590        E2 = self.orthonormal_frame_vector(2, coordinates)
     591        return {1: E1, 2: E2}
     592
     593       
     594    def lie_bracket(self,v,w):
     595        """
     596        Gives the Lie bracket of two vector fields which are given by coordinates with respect to the natural frame
     597         
     598       
     599        INPUT:
     600        v and w are vectors
     601               
     602
     603        OUTPUT:
     604        vector [v,w]
     605
     606
     607        EXAMPLES::
     608
     609            sage: var('u,v')
     610            sage: assume(cos(v)>0)
     611            sage: sphere = parametrized_surface3d([cos(u)*cos(v),sin(u)*cos(v),sin(v)],[u,v],'sphere')
     612            sage: sphere.lie_bracket([u,v],[-v,u])
     613            (0, 0)
     614
     615            sage: EE_int = sphere.orthonormal_frame(coordinates='int')
     616            sage: sphere.lie_bracket(EE_int[1],EE_int[2])
     617            (sin(v)/cos(v)^2, 0)
     618        """
     619        vv = vector([xx for xx in v])
     620        ww = vector([xx for xx in w])
     621        uuu = [self.variables[1],self.variables[2]]
     622        Dvv = matrix([ [ diff(xxx,uu).simplify_full() for uu in uuu ] for xxx in vv])   
     623        Dww = matrix([ [ diff(xxx,uu).simplify_full() for uu in uuu ] for xxx in ww])   
     624        return simplify_vector(vector(Dvv*ww - Dww*vv))
     625
     626
     627
     628    def frame_structure_functions(self,e1,e2):
     629        """
     630        Gives the structure functions $c^k_{ij}$ for a frame field $e_1$, $e_2$, where
     631        $[e_i,e_j] = c^k_{ij}e_k$
     632         
     633       
     634        INPUT:
     635        Frame vectors $e_1$, $e_2$
     636               
     637
     638        OUTPUT:
     639        The dictionary of $c^k_{ij}$.
     640        Warning: the tuple (i,j,k) corresponds to $c^k_{ij}$
     641
     642
     643        EXAMPLES::
     644
     645            sage: var('u,v')
     646            sage: assume(cos(v)>0)
     647            sage: sphere = parametrized_surface3d([cos(u)*cos(v),sin(u)*cos(v),sin(v)],[u,v],'sphere')
     648            sage: sphere.frame_structure_functions([u,v],[-v,u])
     649            {(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}
     650            sage: structfun = sphere.frame_structure_functions(EE_int[1],EE_int[2])
     651            sage: structfun
     652            {(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}
     653            sage: sphere.lie_bracket(EE_int[1],EE_int[2]) - CC[(1,2,1)]*EE_int[1] - CC[(1,2,2)]*EE_int[2]
     654            (0, 0)
     655            """
     656        ee1 = vector([xx for xx in e1])
     657        ee2 = vector([xx for xx in e2])
     658       
     659        lb = simplify_vector(self.lie_bracket(ee1,ee2))
     660        trans_matrix = matrix([[ee1[0],ee2[0]],[ee1[1],ee2[1]]])
     661        zz = simplify_vector(trans_matrix.inverse()*lb)
     662        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}
     663
     664
     665    # jvkersch: cached version of the code to compute second order frame -- original code below.
     666    # Same pattern as always.
     667
     668    @cached_method
     669    def _compute_second_order_frame_element(self, index):
     670        variables = [self.variables[i] for i in index]
     671        ddr_element = vector([diff(f, variables).simplify_full() for f in self.equation])
     672       
     673        return ddr_element
     674
     675    @cached_method
     676    def second_order_natural_frame_new(self):
     677        vectors = {}
     678        for index in self.index_list:
     679            sorted_index = tuple(sorted(index))
     680            vectors[index] = \
     681                self._compute_second_order_frame_element(sorted_index)
     682        return vectors
     683
     684    def second_order_natural_frame_element_new(self, index):
     685        index = tuple(sorted(index))
     686        if index not in self.index_list:
     687            raise ValueError, "Index %s out of bounds." % str(index)
     688        return self._compute_second_order_frame_element(index)
     689
     690
     691    # jvkersch: original code
     692
     693    def second_order_natural_frame(self,index1=0,index2=0):
     694        """
     695        Gives the second derivatives of the equation $\\vec r = \\vec r(u^1,u^2)$ of parametrized surface
     696         
     697       
     698        INPUT:
     699        empty argument, or one of the the following pair of indices: 1,1; 1,2; 2,1; 2,2.
     700               
     701
     702        OUTPUT:
     703        With empty argument the output is the dictionary of $\partial_{ij}\\vec r(u^1,u^2)$.
     704
     705        With given indices the output is the corresponding second derivative of the surface parametric equation,
     706        For example, index1 = 1, index2 = 2 gives $\partial_{12}\\vec r(u^1,u^2)$.
     707
     708
     709        EXAMPLES::
     710
     711            sage: var('u,v')
     712            sage: sphere = parametrized_surface3d([cos(u)*cos(v),sin(u)*cos(v),sin(v)],[u,v],'sphere')
     713            sage: sphere.second_order_natural_frame()
     714            {(1, 2): (sin(u)*sin(v), -sin(v)*cos(u), 0), (1, 1): (-cos(u)*cos(v),
     715            -sin(u)*cos(v), 0), (2, 1): (sin(u)*sin(v), -sin(v)*cos(u), 0), (2, 2):
     716            (-cos(u)*cos(v), -sin(u)*cos(v), -sin(v))}
     717           
     718            sage: sphere.second_order_natural_frame(1,1)
     719            (-cos(u)*cos(v), -sin(u)*cos(v), 0)
     720            sage: sphere.second_order_natural_frame(1,2)
     721            (sin(u)*sin(v), -sin(v)*cos(u), 0)
     722            sage: sphere.second_order_natural_frame(2,2)
     723            (-cos(u)*cos(v), -sin(u)*cos(v), -sin(v)
     724            """
     725        u = self.variables[1]
     726        v = self.variables[2]
     727        indices=[index1,index2]
     728       
     729        if indices == [0,0]:
     730            dr11 = vector([diff(f,u,2).simplify_full() for f in self.equation])
     731            dr12 = vector([diff(f,u,v).simplify_full() for f in self.equation])
     732            dr22 = vector([diff(f,v,2).simplify_full() for f in self.equation])
     733            return {(1,1):dr11,(1,2):dr12,(2,1):dr12,(2,2):dr22}
     734
     735        elif indices == [1,1]:
     736            return vector([diff(f,u,2).simplify_full() for f in self.equation])
     737        elif indices == [1,2]:
     738            return vector([diff(f,u,v).simplify_full() for f in self.equation])
     739        elif indices == [2,1]:
     740            return vector([diff(f,u,v).simplify_full() for f in self.equation])
     741        elif indices == [2,2]:
     742            return vector([diff(f,v,2).simplify_full() for f in self.equation])
     743        else:
     744            return 'The argument is not appropriate. Read doc'
     745
     746
     747
     748    # jvkersch: the following three functions compute the components of the second fundamental form
     749    # in a cached way.  Original code below.
     750
     751    @cached_method
     752    def _compute_second_fundamental_form_coefficient(self, index):
     753        NN = self.normal_vector(normalized=True)
     754        v  = self.second_order_natural_frame_element_new(index)
     755        return (v*NN).simplify_full()
     756       
     757    def second_fundamental_form_coefficient(self, index):
     758        index = tuple(index)
     759        if index not in self.index_list:
     760            raise ValueError, "Index %s out of bounds." % str(index)
     761        return self._compute_second_fundamental_form_coefficient(index)
     762
     763    @cached_method
     764    def second_fundamental_form_coefficients_new(self):
     765        coefficients = {}
     766        for index in self.index_list:
     767            coefficients[index] = \
     768                self._compute_second_fundamental_form_coefficient(index)
     769        return coefficients
     770
     771       
     772    # jvkersch: this is the original function   
     773    def second_fundamental_form_coefficients(self,index1=0,index2=0):
     774        """
     775        This function gives the coefficients $h_{ij}$ of the second fundamental form.
     776       
     777        INPUT:
     778        empty argument, or one of the the following pair of indices: 1,1; 1,2; 2,1; 2,2.
     779
     780        OUTPUT:
     781        with empty argument the output is the dictionary of coefficients of the second fundamental form
     782
     783        with given indices the output is the corresponding coefficient of the second fundamental form,
     784        for example index1 = 1, index2 = 2 gives the coefficient $h_{12}$
     785
     786        EXAMPLES::
     787
     788            sage: var('u,v')
     789            sage: assume(cos(v)>0)
     790            sage: sphere = parametrized_surface3d([cos(u)*cos(v),sin(u)*cos(v),sin(v)],[u,v],'sphere')
     791
     792            sage: sphere.second_fundamental_form_coefficients()
     793            {(1, 2): 0, (1, 1): -cos(v)^2, (2, 1): 0, (2, 2): -1}
     794            sage: sphere.second_fundamental_form_coefficients(1,1)
     795            -cos(v)^2
     796            sage: sphere.second_fundamental_form_coefficients(2,1)
     797            0
     798
     799            sage: sphere.second_fundamental_form_coefficients(3,2)
     800            'The argument is not appropriate. Read doc'
     801        """
     802        NN = self.normal_vector(1)
     803        indices=[index1,index2]
     804
     805        if indices == [0,0]:
     806            ddr = self.second_order_natural_frame()
     807            h11 = (ddr[1,1]*NN).simplify_full()
     808            h12 = (ddr[1,2]*NN).simplify_full()
     809            h21 = h12
     810            h22 = (ddr[2,2]*NN).simplify_full()
     811            return {(1,1):h11,(1,2):h12,(2,1):h21,(2,2):h22}
     812        elif indices == [1,1]:
     813            return (self.second_order_natural_frame(1,1)*NN).simplify_full()
     814        elif indices == [1,2]:
     815            return (self.second_order_natural_frame(1,2)*NN).simplify_full()
     816        elif indices == [2,1]:
     817            return (self.second_order_natural_frame(2,1)*NN).simplify_full()
     818        elif indices == [2,2]:
     819            return (self.second_order_natural_frame(2,2)*NN).simplify_full()
     820        else:
     821            return 'The argument is not appropriate. Read doc'
     822
     823
     824    def second_fundamental_form(self,vector1,vector2):
     825        """
     826        Finds the value of second fundamental form on two vectors
     827
     828        INPUT:
     829        Two vectors $v=(v^1,v^2)$ and $w=(w^1,w^2)$
     830
     831        OUTPUT:
     832        $h_{11} v^1 w^1 + h_{12}(v^1 w^2 + v^2 w^1) + h_{22} v^2 w^2$
     833
     834        EXAMPLES::
     835
     836           sage: var('u,v')
     837           sage: var ('v1,v2,w1,w2')
     838           sage: assume(cos(v) > 0)
     839           sage: sphere = parametrized_surface3d([cos(u)*cos(v),sin(u)*cos(v),sin(v)],[u,v],'sphere')
     840           sage: sphere.second_fundamental_form(vector([v1,v2]),vector([w1,w2]))
     841           -v1*w1*cos(v)^2 - v2*w2
     842           sage: vv = vector([1,2])
     843           sage: sphere.second_fundamental_form(vv,vv)
     844           -cos(v)^2 - 4
     845           sage: sphere.second_fundamental_form([1,1],[2,1])
     846           -2*cos(v)^2 - 1
     847                   
     848        """
     849        hh = self.second_fundamental_form_coefficients()
     850        return sum(hh[ind]*vector1[ind[0]-1]*vector2[ind[1]-1] for ind in self.index_list)
     851
     852
     853    # jvkersch: cached this method
     854
     855    @cached_method
     856    def gauss_curvature(self):
     857        """
     858        Finds the gaussian curvature $K = \\frac{h_{11}h_{22} - h_{12}^2}{g_{11}g_{22} - g_{12}^2}$.
     859
     860        INPUT:
     861        No arguments
     862
     863        OUTPUT:
     864        $K = \\frac{h_{11}h_{22} - h_{12}^2}{g_{11}g_{22} - g_{12}^2}$
     865
     866        EXAMPLES::
     867
     868           sage: var('R')
     869           sage: assume(R>0)
     870           sage: var('u,v')
     871           sage: assume(cos(v)>0)
     872           sage: sphere = parametrized_surface3d([R*cos(u)*cos(v),R*sin(u)*cos(v),R*sin(v)],[u,v],'sphere')
     873           sage: sphere.gauss_curvature()
     874           R^(-2)
     875                   
     876        """
     877        hh = self.second_fundamental_form_coefficients()
     878        return ((hh[(1,1)]*hh[(2,2)]-hh[(1,2)]**2)/self.area_form_squared()).simplify_full()
     879
     880
     881    # jvkersch: cached this method
     882
     883    @cached_method
     884    def mean_curvature(self):
     885        """
     886        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}$.
     887
     888        INPUT:
     889        No arguments
     890
     891        OUTPUT:
     892        $H = \\frac{1}{2}\\frac{g_{22}h_{11} - 2g_{12}h_{12} + g_{11}h_{22}}{g_{11}g_{22} - g_{12}^2}$.
     893
     894        EXAMPLES::
     895           
     896           sage: var('R')
     897           sage: assume(R>0)
     898           sage: var('u,v')
     899           sage: assume(cos(v)>0)
     900           sage: sphere = parametrized_surface3d([R*cos(u)*cos(v),R*sin(u)*cos(v),R*sin(v)],[u,v],'sphere')
     901           sage: sphere.mean_curvature()
     902           -1/R
     903
     904        """
     905        gg = self.first_fundamental_form_coefficients()
     906        hh = self.second_fundamental_form_coefficients()
     907        denom = 2*self.area_form_squared()
     908        enum =(gg[(2,2)]*hh[(1,1)]-2*gg[(1,2)]*hh[(1,2)]+gg[(1,1)]*hh[(2,2)]).simplify_full()
     909        return (enum/denom).simplify_full()
     910
     911
     912    # jvkersch: cached this method
     913   
     914    @cached_method
     915    def principal_curvatures(self):
     916        """
     917        Finds the principal curvatures of the surface
     918
     919        INPUT:
     920        No arguments
     921
     922        OUTPUT:
     923        The dictionary of principal curvatures         
     924
     925        EXAMPLES::
     926                     
     927           sage: var('R')
     928           sage: assume(R>0)
     929           sage: var('u,v')
     930           sage: assume(cos(v)>0)
     931           sage: sphere = parametrized_surface3d([R*cos(u)*cos(v),R*sin(u)*cos(v),R*sin(v)],[u,v],'sphere')
     932           sage: sphere.principal_curvatures()
     933           {1: -1/R, 2: -1/R}
     934
     935           sage: var('u,v')
     936           sage: var('R,r')
     937           sage: assume(R>r,r>0)
     938           sage: torus = parametrized_surface3d([(R+r*cos(v))*cos(u),(R+r*cos(v))*sin(u),r*sin(v)],[u,v],'torus')
     939           sage: torus.principal_curvatures()
     940           {1: -cos(v)/(r*cos(v) + R), 2: -1/r}
     941
     942        """
     943
     944        from sage.symbolic.assumptions import assume
     945        from sage.symbolic.relation import solve
     946        from sage.calculus.var import var
     947
     948        KK = self.gauss_curvature()
     949        HH = self.mean_curvature()
     950
     951        # jvkersch: when this assumption is uncommented, Sage raises an error stating that the assumption
     952        # is redundant... Can we safely omit this, based on some geometric reasoning?
     953       
     954        # assume(HH**2-KK>=0)
     955
     956        x = var('x')
     957        sol = solve(x**2 -2*HH*x + KK == 0,x)
     958       
     959        #k1=var('k1')
     960        #k2=var('k2')
     961
     962        # 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.
     963
     964        k1 = (x.substitute(sol[0])).simplify_full()
     965        if len(sol) == 1:
     966            k2 = k1
     967        else:
     968            k2 = (x.substitute(sol[1])).simplify_full()
     969       
     970        #k1 = (x.substitute(sol[0])).simplify_full()
     971        #k2 = (x.substitute(sol[1])).simplify_full()
     972
     973        # 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]
     974
     975        return {1:k1,2: k2}
     976
     977
     978###############################################
     979#
     980#  jvkersch: below I didn't add any code to cache the member functions
     981#  except for the numerical code that computes the geodesics.
     982#
     983###############################################
     984
     985
     986    def principal_directions(self):
     987        """
     988        Finds the principal curvatures and principal directions of the surface
     989
     990        INPUT:
     991        No arguments
     992
     993        OUTPUT:
     994        The dictionary of lists [a principal curvature, the corresponding principal direction]         
     995
     996        If principal curvatures coincide, gives the warning that the surface is a sphere.
     997
     998        EXAMPLES::
     999
     1000           sage: var('u,v')
     1001           sage: var('R,r')
     1002           sage: assume(R>r,r>0)
     1003           sage: torus = parametrized_surface3d([(R+r*cos(v))*cos(u),(R+r*cos(v))*sin(u),r*sin(v)],[u,v],'torus')
     1004           sage: pdd = torus.principal_directions()
     1005           sage: pdd[1]
     1006           [-cos(v)/(r*cos(v) + R), (1, 0)]
     1007           sage: pdd[2]
     1008           [-1/r, (0, -(R*r*cos(v) + R^2)/r)]
     1009
     1010           sage: var('RR')
     1011           sage: assume(RR>0)
     1012           sage: var('u,v')
     1013           sage: assume(cos(v)>0)
     1014           sage: sphere = parametrized_surface3d([RR*cos(u)*cos(v),RR*sin(u)*cos(v),RR*sin(v)],[u,v],'sphere')
     1015           sage: sphere.principal_directions()
     1016           'This is a sphere, so any direction is principal'
     1017
     1018           sage: var('aa')
     1019           sage: assume(aa>0)
     1020           sage: catenoid = parametrized_surface3d([aa*cosh(v)*cos(u),aa*cosh(v)*sin(u),v],[u,v],'catenoid')
     1021           sage: pd = catenoid.principal_directions()
     1022           sage: pd[1][1]
     1023           (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))
     1024           sage: pd[2][1]
     1025           (1, 0)
     1026           sage: pd[1][1]*pd[2][1]
     1027           0
     1028        """
     1029        gg = self.first_fundamental_form_coefficients()
     1030        hh = self.second_fundamental_form_coefficients()
     1031        kk = self.principal_curvatures()
     1032        if kk[1]==kk[2]:
     1033            return "This is a sphere, so any direction is principal"
     1034        pd1 = simplify_vector([hh[(1,2)]-kk[1]*gg[(1,2)],-hh[(1,1)]+kk[1]*gg[(1,1)]])
     1035        if pd1==vector([0,0]):
     1036           pd1 = vector([1,0])
     1037        pd2 = simplify_vector([hh[(1,2)]-kk[2]*gg[(1,2)],-hh[(1,1)]+kk[2]*gg[(1,1)]])
     1038        if pd2==vector([0,0]):
     1039           pd2 = vector([1,0])
     1040        return {1:[kk[1],pd1],2:[kk[2],pd2]}
     1041       
     1042    def connection_coefficients(self):
     1043        """
     1044        Finds the connection coefficients of the surface
     1045
     1046        INPUT:
     1047        No arguments
     1048
     1049        OUTPUT:
     1050        The dictionary of connection coefficients.
     1051
     1052        Warning: the triple $(i,j,k)$ corresponds to $\Gamma^k_{ij}$. 
     1053
     1054        EXAMPLES::
     1055                                 
     1056           sage: var('r')
     1057           sage: assume(r > 0)
     1058           sage: var('u,v')
     1059           sage: assume(cos(v)>0)
     1060           sage: sphere = parametrized_surface3d([r*cos(u)*cos(v),r*sin(u)*cos(v),r*sin(v)],[u,v],'sphere')
     1061           sage: sphere.connection_coefficients()
     1062           {(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}
     1063           
     1064        """
     1065        x = self.variables
     1066        gg = self.first_fundamental_form_coefficients()
     1067        gi = self.first_fundamental_form_inverse_coefficients()
     1068
     1069        dg = {}
     1070        for kkk in self.index_list_3:
     1071            dg[kkk]=gg[(kkk[1],kkk[2])].differentiate(x[kkk[0]]).simplify_full()
     1072        structfun={}
     1073
     1074        for kkk in self.index_list_3:
     1075            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()
     1076        return structfun
     1077
     1078
     1079    def geodesics_numerical(self,p0,v0,tinterval):
     1080        """
     1081        This method gives the numerical solution for the geodesic equations 
     1082
     1083        INPUT:
     1084        p0 is the list of the coordinates of the initial point
     1085       
     1086        v0 is the  list of the coordinates of the initial vector
     1087
     1088        tinterval is the list [a,b,M], where (a,b) is the domain of the geodesic, M is the number of division points
     1089
     1090        OUTPUT:
     1091        The list consisting of lists [t, [u1(t),u2(t)], [v1(t),v2(t)], [x1(t),x2(t),x3(t)]], where
     1092
     1093        t is a subdivision point
     1094
     1095        [u1(t),u2(t)] is the list of coordinates of the geodesic point
     1096         
     1097        [v1(t),v2(t)] is the list of coordinates of the vector tanget to the geodesic 
     1098
     1099        [x1(t),x2(t),x3(t)] is the list of coordinates of the geodesic point in the three-dimensional space
     1100
     1101
     1102        EXAMPLES::
     1103
     1104           sage: var('p,q')
     1105           sage: v = [p,q]
     1106           sage: assume(cos(q)>0)
     1107           sage: sphere = parametrized_surface3d([cos(q)*cos(p),sin(q)*cos(p),sin(p)],v,'sphere')
     1108           sage: gg_array = sphere.geodesics_numerical([0,0],[1,1],[0,2*pi,5])
     1109           sage: gg_array[0]
     1110           [0.0, [0.0, 0.0], [1.0, 1.0], (1, 0, 0)]
     1111           sage: gg_array[1]
     1112           [1.2566370614359172, [0.76440104189216407, 1.8586223516062499], [-0.2838683571264714, 1.9194187087799912], (-0.204895333443, 0.692104654602, 0.692104796553)]
     1113                                           
     1114         
     1115        """
     1116
     1117        # jvkersch: rewrote a few things in this routine to make it pure python, i.e. without the sage extensions.
     1118
     1119        #from sage.calculus.var import var
     1120
     1121        u1 = self.variables[1]
     1122        u2 = self.variables[2]
     1123        #uu1, uu2 = var('uu1,uu2')
     1124
     1125        # jvkersch: hack for now
     1126        #def f(ux, vx):
     1127        #    return [self.equation[0].subs({u1: ux, u2: vx}),
     1128        #            self.equation[1].subs({u1: ux, u2: vx}),
     1129        #            self.equation[2].subs({u1: ux, u2: vx})]
     1130
     1131        C111 = self.connection_coefficients()[(1,1,1)]
     1132        C121 = self.connection_coefficients()[(1,2,1)]
     1133        C221 = self.connection_coefficients()[(2,2,1)]
     1134        C112 = self.connection_coefficients()[(1,1,2)]
     1135        C122 = self.connection_coefficients()[(1,2,2)]
     1136        C222 = self.connection_coefficients()[(2,2,2)]
     1137 
     1138        def ode_system(unknown_functions,t):
     1139            uu1,uu2,v1,v2 = unknown_functions
     1140
     1141            # 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.
     1142            c111 = C111.subs({u1: uu1, u2: uu2})
     1143            c121 = C121.subs({u1: uu1, u2: uu2})
     1144            c221 = C221.subs({u1: uu1, u2: uu2})
     1145            c112 = C112.subs({u1: uu1, u2: uu2})
     1146            c122 = C122.subs({u1: uu1, u2: uu2})
     1147            c222 = C222.subs({u1: uu1, u2: uu2})
     1148
     1149            du1 = v1
     1150            du2 = v2
     1151            dv1 = - c111*v1**2 - 2*c121*v1*v2 - c221*v2**2
     1152            dv2 = - c112*v1**2 - 2*c122*v1*v2 - c222*v2**2
     1153            return float(du1), float(du2), float(dv1), float(dv2)
     1154
     1155
     1156        step = float((tinterval[1]-tinterval[0])/tinterval[2])
     1157        ttt =  float(tinterval[0])
     1158        tarray = [ ttt ]
     1159        for counter in range(1,tinterval[2]+1):
     1160           ttt = ttt + step
     1161           tarray = tarray + [ ttt ]
     1162
     1163        # import ode solving routine
     1164        import scipy.integrate
     1165
     1166        initial_data = (p0[0],p0[1],v0[0],v0[1])
     1167        solution = scipy.integrate.odeint(ode_system,initial_data,tarray)
     1168        return [[ tarray[counter],[solution[counter,0], solution[counter,1]], [solution[counter,2], solution[counter,3]], self.eq_callable(solution[counter,0],solution[counter,1])]   for counter in range(0,len(tarray))]
     1169
     1170
     1171    # jvkersch: this private method creates an ode_solver object, which can be used
     1172    # to integrate the geodesic equations numerically.
     1173   
     1174    @cached_method
     1175    def _create_ode_system(self):
     1176        from sage.ext.fast_eval import fast_float
     1177        from sage.calculus.var import var
     1178        from sage.gsl.ode import ode_solver
     1179
     1180        u1 = self.variables[1]
     1181        u2 = self.variables[2]
     1182        v1, v2 = var('v1, v2')
     1183
     1184        C = self.connection_coefficients()
     1185       
     1186        dv1 = - C[(1,1,1)]*v1**2 - 2*C[(1,2,1)]*v1*v2 - C[(2,2,1)]*v2**2
     1187        dv2 = - C[(1,1,2)]*v1**2 - 2*C[(1,2,2)]*v1*v2 - C[(2,2,2)]*v2**2
     1188        fun1 = fast_float(dv1, str(u1), str(u2), str(v1), str(v2))
     1189        fun2 = fast_float(dv2, str(u1), str(u2), str(v1), str(v2))
     1190
     1191        geodesic_ode = ode_solver()
     1192        geodesic_ode.function = \
     1193                              lambda t, (u1, u2, v1, v2) : \
     1194                              [v1, v2, fun1(u1, u2, v1, v2), fun2(u1, u2, v1, v2)]
     1195        return geodesic_ode
     1196
     1197
     1198    # jvkersch: integrate the geodesic equations numerically
     1199       
     1200    def geodesics_numerical_fast(self, p0, v0, tinterval):
     1201        u1 = self.variables[1]
     1202        u2 = self.variables[2]
     1203
     1204        solver = self._create_ode_system()
     1205           
     1206        t_interval, n = tinterval[0:2], tinterval[2]
     1207        solver.y_0 = [p0[0], p0[1], v0[0], v0[1]]
     1208        solver.ode_solve(t_span=t_interval, num_points=n)
     1209
     1210        parsed_solution = \
     1211          [[vec[0], vec[1][0:2], vec[1][2:], self.eq_callable(vec[1][0], vec[1][1])] \
     1212           for vec in solver.solution]
     1213
     1214        return parsed_solution
     1215
     1216
     1217    def parallel_translation_numerical(self,curve,t,v0,tinterval):
     1218        """
     1219        This method gives the numerical solution to the equation of parallel translation of a vector
     1220
     1221        INPUT:
     1222        curve equation = list of functions which determine the curve wrt the local coordinate
     1223
     1224        t - curve parameter
     1225
     1226        v0 - initial vector
     1227
     1228        tinterval = [a,b,N], (a,b) is the domain of the curve, N is the number of subdivision points
     1229
     1230        OUTPUT:
     1231        The list consisting of lists [t, [v1(t),v2(t)]], where
     1232
     1233        t is a subdivision point
     1234
     1235        [v1(t),v2(t)] is the list of coordinates of the vector translated parallely along the curve
     1236
     1237        EXAMPLES::
     1238
     1239           sage: var('p,q')
     1240           sage: v = [p,q]
     1241           sage: assume(cos(q)>0)
     1242           sage: sphere = parametrized_surface3d([cos(q)*cos(p),sin(q)*cos(p),sin(p)],v,'sphere')
     1243           sage: var('ss')
     1244           sage: vv_array = sphere.parallel_translation_numerical([ss,ss],ss,[1,1],[0,pi/4,20])
     1245           sage: vv_array[0]
     1246           [0.0, [1.0, 1.0]]
     1247           sage: vv_array[5]
     1248           [0.19634954084936207, [0.98060182522638917, 1.0389930474425113]]
     1249
     1250        """
     1251        u1 = self.variables[1]
     1252        u2 = self.variables[2]
     1253        #var('uu1,uu2')
     1254       
     1255        #f(u1,u2) = self.equation
     1256
     1257        def f(ux, vx):
     1258            return [self.equation[0].subs({u1: ux, u2: vx}),
     1259                    self.equation[1].subs({u1: ux, u2: vx}),
     1260                    self.equation[2].subs({u1: ux, u2: vx})]
     1261
     1262        #from sage.calculus.var import var
     1263        #t = var('t')
     1264        tt = t
     1265
     1266        C111 = self.connection_coefficients()[(1,1,1)]
     1267        C121 = self.connection_coefficients()[(1,2,1)]
     1268        C221 = self.connection_coefficients()[(2,2,1)]
     1269        C112 = self.connection_coefficients()[(1,1,2)]
     1270        C122 = self.connection_coefficients()[(1,2,2)]
     1271        C222 = self.connection_coefficients()[(2,2,2)]
     1272
     1273        du1=diff(curve[0],tt)
     1274        du2=diff(curve[1],tt)
     1275        uu1=curve[0]
     1276        uu2=curve[1]
     1277
     1278
     1279        def ode_system(unknown_functions,t1):
     1280            du1p = du1.subs(tt==t1)
     1281            du2p = du2.subs(tt==t1)
     1282            uu1p = uu1.subs(tt==t1)
     1283            uu2p = uu2.subs(tt==t1)
     1284
     1285            c111 = C111.subs({u1: uu1p, u2: uu2p})
     1286            c121 = C121.subs({u1: uu1p, u2: uu2p})
     1287            c221 = C221.subs({u1: uu1p, u2: uu2p})
     1288            c112 = C112.subs({u1: uu1p, u2: uu2p})
     1289            c122 = C122.subs({u1: uu1p, u2: uu2p})
     1290            c222 = C222.subs({u1: uu1p, u2: uu2p})
     1291
     1292            #print c111, c121, c221, c112, c122, c222
     1293           
     1294            v1,v2 = unknown_functions
     1295            f1 = - c111*du1p*v1 - c121*(du1p*v2 + du2p*v1) - c221*du2p*v2
     1296            f2 = - c112*du1p*v1 - c122*(du1p*v2 + du2p*v1) - c222*du2p*v2
     1297            dv1 = f1.subs(tt==t1)
     1298            dv2 = f2.subs(tt==t1)
     1299
     1300            #print dv1, dv2
     1301           
     1302            return float(dv1), float(dv2)
     1303
     1304        step = float((tinterval[1]-tinterval[0])/tinterval[2])
     1305        ttt =  float(tinterval[0])
     1306        tarray = [ ttt ]
     1307        for counter in range(1,tinterval[2]+1):
     1308           ttt = ttt + step
     1309           tarray = tarray + [ ttt ]
     1310
     1311
     1312        # import ode solving routine
     1313        import scipy.integrate
     1314
     1315        initial_data = (v0[0],v0[1])
     1316        solution = scipy.integrate.odeint(ode_system,initial_data,tarray)
     1317        return [[ tarray[counter],[solution[counter,0], solution[counter,1]]]   for counter in range(0,len(tarray))]
     1318
     1319
  • new file sage/riemann/vector_functions.py

    diff -r f667e86a25fb -r 38ffdf195219 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
     7#def scalar_product(a,b):
     8#    return(a[0]*b[0]+a[1]*b[1]+a[2]*b[2])
     9#def 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 f667e86a25fb -r 38ffdf195219 setup.py
    a b  
    985985                     'sage.quadratic_forms',
    986986                     'sage.quadratic_forms.genera',
    987987
     988                     'sage.riemann',
     989
    988990                     'sage.rings',
    989991                     'sage.rings.finite_rings',
    990992                     'sage.rings.number_field',