examples/relativity.py

   1 """
   2 This example calculates the Ricci tensor from the metric and does this
   3 on the example of Schwarzschild solution.
   4 """
   5 import sys
   6 sys.path.append(".")
   7 sys.path.append("..")
   8
   9 from sympy import exp, Symbol, sin, Rational, Derivative, dsolve, Function, Matrix
  10
  11 def grad(f,X):
  12     a=[]
  13     for x in X:
  14         a.append( f.diff(x) )
  15     return a
  16
  17 def d(m,x):
  18     return grad(m[0,0],x)
  19
  20 class MT(object):
  21     def __init__(self,m):
  22         self.gdd=m
  23         self.guu=m.inv()
  24
  25     def __str__(self):
  26         return "g_dd =\n" + str(self.gdd)
  27
  28     def dd(self,i,j):
  29         return self.gdd[i,j]
  30
  31     def uu(self,i,j):
  32         return self.guu[i,j]
  33
  34 class G(object):
  35     def __init__(self,g,x):
  36         self.g = g
  37         self.x = x
  38
  39     def udd(self,i,k,l):
  40         g=self.g
  41         x=self.x
  42         r=0
  43         for m in [0,1,2,3]:
  44             r+=g.uu(i,m)/2 * (g.dd(m,k).diff(x[l])+g.dd(m,l).diff(x[k]) \
  45                     - g.dd(k,l).diff(x[m]))
  46         return r
  47
  48 class Riemann(object):
  49     def __init__(self,G,x):
  50         self.G = G
  51         self.x = x
  52
  53     def uddd(self,rho,sigma,mu,nu):
  54         G=self.G
  55         x=self.x
  56         r=G.udd(rho,nu,sigma).diff(x[mu])-G.udd(rho,mu,sigma).diff(x[nu])
  57         for lam in [0,1,2,3]:
  58             r+=G.udd(rho,mu,lam)*G.udd(lam,nu,sigma) \
  59                 -G.udd(rho,nu,lam)*G.udd(lam,mu,sigma)
  60         return r
  61
  62 class Ricci(object):
  63     def __init__(self,R,x):
  64         self.R = R
  65         self.x = x
  66         self.g = R.G.g
  67
  68     def dd(self,mu,nu):
  69         R=self.R
  70         x=self.x
  71         r=0
  72         for lam in [0,1,2,3]:
  73             r+=R.uddd(lam,mu,lam,nu)
  74         return r
  75
  76     def ud(self,mu,nu):
  77         r=0
  78         for lam in [0,1,2,3]:
  79             r+=self.g.uu(mu,lam)*self.dd(lam,nu)
  80         return r.expand()
  81
  82 def curvature(Rmn):
  83     return Rmn.ud(0,0)+Rmn.ud(1,1)+Rmn.ud(2,2)+Rmn.ud(3,3)
  84
  85 #class nu(Function):
  86 #    def getname(self):
  87 #        return r"\nu"
  88 #        return r"nu"
  89
  90 #class lam(Function):
  91 #    def getname(self):
  92 #        return r"\lambda"
  93 #        return r"lambda"
  94 nu = Function("nu")
  95 lam = Function("lam")
  96
  97 t=Symbol("t")
  98 r=Symbol("r")
  99 theta=Symbol(r"\theta")
 100 phi=Symbol(r"\phi")
 101
 102 #general, spherically symmetric metric
 103 gdd=Matrix((
 104     (-exp(nu(r)),0,0,0),
 105     (0, exp(lam(r)), 0, 0),
 106     (0, 0, r**2, 0),
 107     (0, 0, 0, r**2*sin(theta)**2)
 108     ))
 109 #spherical - flat
 110 #gdd=Matrix((
 111 #    (-1, 0, 0, 0),
 112 #    (0, 1, 0, 0),
 113 #    (0, 0, r**2, 0),
 114 #    (0, 0, 0, r**2*sin(theta)**2)
 115 #    ))
 116 #polar - flat
 117 #gdd=Matrix((
 118 #    (-1, 0, 0, 0),
 119 #    (0, 1, 0, 0),
 120 #    (0, 0, 1, 0),
 121 #    (0, 0, 0, r**2)
 122 #    ))
 123 #polar - on the sphere, on the north pole
 124 #gdd=Matrix((
 125 #    (-1, 0, 0, 0),
 126 #    (0, 1, 0, 0),
 127 #    (0, 0, r**2*sin(theta)**2, 0),
 128 #    (0, 0, 0, r**2)
 129 #    ))
 130 g=MT(gdd)
 131 X=(t,r,theta,phi)
 132 Gamma=G(g,X)
 133 Rmn=Ricci(Riemann(Gamma,X),X)
 134
 135 def main():
 136
 137     #print g
 138     print "-"*40
 139     print "Christoffel symbols:"
 140     print Gamma.udd(0,1,0)
 141     print Gamma.udd(0,0,1)
 142     print
 143     print Gamma.udd(1,0,0)
 144     print Gamma.udd(1,1,1)
 145     print Gamma.udd(1,2,2)
 146     print Gamma.udd(1,3,3)
 147     print
 148     print Gamma.udd(2,2,1)
 149     print Gamma.udd(2,1,2)
 150     print Gamma.udd(2,3,3)
 151     print
 152     print Gamma.udd(3,2,3)
 153     print Gamma.udd(3,3,2)
 154     print Gamma.udd(3,1,3)
 155     print Gamma.udd(3,3,1)
 156     print "-"*40
 157     print "Ricci tensor:"
 158     print Rmn.dd(0,0)
 159     e =  Rmn.dd(1,1)
 160     print e
 161     print Rmn.dd(2,2)
 162     print Rmn.dd(3,3)
 163     #print
 164     #print "scalar curvature:"
 165     #print curvature(Rmn)
 166     print "-"*40
 167     print "solve the Einstein's equations:"
 168     e = e.subs(nu(r), -lam(r))
 169     l =  dsolve(e, [lam(r)])
 170     print lam(r)," = ",l
 171     metric = gdd.subs(lam(r), l).subs(nu(r),-l)#.combine()
 172     print "metric:"
 173     print metric
 174
 175 if __name__ == "__main__":
 176     main()