Ticket #10132: parametrized_surface3d.py

File parametrized_surface3d.py, 50.3 KB (added by mikarm, 10 years ago)

version 2 of parametrized_surface3d.py with comments

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