Pruefung3.mw

Allgemeine Relativitätstheorie mit dem Computer 

General Theory of Relativity on the Computer 

Vorlesung gehalten an der J.W.Goethe-Universität in Frankfurt am Main (Sommersemester 2016)  

von Dr.phil.nat. Dr.rer.pol. Matthias Hanauske 

Frankfurt am Main 11.04.2016 

 

Erster Vorlesungsteil: Allgemeine Relativitätstheorie mit Maple  

Probe Prüfung 

Probe Prüfung 

Aufgabe 3. Maple Worksheet Vorlesung3.mw und Definition ds^2 aus Vorlesung2.mw verwendet. 

Bewegung eines Probekörpers um ein schwarzes Loch in der Ebene 

Im folgenden wird die Geodätengleichung in vorgegebener Schwarzschild Raumzeit betrachtet. Die Geodätengleichung beschreibt wie sich ein Probekörper (Masse = 0) im Raum bewegt und sagt voraus, dass diese Bewegung sich stehts entlang der kürzesten Kurve, in der durch die Metrik beschriebenen gekrümmten Raumzeit, vollzieht. 

Die Geodätengleichung und das effektive Potential V(r) 

> restart:
with( tensor ):
with(plots):
with(plottools):
 

Definition der kovarianten Raumzeit-Metrik eines schwarzen Lochs der Masse M in Schwarzschildkoordinaten: 

> coord := [t, r, theta, phi]:
g_compts := array(symmetric,sparse, 1..4, 1..4):
g_compts[1,1] := 1-2*M/r: g_compts[2,2] := -1/g_compts[1,1]:
g_compts[3,3] := -r^2:    g_compts[4,4] := -r^2*sin(theta)^2:
g := create( [-1,-1], eval(g_compts));
 

`:=`(g, table([index_char = [-1, -1], compts = Matrix(%id = 42830864)]))
`:=`(g, table([index_char = [-1, -1], compts = Matrix(%id = 42830864)]))
(2.1.1)
 

Berechnung der kontravarianten Metrik und der Christoffel Symbole: 

> ginv := invert( g, 'detg' ):
D1g := d1metric( g, coord ):  
D2g := d2metric( D1g, coord ):
Cf1 := Christoffel1( D1g ):
Cf2:= Christoffel2( ginv, Cf1 ):
 

Berechnung der Geodätengleichung als Funktion des affinen Parameters lambda. Die Geodätengleichung ist ein System gekoppelter Differentialgleichungen:   

> eqns:=geodesic_eqns( coord, lambda, Cf2 );
 

{`+`(diff(diff(t(lambda), lambda), lambda), `-`(`/`(`*`(2, `*`(M, `*`(diff(t(lambda), lambda), `*`(diff(r(lambda), lambda))))), `*`(r, `*`(`+`(`-`(r), `*`(2, `*`(M)))))))) = 0, `+`(diff(diff(phi(lambd...
{`+`(diff(diff(t(lambda), lambda), lambda), `-`(`/`(`*`(2, `*`(M, `*`(diff(t(lambda), lambda), `*`(diff(r(lambda), lambda))))), `*`(r, `*`(`+`(`-`(r), `*`(2, `*`(M)))))))) = 0, `+`(diff(diff(phi(lambd...
{`+`(diff(diff(t(lambda), lambda), lambda), `-`(`/`(`*`(2, `*`(M, `*`(diff(t(lambda), lambda), `*`(diff(r(lambda), lambda))))), `*`(r, `*`(`+`(`-`(r), `*`(2, `*`(M)))))))) = 0, `+`(diff(diff(phi(lambd...
{`+`(diff(diff(t(lambda), lambda), lambda), `-`(`/`(`*`(2, `*`(M, `*`(diff(t(lambda), lambda), `*`(diff(r(lambda), lambda))))), `*`(r, `*`(`+`(`-`(r), `*`(2, `*`(M)))))))) = 0, `+`(diff(diff(phi(lambd...
{`+`(diff(diff(t(lambda), lambda), lambda), `-`(`/`(`*`(2, `*`(M, `*`(diff(t(lambda), lambda), `*`(diff(r(lambda), lambda))))), `*`(r, `*`(`+`(`-`(r), `*`(2, `*`(M)))))))) = 0, `+`(diff(diff(phi(lambd...
(2.1.2)
 

Wir lassen nur ebene Bewegungen zu (theta=Pi/2, dtheta=0) und setzen M=1:Hier wurde das Maple-File verändert 

> setM:=0.8:
eq1:=subs({r=r(lambda),sin(theta)=1,cos(theta)=0,diff(theta(lambda),lambda)=0,theta=Pi/2,M=setM},eqns[1]):
eq2:=subs({r=r(lambda),sin(theta)=1,cos(theta)=0,diff(theta(lambda),lambda)=0,theta=Pi/2,M=setM},eqns[2]):
eq3:=subs({r=r(lambda),sin(theta)=1,cos(theta)=0,diff(theta(lambda),lambda)=0,theta=Pi/2,M=setM},eqns[3]):
eq4:=subs({r=r(lambda),sin(theta)=1,cos(theta)=0,diff(theta(lambda),lambda)=0,theta=Pi/2,M=setM},eqns[4]):
eq1:=subs({r(lambda)(lambda)=r(lambda)},eq1);
eq2:=subs({r(lambda)(lambda)=r(lambda)},eq2);
eq3:=subs({r(lambda)(lambda)=r(lambda)},eq3);
eq4:=subs({r(lambda)(lambda)=r(lambda)},eq4);
 

 

 

 

`+`(diff(diff(t(lambda), lambda), lambda), `-`(`/`(`*`(1.6, `*`(diff(t(lambda), lambda), `*`(diff(r(lambda), lambda)))), `*`(r(lambda), `*`(`+`(`-`(r(lambda)), 1.6)))))) = 0
`+`(diff(diff(phi(lambda), lambda), lambda), `/`(`*`(2, `*`(diff(r(lambda), lambda), `*`(diff(phi(lambda), lambda)))), `*`(r(lambda)))) = 0
0 = 0
`+`(diff(diff(r(lambda), lambda), lambda), `-`(`/`(`*`(.8, `*`(`+`(`-`(r(lambda)), 1.6), `*`(`^`(diff(t(lambda), lambda), 2)))), `*`(`^`(r(lambda), 3)))), `/`(`*`(.8, `*`(`^`(diff(r(lambda), lambda), ...
`+`(diff(diff(r(lambda), lambda), lambda), `-`(`/`(`*`(.8, `*`(`+`(`-`(r(lambda)), 1.6), `*`(`^`(diff(t(lambda), lambda), 2)))), `*`(`^`(r(lambda), 3)))), `/`(`*`(.8, `*`(`^`(diff(r(lambda), lambda), ...
(2.1.3)
 

In Abhängigkeit von den Anfangswerten können unterschiedliche Bahnen der Bewegung entstehen. Wir wählen zunächst als Beispiel die Anfangswerte einer geschlossenen Bahn:  

Zur Zeit t=0 sei der fallende Körper bei einem Radius von r=10=5*(Schwarzschildradius), die Anfangsgeschwindigkeit des Körpers sei in radialer Richtung dr=0 und in phi-Richtung dphi=0.036. Wir beschreiben die Bewegung aus der Sichtweise eines im Unendlichen ruhenden Beobachters. Hier wurde das Maple-File verändert, ds^2 aus Vorlesung2.mw 

Bemerkung: Der Anfangswert dt0 ergibt sich hierbei aus dem infinitesimalen Weglängenelements ds²=1 eines massiven Probekörpers: 

> dx:=create([1], array([dt,dr,dtheta,dphi])):
ds2:=get_compts(prod(dx,lower(g,dx,1),[1,1])):
ds2:=collect( simplify(ds2), [dt,dr,dtheta,dphi]);
Gl1:=subs(sin(theta)=1,cos(theta)=0,dtheta=0,dr=0,ds2=1);
dt=solve(Gl1,dt);
 

 

 

`+`(`-`(`/`(`*`(`+`(`*`(`^`(r, 2)), `-`(`*`(4, `*`(r, `*`(M)))), `*`(4, `*`(`^`(M, 2)))), `*`(`^`(dt, 2))), `*`(r, `*`(`+`(`-`(r), `*`(2, `*`(M))))))), `/`(`*`(`^`(dr, 2), `*`(r)), `*`(`+`(`-`(r), `*`...
`+`(`-`(`/`(`*`(`+`(`*`(`^`(r, 2)), `-`(`*`(4, `*`(r, `*`(M)))), `*`(4, `*`(`^`(M, 2)))), `*`(`^`(dt, 2))), `*`(r, `*`(`+`(`-`(r), `*`(2, `*`(M))))))), `/`(`*`(`^`(dr, 2), `*`(r)), `*`(`+`(`-`(r), `*`...
`+`(`-`(`/`(`*`(`+`(`*`(`^`(r, 2)), `-`(`*`(4, `*`(r, `*`(M)))), `*`(4, `*`(`^`(M, 2)))), `*`(`^`(dt, 2))), `*`(r, `*`(`+`(`-`(r), `*`(2, `*`(M))))))), `-`(`/`(`*`(`+`(`-`(`*`(`^`(r, 4))), `*`(2, `*`(...
dt = (`/`(`*`(`^`(`+`(`-`(`*`(`+`(`-`(r), `*`(2, `*`(M))), `*`(r, `*`(`+`(1, `*`(`^`(dphi, 2), `*`(`^`(r, 2))))))))), `/`(1, 2))), `*`(`+`(`-`(r), `*`(2, `*`(M))))), `+`(`-`(`/`(`*`(`^`(`+`(`-`(`*`(`+... (2.1.4)
 

 

> r0:=7:
t0:=0:
phi0:=0:
theta0:=Pi/2:
dr0:=0:
setl:=3.5;
dphi0:=setl/r0^2:
dt0:=evalf(-sqrt(-(-r0+2*setM)*r0*(1+dphi0^2*r0^2))/(-r0+2*setM));
dtheta0:=0:
 

 

3.5
1.272937693 (2.1.5)
 

In der Literatur wird die Bewegung eines Probekörpers um ein schwarzes Loch mittels eines definierten, effektiven Potentials illustriert: 

1. General relativity : An introduction for physicists by M. P. Hobson, G. P. Efstathiou and A. N. Lasenby  

2. Gravity : An introduction to Einstein's general relativity by James B. Hartle 

3. Allgemeine Relativitätstheorie by Torsten Fließbach  

4. Relativistic hydrodynamics by Luciano Rezzolla and Olindo Zanotti  

Dieses Potential hängt von dem, bei der Bewegung erhaltenen Drehimpuls pro Masse m ab. Die im Zentralfeld möglichen Bewegungen werden mittels zweier erhaltener Größen (l: Drehimpuls pro Masse m und E: Energie pro Masse) charakteriesiert. Die folgende Abbildung zeigt (in der Nomenklatur des 4.Buches (Diagramm links) bzw. des 1. bis 3. Buches (Diagramm rechts) ) das effektive Potential als Funktion des Radius bei den obigen gewählten Anfangswerten:Hier wurde das Maple-File verändert 

> setl:=r0^2*dphi0;
setE:=dt0*evalf((1-2*setM/r0));
VeffFb:=(r,M,l)->-M/r+l^2/(2*r^2)-M*l^2/r^3;
VeffRez:=(r,M,l)->sqrt(1-2*M/r)*sqrt(1+l^2/r^2);
Pot:=plot(VeffRez(r,setM,setl),r=2.1..30,title="Effektives Potential V(r) (Definition nach Referenz 4)"):
Pot1:=plot(VeffFb(r,setM,setl),r=2.1..30,title="Effektives Potential V(r) (Definition nach Referenz 1-3)"):
B:=pointplot({[r0, VeffRez(r0,setM,setl)]}, symbol=circle,symbolsize=20):
B1:=pointplot({[r0, VeffFb(r0,setM,setl)]}, symbol=circle,symbolsize=20):
display(Matrix(1,2,[display(Pot,B),display(Pot1,B1)]));
 

 

 

 

 

3.500000000
.9819805060
proc (r, M, l) options operator, arrow; `+`(`-`(`/`(`*`(M), `*`(r))), `/`(`*`(`/`(1, 2), `*`(`^`(l, 2))), `*`(`^`(r, 2))), `-`(`/`(`*`(M, `*`(`^`(l, 2))), `*`(`^`(r, 3))))) end proc
proc (r, M, l) options operator, arrow; `*`(sqrt(`+`(1, `-`(`/`(`*`(2, `*`(M)), `*`(r))))), `*`(sqrt(`+`(1, `/`(`*`(`^`(l, 2)), `*`(`^`(r, 2))))))) end proc

Plot_2d Plot_2d

 

Numerisches Lösen der Geodätengleichung: 

> Loes:=dsolve({eq1,eq2,eq4,t(0)=t0,r(0)=r0,D(r)(0)=dr0,D(t)(0)=dt0,phi(0)=phi0,D(phi)(0)=dphi0},{r(lambda),t(lambda),phi(lambda)},type=numeric,output=listprocedure);
 

[lambda = proc (lambda) local _res, _dat, _solnproc, _xout, _ndsol, _pars, _i; option `Copyright (c) 2000 by Waterloo Maple Inc. All rights reserved.`; if `<`(1, nargs) then error
[lambda = proc (lambda) local _res, _dat, _solnproc, _xout, _ndsol, _pars, _i; option `Copyright (c) 2000 by Waterloo Maple Inc. All rights reserved.`; if `<`(1, nargs) then error
[lambda = proc (lambda) local _res, _dat, _solnproc, _xout, _ndsol, _pars, _i; option `Copyright (c) 2000 by Waterloo Maple Inc. All rights reserved.`; if `<`(1, nargs) then error
[lambda = proc (lambda) local _res, _dat, _solnproc, _xout, _ndsol, _pars, _i; option `Copyright (c) 2000 by Waterloo Maple Inc. All rights reserved.`; if `<`(1, nargs) then error
[lambda = proc (lambda) local _res, _dat, _solnproc, _xout, _ndsol, _pars, _i; option `Copyright (c) 2000 by Waterloo Maple Inc. All rights reserved.`; if `<`(1, nargs) then error
[lambda = proc (lambda) local _res, _dat, _solnproc, _xout, _ndsol, _pars, _i; option `Copyright (c) 2000 by Waterloo Maple Inc. All rights reserved.`; if `<`(1, nargs) then error
(2.1.6)
 

Grafische Veranschaulichung der Lösung:Hier wurde das Maple-File verändert 

> lend:=1200:
Plot1:=odeplot(Loes,[lambda,phi(lambda)],0..lend,numpoints=200,color=blue,thickness=2,title="Koordinatenzeit t vs affiner Parameter lambda"):
Plot2:=odeplot(Loes,[lambda,r(lambda)],0..lend,numpoints=200,color=blue,thickness=2,title="radius vs affiner Parameter lambda"):
Plot3:=odeplot(Loes,[r(lambda),t(lambda)],0..lend,numpoints=700,color=blue,thickness=2,title="Koordinatenzeit t vs radius"):
display(Matrix(1,3,[Plot1,Plot2,Plot3]));
 

Plot_2d Plot_2d Plot_2d

 

> phi=evalf(rhs(Loes(760)[2])-2*Pi);
 

phi = 2.670920897 (2.1.7)
 

> frames:=300:
lend:=1200:
BH:=display(disk([0,0],2,color=black)):
for i from 0 by 1 to frames do
Koerper[i]:=display(disk([rhs(Loes(i*lend/frames)[4])*cos(rhs(Loes(i*lend/frames)[2])),rhs(Loes(i*lend/frames)[4])*sin(rhs(Loes(i*lend/frames)[2]))],0.4,color=blue)):
Ani[i]:=display({Koerper[i],BH});
od:
 

> display([seq(Ani[i],i=0..frames)],insequence=true,scaling=constrained);
 

Plot_2d
 

Während der Bewegung erhaltenen Größen (l: Drehimpuls pro Masse m, E: Energie pro Masse und Weglängenelement ds²): Hier wurde das Maple-File verändert 

> Plot4:=odeplot(Loes,[lambda,r(lambda)^2*diff(phi(lambda), lambda)],0..lend,numpoints=700,color=blue,thickness=2,title="Dehimpuls l vs affiner Parameter lambda"):
Plot5:=odeplot(Loes,[lambda,(1-2*setM/r(lambda))*diff(t(lambda), lambda)],0..lend,numpoints=700,color=blue,thickness=2,title="Energie vs affiner Parameter lambda"):
Plot6:=odeplot(Loes,[lambda,(1-2*setM/r(lambda))*(diff(t(lambda), lambda))^2 - 1/(1-2*setM/r(lambda))*(diff(r(lambda), lambda))^2 - r(lambda)^2*(diff(phi(lambda), lambda))^2],0..lend,numpoints=700,color=blue,thickness=2,title="Weglängenelement ds² vs affiner Parameter lambda"):
display(Matrix(1,3,[Plot4,Plot5,Plot6]));
 

Plot_2d Plot_2d Plot_2d

 

Animation der Bewegung im effektiven Potential:Hier wurde das Maple-File verändert 

> frames:=300:
lend:=1200:
Pot:=plot(VeffRez(r,setM,setl),r=2.85..40):
for i from 0 by 1 to frames do
Koerper[i]:=pointplot({[rhs(Loes(i*lend/frames)[4]), setE]}, symbol=solidcircle,symbolsize=20,color=blue):
Ani1[i]:=display({Pot,Koerper[i]});
od:
 

> display([seq(Ani1[i],i=0..frames)],insequence=true):
 

> Animat1:=display([seq(Ani[i],i=0..frames)],insequence=true,scaling=constrained):
Animat2:=display([seq(Ani1[i],i=0..frames)],insequence=true,title="Bewegung im effektives Potential V(r) (Def. nach Ref. 4)"):
display(Array([Animat1,Animat2]));
 

Plot_2d Plot_2d

 

>