Ticket #19978: Rn.py

File Rn.py, 15.8 KB (added by ernestyalumni, 6 years ago)

Rn.py

Line 
1## Rn.py
2## Sage Math Sagemanifolds implementation of Euclidean R^n as classes
3## namespace or names follow closely the Tutorial pdf on sagemanifolds webpage:
4## http://sagemanifolds.obspm.fr/examples/pdf/SM_tutorial.pdf
5############################################################################
6## Copyleft 2015, Ernest Yeung <ernestyalumni@gmail.com>                           
7## 20160109
8##                                                                               
9## This program, along with all its code, is free software; you can redistribute
10## it and/or modify it under the terms of the GNU General Public License as
11## published by the Free Software Foundation; either version 2 of the License, or   
12## (at your option) any later version.                                       
13##     
14## linkedin     : ernestyalumni                                                   
15## wordpress    : ernestyalumni                                                   
16############################################################################
17t = var('t')
18assume(t,"real")
19
20class R1(object):
21      def __init__(self):
22          self.M = Manifold(2,'R1',r'\mathbb{R}',start_index=1)
23          self.cart_ch = self.M.chart('x')
24         
25
26class R2(object):
27      def __init__(self):
28          self.M = Manifold(2,'R2',r'\mathbb{R}^2',start_index=1)
29          self.cart_ch = self.M.chart('x y')
30          self.U       = self.M.open_subset('U',
31                                coord_def={self.cart_ch: 
32                                        (self.cart_ch[1]<0, self.cart_ch[2]!=0)})         # cf. http://sagemanifolds.obspm.fr/examples/pdf/SM_tutorial.pdf "Introducing a second chart on the manifold" the condition AND written with [] instead of ()
33
34          self.sph_ch  = self.U.chart(r'r:(0,+oo) ph:(0,2*pi):\phi')
35          self.cart_ch_U = self.cart_ch.restrict(self.U)
36          self.transit_sph_to_cart = self.sph_ch.transition_map(self.cart_ch_U, 
37                                   [self.sph_ch[1]*cos(self.sph_ch[2]), 
38                                    self.sph_ch[1]*sin(self.sph_ch[2])])
39          Eucnorm = sqrt( sum([self.cart_ch_U[i[0]]**2 for i in self.M.index_generator(1)]) )
40          self.transit_sph_to_cart.set_inverse( Eucnorm, 
41                                                atan2(self.cart_ch_U[2],self.cart_ch_U[1]) )
42      def equip_metric(self):
43          self.g = self.M.riemannian_metric('g')
44          for i in self.M.index_generator(1):
45              self.g[i[0],i[0]] = 1 
46
47      def make_orthon_frames(self,ch):
48          """
49          make_orthon_frames=make_orthon_frames(self,ch)
50          This method creates a change-of-basis matrix for orthonormal coordinates as
51               to_orthonormal
52          and a new orthonormal frame from the input of a (spherical coordinates) chart, ch
53
54          PARAMETERS (INPUTS):
55          ch - <a sagemanifolds chart>
56
57          OUTPUT
58          to_orthonormal - a change-of-basis matrix
59          eo - new orthonormal frame
60         
61          EXAMPLES of USAGE:
62          R2eg = R2()     
63          R2.equip_metric()
64          R2.make_orthon_frames(R2.sph_ch)       
65
66          """
67          try:
68              to_orthonormal = ch.domain().automorphism_field()
69              for i,j in self.M.index_generator(2):
70                  if self.g[ch.frame(),i,j,ch]!= 0:
71                     g_ij = self.g[ch.frame(),i,j,ch]
72                     to_orthonormal[ch.frame(),i,j,ch] = 1/sqrt(g_ij)
73              esph = ch.frame().new_frame(to_orthonormal,'e')               
74
75              # cf. https://github.com/sagemanifolds/SageManifolds/blob/master/Worksheets/v0.9/SM_Cartesian_spherical-3D.ipynb for explanation on the change of frame
76              Jacobian_sph_to_cart = ch.domain().change_of_frame(ch.frame(), ch.domain().default_chart().frame() )
77              Jacobian_och_to_sph  = ch.domain().change_of_frame(esph, ch.frame() )
78
79              ch.domain().set_change_of_frame(ch.domain().default_frame(), esph, 
80                                                Jacobian_och_to_sph.inverse()*Jacobian_sph_to_cart.inverse())
81
82              ch.domain().set_change_of_frame( esph, ch.domain().default_frame(),
83                                               Jacobian_sph_to_cart*Jacobian_och_to_sph )
84
85              return to_orthonormal, esph, [Jacobian_sph_to_cart, Jacobian_och_to_sph]
86
87          except AttributeError:
88              print "Equip the manifold with a metric g by doing the method .equip_metric() first!"
89
90
91class R3(object):
92      def __init__(self):
93          self.M = Manifold(3,'R3',r'\mathbb{R}^3',start_index=1)
94          self.cart_ch = self.M.chart('x y z')
95         
96          self.U = self.M.open_subset('U',coord_def={self.cart_ch: (self.cart_ch[1]<0, self.cart_ch[2]!=0)})
97          self.cart_ch_U = self.cart_ch.restrict(self.U)
98         
99          self.sph_ch  = self.U.chart(r'rh:(0,+oo):\rho th:(0,pi):\theta ph:(0,2*pi):\phi')
100          rh,th,ph = [self.sph_ch[i[0]] for i in self.M.index_generator(1)] 
101          self.transit_sph_to_cart = self.sph_ch.transition_map(self.cart_ch_U,
102                                        [rh*sin(th)*cos(ph),rh*sin(th)*sin(ph),rh*cos(th)])
103          Sphnorm = sqrt(sum([self.cart_ch_U[i[0]]**2 for i in self.M.index_generator(1)]))
104          self.transit_sph_to_cart.set_inverse( Sphnorm,atan2( sqrt( sum([ self.cart_ch_U[i]**2 for i in range(1,3)])), self.cart_ch_U[3]), atan2( self.cart_ch_U[2],self.cart_ch_U[1]) )
105
106         
107          self.cyl_ch = self.U.chart(r'r:(0,+oo) phi:(0,2*pi):\phi zc')
108          r,phi,zc = [self.cyl_ch[i[0]] for i in self.M.index_generator(1)]
109          self.transit_cyl_to_cart = self.cyl_ch.transition_map(self.cart_ch_U,
110                                        [r*cos(phi),r*sin(phi),zc])
111          self.transit_cyl_to_cart.set_inverse( sqrt( self.cart_ch_U[1]**2 + self.cart_ch_U[2]**2 ) , atan2( self.cart_ch_U[2],self.cart_ch_U[1]), self.cart_ch_U[3] )
112         
113      def equip_metric(self):
114          self.g = self.M.riemannian_metric('g')
115          for i in self.M.index_generator(1):
116              self.g[i[0],i[0]] = 1
117
118      def make_orthon_frames(self,ch):
119          """
120          make_orthon_frames=make_orthon_frames(self,ch)
121          This method creates a change-of-basis matrix for orthonormal coordinates as
122               to_orthonormal
123          and a new orthonormal frame from the input of a (spherical coordinates) chart, ch
124
125          PARAMETERS (INPUTS):
126          ch - <a sagemanifolds chart>
127
128          OUTPUT
129          to_orthonormal - a change-of-basis matrix
130          eo - new orthonormal frame
131         
132          EXAMPLES of USAGE:
133          R3eg = R3()     
134          R3.equip_metric()
135          R3.make_orthon_frames(R3.sph_ch)       
136
137          """
138          try:
139              to_orthonormal = ch.domain().automorphism_field()
140              for i,j in self.M.index_generator(2):
141                  if self.g[ch.frame(),i,j,ch]!= 0:
142                     g_ij = self.g[ch.frame(),i,j,ch]
143                     to_orthonormal[ch.frame(),i,j,ch] = 1/sqrt(g_ij)
144              eoch = ch.frame().new_frame(to_orthonormal,'e')               
145             
146              # cf. https://github.com/sagemanifolds/SageManifolds/blob/master/Worksheets/v0.9/SM_Cartesian_spherical-3D.ipynb for explanation on the change of frame
147              Jacobian_ch_to_cart = ch.domain().change_of_frame(ch.frame(), ch.domain().default_chart().frame() )
148              Jacobian_och_to_ch  = ch.domain().change_of_frame(eoch, ch.frame() )
149
150              ch.domain().set_change_of_frame(ch.domain().default_frame(), eoch, 
151                                                Jacobian_och_to_ch.inverse()*Jacobian_ch_to_cart.inverse())
152
153              ch.domain().set_change_of_frame( eoch, ch.domain().default_frame(),
154                                               Jacobian_ch_to_cart*Jacobian_och_to_ch )
155
156
157              return to_orthonormal, eoch, [Jacobian_ch_to_cart, Jacobian_och_to_ch]
158
159          except AttributeError:
160              print "Equip the manifold with a metric g by doing the method .equip_metric() first!"
161
162
163class Rn(object):
164      def __init__(self,n):
165          assert n>0
166          if n == 2:
167             print "Use the class R2"
168          elif n == 3:
169             print "Use the class R3"
170          else:
171             self.M = Manifold(n,'R'+str(n),r'\mathbb{R}^'+str(n),start_index=1)
172             self.cart_ch = self.M.chart(r" ".join([r"x"+str(i) for i in range(1,n+1)]))
173             xis = [self.cart_ch[i[0]] for i in self.M.index_generator(1)]
174
175             self.U = self.M.open_subset('U',coord_def={self.cart_ch:(xis[0]<0,xis[1]!=0)})
176             self.cart_ch_U = self.cart_ch.restrict(self.U)
177
178             # spherical coordinates         
179             self.sph_ch = self.U.chart(r'rh:(0,+oo):\rho '+r" ".join([r"th"+str(i)+r":(0,pi)" for i in range(1,n+1-2)])+r' ph:(0,2*pi):\phi')
180             sphs = [self.sph_ch[i[0]] for i in self.M.index_generator(1)]
181             self.transit_sph_to_cart = self.sph_ch.transition_map(self.cart_ch_U,
182                                      [sphs[0]*prod([sin(sphs[i]) for i in range(1,n-1)])*cos(sphs[-1]), sphs[0]*prod([sin(sphs[i]) for i in range(1,n-1)])*sin(sphs[-1])]+
183                                       [sphs[0]*prod([sin(sphs[i]) for i in range(1,j)])*cos(sphs[j]) for j in range(n-2,1,-1)]+[sphs[0]*cos(sphs[1]),])
184             gen_transit_list_sph = [ sqrt(sum([ xis[i]**2 for i in range(len(xis))])),] + [atan2( sqrt( sum([xis[i]**2 for i in range(j)])),xis[j]) for j in range(n-1,1,-1)]+[atan2(xis[1],xis[0]),] 
185             self.transit_sph_to_cart.set_inverse(*gen_transit_list_sph)
186
187             # cylindrical coordinates
188             self.cyl_ch = self.U.chart(r'r:(0,+oo) '+r" ".join([r"the"+str(i)+r":(0,pi)" for i in range(1,n+1-3)])+r' phi:(0,2*pi):\varphi z')
189             cyls = [self.cyl_ch[i[0]] for i in self.M.index_generator(1)]
190             self.transit_cyl_to_cart = self.cyl_ch.transition_map(self.cart_ch_U, [cyls[0]*prod([sin(cyls[i]) for i in range(1,n-2)])*cos(cyls[-2]), cyls[0]*prod([sin(cyls[i]) for i in range(1,n-2)])*sin(cyls[-2])]+ [cyls[0]*prod([sin(cyls[i]) for i in range(1,j)])*cos(cyls[j]) for j in range(n-3,1,-1)]+[cyls[0]*cos(cyls[1]),cyls[-1]] ) 
191             gen_transit_list_cyl = [ sqrt(sum([ xis[i]**2 for i in range(len(xis)-1)])),]+[atan2(sqrt( sum([xis[i]**2 for i in range(j)])),xis[j]) for j in range(n-2,1,-1)]+[atan2(xis[1],xis[0]),xis[n-1]]
192             self.transit_cyl_to_cart.set_inverse(*gen_transit_list_cyl)
193
194      def equip_metric(self):
195          self.g=self.M.riemannian_metric('g')
196          for i in self.M.index_generator(1):
197              self.g[i[0],i[0]]=1           
198
199      def make_orthon_frames(self,ch):
200          """
201          make_orthon_frames=make_orthon_frames(self,ch)
202          This method creates a change-of-basis matrix for orthonormal coordinates as
203               to_orthonormal
204          and a new orthonormal frame from the input of a (spherical coordinates) chart, ch
205
206          PARAMETERS (INPUTS):
207          ch - <a sagemanifolds chart>
208
209          OUTPUT
210          to_orthonormal - a change-of-basis matrix
211          eoch - new orthonormal frame
212         
213          EXAMPLES of USAGE:
214          R4 = Rn(4)     
215          R4.equip_metric()
216          R4.make_orthon_frames(R4.sph_ch)       
217
218          """
219          try:
220              to_orthonormal = ch.domain().automorphism_field()
221              for i,j in self.M.index_generator(2):
222                  if self.g[ch.frame(),i,j,ch]!= 0:
223                     g_ij = self.g[ch.frame(),i,j,ch]
224                     to_orthonormal[ch.frame(),i,j,ch] = 1/sqrt(g_ij)
225              eoch = ch.frame().new_frame(to_orthonormal,'e')               
226             
227              # cf. https://github.com/sagemanifolds/SageManifolds/blob/master/Worksheets/v0.9/SM_Cartesian_spherical-3D.ipynb for explanation on the change of frame
228              Jacobian_ch_to_cart = ch.domain().change_of_frame(ch.frame(), ch.domain().default_chart().frame() )
229              Jacobian_och_to_ch  = ch.domain().change_of_frame(eoch, ch.frame() )
230
231              ch.domain().set_change_of_frame(ch.domain().default_frame(), eoch, 
232                                                Jacobian_och_to_ch.inverse()*Jacobian_ch_to_cart.inverse())
233
234              ch.domain().set_change_of_frame( eoch, ch.domain().default_frame(),
235                                               Jacobian_ch_to_cart*Jacobian_och_to_ch )
236
237
238              return to_orthonormal, eoch, [Jacobian_ch_to_cart,Jacobian_och_to_ch]
239
240          except AttributeError:
241              print "Equip the manifold with a metric g by doing the method .equip_metric() first!"
242     
243         
244def make_pt(ch):
245    """
246    make_pt = make_pt(ch)
247    INPUT
248    ch = sagemanifold chart
249
250    EXAMPLES of USAGE
251    p = make_pt(R3.cart_ch)
252    """
253    coords = ch[:]
254    farglst = ['p',]+list(coords)
255    p = ch.scalar_field( function(*farglst) )
256    return p
257
258
259def make_u(ch):
260    """
261    make_u = make_u(ch)
262    make_u creates a time-INDEPENDENT velocity vector field
263   
264    INPUT
265    ch = sage manifold chart
266
267    EXAMPLEs of USAGE:
268    R2 = Rd(2)
269    u2 = make_u(R2.X_U)
270   
271    R3 = Rd(3)
272    u3 = make_u(R3.X_U)
273
274    u3[1].expr().diff(t) # 0 ; this demonstrates that this velocity vector is time-INDEPENDENT
275    """                           
276    n_0 = ch.domain().manifold().dim()
277    # ucomplst components list of u                               
278    ucomplst = []
279    for i in ch.domain().manifold().index_generator(1):
280        farglst = ['u'+str(i[0]),] + list(ch[:]) 
281        ucomplst.append( function( *farglst ) )
282    u = ch.domain().vector_field()
283    u[ch.frame(),:,ch] = ucomplst
284    return u
285
286
287def make_ut(ch):
288    """
289    make_ut = make_ut(ch)
290   
291    INPUT
292    ch = sage manifold chart
293
294    EXAMPLEs of USAGE:
295    R2 = Rd(2)
296    ut2 = make_ut(R2.X_U)
297   
298    R3 = Rd(3)
299    ut3 = make_ut(R3.X_U)
300    """                           
301    n_0 = ch.domain().manifold().dim()
302    # ucomplst components list of u                               
303    ucomplst = []
304    for i in ch.domain().manifold().index_generator(1):
305        farglst = ['u'+str(i[0]),] + [t,] + list(ch[:]) 
306        ucomplst.append( function( *farglst ) )
307    u = ch.domain().vector_field()
308    u[ch.frame(),:,ch] = ucomplst
309    return u
310
311
312def make_material_der(u, ch):
313    """
314    make_material_der = make_material_der(u,ch)
315
316    EXAMPLES of USAGE:
317    R3 = Rd(3)
318    u3t = make_ut(R3.X_U)
319    udu = make_material_der(u3t, R3.X_U)
320    """
321    uedcomp = []
322    for ui in u[ch.frame(),:,ch]:
323        uidict = dict( [(ch,ui),])
324        uedcomp.append( u( ch.domain().scalar_field( uidict ) ))
325    X = sum( [ uedcomp[i[0]-1]*ch.frame()[i[0]] for i in ch.domain().manifold().index_generator(1) ] )
326    return X
327
328
329def div(u,g):
330    """
331    div = div(u,g)
332    Return the divergence of vector field u \in \mathfrak{X}(M), given the metric g for the manifold M
333    """
334    uflat = g['_ij']*u['^j']
335    return xder( uflat.hodge_star(g) )
336
337def grad(p,g):
338    """
339    grad = grad(p,g)
340
341    EXAMPLE of USAGE
342    R3 = Rd(3)
343    p = make_pt(R3.M)
344    grad(p,R3.g)
345    """
346    dp = xder( p )
347    gradp = g.inverse()['^ij']*dp['_j']
348    return gradp
349
350def curl(u,g):
351    """
352    curl = curl(u,g)
353    Return the curl of vector field u \in \mathfrak{X}(M), given the metric g for the manifold M
354    """
355    uflat = g['_ij']*u['^j']
356    duflat = xder( uflat )
357    return duflat.hodge_star(g)
358
359def buildrho(ch):
360    """
361    buildrho = buildrho(ch)
362    build a time-dependent $\rho$ the mass density, as a scalar function on a chart of a manifold
363   
364    EXAMPLE of USAGE:
365    R2=Rd(2)
366    rho2=buildrho(R2.X_U)
367    """
368    n_0 = ch.domain().manifold().dim()
369    variables = [t,]+[ch[i] for i in range(1,n_0+1)]
370    rho = ch.domain().scalar_field(function('rho',*variables),name='rho',latex_name=r'\rho' )
371    return rho
372
373
374
375##############################
376## Usage Examples
377##############################
378"""
379R2eg = R2()
380R2eg.transit_sph_to_cart.display()
381R2eg.equip_metric()
382R2eg.g.display(R2eg.sph_ch.frame(),R2eg.sph_ch)
383
384to_orthonormal2, e2, Jacobians2 = R2eg.make_orthon_frames(R2eg.sph_ch)
385to_orthonormal2.display(R2eg.sph_ch.frame(),R2eg.sph_ch)
386e2[1].display( R2eg.sph_ch.frame(), R2eg.sph_ch)
387e2[2].display( R2eg.sph_ch.frame(), R2eg.sph_ch)
388Jacobians2[0].inverse()[:,R2eg.sph_ch]
389Jacobians2[1].inverse()[R2eg.sph_ch.frame(),:,R2eg.sph_ch]
390
391R3eg = R3()
392R3eg.transit_sph_to_cart.display()
393R3eg.transit_cyl_to_cart.display()
394R3eg.equip_metric()
395R3eg.g.display(R3eg.sph_ch.frame(),R3eg.sph_ch)
396R3eg.g.display(R3eg.cyl_ch.frame(),R3eg.cyl_ch)
397
398to_orthonormal3sph, e3sph, Jacobians3sph = R3eg.make_orthon_frames(R3eg.sph_ch)
399to_orthonormal3cyl, e3cyl, Jacobians3cyl = R3eg.make_orthon_frames(R3eg.cyl_ch)
400to_orthonormal3sph.display(R3eg.sph_ch.frame(),R3eg.sph_ch)
401to_orthonormal3cyl.display(R3eg.cyl_ch.frame(),R3eg.cyl_ch)
402
403for i in range(1,3+1): 
404    e3sph[i].display( R3eg.sph_ch.frame(), R3eg.sph_ch )
405for i in range(1,3+1):
406    e3cyl[i].display( R3eg.cyl_ch.frame(), R3eg.cyl_ch )
407Jacobians3sph[0].inverse()[:,R3eg.sph_ch]
408Jacobians3cyl[0].inverse()[:,R3eg.cyl_ch]
409
410R4 = Rn(4)
411R4.transit_sph_to_cart.display()
412R4.transit_cyl_to_cart.display()
413R4.equip_metric()
414R4.g.display(R4.sph_ch.frame(),R4.sph_ch)
415R4.g.display(R4.cyl_ch.frame(),R4.cyl_ch)
416
417to_orthonormal4sph, e4sph, Jacobians4sph = R4.make_orthon_frames(R4.sph_ch)
418to_orthonormal4cyl, e4cyl, Jacobians4cyl = R4.make_orthon_frames(R4.cyl_ch)
419to_orthonormal4sph.display(R4.sph_ch.frame(),R4.sph_ch)
420to_orthonormal4cyl.display(R4.cyl_ch.frame(),R4.cyl_ch)
421
422for i in range(1,4+1): 
423    e4sph[i].display( R4.sph_ch.frame(), R4.sph_ch )
424for i in range(1,4+1):
425    e4cyl[i].display( R4.cyl_ch.frame(), R4.cyl_ch )
426Jacobians4sph[0].inverse()[:,R4.sph_ch]
427Jacobians4cyl[0].inverse()[:,R4.cyl_ch]
428
429"""