Opsta Teorija Relativnostich18mws

> restart: # 18.1 > with(tensor): > coords := [t, x, y, z]: g := array(symmetric, sparse, 1..4, 1..4): g[1,1] := -c^2 - 2*Phi(x,y,z): g[2,2] := 1 2*Phi(x,y,z)/c^2: g[3,3] := g[2,2]: g[4,4] := g[2,2]: metric := create([-1,-1], eval(g));

> tensorsGR(coords, metric, contra_metric, det_met, C1, C2, Rm, Rc, R, GEin, C): > rGEin := raise(contra_metric, GEin, 1): > comp_rGEin := get_compts(rGEin): > comp_rGEin[1,1];

> limit(comp_rGEin[1,1]*c^2, c=infinity);

> > restart: # 18.2 > Eq1 := G*M*m/r^2 = m*w^2*r; > Soln1 := solve(Eq1, r);

1

> rsat := Soln1[1]; > G:=6.673e-11; c:=299792458; M:=5.974e24; rearth:=6.38e6;

> Tsat := 12*60*60; Tearth := 24*60*60; > wsat := 2*Pi/Tsat: > rsat := eval(rsat,w=wsat): > rsat := evalf(rsat); > vsat := wsat*rsat: > vsat := evalf(vsat); > vearth := 2*Pi*rearth/Tearth: > vearth := evalf(vearth); > Digits := 20: > tdf := (sqrt((1 - 2*G*M/(rsat*c^2)) - (vsat/c)^2) sqrt((1 - 2*G*M/(rearth*c^2)) - (vearth/c)^2))*Tearth; > > restart: # 18.3 > with(tensor): > coords:=[t, r, theta,phi]: g := array(symmetric, sparse, 1..4, 1..4): g[1,1] := -(c^2 - 2*G*m/r): g[2,2] := -c^2/g[1,1]: g[3,3] := r^2: g[4,4] := r^2*sin(theta)^2: metric := create([-1,-1], eval(g));

> tensorsGR(coords, metric, contra_metric, det_met, C1, C2, Rm, Rc, R, GE, C): > displayGR(Christoffel2, C2);

2

> Eqns := geodesic_eqns(coords, s, C2 );

> > restart: # 18.4 > f := 1/2*(-(c^2 - 2*G*M/r(s))*diff(t(s),s)^2 + 1/(1 2*G*M/(c^2*r(s)))*diff(r(s),s)^2 + r(s)^2*diff(phi(s),s)^2);

3

> f1 := subs({t(s)=var1, diff(t(s),s)=var2, r(s)=var3, diff(r(s),s)=var4, phi(s)=var5, diff(phi(s),s)=var6}, f): > Epr11 := diff(f1, var2): > Epr12 := diff(f1, var1): > Epr13 := subs({var1=t(s), var2=diff(t(s),s), var3=r(s), var5=phi(s), var4=diff(r(s),s), var6=diff(phi(s),s)}, Epr11): > Eq14 := Epr13 = -a/c; > Epr21 := diff(f1, var4): > Epr22 := diff(f1, var3): > Epr23 := subs({var1=t(s), var2=diff(t(s),s), var3=r(s), var5=phi(s), var4=diff(r(s),s), var6=diff(phi(s),s)}, Epr21): > Epr24 := subs({var1=t(s), var2=diff(t(s),s), var3=r(s), var5=phi(s), var4=diff(r(s),s), var6=diff(phi(s),s)}, Epr22): > Epr25 := diff(Epr23, s): > Eq26 := Epr25 - Epr24 = 0;

> Epr31 := diff(f1, var6): > Epr32 := diff(f1, var5): > Epr33 := subs({var1=t(s), var2=diff(t(s),s), var3=r(s), var5=phi(s), var4=diff(r(s),s), var6=diff(phi(s),s)}, Epr31): > Eq34 := Epr33 = h/c; > Eq53 := isolate(Eq14, diff(t(s),s));

> Eq54 := isolate(Eq34, diff(phi(s),s));

4

> Eq55 := eval(Eq26, {Eq53, Eq54});

> M := 1; G := 1; c := 1;

> Eq77 := r(0) = 11; > Eq78 := D(r)(0) = 0; > Eq79 := phi(0) = 0; > Eq80 := D(phi)(0) = 0.0295; > h := eval(lhs(Eq34), {r(s)=rhs(Eq77), diff(phi(s),s)=rhs(Eq80)}); > a := sqrt((1 - 2*M/rhs(Eq77))*(1 + h^2/rhs(Eq77)^2)); > Vsq := sqrt((1 - 2*M/x)*(1 + h^2/x^2)): > plot([Vsq, a], x=2..25, a-0.01*a..a+0.01*a, legend=["effective potential", "energy"]);

5

> ini := Eq77, Eq78, Eq79; > Eq91 := dsolve({Eq54, Eq55, ini}, {r(s), phi(s)}, numeric, output=listprocedure);

> with(plots): > polarplot([rhs(Eq91(s)[3]), rhs(Eq91(s)[2]), s=0..500], scaling=constrained);

6

> p1 := plot3d([x*cos(t), x*sin(t), sqrt(8*(x-2))], x=2..14, t=0..2*Pi, style=hidden): > p2 := spacecurve([rhs(Eq91(s)[3]), rhs(Eq91(s)[2]), sqrt(8*(rhs(Eq91(s)[3])-2)), s=0..500], coords=cylindrical, color=black, numpoints=400): > display([p1, p2]);

7

> animate(polarplot, [[rhs(Eq91(s)[3]), rhs(Eq91(s)[2]), s=0..t]], t=0..500, scaling=constrained);

8

> p3 := animate(spacecurve, [[rhs(Eq91(s)[3]), rhs(Eq91(s)[2]), sqrt(8*(rhs(Eq91(s)[3])-2))], s=0..t, coords=cylindrical, numpoints=400, color=black], t=0..500): > display([p1, p3]);

9

> > restart: # 18.5 > with(tensor): > coords := [t, r, theta, phi]: g := array(symmetric, sparse, 1..4, 1..4): g[1,1] := -c^2: g[2,2] := (a(t))^2/(1 - k*r^2): g[3,3] := (a(t)*r)^2: g[4,4] := (a(t)*r*sin(theta))^2: metric := create([-1,-1], eval(g));

> tensorsGR(coords, metric, contra_metric, det_met, C1, C2, Rm, Rc, R, GEin, C): > displayGR(Einstein, GEin);

10

> rGEin := raise(contra_metric, GEin, 1);

> restart: # 18.6 > Eq1 := diff(a(t),t)^2 = 1/a(t); > Eq2 := dsolve({ Eq1, a(0)=0}, a(t));

> > restart: # 18.7 > Eq1 := diff(a(t),t)^2 + 1 = 1/a(t);

11