Relativistische Umlaufbahnen und Schwarze Löcher Die aufregende Geschichte eines Partikels (Photon), welches mit dem Schwarzen Loch Samba tanzt Wilfried Aufbauend auf der Lagrange DGL und deren numerischer Lösung werden Trajektorien von Körpern sowohl mit der klassischen Newton Theorie als auch mit der ART inklusive Schwarzschild und Kerr Metrik gerechnet. Die Lagrange DGL "L" der Kerr'schen Lösung der ART beschreibt ein spärischens Objekt der Masse M mit einem Drehwinkel-Moment (Spin Moment) J: L = -(1-2*M/r)*diff(t(tau),tau)^2/2-2*a*M/r*diff(t(tau),tau)*diff(phi(tau),tau)+r^2/(2*Delta)*diff(r(tau),tau)^2+(r^2+a^2+2*M*a^2/r)*diff(phi(tau),tau)^2/2; mit: a = J/M; and Delta = r^2-2*M*r+a^2; . Ist a = 0. heißt das, dieses Objekt rotiert nicht. Die Kerr Lösung mutiert zur Schwarzschild Lösung. Für c = unendlich (infinity) und t = tau reduziert sich die Schwarzschildlösung auf die Newton Lösung der klassischen Mechanik. Das heißt als Differenz zwischen der kinetischen (Bewegungs-) und potentiellen (Lage oder Ruhe) Energie (das ist die Aussage der Lagrange Gleichung) Die Gleichung der bewegung erfolgt aus einer Vereinfachung der Euler-Lagrange Beziehung: d/dt; diff(L,`q'`[i]); - diff(L,q[i]); = 0 wobei q eine Koordinate darstellt und q' deren Differential zur Zeit > restart: with(plots): with(plottools): Warning, the name changecoords has been redefined Warning, the name arrow has been redefined Newtonian mechanics Zunächst Definitionen: kinetische und potentielle Energie sowie Lagrange Gleichung > T := 1/2*(diff(r(t),t)^2 + r(t)^2*diff(phi(t),t)^2); 2 2 1 /d \ 1 2 /d \ T := - |-- r(t)| + - r(t) |-- phi(t)| 2 \dt / 2 \dt / > V := -G*M/r(t); G M V := - ---- r(t) > L := T - V; 2 2 1 /d \ 1 2 /d \ G M L := - |-- r(t)| + - r(t) |-- phi(t)| + ---- 2 \dt / 2 \dt / r(t) > L1 := subs({r(t)=var1, diff(r(t),t)=var2, phi(t)=var3, diff(phi(t),t)=var4}, L): Euler Lagrange Beziehung für die Newton-Lagrange Gleichungen > E11 := diff(L1, var4): > E12 := diff(L1, var3): > E13 := subs({var1=r(t),var2=diff(r(t),t),var3=phi(t),var4=diff(phi(t),t)},E11): > Eq14 := E13 = l; 2 /d \ Eq14 := r(t) |-- phi(t)| = l \dt / > E21 := diff(L1, var2): > E22 := diff(L1, var1): > E23 := subs({var1=r(t), var2=diff(r(t),t), var3=phi(t), var4=diff(phi(t),t)}, E21): > E24 := subs({var1=r(t), var2=diff(r(t),t), var3=phi(t), var4=diff(phi(t),t)}, E22): > E25 := diff(E23,t): > Eq26 := E25 - E24 = 0; / 2 \ 2 |d | /d \ G M Eq26 := |--- r(t)| - r(t) |-- phi(t)| + ----- = 0 | 2 | \dt / 2 \dt / r(t) Eq31 und Eq32 sind entkoppelte DGLs aus Eq14 . > Eq31 := isolate(Eq14, diff(phi(t),t)); d l Eq31 := -- phi(t) = ----- dt 2 r(t) > Eq32 := eval(Eq26, Eq31); / 2 \ 2 |d | l G M Eq32 := |--- r(t)| - ----- + ----- = 0 | 2 | 3 2 \dt / r(t) r(t) Die init Werte werden behalten, die DGLs werden numerisch gelöst und es wird ein schön Bildchen gemalt G und M werden normiert zu 1 > G:=1; M:=1; G := 1 M := 1 > Eq41 := r(0) = 26; Eq41 := r(0) = 26 > Eq42 := D(r)(0) = 0; Eq42 := D(r)(0) = 0 > Eq43 := phi(0) = 0; Eq43 := phi(0) = 0 > Eq44 := D(phi)(0) = 0.0071; Eq44 := D(phi)(0) = 0.0071 > En := eval(T + V, {r(t)=rhs(Eq41), diff(r(t),t)=rhs(Eq42), diff(phi(t),t)=rhs(Eq44)}); En := -0.02142295846 > l := eval(lhs(Eq14), {r(t)=rhs(Eq41), diff(phi(t),t)=rhs(Eq44)}); l := 4.7996 Die Exzentrizität wird nicht benötigt, wird aber als Referenz auzfgeführt > epsilon := sqrt(1 + 2*En*l^2/((G*M)^2)); epsilon := 0.1139938402 > ini1 := Eq41, Eq42, Eq43; ini1 := r(0) = 26, D(r)(0) = 0, phi(0) = 0 > Eq51:=dsolve({Eq31, Eq32, ini1}, {r(t), phi(t)}, numeric, output=listprocedure); [ Eq51 := [t = proc(t) ... end proc , phi(t) = proc(t) ... end proc , [ r(t) = proc(t) ... end proc , d ] -- r(t) = proc(t) ... end proc ] dt ] > polarplot([rhs(Eq51(t)[3]), rhs(Eq51(t)[2]), t=0..720], scaling=constrained, legend="Newtonian"); > pnewton := polarplot([rhs(Eq51(t)[3]), rhs(Eq51(t)[2]), t=0..720], scaling=constrained, color=black, legend="Newtonian", linestyle=3, thickness=2): Schwarzschild Metrik > G := 'G': M := 'M': Die Berechnung der Schwarzschild metrik ist nahezu identisch mit der Newton Mechanik Wir nutzen die Lagrange Beziehung: > f := 1/2*(-(c^2 - 2*G*M/r(tau))*diff(t(tau),tau)^2 > + 1/(1 - 2*G*M/(c^2*r(tau)))*diff(r(tau),tau)^2 + r(tau)^2*diff(phi(tau),tau)^2); 2 / d \ 2 |---- r(tau)| 1 / 2 2 G M \ / d \ \dtau / f := - - |c - ------| |---- t(tau)| + ----------------- 2 \ r(tau)/ \dtau / / 2 G M \ 2 |1 - ---------| | 2 | \ c r(tau)/ 2 1 2 / d \ + - r(tau) |---- phi(tau)| 2 \dtau / > f1 := subs({t(tau)=var1, diff(t(tau),tau)=var2, r(tau)=var3, diff(r(tau),tau)=var4, > phi(tau)=var5, diff(phi(tau),tau)=var6}, f): Die Eulker-Lagrange beziehungen werden gerechnet (E11 through Eq34) > E11 := diff(f1, var2): > E12 := diff(f1, var1): > E13 := subs({var1=t(tau), var2=diff(t(tau),tau), var3=r(tau), var5=phi(tau), > var4=diff(r(tau),tau), var6=diff(phi(tau),tau)}, E11): > Eq14 := E13 = -a; / 2 2 G M \ / d \ Eq14 := -|c - ------| |---- t(tau)| = -a \ r(tau)/ \dtau / > E21 := diff(f1, var4): > E22 := diff(f1, var3): > E23 := subs({var1=t(tau), var2=diff(t(tau),tau), var3=r(tau), var5=phi(tau), > var4=diff(r(tau),tau), var6=diff(phi(tau),tau)}, E21): > E24 := subs({var1=t(tau), var2=diff(t(tau),tau), var3=r(tau), var5=phi(tau), > var4=diff(r(tau),tau), var6=diff(phi(tau),tau)}, E22): > E25 := diff(E23, tau): > Eq26 := E25 - E24 = 0; 2 2 d 2 / d \ ----- r(tau) / d \ |---- r(tau)| G M 2 G M |---- t(tau)| \dtau / dtau \dtau / Eq26 := - --------------------------- + ------------- + ------------------ 2 2 G M 2 / 2 G M \ 2 2 1 - --------- r(tau) |1 - ---------| c r(tau) 2 | 2 | c r(tau) \ c r(tau)/ 2 / d \ - r(tau) |---- phi(tau)| = 0 \dtau / > E31 := diff(f1, var6): > E32 := diff(f1, var5): > E33 := subs({var1=t(tau), var2=diff(t(tau),tau), var3=r(tau), var5=phi(tau), > var4=diff(r(tau),tau), var6=diff(phi(tau),tau)}, E31): > Eq34 := E33 = h; 2 / d \ Eq34 := r(tau) |---- phi(tau)| = h \dtau / Die DGL werden in den folgenden 3 DGLs aufgerdröselt > Eq53 := isolate(Eq14, diff(t(tau),tau)); d a Eq53 := ---- t(tau) = - ------------ dtau 2 2 G M -c + ------ r(tau) > Eq54 := isolate(Eq34, diff(phi(tau),tau)); d h Eq54 := ---- phi(tau) = ------- dtau 2 r(tau) > Eq55 := subs({Eq53, Eq54}, Eq26); 2 2 d / d \ ----- r(tau) |---- r(tau)| G M 2 \dtau / dtau Eq55 := - --------------------------- + ------------- 2 2 G M / 2 G M \ 2 2 1 - --------- |1 - ---------| c r(tau) 2 | 2 | c r(tau) \ c r(tau)/ 2 2 G M a h + ----------------------- - ------- = 0 2 3 2 / 2 2 G M \ r(tau) r(tau) |-c + ------| \ r(tau)/ Die Init Werte werden beibehalten, wir lösen wieder numerisch und bedienen uns der Normierung von M =1 G = 1 und c = 1 Dann malen wir ein schön Bildchen > M := 1; G := 1; c := 1; M := 1 G := 1 c := 1 > Eq77 := r(0) = 26; Eq77 := r(0) = 26 > Eq78 := D(r)(0) = 0; Eq78 := D(r)(0) = 0 > Eq79 := phi(0) = 0; Eq79 := phi(0) = 0 > Eq80 := D(phi)(0) = 0.0071; Eq80 := D(phi)(0) = 0.0071 > h := eval(lhs(Eq34), {r(tau)=rhs(Eq77), diff(phi(tau),tau)=rhs(Eq80)}); h := 4.7996 > a := sqrt((1 - 2*M/rhs(Eq77))*(1 + h^2/rhs(Eq77)^2)); a := 0.9770019258 Es wird jetzt das effektive Potential mit seinem Umkehrpunkt visualisiert > Vsq := sqrt((1 - 2*M/r)*(1 + h^2/r^2)): > plot([Vsq, a], r=2..40, 0.973..1.12, legend=["Effectives Potential", "Energie"]); > fsolve(Vsq=a, r, 10..20); 15.46812439 > fsolve(Vsq=a, r, 20..30); 25.99999991 > ini := Eq77, Eq78, Eq79; ini := r(0) = 26, D(r)(0) = 0, phi(0) = 0 > Eq91 := dsolve({Eq54, Eq55, ini}, {r(tau), phi(tau)}, numeric, > output=listprocedure); [ Eq91 := [tau = proc(tau) ... end proc , [ phi(tau) = proc(tau) ... end proc , r(tau) = proc(tau) ... end proc , d ] ---- r(tau) = proc(tau) ... end proc ] dtau ] > psch1 := polarplot([rhs(Eq91(tau)[3]), rhs(Eq91(tau)[2]), tau=0..2300], > scaling=constrained, axesfont=[TIMES, ROMAN, 12], legend="Schwarzschild"): > psch2 := disk([15.47*cos(3.7797), 15.47*sin(3.7797)], 0.5, color=black): > psch3 := disk([15.47*cos(11.3399), 15.47*sin(11.3399)], 0.5, color=black): > psch4 := disk([15.47*cos(18.9008), 15.47*sin(18.9008)], 0.5, color=black): > psch5 := disk([15.47*cos(26.4221), 15.47*sin(26.4221)], 0.5, color=black): > pns := disk([0,0],5.8, color=orange): Nun die ganze Schönheit der Bahnkurvenschar um einen Massekörper > display([psch1, psch2, psch3, psch4, psch5, pns, pnewton]); > pschwarz := polarplot([rhs(Eq91(tau)[3]), rhs(Eq91(tau)[2]), tau=0..1000], > scaling=constrained, color=blue, legend=Schwarzschild): > pschwarz3d := spacecurve([rhs(Eq91(s)[3]), rhs(Eq91(s)[2]), > sqrt(8*(rhs(Eq91(s)[3])-2)), s=0..2500], > coords=cylindrical, color=black, numpoints=400): Kerr Metrik > M := 'M': l := 'l': The same procedure as every year: Das gleiche mit der Kerr Metrik. Das einzige was wir ändern ist die lagrange beziehung > a := J/M; Delta := r(tau)^2 - 2*M*r(tau) + a^2; J a := - M 2 2 J Delta := r(tau) - 2 M r(tau) + -- 2 M > f := 1/2*(-(1 - 2*M/r(tau))*diff(t(tau),tau)^2 > + r(tau)^2/(Delta)*diff(r(tau),tau)^2 + (r(tau)^2 + a^2 + 2*M*a^2/r(tau))*diff(phi(tau),tau)^2-4*a*M/r(tau)*diff(t(tau), tau)*diff(phi(tau), tau)); 2 2 / d \ 2 r(tau) |---- r(tau)| 1 / 2 M \ / d \ \dtau / f := - - |1 - ------| |---- t(tau)| + ----------------------------- 2 \ r(tau)/ \dtau / / 2\ | 2 J | 2 |r(tau) - 2 M r(tau) + --| | 2| \ M / / 2 2 \ 2 1 | 2 J 2 J | / d \ + - |r(tau) + -- + --------| |---- phi(tau)| 2 | 2 M r(tau)| \dtau / \ M / / d \ / d \ 2 J |---- t(tau)| |---- phi(tau)| \dtau / \dtau / - --------------------------------- r(tau) > f1 := subs({t(tau)=var1, diff(t(tau),tau)=var2, r(tau)=var3, diff(r(tau),tau)=var4, > phi(tau)=var5, diff(phi(tau),tau)=var6}, f): Entwicklung der Euler-Lagrange Gleichungen E11 bis E34 > E11 := diff(f1, var2): > E12 := diff(f1, var1): > E13 := subs({var1=t(tau), var2=diff(t(tau),tau), var3=r(tau), var5=phi(tau), > var4=diff(r(tau),tau), var6=diff(phi(tau),tau)}, E11): Da passiert immer die Änderung: erst war es 1, dann -a und jetzt -E > Eq14 := E13 = -E; / d \ 2 J |---- phi(tau)| / 2 M \ / d \ \dtau / Eq14 := -|1 - ------| |---- t(tau)| - ------------------- = -E \ r(tau)/ \dtau / r(tau) > E21 := diff(f1, var4): > E22 := diff(f1, var3): > E23 := subs({var1=t(tau), var2=diff(t(tau),tau), var3=r(tau), var5=phi(tau), > var4=diff(r(tau),tau), var6=diff(phi(tau),tau)}, E21): > E24 := subs({var1=t(tau), var2=diff(t(tau),tau), var3=r(tau), var5=phi(tau), > var4=diff(r(tau),tau), var6=diff(phi(tau),tau)}, E22): > E25 := diff(E23, tau): > Eq26 := E25 - E24 = 0; 2 / d \ r(tau) |---- r(tau)| \dtau / Eq26 := ------------------------- 2 2 J r(tau) - 2 M r(tau) + -- 2 M 1 / 2 / d \ / / d - ---------------------------- |r(tau) |---- r(tau)| |2 r(tau) |---- r(tau 2 \ \dtau / \ \dtau / 2\ | 2 J | |r(tau) - 2 M r(tau) + --| | 2| \ M / / 2 \ 2 | d | 2 r(tau) |----- r(tau)| / d \ | 2 | M |---- t(tau)| \ / d \\\ \dtau / \dtau / )| - 2 M |---- r(tau)||| + ------------------------- + ---------------- / \dtau /// 2 2 2 J r(tau) r(tau) - 2 M r(tau) + -- 2 M 2 2 / d \ r(tau) |---- r(tau)| (2 r(tau) - 2 M) \dtau / + --------------------------------------- 2 / 2\ | 2 J | 2 |r(tau) - 2 M r(tau) + --| | 2| \ M / / 2 \ 2 1 | 2 J | / d \ - - |2 r(tau) - ---------| |---- phi(tau)| 2 | 2| \dtau / \ M r(tau) / / d \ / d \ 2 J |---- t(tau)| |---- phi(tau)| \dtau / \dtau / - --------------------------------- = 0 2 r(tau) > E31 := diff(f1, var6): > E32 := diff(f1, var5): > E33 := subs({var1=t(tau), var2=diff(t(tau),tau), var3=r(tau), var5=phi(tau), > var4=diff(r(tau),tau), var6=diff(phi(tau),tau)}, E31): > Eq34 := E33 = l; / d \ / 2 2 \ 2 J |---- t(tau)| | 2 J 2 J | / d \ \dtau / Eq34 := |r(tau) + -- + --------| |---- phi(tau)| - ----------------- = l | 2 M r(tau)| \dtau / r(tau) \ M / A similar decoupling manipulation (Eq53 through Eq55) > Eq53 := isolate(Eq14, diff(t(tau),tau)); / d \ 2 J |---- phi(tau)| \dtau / -E + ------------------- d r(tau) Eq53 := ---- t(tau) = ------------------------ dtau 2 M -1 + ------ r(tau) > Eq53p := subs(Eq53, Eq34); / 2 2 \ | 2 J 2 J | / d \ Eq53p := |r(tau) + -- + --------| |---- phi(tau)| | 2 M r(tau)| \dtau / \ M / / / d \\ | 2 J |---- phi(tau)|| | \dtau /| 2 J |-E + -------------------| \ r(tau) / - ------------------------------ = l / 2 M \ r(tau) |-1 + ------| \ r(tau)/ > Eq54 := isolate(Eq53p, diff(phi(tau),tau)); 2 2 d -l M (-r(tau) + 2 M) + 2 J M E Eq54 := ---- phi(tau) = ------------------------------------- dtau 3 2 2 3 2 r(tau) M - 2 r(tau) M + J r(tau) > Eq54p := subs(Eq54, Eq53); / 2 2 \ 2 J \-l M (-r(tau) + 2 M) + 2 J M E/ -E + ---------------------------------------------- / 3 2 2 3 2 \ d r(tau) \r(tau) M - 2 r(tau) M + J r(tau)/ Eq54p := ---- t(tau) = --------------------------------------------------- dtau 2 M -1 + ------ r(tau) > Eq55 := subs({Eq53, Eq54}, Eq26); 2 / d \ r(tau) |---- r(tau)| \dtau / Eq55 := ------------------------- 2 2 J r(tau) - 2 M r(tau) + -- 2 M 1 / 2 / d \ / / d - ---------------------------- |r(tau) |---- r(tau)| |2 r(tau) |---- r(tau 2 \ \dtau / \ \dtau / 2\ | 2 J | |r(tau) - 2 M r(tau) + --| | 2| \ M / / 2 \ 2 | d | r(tau) |----- r(tau)| | 2 | \ / d \\\ \dtau / )| - 2 M |---- r(tau)||| + ------------------------- / \dtau /// 2 2 J r(tau) - 2 M r(tau) + -- 2 M 2 / / d \\ | 2 J |---- phi(tau)|| 2 | \dtau /| 2 / d \ M |-E + -------------------| r(tau) |---- r(tau)| (2 r(tau) - 2 M) \ r(tau) / \dtau / + ----------------------------- + --------------------------------------- 2 2 2 / 2 M \ / 2\ r(tau) |-1 + ------| | 2 J | \ r(tau)/ 2 |r(tau) - 2 M r(tau) + --| | 2| \ M / / 2 \ 2 | 2 J | / 2 2 \ |2 r(tau) - ---------| \-l M (-r(tau) + 2 M) + 2 J M E/ | 2| \ M r(tau) / - ---------------------------------------------------------- 2 / 3 2 2 3 2 \ 2 \r(tau) M - 2 r(tau) M + J r(tau)/ / / | | 1 | | - ------------------------------------------------------------- |2 J |-E 2 / 2 M \ / 3 2 2 3 2 \ \ \ r(tau) |-1 + ------| \r(tau) M - 2 r(tau) M + J r(tau)/ \ r(tau)/ / d \\ \ 2 J |---- phi(tau)|| | \dtau /| / 2 2 \| + -------------------| \-l M (-r(tau) + 2 M) + 2 J M E/| = 0 r(tau) / / Lösung der DGL mit den angenommenen init-Bedingungn > M := 1; J := 0.37; M := 1 J := 0.37 > l := 4.7996; E := 0.9772; l := 4.7996 E := 0.9772 > Eq77 := r(0) = 26; Eq77 := r(0) = 26 > Eq78 := D(r)(0) = 0; Eq78 := D(r)(0) = 0 > Eq79 := phi(0) = 0; Eq79 := phi(0) = 0 > Eq80 := D(phi)(0) = eval(rhs(Eq54), {r(tau)=rhs(Eq77)}); Eq80 := D(phi)(0) = 0.007143004388 > ini := Eq77, Eq78, Eq79; ini := r(0) = 26, D(r)(0) = 0, phi(0) = 0 > Eq91 := dsolve({Eq54, Eq55, ini}, {r(tau), phi(tau)}, numeric, > output=listprocedure); [ Eq91 := [tau = proc(tau) ... end proc , [ phi(tau) = proc(tau) ... end proc , r(tau) = proc(tau) ... end proc , d ] ---- r(tau) = proc(tau) ... end proc ] dtau ] > polarplot([rhs(Eq91(tau)[3]), rhs(Eq91(tau)[2]), tau=0..2500], > scaling=constrained, color=red, legend="Kerr direct"); > pkerrd := polarplot([rhs(Eq91(tau)[3]), rhs(Eq91(tau)[2]), tau=0..1000], > scaling=constrained, color=red, legend="Kerr direct"): > M := 1; J := -0.37; M := 1 J := -.37 > l := 4.7996; E := 0.976797; l := 4.7996 E := 0.976797 > Eq77 := r(0) = 26; Eq77 := r(0) = 26 > Eq78 := D(r)(0) = 0; Eq78 := D(r)(0) = 0 > Eq79 := phi(0) = 0; Eq79 := phi(0) = 0 > Eq80 := D(phi)(0) = eval(rhs(Eq54), {r(tau)=rhs(Eq77)}); Eq80 := D(phi)(0) = 0.007053899319 > ini := Eq77, Eq78, Eq79; ini := r(0) = 26, D(r)(0) = 0, phi(0) = 0 > Eq91 := dsolve({Eq54, Eq55, ini}, {r(tau), phi(tau)}, numeric, > output=listprocedure); [ Eq91 := [tau = proc(tau) ... end proc , [ phi(tau) = proc(tau) ... end proc , r(tau) = proc(tau) ... end proc , d ] ---- r(tau) = proc(tau) ... end proc ] dtau ] > polarplot([rhs(Eq91(tau)[3]), rhs(Eq91(tau)[2]), tau=0..2500], > scaling=constrained, legend="Kerr retro"); > pkerro := polarplot([rhs(Eq91(tau)[3]), rhs(Eq91(tau)[2]), tau=0..1000], > scaling=constrained, color=green, legend="Kerr retro", axesfont=[TIMES, ROMAN, 12]): > pkerr1 := disk([16.31*cos(3.63467), 16.31*sin(3.63467)], 0.5, color=black): > pkerr2 := disk([14.53*cos(3.923), 14.53*sin(3.923)], 0.5, color=black): > pkerr3 := disk([16.31*cos(10.9765), 16.31*sin(10.9765)], 0.5, color=black): > pkerr4 := disk([14.53*cos(11.7778), 14.53*sin(11.778)], 0.5, color=black): > display([pnewton,pschwarz,pkerrd,pkerro,pkerr1,pkerr2,psch2,pkerr3,pkerr4,psch3,pns]); Das "Gummibett" Die Gummibett gleichung, siehe: Misner, Thorne, and Wheeler, Gravitation, p. 613. Der Stern wird mit gleichförmig verteilter Masse angenommen > M := 'M': m := rp^3/R^3*M; 3 rp M m := ----- 3 R > zin := int(sqrt(1/(1 - 2*m/rp) -1), rp=2*M..r) + C; / (1/2) \ |/ 1 \ | zin := int||----------- - 1| , rp = (2 M .. r)| + C || 2 | | || 2 rp M | | ||1 - ------- | | || 3 | | \\ R / / > zout := int(sqrt(1/(1 - 2*M/rp) -1), rp=2*M..r); (1/2) (1/2) / M \ zout := 2 2 |-------| (r - 2 M) \r - 2 M/ > Eq201 := eval(zin, r=R) = eval(zout, r=R): > C := solve(Eq201, C): > pin := plot3d([r*cos(phi), r*sin(phi), eval(zin, {M=1, R=5.8})], phi=0..2*Pi, r=0..5.8, grid=[25,6], orientation=[88,135], shading=zhue): > pout := plot3d([r*cos(phi), r*sin(phi), eval(zout, {M=1, R=5.8})], phi=0..2*Pi, r=5.8..27, grid=[25,12]): Die Trajektorie (Bahnkurve) eird in dieser Gummibett Darstellung eingemalt Diesen Federball kann man drehen und wenden und sich die Kurven betrachten. > display([pin, pout, pschwarz3d], style=hidden); 1. Beispiel Ein Partikel mit einem 0 Grad Moment wird von einem Schwarzen Loch "Kerr Typ" eingefangen max. Spin-Winkel-Moment J; (a = MWink. > restart: > a := J/M; Delta := r(tau)^2 - 2*M*r(tau) + a^2; J a := - M 2 2 J Delta := r(tau) - 2 M r(tau) + -- 2 M > f := 1/2*(-(1 - 2*M/r(tau))*diff(t(tau),tau)^2 > + r(tau)^2/(Delta)*diff(r(tau),tau)^2 + (r(tau)^2 + a^2 + 2*M*a^2/r(tau))*diff(phi(tau),tau)^2-4*a*M/r(tau)*diff(t(tau), tau)*diff(phi(tau), tau)); 2 2 / d \ 2 r(tau) |---- r(tau)| 1 / 2 M \ / d \ \dtau / f := - - |1 - ------| |---- t(tau)| + ----------------------------- 2 \ r(tau)/ \dtau / / 2\ | 2 J | 2 |r(tau) - 2 M r(tau) + --| | 2| \ M / / 2 2 \ 2 1 | 2 J 2 J | / d \ + - |r(tau) + -- + --------| |---- phi(tau)| 2 | 2 M r(tau)| \dtau / \ M / / d \ / d \ 2 J |---- t(tau)| |---- phi(tau)| \dtau / \dtau / - --------------------------------- r(tau) > f1 := subs({t(tau)=var1, diff(t(tau),tau)=var2, r(tau)=var3, diff(r(tau),tau)=var4, > phi(tau)=var5, diff(phi(tau),tau)=var6}, f): > E11 := diff(f1, var2): > E12 := diff(f1, var1): > E13 := subs({var1=t(tau), var2=diff(t(tau),tau), var3=r(tau), var5=phi(tau), > var4=diff(r(tau),tau), var6=diff(phi(tau),tau)}, E11): > Eq14 := E13 = -E; / d \ 2 J |---- phi(tau)| / 2 M \ / d \ \dtau / Eq14 := -|1 - ------| |---- t(tau)| - ------------------- = -E \ r(tau)/ \dtau / r(tau) > E21 := diff(f1, var4): > E22 := diff(f1, var3): > E23 := subs({var1=t(tau), var2=diff(t(tau),tau), var3=r(tau), var5=phi(tau), > var4=diff(r(tau),tau), var6=diff(phi(tau),tau)}, E21): > E24 := subs({var1=t(tau), var2=diff(t(tau),tau), var3=r(tau), var5=phi(tau), > var4=diff(r(tau),tau), var6=diff(phi(tau),tau)}, E22): > E25 := diff(E23, tau): > Eq26 := E25 - E24 = 0; 2 / d \ r(tau) |---- r(tau)| \dtau / Eq26 := ------------------------- 2 2 J r(tau) - 2 M r(tau) + -- 2 M 1 / 2 / d \ / / d - ---------------------------- |r(tau) |---- r(tau)| |2 r(tau) |---- r(tau 2 \ \dtau / \ \dtau / 2\ | 2 J | |r(tau) - 2 M r(tau) + --| | 2| \ M / / 2 \ 2 | d | 2 r(tau) |----- r(tau)| / d \ | 2 | M |---- t(tau)| \ / d \\\ \dtau / \dtau / )| - 2 M |---- r(tau)||| + ------------------------- + ---------------- / \dtau /// 2 2 2 J r(tau) r(tau) - 2 M r(tau) + -- 2 M 2 2 / d \ r(tau) |---- r(tau)| (2 r(tau) - 2 M) \dtau / + --------------------------------------- 2 / 2\ | 2 J | 2 |r(tau) - 2 M r(tau) + --| | 2| \ M / / 2 \ 2 1 | 2 J | / d \ - - |2 r(tau) - ---------| |---- phi(tau)| 2 | 2| \dtau / \ M r(tau) / / d \ / d \ 2 J |---- t(tau)| |---- phi(tau)| \dtau / \dtau / - --------------------------------- = 0 2 r(tau) > E31 := diff(f1, var6): > E32 := diff(f1, var5): > E33 := subs({var1=t(tau), var2=diff(t(tau),tau), var3=r(tau), var5=phi(tau), > var4=diff(r(tau),tau), var6=diff(phi(tau),tau)}, E31): > Eq34 := E33 = l; / d \ / 2 2 \ 2 J |---- t(tau)| | 2 J 2 J | / d \ \dtau / Eq34 := |r(tau) + -- + --------| |---- phi(tau)| - ----------------- = l | 2 M r(tau)| \dtau / r(tau) \ M / > Eq53 := isolate(Eq14, diff(t(tau),tau)); / d \ 2 J |---- phi(tau)| \dtau / -E + ------------------- d r(tau) Eq53 := ---- t(tau) = ------------------------ dtau 2 M -1 + ------ r(tau) > Eq53p := subs(Eq53, Eq34); / 2 2 \ | 2 J 2 J | / d \ Eq53p := |r(tau) + -- + --------| |---- phi(tau)| | 2 M r(tau)| \dtau / \ M / / / d \\ | 2 J |---- phi(tau)|| | \dtau /| 2 J |-E + -------------------| \ r(tau) / - ------------------------------ = l / 2 M \ r(tau) |-1 + ------| \ r(tau)/ > Eq54 := isolate(Eq53p, diff(phi(tau),tau)); 2 2 d l M (-r(tau) + 2 M) - 2 J M E Eq54 := ---- phi(tau) = -------------------------------------- dtau 3 2 2 3 2 -r(tau) M + 2 r(tau) M - J r(tau) > Eq54p := subs(Eq54, Eq53); / 2 2 \ 2 J \l M (-r(tau) + 2 M) - 2 J M E/ -E + ----------------------------------------------- / 3 2 2 3 2 \ d r(tau) \-r(tau) M + 2 r(tau) M - J r(tau)/ Eq54p := ---- t(tau) = ---------------------------------------------------- dtau 2 M -1 + ------ r(tau) > Eq55 := subs({Eq53, Eq54}, Eq26); 2 / d \ r(tau) |---- r(tau)| \dtau / Eq55 := ------------------------- 2 2 J r(tau) - 2 M r(tau) + -- 2 M 1 / 2 / d \ / / d - ---------------------------- |r(tau) |---- r(tau)| |2 r(tau) |---- r(tau 2 \ \dtau / \ \dtau / 2\ | 2 J | |r(tau) - 2 M r(tau) + --| | 2| \ M / / 2 \ 2 | d | r(tau) |----- r(tau)| | 2 | \ / d \\\ \dtau / )| - 2 M |---- r(tau)||| + ------------------------- / \dtau /// 2 2 J r(tau) - 2 M r(tau) + -- 2 M 2 / / d \\ | 2 J |---- phi(tau)|| 2 | \dtau /| 2 / d \ M |-E + -------------------| r(tau) |---- r(tau)| (2 r(tau) - 2 M) \ r(tau) / \dtau / + ----------------------------- + --------------------------------------- 2 2 2 / 2 M \ / 2\ r(tau) |-1 + ------| | 2 J | \ r(tau)/ 2 |r(tau) - 2 M r(tau) + --| | 2| \ M / / 2 \ 2 | 2 J | / 2 2 \ |2 r(tau) - ---------| \l M (-r(tau) + 2 M) - 2 J M E/ | 2| \ M r(tau) / - --------------------------------------------------------- 2 / 3 2 2 3 2 \ 2 \-r(tau) M + 2 r(tau) M - J r(tau)/ / / d \\ | 2 J |---- phi(tau)|| | \dtau /| / 2 2 \ 2 J |-E + -------------------| \l M (-r(tau) + 2 M) - 2 J M E/ \ r(tau) / - ---------------------------------------------------------------- = 0 2 / 2 M \ / 3 2 2 3 2 \ r(tau) |-1 + ------| \-r(tau) M + 2 r(tau) M - J r(tau)/ \ r(tau)/ > M := 1; J := 1; M := 1 J := 1 > l := 0; E := 1; l := 0 E := 1 > Eq77 := r(0) = 5; Eq77 := r(0) = 5 > Eq78 := D(r)(0) = -0.64; Eq78 := D(r)(0) = -.64 > Eq79 := phi(0) = 0; Eq79 := phi(0) = 0 > Eq80 := D(phi)(0) = eval(rhs(Eq54), {r(tau)=rhs(Eq77)}); 1 Eq80 := D(phi)(0) = -- 40 > ini := Eq77, Eq78, Eq79; ini := r(0) = 5, D(r)(0) = -.64, phi(0) = 0 > Eq91 := dsolve({Eq54, Eq55, ini}, {r(tau), phi(tau)}, numeric, > output=listprocedure); [ Eq91 := [tau = proc(tau) ... end proc , [ phi(tau) = proc(tau) ... end proc , r(tau) = proc(tau) ... end proc , d ] ---- r(tau) = proc(tau) ... end proc ] dtau ] > with(plots): Warning, the name changecoords has been redefined > polarplot([rhs(Eq91(tau)[3]), rhs(Eq91(tau)[2]), tau=0..4.49], > scaling=constrained, numpoints=2400); 2. Beispile Hier wird das Verhalten eines Partikels nahe des instabilen Orbits demonstriert. Das Partikelchen umkreist den Körper einige male und wird dann rausgeschleudert Klar, das werden wir auch animieren. > restart: > 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); 2 / d \ 2 |-- r(s)| 2 1 / 2 2 G M\ /d \ \ds / 1 2 /d \ f := - - |c - -----| |-- t(s)| + --------------- + - r(s) |-- phi(s)| 2 \ r(s) / \ds / / 2 G M \ 2 \ds / 2 |1 - -------| | 2 | \ c r(s)/ > 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): > E11 := diff(f1, var2): > E12 := diff(f1, var1): > E13 := 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)}, E11): > Eq14 := E13 = -a/c; / 2 2 G M\ /d \ a Eq14 := -|c - -----| |-- t(s)| = - - \ r(s) / \ds / c > E21 := diff(f1, var4): > E22 := diff(f1, var3): > E23 := 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)}, E21): > E24 := 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)}, E22): > E25 := diff(E23, s): > Eq26 := E25 - E24 = 0; 2 2 d 2 / d \ --- r(s) /d \ |-- r(s)| G M 2 G M |-- t(s)| \ds / ds \ds / Eq26 := - ----------------------- + ----------- + -------------- 2 2 G M 2 / 2 G M \ 2 2 1 - ------- r(s) |1 - -------| c r(s) 2 | 2 | c r(s) \ c r(s)/ 2 /d \ - r(s) |-- phi(s)| = 0 \ds / > E31 := diff(f1, var6): > E32 := diff(f1, var5): > E33 := 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)}, E31): > Eq34 := E33 = h/c; 2 /d \ h Eq34 := r(s) |-- phi(s)| = - \ds / c > Eq53 := isolate(Eq14, diff(t(s),s)); d a Eq53 := -- t(s) = - --------------- ds / 2 2 G M\ c |-c + -----| \ r(s) / > Eq54 := isolate(Eq34, diff(phi(s),s)); d h Eq54 := -- phi(s) = ------- ds 2 c r(s) > Eq55 := eval(Eq26, {Eq53, Eq54}); 2 2 d / d \ --- r(s) |-- r(s)| G M 2 2 \ds / ds G M a Eq55 := - ----------------------- + ----------- + ----------------------- 2 2 G M 2 / 2 G M \ 2 2 1 - ------- 2 2 / 2 2 G M\ |1 - -------| c r(s) 2 r(s) c |-c + -----| | 2 | c r(s) \ r(s) / \ c r(s)/ 2 h - -------- = 0 3 2 r(s) c > M := 1; G := 1; c := 1; M := 1 G := 1 c := 1 > Eq77 := r(0) = 3.76777; Eq77 := r(0) = 3.76777 > Eq78 := D(r)(0) = 0; Eq78 := D(r)(0) = 0 > Eq79 := phi(0) = 0; Eq79 := phi(0) = 0 > Eq80 := D(phi)(0) = 0.30290054; Eq80 := D(phi)(0) = 0.30290054 > h := eval(lhs(Eq34), {r(s)=rhs(Eq77), diff(phi(s),s)=rhs(Eq80)}); h := 4.300003560 > a := sqrt((1 - 2*M/rhs(Eq77))*(1 + h^2/rhs(Eq77)^2)); a := 1.039364786 > Vsq := sqrt((1 - 2*M/x)*(1 + h^2/x^2)): > plot([Vsq, a], x=2..25, a-0.1*a..a+0.1*a, legend=["Effectives Potential", "Energie"]); > ini := Eq77, Eq78, Eq79; ini := r(0) = 3.76777, D(r)(0) = 0, phi(0) = 0 > Eq91 := dsolve({Eq54, Eq55, ini}, {r(s), phi(s)}, numeric, output=listprocedure); [ Eq91 := [s = proc(s) ... end proc , phi(s) = proc(s) ... end proc , [ r(s) = proc(s) ... end proc , d ] -- r(s) = proc(s) ... end proc ] ds ] > with(plots): Warning, the name changecoords has been redefined > polarplot([rhs(Eq91(s)[3]), rhs(Eq91(s)[2]), s=0..100], scaling=constrained); > 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..100], coords=cylindrical, color=black, numpoints=400): > display([p1, p2]); > N := 75; nstep := 2; N := 75 nstep := 2 > for i from 0 by nstep to N*nstep do > pl[i] := polarplot([rhs(Eq91(s)[3]), rhs(Eq91(s)[2]), s=i..i+nstep], numpoints=5): > pl3d[i] := spacecurve([rhs(Eq91(s)[3]), rhs(Eq91(s)[2]), sqrt(8*(rhs(Eq91(s)[3])-2)), s=i..i+nstep], coords=cylindrical, numpoints=5, color=black): > plsum[i] := display(seq(pl[j*nstep], j=0..i/nstep)): > pl3dsum[i] := display({p1, seq(pl3d[j*nstep], j=0..i/nstep)}): > end do: > display([seq(plsum[i*nstep], i=0..N)], insequence=true, scaling=constrained); > display([seq(pl3dsum[i*nstep], i=0..N)], insequence=true); 3. Beispiel In einer Veröffentlichung von M. Johnston and R. Ruffini, "Generalized Wilkins effect and selected orbits in a Kerr-Newman geometry," Phys. Rev. D 10, 2324-2329 (1974) stehen einige interessante Dinge. Dieses Beispiel nutzt die Randwertbedingung der Fig. 6 dieses papers. Darin wird ein Partikel Orbit gegenüber einem maximal rotierenden Schwarzen Loch beschrieben. Dieses Beispiel zeigt einige interessante Details, welche die Autoren -sicher auf Grund mangelnder rechenleistung damals!!- noch nicht erkannt haben. > restart: > a := J/M; Delta := r(tau)^2 - 2*M*r(tau) + a^2; J a := - M 2 2 J Delta := r(tau) - 2 M r(tau) + -- 2 M > f := 1/2*(-(1 - 2*M/r(tau))*diff(t(tau),tau)^2 > + r(tau)^2/(Delta)*diff(r(tau),tau)^2 + (r(tau)^2 + a^2 + 2*M*a^2/r(tau))*diff(phi(tau),tau)^2-4*a*M/r(tau)*diff(t(tau), tau)*diff(phi(tau), tau)); 2 2 / d \ 2 r(tau) |---- r(tau)| 1 / 2 M \ / d \ \dtau / f := - - |1 - ------| |---- t(tau)| + ----------------------------- 2 \ r(tau)/ \dtau / / 2\ | 2 J | 2 |r(tau) - 2 M r(tau) + --| | 2| \ M / / 2 2 \ 2 1 | 2 J 2 J | / d \ + - |r(tau) + -- + --------| |---- phi(tau)| 2 | 2 M r(tau)| \dtau / \ M / / d \ / d \ 2 J |---- t(tau)| |---- phi(tau)| \dtau / \dtau / - --------------------------------- r(tau) > f1 := subs({t(tau)=var1, diff(t(tau),tau)=var2, r(tau)=var3, diff(r(tau),tau)=var4, > phi(tau)=var5, diff(phi(tau),tau)=var6}, f): > E11 := diff(f1, var2): > E12 := diff(f1, var1): > E13 := subs({var1=t(tau), var2=diff(t(tau),tau), var3=r(tau), var5=phi(tau), > var4=diff(r(tau),tau), var6=diff(phi(tau),tau)}, E11): > Eq14 := E13 = -E; / d \ 2 J |---- phi(tau)| / 2 M \ / d \ \dtau / Eq14 := -|1 - ------| |---- t(tau)| - ------------------- = -E \ r(tau)/ \dtau / r(tau) > E21 := diff(f1, var4): > E22 := diff(f1, var3): > E23 := subs({var1=t(tau), var2=diff(t(tau),tau), var3=r(tau), var5=phi(tau), > var4=diff(r(tau),tau), var6=diff(phi(tau),tau)}, E21): > E24 := subs({var1=t(tau), var2=diff(t(tau),tau), var3=r(tau), var5=phi(tau), > var4=diff(r(tau),tau), var6=diff(phi(tau),tau)}, E22): > E25 := diff(E23, tau): > Eq26 := E25 - E24 = 0; 2 / d \ r(tau) |---- r(tau)| \dtau / Eq26 := ------------------------- 2 2 J r(tau) - 2 M r(tau) + -- 2 M 1 / 2 / d \ / / d - ---------------------------- |r(tau) |---- r(tau)| |2 r(tau) |---- r(tau 2 \ \dtau / \ \dtau / 2\ | 2 J | |r(tau) - 2 M r(tau) + --| | 2| \ M / / 2 \ 2 | d | 2 r(tau) |----- r(tau)| / d \ | 2 | M |---- t(tau)| \ / d \\\ \dtau / \dtau / )| - 2 M |---- r(tau)||| + ------------------------- + ---------------- / \dtau /// 2 2 2 J r(tau) r(tau) - 2 M r(tau) + -- 2 M 2 2 / d \ r(tau) |---- r(tau)| (2 r(tau) - 2 M) \dtau / + --------------------------------------- 2 / 2\ | 2 J | 2 |r(tau) - 2 M r(tau) + --| | 2| \ M / / 2 \ 2 1 | 2 J | / d \ - - |2 r(tau) - ---------| |---- phi(tau)| 2 | 2| \dtau / \ M r(tau) / / d \ / d \ 2 J |---- t(tau)| |---- phi(tau)| \dtau / \dtau / - --------------------------------- = 0 2 r(tau) > E31 := diff(f1, var6): > E32 := diff(f1, var5): > E33 := subs({var1=t(tau), var2=diff(t(tau),tau), var3=r(tau), var5=phi(tau), > var4=diff(r(tau),tau), var6=diff(phi(tau),tau)}, E31): > Eq34 := E33 = l; / d \ / 2 2 \ 2 J |---- t(tau)| | 2 J 2 J | / d \ \dtau / Eq34 := |r(tau) + -- + --------| |---- phi(tau)| - ----------------- = l | 2 M r(tau)| \dtau / r(tau) \ M / > Eq53 := isolate(Eq14, diff(t(tau),tau)); / d \ 2 J |---- phi(tau)| \dtau / -E + ------------------- d r(tau) Eq53 := ---- t(tau) = ------------------------ dtau 2 M -1 + ------ r(tau) > Eq53p := subs(Eq53, Eq34); / 2 2 \ | 2 J 2 J | / d \ Eq53p := |r(tau) + -- + --------| |---- phi(tau)| | 2 M r(tau)| \dtau / \ M / / / d \\ | 2 J |---- phi(tau)|| | \dtau /| 2 J |-E + -------------------| \ r(tau) / - ------------------------------ = l / 2 M \ r(tau) |-1 + ------| \ r(tau)/ > Eq54 := isolate(Eq53p, diff(phi(tau),tau)); 2 2 d l M (r(tau) - 2 M) + 2 J M E Eq54 := ---- phi(tau) = ------------------------------------- dtau 3 2 2 3 2 r(tau) M - 2 r(tau) M + J r(tau) > Eq54p := subs(Eq54, Eq53); / 2 2 \ 2 J \l M (r(tau) - 2 M) + 2 J M E/ -E + ---------------------------------------------- / _________________ Dis Symmetrie ist der entscheidende Ansatz Dinge zu verstehen: -rot E - dB / (c dt) = (4 pi k ) / c rot B - dE/ / (c dt) = (4 pi j ) / c div B = 4 pi rho_m div E = 4 pi rho_e