Zentralmasse

[Orbit]

[Nordlicht]

[quote="Barney"]


Die homogene Scheibe hat dagegen nur eine Rotationssymmetrie. Bei der Integration der Feldstärke fehlt mir deshalb in der Nähe der Scheibe die z-Abhängigkeit. Von der Anschauung her ist diese Abhängigkeit erst bei einer großen Entfernung zur Scheibe zu vernachlässigen und deshalb bin ich bei den Formeln für die Scheibe auch eher skeptisch, ob da das Stabilitätskriterium wirklich erfüllt ist.
mfg[/quote]

Der äußere Anteil bringt an sich nur eine Konstante, die beim Differenzieren
wieder verschwindet. Bleibe natürlich auch selbst skeptisch, vielleicht hat jemand
noch eine zündende Idee.

Nordlicht

[Nordlicht]

Eine Näherungslösung für das Gravitationspotential einer flachen Scheibe wäre in Zylinderkoordinaten sinnvoll, da hier die Symmetrie voll genutzt werden kann.
Für V = f, r = x und z = y wäre die partielle Differentialgleichung dann vom Typ

1/x df/dx + d^2f/dx^2 + dˆ2f/dy^2 = k , (d= partiell)

wobei f(x,y) für x und y gegen unendlich nach Null konvergieren soll.
Die Konstante k = 0 für x > R und y > d, sonst k = 4πGρ.

So einfach die Gleichung aussieht, dürfte selbst ein komfortables Programm damit
überfordert sein. Vielleicht hat jemand die Möglichkeit und kann einen Versuch wagen?

[Barney]

Hallo Nordlicht,

will man die angegebene partielle DGL numerisch auf einem begrenzten Gebiet lösen, betrachtet man den z.B. folgenden Definitionsbereich: 0 < x < 10 R und 0 < y < 10 R. Der Definitionsbereich muss dabei nur quadratisch sein. Die Grenze bei 10 R ist willkürlich, aber einigermaßen vernünftig gewählt. Über diesen Definitionsbereich kann dann ein quadratisches Gitter gemäß:
$x=j \cdot \frac{10 R}{N}$ mit j = 0,1, ... , N
$y = l \cdot \frac{10R}{N}$ mit l = 0,1, ... , N
gelegt werden. Die Funktion f und auch die Konstante k wird damit zu einer (N+1)x(N+1)-Matrix und man hat dann die unbekannten Werte $f_{j,l}$, sowie die numerisch bekannten Werte $k_{j,l}$.

Die unbekannten $f_{j,l}$ können über das folgende Gleichungssystem berechnet werden:
$f_{j+1,l}\quad + \quad f_{j-1,l}\quad + \quad f_{j,l+1}\quad + \quad f_{j,l-1}\quad - \quad 4 f_{j,l}\quad + \quad \frac{1}{2j}\left(f_{j+1,l}\quad - \quad f_{j-1,l} \right) = \frac{100 R^2}{N^2} k_{j,l} $ Gleichung 1
mit j = 1,2,...,N-1 und l = 1,2,...,N-1. Kennt man die Werte von f an den Rändern, also für j=0, j=N, l=0, l=N, so reduziert sich die Anzahl der Gleichungen auf (N-1)^2 mit (N-1)^2 Unbekannten. Für j=N und l=N setzt man $f(x,y) = -\frac{GM}{\sqrt{x^2+y^2}}$, mit M: Gesamtmasse der Scheibe.

Aufgrund der Achsensymmetrie und der Stetigkeit der ersten partiellen Ableitungen an den Koordinatenachsen, gilt zusätzlich $f_{0,l} = f_{1,l}$ für alle l = 0,1,...,N-1 bzw $f_{j,0} = f_{j,1}$ für alle j = 0,1,...,N-1. Damit hat man N^2 Gleichungen für N^2 Unbekannte. Zuletzt wird noch der Ausdruck $\frac{100 R^2}{N^2} k_{j,l}$ zu $\frac{200}{N^2}\frac{GM}{d}$, bzw. 0 (je nach j und l) umgeformt, damit alle physikalischen Einheiten herausfallen.

[Barney]

Numerische Lösung von Gleichung 1:

f_j,l ist, wie gesagt, eine NxN Matrix. Werden jedoch sämtliche Werte von f als Kette hintereinander geschrieben ergibt sich ein Vektor mit N^2 Komponenten und Gleichung 1 läßt sich sehr kompakt als Matrixgleichung anschreiben:
Af = k Gleichung 2

Die beiden Indizes j und l werden dabei über die Gleichung i=Nj+l zu dem neuen Index i zusammengefaßt. Dieser Index durchläuft dann die Werte 0,1, ...., N^2-1. Für solche Matrixgleichungen gibt es effektive Lösungsalgorithmen, wie LU-Zerlegung oder das Gauss´sche Eliminationsverfahren, allerdings ergibt sich die Anzahl der Zahlenwerte der Matrix A zu N^4. Aus diesem Grunde kann man hier über die Invertierung der Matrix A nur relativ grobe Näherungen mit beispielsweise N=50 ausrechnen. Für N=50 ergibt sich bereits ein Speicherbedarf von 48 MB.
MfG

[Barney]

Eine numerisch weniger anspruchsvolle Möglichkeit zur Berechnung von \(\Phi\) ergibt sich über die Verwendung von Kugelflächenfunktionen. Aufgrund der Zylindersymmetrie und Endlichkeit des gesuchten Potentials kann man o.B.d.A von folgendendem Ansatz ausgehen:
\(\Phi = \sum_{l=0}^{\infty}R_l(r)\cdot Y_l^0\) mit l=0,1,2,...
Aufgrund der Zylindersymmetrie ist der zweite Parameter m der Kugelflächenfunktionen immer gleich Null und wird deswegen im folgenden auch nicht mehr angeschrieben. Nach einigen kurzen Umformungen schreibt sich die Poisson-Gleichung unseres Problems dann zu:
\(\sum_{l=0}^{\infty}Y_l\frac{1}{r}\frac{\partial ^2}{\partial r^2}rR_l - \sum_{l=0}^{\infty}\frac{l(l+1)}{r^2}Y_lR_l = k(r,\vartheta)\)
Aufgrund der Vollständigkeit der Kugelflächenfunktionen bezüglich l, kann diese Gleichung zu
$\frac{d^2 R_l}{dr^2} \quad + \quad \frac{2}{r}\frac{dR_l}{dr}\quad - \quad \frac{l(l+1)}{r^2}R_l \quad = \quad \int_0^{\pi}k(r,\vartheta)\sin(\vartheta)Y_l d\vartheta$
umgeschrieben werden. Jetzt muss man nur noch das Integral rechts auswerten und bekommt damit eine gewöhnliche DGL zweiter Ordnung mit dem Parameter l und der Variablen r. Diese gewöhnliche DGL kann numerisch relativ leicht gelöst werden (z.B. mit Runge-Kutta) und man hat dann sofort entsprechende Näherungen für \(\Phi\).

Die Berechnung des Integrals ist im Falle der homogenen Scheibe wiederum eine ganz interessante Aufgabe. Eine Zusammenfassung der Lösung schreibe ich in den nächsten Beitrag.

[Orbit]

Barney
Ich kann Deinen Ausführungen leider nicht folgen. Und so wird es wohl auch dem Thread-Eröffner 'Nordlicht' gehen. Von den zur Zeit noch mit diskutierenden Usern, ist 'zeitgenosse' wohl der einzige, der noch mit steigt, und den scheint das Thema nicht sonderlich zu interessieren. Deshalb führst Du hier einen Monolog.

Orbit

[Barney]

[Barney]

hier also noch die Berechnung des Integrals und Fortsetzung zur Lösung der Poissongleichung für die homogene Scheibe:

Um das Integral \(\int_0^{\pi}k(r,\vartheta)\sin(\vartheta)Y_l d\vartheta\) zu berechnen, macht man sich zuerst klar, dass r lediglich das Verhalten der Funktion k steuert. Die Funktion k nimmt in Abhängigkeit von r und \(\vartheta\) nur zwei verschieden Werte an, also entweder Null oder die Konstante \(4\pi G\rho\). Bei konstantem r und variablem \(\vartheta\) enstspricht der Definitionsbereich (Input) von k gerade einem Halbkreis mit Radius r. Dieser Halbkreis liegt mit steigendem r entweder voll innerhalb der homogenen Scheibe, schneidet diese in zwei bzw. vier Punkten oder liegt komplett außerhalb der Scheibe. Genau mit Hilfe dieser Veranschaulichung läßt sich auch das Integral berechnen. Dazu ersetze ich zuerst einmal die speziellen Kugelflächenfunktionen Y_l (m=0) durch die normalen Legendrepolynome gemäß:

\[ \int_0^{\pi}k(r,\vartheta)\sin(\vartheta)Y_l d\vartheta
= \int_0^{\pi}k_1(r,\vartheta)\sin(\vartheta)P_l(\cos(\vartheta)) d\vartheta \]
mit
\[ k_1(r,\vartheta) := \sqrt{\frac{2l+1}{4\pi}}k(r,\vartheta). \]

Aufgrund der Symmetrieeigenschaften der Legendre-Polynome ergibt sich zuerst $ \int_0^{\pi}k_1(r,\vartheta)\sin(\vartheta)P_l d\vartheta = 0 $ falls l=1,3,5,7,... also l ungerade. Für ungerade l kann man deswegen \(R_l\) unmittelbar angeben. Dazu erinnere man sich, dass man zu jeder Lösung der Poissongleichung auch eine Lösung der Laplace-Gleichung (=homogene Poisson-Gleichung) addieren darf. In unserem Fall ergibt jedoch das Grenzwertverhalten \(\lim_{r\rightarrow \infty} \Phi = 0\), dass dieser zusätzliche Anteil der homogenen Gleichung insgesamt gleich Null ist und deswegen gilt \(R_l = 0\) falls l ungerade. Mit der Substitution \(x = \cos(\vartheta)\), einer kleinen Skizze und den trigonometrischen Funktionen (sin(x), cos(x)) erhält man für gerade l das folgende Ergebnis:

\[ \int_0^{\pi}k_1(r,\vartheta)\sin(\vartheta)P_l d\vartheta = \left\{ \begin{array}{ll} 2k_1\int_{0}^{1}P_l(x)dx & \mbox{falls}\quad r \leq d/2 \\ 2k_1\int_{0}^{\frac{d}{2r}}P_l(x)dx & \mbox{falls}\quad d/2 < r \leq R \\ 2k_1 \int_{\sqrt{1-\frac{R^2}{r^2}}}^{\frac{d}{2r}}P_l(x)dx & \mbox{falls}\quad R < r < \sqrt{R^2 + d^2/4} \\ 0 & \mbox{falls}\quad r \geq \sqrt{R^2 + d^2/4}\end{array}\right. \]
mit
\[ k_1 := \sqrt{(2l+1)4\pi}G\rho. \]

Somit bekommt man sehr gute Näherungslösungen, wenn man die \(R_l\) numerisch für \(l=0,2,4, \ldots ,l_{max}\) berechnet und \(l_{max}\) muss dabei nicht allzu groß gewählt werden. Der Koeffizient vor den Legendre-Polynomen $\frac{1}{2^l l!}$ zeigt, dass man bei l=14 abbrechen kann. Die rechte Seite der gewöhnlichen Differentialgleichung für \(R_l\) ist dann immer kleiner als 1e-15. Also kann man \(R_{14} = 0\), \(R_{16} = 0\) usw. setzen.

Wesentlich für die Bestimmung von \(\Phi\) sind also nur die Funktionen \(R_0\), \(R_2\) bis \(R_{12}\)! Die numerische Berechnung dieser Funktionen ist deutlich weniger aufwändig als das Invertieren einer 2500x2500 Matrix (Ansatz mit 50x50 Gitterpunkten) und das Ergebnis dürfte auch um einiges genauer sein.
MfG

[Barney]

NACHTRAG: Die Abschätzung, ab wann die Berechnung der R_l abgebrochen werden kann, muss ich mir noch genauer ansehen. Teilnehmer babau hat mich darauf aufmerksam gemacht, dass es in der Wikipedia eine geschlossene Formel für die Legendre-Polynome gibt. Damit kann man die rechte Seite der DGL (mit Parameter l) so auswerten, dass man \(l_{max}\) voraussichtlich gut ermitteln kann.
MfG

[Barney]

NACHTRAG1: Mit dem oben gesagten ergibt sich noch eine interessante Beziehung:
$\int_0^1 P_l(x)dx = 0 \quad \mbox{falls l gerade und}\quad l \neq 0$.
Den Beweis dazu sehe ich momentan nicht, aber die Beziehung gilt bis einschließlich l=12.
Der dritte Abschnitt von oben, rechts neben der geschweiften Klammer (s. oben) kann in einer ersten Rechnung ebenfalls gleich Null gesetzt werden. Man erhält dann das Potential einer homogenen, symmetrischen Kugelschicht.

[babau]

[Barney]

@ babau,

vielen Dank für den eleganten Beweis.

@ All:

Das Konvergenzverhalten des zweiten Integrals auf der rechten Seite der gewöhnlichen DGL ist für große l leider nicht ganz so gut wie oben angenommen. Für l=12 und $\frac{d}{2r} \in (0,1)$ liegt das Maximum bei \(\approx \pm \frac{GM}{R^2d} \cdot 0.2\). Für eine brauchbare Genauigkeit muss die gewöhnliche DGL also auch für l > 12 gelöst werden. Die gewöhnliche DGL läßt sich übrigens gemäß T. Fließbach, "Elektrodynamik" durch die neuen Funktionen $A_l := \frac{R_l}{r}$ weiter vereinfachen. Durch Einsetzen kann man die folgende Identität beweisen:
$\frac{1}{r}\frac{d^2A_l}{dr^2}=\frac{d^2R_l}{dr^2}+\frac{2}{r}\frac{dR_l}{dr}$.
Damit vereinfacht sich die linke Seite der gewöhnlichen DGL wesentlich und kann für l=0 analytisch gelöst werden. Für l=2,4,6,... kann die DGL abschnittsweise gelöst werden. Die rechte Seite der DGL ist für r kleiner d/2 und r > R (l=2,4,6,...) gleich Null. Die allgemeine Lösung für \(A_l\) kann damit aus dem Buch von T. Fließbach übernommen werden:
$A_l(r) = a_lr^{l+1}+\frac{b_l}{r^l}$.
Damit das Potential keine mathematischen Singularitäten enthält (was aus physikalischen Gründen ausgeschlossen werden kann) lautet die Lösung für r kleiner d/2 $A_l = a_l r^{l+1}$ und für r > R, $A_l = \frac{b_l}{r^l}$. Die beiden Integrationskonstanten \(a_l\) und \(b_l\) können über die Bedingung nach Stetigkeit und Differenzierbarkeit der \(A_l\) bei r=d/2 und r = R zusammen mit einer numerischen Lösung im Bereich d/2 < r < R ebenfalls numerisch berechnet werden. Als erstes Ergebnis zeigt sich nun, dass \(R_l(0)=0\) gilt, falls l=2,4,6,.... Dies hat wiederum zur Folge, dass das Potential bei r=0 nur von \(R_0(0)\) festgelegt wird. Da sich \(R_0\) komplett analytisch berechnen läßt, ergibt sich demnach

$\Phi (0) = -\frac{GM}{24\pi R^3}\left( 12 R^2-6Rd+d^2\right)$ für eine homogene, symmetrische Kugelschicht mit Radius R und Dicke d.

Damit hat man eine gute Vergleichsmöglichkeit für eine komplett numerische Lösung und in diesem Sinne betrachte ich die Berechnung des Potentials vorerst als gelöst an.

[Barney]

Für die weitere Untersuchung der Stabilität homogener Gaswolken kann man sich noch die Navier-Stokes-Gleichung mit externen Kräften (Gravitation) ansehen. Über die zugehörige(n), partielle(n) DGL sollten noch weitere, bzw. präzisere Aussagen über die Dynamik (also die zeitliche Entwicklung) homogener Gaswolken gemacht werden können.
MfG

[Barney]