Bestimmung von Frequenzänderungen
▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀
periodischer Funktionen aus
▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀
Messungen phasenbekannter
▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀
Zeitpunkte
▀▀▀▀▀▀▀▀▀▀
Rolf Schröder, 1999-06-07
Vielfältig treten in der Natur schwingende bzw. rotierende
Systeme auf, die uns durch die Exaktheit, mit der sie scheinbar
ihre Frequenz beibehalten, überraschen. Denken wir nur an die
Genauigkeit, mit der die Jahreslänge oder gar die Tageslänge
eingehalten wird! Dieses Verhalten hat zum Teil in Form von
Erhaltungssätzen, wie etwa dem der Drehimpulserhaltung, seinen
Niederschlag gefunden.
Aber auch andere Vorgänge laufen mit außerordentlich hoher
Präzision ab, wie uns verschiedenste Uhrenarten demonstrieren:
angefangen mit der Pendeluhr über Uhren mit Unruhe bis zur
Quarzuhr. Die schärfsten uns bekannten Resonanzen aber treten im
Bereich der Atom- und Kernphysik auf.
All die hier beispielhaft aufgeführten Systeme zeigen trotz ihrer
vom Prinzip her diskret definierten Eigenfrequenz kleine, aber
dennoch meßbare Frequenzänderungen. Diese haben ihre Ursache
einerseits in äußeren Störungen, wie z. B. die Störung der
Erdbahn durch Jupiter, die andererseits aber auch auf eine
Änderung der inneren Eigenschaften zurückgehen, wie z. B. die
Veränderungen im Kristallgitter bei Quarzuhren. Ja selbst bei
Atomuhren, bei denen die Genauigkeit letztlich durch die
Heisenbergsche Unschärferelation gegeben ist, wird die theoretisch
mögliche relative Genauigkeit von ca. 3*10e-25 (Wasserstoffmaser)
nicht erreicht, da es im Labor unmöglich ist, äußere Störeinflüsse
auf das Wasserstoffatom völlig auszuschalten. Praktisch erreicht
man heute aber eine relative Genauigkeit von etwa 10e-15.
Die Frequenzen schwingender Systeme sind durch gewisse Eigen-
schaften, wie Geometrie, Massenverteilung, innere Kräfte
(Zustandsgleichung!) usw., bestimmt. Man kann die zeitlichen
Veränderungen solcher Größen durch Frequenzmessungen mit großer
Genauigkeit ermitteln. Noch höhere Genauigkeiten werden erzielt,
wenn statt einfacher Frequenzmessungen die Zeitpunkte des
Phasenablaufs der Schwingung registriert werden.
Diese Methode wollen wir im folgenden ausführlich darstellen. Wir
beschränken uns dabei auf die Messung ungedämpft schwingender
Systeme, deren Schwingungsablauf sich durch periodische Funktionen
beschreiben läßt.
Gegeben sei eine über dem Einheitsintervall periodische Funktion
h, so daß
h(E) = h(E+1) für beliebige reelle E , (1)
mit E als Phase. Wir wollen den Begriff der periodischen Funktion
auch für Funktionen verwenden, die im strengen Sinne nicht die
Bedingung (1) erfüllen. Wir nennen eine zeitabhängige Funktion
f(t) periodisch in der Phase, wenn es eine Funktion h nach (1)
gibt und dazu eine streng monotone stetige Phasenfunktion E(t)
gibt mit dE/dt ≥ k > 0, k eine Konstante, so daß die Funktion f(t)
beschrieben werden kann durch:
f(t) = h(E(t)) . (2)
Haben wir z.B. eine in der Zeit periodische Funktion mit einer
konstanten Frequenz Oo , so gilt für die Phase:
E(t) = Oo∙(t-to) + Eo . (3)
Ist die Frequenz zeitabhängig, O = O(t), so können wir leider
nicht in (3) Oo∙(t-to) durch O(t) ersetzen. Es würde nämlich
bedeuten, daß die Phase E zu einem bestimmten Zeitpunkt t1 nur von
der Frequenz zu diesem Zeitpunkt abhinge. Das kann aber nicht
sein, da Änderungen der Frequenz vor diesem Zeitpunkt ebenfalls
die Phase zur Zeit t1 beeinflussen müssen! Richtig ist, daß
dE = O(t)∙dt , (4)
und daß sich daher E durch Integration über die Zeit ergibt:
t
⌠
E(t) = │ O(t')dt' + Eo (5)
⌡
to
mit der Phase Eo = E(to) als Integrationskonstante zum
(willkürlich vorgegebenen) Zeitpunkt to.
Unsere Aufgabe besteht nun darin, aus N Zeitpunktmessungen ti
(z. B. die Zeiten bestimmter Kurvenpunkte, wie Maxima, Minima
usw.) der Funktion f(t) zu bekannten Phasen Ei die
Frequenzfunktion O(t) zu bestimmen.
Dazu wird O(t) zweckmäßig parametrisiert. Die Art der Parametri-
sierung wird wesentlich bestimmt durch die physikalischen Eigen-
schaften des selbsterregt schwingenden Systems. Wir wollen hier
Systeme behandeln, bei denen mit zwei Arten von Frequenzänderun-
gen zu rechnen ist: Erstens erwarten wir, daß das System langsame
Änderungen (wie Entwicklungseffekte oder Ermüdungserscheinungen
usw.) aufweist, einen Trend. Zweitens vermuten wir periodische
Einflüsse im Sinne einer Frequenzmodulation (wie z. B. die
Modulation der Erdrotation mit der Jahreslänge).
Den Trend beschreiben wir durch ein Polynom bis zum Grad NP, die
Frequenzmodulation beschreiben wir durch eine Fourierreihe bis zur
Ordnung NF. Für O(t) erhalten wir daher:
(n) n
NP Oo∙(t-to)
O(t) = Σ ──────────
n=o n!
(6)
NF
+ Σ Am∙cos(2π∙m(t-to)/Po)+Bm∙sin(2π∙m(t-to)/Po) .
m=1
Als Parameter haben wir die Zeitableitungen des Trends im
Zeitpunkt to von nullter bis NP-ter Ordnung, die Fourierkoeffizi-
enten für die Glieder erster bis NF-ter Ordnung und die Grund-
periode Po der Modulation. Mit (6) erhält man so aus (5):
(n) (n+1)
NP Oo∙(t-to)
E(t) = Σ ────────────── + Eo
n=o (n+1)!
NF
+ Σ Am∙Po/(2π∙m)∙sin(2π∙m(t-to)/Po) (7)
m=1
+ Bm∙Po/(2π∙m)∙(1-cos(2π∙m(t-to)/Po)) .
Setzen wir die N Messungen ti zu den Phasen Ei in (7) ein, so
haben wir ein Gleichungssystem mit N Gleichungen und
M = NP + 3 + 2∙NF Unbekannten. Anzustreben ist eine Lösung nach
der Methode der kleinsten Quadrate, da diese biasfreie Schätzungen
der Parameter mit minimalen Varianzen liefert.
Zur direkten Lösung des Gleichungssystems (7) müßten die Messungen
ti explizit auftreten, und alle Unbekannten müßten linear im
System vorkommen. Leider sind die Meßwerte ti implizit in (7)
vorhanden, und die Unbekannte Po ist nichtlinear in den Gleichun-
gen enthalten. Man könnte das System zu lösen versuchen, indem
man es nach den Meßwerten und den Unbekannten linearisiert, um
durch Iteration die Lösung zu finden. Dann hat man das Problem,
geeignete Anfangswerte zu finden.
Glücklicherweise ist es aber gar nicht nötig, das System (7) nach
den Meßwerten ti zu linearisieren. Wir können nämlich für unsere
Fälle von der allgemein gut erfüllten Annahme ausgehen, daß die
Frequenz O(t) relativ kleine Schwankungen um eine feste mittlere
Frequenz Oo ausführt, die dem Betrage nach eine gewisse Schranke
C, die wesentlich kleiner als eins (und größer Null) ist, nicht
übersteigt:
│O(t)-Oo│
───────── < C << 1 . (8)
Oo
Aufgrund der daher nahezu linearen Beziehung zwischen den ti und
den Ei = E(ti) können statt der eigentlichen Meßwerte ti die
Phasen Ei als Meßwerte angesehen werden, und die ti sind als fest
gegeben anzusehen. Mit Hilfe von (4) können wir aus den
Meßfehlern σ(ti) der ti dann auch einen entsprechenden Fehler
σ(Ei) der Ei angeben:
σ(Ei) = Oo∙σ(ti) . (9)
Auch die Gewichte Pi zur Lösung des Gleichungssystems (7) nach der
Methode der kleinsten Quadrate mit den Phasen Ei als Messungen
sind dann bekannt, sie sind proportional zu 1/σ(ti)², da
2 2
Pi = 1/σ(Ei) = 1/(Oo∙σ(ti)) . (10)
Wäre die Grundperiode Po der Modulation fest gegeben, so könnte
man (7) jetzt direkt lösen. Tatsächlich ist aber Po unbekannt.
Eine bezüglich Po iterative Lösung ist nur möglich, wenn bereits
ein hinreichend guter Näherungswert bekannt ist. Wir diskutieren
daher zunächst andere Lösungswege, und sei es, um wenigsten einen
guten Anfangswert Poo für eine anschließende Iteration nach Po zu
finden. Eine Möglichkeit der Lösung des Gleichungssystems (7)
besteht darin, daß wir die Anzahl der Fourierkoeffizienten soweit
erhöhen, bis die Zahl der insgesamt zu bestimmenden Unbekannten
gleich der Zahl der Gleichungen ist. Im Falle äquidistanter
Zeitpunkte ti hätten wir dann eine klassische Fourieranalyse (hier
jedoch mit zusätzlichen Polynomgliedern). Die Periode Po wäre
dann durch
Po = 2∙NF∙Dt (Dt: Abstand der Zeitpunkte ti) (11)
festgelegt und wir hätten nicht die Qual der Wahl von Po. Dieses
Vorgehen würde sogar für nichtäquidistante ti funktionieren (mit
Dt = (tN-t1)/(N-1) als mittlerem Zeitabstand), solange nur die in
(7) auftretenden Funktionen bezüglich der Verteilung der ti linear
unabhängig blieben.
Dieses Verfahren wäre zweckmäßig, wenn in einer größeren
Datenmenge unbekannte Frequenzen gesucht werden und keine
Signifikanzaussage, insbesondere über die Polynomglieder, benötigt
wird.
Ist man jedoch an einer physikalisch relevanten Aussage der
Polynomglieder interessiert, so müssen die Varianzen der
Ei = E(ti), i=1,N, bekannt sein, so daß mit Hilfe der über (7)
durch Linearisierung nach den Parametern gegebenen Matrix die
Varianzen der Parameter (und damit deren Signifikanzen) bestimmt
werden können.
Sind jedoch die Varianzen der Ei unbekannt, so können diese
dennoch mit Hilfe einer Ausgleichsrechnung nach der Methode der
kleinsten Quadrate geschätzt werden, allerdings muß dazu die Zahl
der zu bestimmenden Parameter beträchtlich verringert werden, da
sonst die Signifikanzaussage (in Gestalt der mittleren Fehler)
wenig verläßlich ist. Das bedeutet aber auch eine beträchtliche
Verringerung der Zahl der Fourierparameter, so daß nicht annähernd
alle über das Meßzeitintervall N*Dt möglichen orthogonalen
Fourierkomponenten zugleich bestimmbar sind. Besonders schwer-
wiegend wird dieses Problem, wenn, wie leider häufig der Fall, die
Zahl der Meßwerte nicht sehr groß ist.
Im allgemeinen werden wir daher anders vorgehen müssen. Wir haben
zwei durch die Länge des Meßzeitraumes gegebene grundsätzliche
Schranken für den Bereich der überhaupt bestimmbaren
Modulationsfrequenzen: Einerseits können schwerlich Frequenzen
ermittelt werden, deren Perioden größer als N*Dt sind, da dann die
Korrelation zwischen Polynom- und Fourierkomponenten gegen eins
strebt. Andererseits besagt das Abtasttheorem, daß Perioden
kleiner als etwa 2*Dt nicht bestimmbar sind. Wir brauchen also
nur nach Perioden Po zu suchen, die die folgende Bedingung (12)
erfüllen:
N∙Dt ≥ Po ≥ 2∙Dt∙NF . (12)
Wir lösen nun unter Vorgabe von Polynomgrad NP und Fourierordnung
NF das Gleichungssystem (7) mit der Grundperiode Po als unabhängig
gegebener Variablen. Wir erhalten so die Schätzung der Varianz
(bzw. des mittleren Fehlers) des Meßwertes mit Einheitsgewicht
als Funktion von Po. Beim absoluten Minimum dieser Funktion haben
wir dann die wahrscheinlichste Lösung gefunden. Ob die Modulation
mit der Grundperiode Po tatsächlich signifikant ist, läßt sich mit
einer anschließenden Ausgleichsrechnung unter Einbeziehung von Po
als Unbekannter aus den Varianzen der Fourierkomponenten
ermitteln.
* * *
An dieser Stelle muß der gewissenhafte Statistiker warnend den
Zeigefinger erheben: Indem für eine quasi dichte Periodenfolge
Ausgleichungen durchgeführt werden, zeigt die damit berechnete
Varianz als Funktion von Po einen um die tatsächliche Varianz
schwankenden Verlauf, und zwar auch im Falle nicht vorhandener
periodischer Glieder in den Daten! D.h., der Erwartungswert
dieser Funktion ist zwar gleich der Varianz der Meßwerte, das
Minimum dieser Funktion ist jedoch systematisch kleiner als die
Varianz der Meßwerte. Daher werden Perioden vorgetäuscht. Der
Elimination solcher nicht signifikanter Perioden werden wir uns
noch später zuwenden.
* * *
Praktisch führen wir für eine Folge äquidistanter Modulations-
grundfrequenzen Oo = 1/Po innerhalb der durch (12) definierten
Grenzen Ausgleichslösungen nach der Methode der kleinsten Quadrate
durch, und entnehmen der Lösung mit kleinster Varianz der Meßwerte
einen Anfangswert Poo für die anschließende iterative Ausgleichs-
lösung nach allen Parametern (für die (7) nur nach Po linearisiert
wurde). Das Ganze muß für eine Reihe von Polynomgraden NP und
Fourierordnungen NF durchgeführt werden. Man wählt diejenige
Lösung aus, für die einerseits die Zahl M der Parameter möglichst
klein ist, und für die andererseits keine signifikante
Verkleinerung der Varianz des Meßwertes mit Einheitsgewicht bei
Vergrößerung von NP bzw. NF mehr zu erzielen ist. Wir können
dann annehmen, das die Frequenz O(t) im Rahmen der Meßfehler nach
(6) mit entsprechenden NP,NF beschrieben wird.
Da wir die Zahl M der Unbekannten immer merklich kleiner als die
Zahl N der Gleichungen haben müssen, muß insbesondere auch NF
relativ klein sein, d.h., die Modulation muß sich mit wenigen
Fouriergliedern darstellen lassen.
Tatsächlich ist das praktisch immer der Fall. Hätte man z.B.
eine rechteckförmige Modulation, so müßte diese in (7) durch eine
unendliche Fourierreihe einer Sägezahnfunktion dargestellt werden.
Die Amplitudenbeträge der Fourierreihe in (7) nehmen mit
wachsendem Index m aber erheblich schneller ab als die in (6), da
(7) durch Integration von (6) erzeugt wurde. Es bedeutet für uns,
daß wir im Rahmen der durch die Meßfehler gegebenen Signifikanz-
grenzen praktisch immer mit wenigen Fouriergliedern auskommen.
Pessimistisch formuliert bedeutet es aber auch, daß sich die
Frequenzmodulation grundsätzlich nur bis auf wenige Fourierglieder
signifikant bestimmen läßt (und damit übrigens auch die zuvor
erwähnte "große" Fourieranalyse unpraktikabel ist).
Solange in den Daten keine Modulation außerhalb der Bedingung (12)
enthalten ist, und solange die Grundkomponente der Fourierreihe
nicht zu stark mit dem Polynom korreliert ist, sind keine Probleme
bei der Reduktion zu erwarten.
Treten jedoch Modulationsperioden in den Daten größer als Po auf,
so können sie durch das Polynom dargestellt werden, insbesondere,
wenn sie nur wenige Fourierkomponenten enthalten. Damit stellen
aber die Polynomglieder nicht mehr den Trend im obigen Sinne
unverfälscht dar, eine Trennung von Modulation und Trend ist dann
prinzipiell unmöglich. Das ist übrigens auch der Fall, wenn
Fourierkomponenten durch die Wahl eines zu hohen Polynomgrades NP
immer besser approximierbar (bzw. mit Polynomkomponenten immer
stärker korreliert) sind.
Treten andererseits Modulationsanteile mit P < 2∙Dt in den Daten
auf, so wirken diese für die Ausgleichung wie eine zusätzliche
Streuung der Meßwerte, solange die Meßwerte ti über die Phasen der
Modulationsanteile statistisch gleichverteilt sind. Dann ist auch
eine systematische Verfälschung des Trends nicht zu erwarten.
-------
Für die Beobachtung weiterer phasenbekannter Zeitpunkte ist es von
Interesse, diese mit möglichst großer Sicherheit im voraus zu
berechnen. Wir wollen also für bekannte Parameter A der
Frequenzfunktion O(t) und der zum Zeitpunkt to bekannten Phase Eo
für eine beliebige Phase E den zugehörigen Zeitpunkt t = t(E) mit
Varianz σ(t) bestimmen. Aufgrund der Eigenschaft (8) können wir
uns eine erste Näherung t1 nach
t1 = to + (E-Eo)/Oo (13)
verschaffen. Sodann linearisieren wir E(t) an der Stelle t1,
E = E(t1) + O(t1)∙(t-t1) + Rest (14)
mit O(t1) = dE/dt(t1), so daß man für die nächst bessere Näherung
t2 erhält:
t2 = t1 + (E-E(t1))/O(t1) . (15)
Wir haben so eine Iterationsvorschrift erhalten, und da wegen (8)
die Ableitung der Funktion E(t) nach der Zeit nahezu konstant ist,
ist die Lösung t nach wenigen Schritten im Rahmen der
Rechengenauigkeit erreicht.
Es bleibt noch die Varianz von t zu bestimmen. Da uns von der
Ausgleichslösung des Gleichungssystems (7) nach der Methode der
kleinsten Quadrate auch die Kovarianzmatrix CA der A und von Eo
bekannt ist, kann man über (15) die Varianz σ(t) der Zeit t
erhalten. Dazu betrachten wir im Ausdruck (15) die Auswirkung
einer Variation DA aller Parameter A (im folgenden denken wir uns
Eo mit in den A enthalten) auf t2. Wir können die daraus
resultierende Variation Dt2 von t2 durch die Matrizengleichung
(16) beschreiben:
Dt2 = DTDA˜͂∙DA . (16)
Dabei ist DTDA˜͂ der transponierte Spaltenvektor der partiellen
Ableitungen von t2 nach den Parametern A, und DA ist der
Spaltenvektor der Variation von A. Dann erhält man nach bekannten
Regeln der Statistik die Kovarianzmatrix Ct von t2 (die in unserem
Fall ein Skalar ist) aus (16) nach:
Ct = DTDA˜͂∙CA∙DTDA . (17)
Da im Konvergenzfall t2 in t übergeht, erhält man für DTDA:
2
DTDA = -1/O(t) ∙DEDA(t) (18)
mit DEDA(t) als Spaltenvektor der partiellen Ableitungen von E(t)
nach den Parametern A zur Zeit t.
Anmerkung:
Hat man bereits to so gewählt, daß E(to) = Eo in der Nähe von Null
liegt, so kann man mit sehr hoher Genauigkeit nach (13) den
Zeitpunkt t(E = 0) angeben: t(E = 0) = -Eo/Oo. Da in der
Umgebung von E = 0 t(E) praktisch nur von den fehlerbehafteten
Größen Eo,Oo nach (13) abhängt, ist zur Berechnung der Varianz
Ct(E = 0) nur die Kovarianzmatrix CEO der beiden Größen Eo,Oo
nötig. Mit (dt/dEo,dt/dOo) = DtDEO˜͂ (von (13)) gilt daher
ebenfalls mit sehr großer Genauigkeit:
Ct(E = 0) = DtDEO˜͂∙CEO∙DTDEO.
Zur Beschreibung der Änderung des physikalischen Systems haben wir
bisher die Frequenzfunktion O(t) benutzt. Für manche Zwecke - zum
Teil auch aus Gewohnheit - ist es wünschenswert, statt der
Frequenz die Periode P(t) zu benutzen. Speziell am Trend- bzw.
Polynomanteil der Periodendarstellung ist man interessiert. Wir
können eine umkehrbar eindeutige Transformation der einander
entsprechenden Parameter von O(t) und P(t) einschließlich
Kovarianzmatrix exakt angeben, d.h., sind von O(t) G Komponenten
der Taylorreihe bekannt, so können umkehrbar eindeutig G
Komponenten der Taylorreihe der Funktion P(t) = 1/O(t) bestimmt
werden; allerdings muß das konstante Glied der Taylorreihe von
O(t) ungleich Null sein.
Wer nicht an rechentechnischen Einzelheiten interessiert ist, der
kann den folgenden Abschnitt auch übergehen.
Wir gehen von der Grundbeziehung zwischen O(t) und P(t) aus:
O(t)∙P(t) = 1 . (19)
Für die n-te Ableitung von (19) erhalten wir für einen Polynomteil
in O(t) mit G Komponenten, d.h. mit nullter bis (G-1)-ter
Zeitableitung:
(n) n (n-i) (i)
[O(t)∙P(t)] = Σ (n÷i) ∙ O(t) ∙ P(t) = 0 für n=1, ..,G-1 .(20)
i=0
( (n÷i) lies: n über i )
Wir definieren rechnerkonform:
(i-1) (i-1)
Pi = Po und Oi = Oo ; i=1,... ,G , (21)
der Index "o" an den Größen P bzw. O bedeutet, daß sie für den
Zeitpunkt to gelten. Die Ausdrücke (19) und (20) lauten damit:
O1 ∙ P1 = 1 und
n (22)
Σ ((n-1)÷(i-1)) ∙ On-i+1 ∙ Pi = 0 für n=2,... ,G .
i=1
Nun können wir (22) rekursiv lösen:
P1 = 1/O1 und
n-1 (23)
Pn = -(1/O1) ∙ Σ ((n-1)÷(i-1)) ∙ On-i+1 ∙ Pi für n=2,... ,G .
i=1
Wir wollen (23) nochmals aufschreiben, jedoch mit Betonung der
funktionellen Verknüpfungen:
P1 = P1(O1) = 1/O1 für i=1 und
(24)
Pi = Pi(O1,... ,Oi;P1,... ,Pi-1) für i>1 .
Zur Berechnung der Kovarianzmatrix der Pi benötigen wir die
totalen Ableitungen der Pi nach den Oj, wobei die totale bzw.
partielle Ableitung mit D bzw. d geschrieben wird:
i-1
DPi/DOj = dPi/dOj + Σ dPi/dPr ∙ dPr/dOj
r=1 (25)
0
für j <= i und 0 sonst und Σ ... = 0 .
r=1
Die Fehlertransformation erfolgt über:
i
DPi = Σ DPi/DOj∙DOj
j=1
i i-1
= Σ { dPi/dOj + Σ (dPi/dPr)∙(dPr/dOj) }∙DOj (26)
j=1 r=1
für i = 1,... ,G ,
in Matrizenschreibweise (man beachte, daß dPr/dOj = 0 für r < j):
DP = T∙DO mit
Tij = DPi/DOj
i-1
= dPi/dOj + Σ (dPi/dPr)∙(dPr/dOj) für j <= i (27)
r=j
i-1
und 0 sonst und Σ ... = 0 .
r=i
T ist eine untere Dreiecksmatrix. Da nun aus (23) folgt:
dPi/dOj = -P1∙((i-1)÷(i-j))∙Pi-j+1 für j=1,... ,i und
(28)
dPi/dPr = -P1∙((i-1)÷(r-1))∙Oi-r+1 für r < i
können wir (27) auch so schreiben:
i-1
Tij = -P1∙((i-1)÷(j-1))∙Pi-j+1 - P1∙Σ ((i-1)÷(r-1))∙Oi-r+1∙Trj
r=j
(29)
0
für j <= i und Σ... = 0 (und ((i-1)÷(i-j))=((i-1)÷(j-1)) ).
0
Mit der Formel (29) läßt sich T rekursiv aufbauen. Da von der
Ausgleichslösung auch die Kovarianzmatrix CO der Oi, i=1,... ,G
bekannt ist, kann nun die Kovarianzmatrix der Pj, j=1,... ,G
berechnet werden:
CP = T∙CO∙T˜͂ . (30)
Man kommt bei geschickter Programmierung mit nur einem Vektor der
Länge G und einer GxG-Matrix aus, um aus dem O-Vektor und der
Kovarianzmatrix CO den P-Vektor und die Kovarianzmatrix CP zu
berechnen. Dabei gehen allerdings die ursprünglichen O, CO
verloren.
* * *
In diesem Abschnitt wollen wir zunächst einen Aspekt bei der
Untersuchung von Periodenänderungen bei δ-Cep-Sternen aufgreifen:
wir wollen uns mit der Ursache des periodischen Gliedes in (6)
befassen, und zwar speziell mit der Entscheidung, ob dieses durch
die Bewegung des Cepheiden als Mitglied eines Doppelsternes
verursacht sein kann.
Zur Abschätzung des Problems nehmen wir eine Kreisbahn an und
werten entsprechend auch nur das erste Glied der Fourierreihe von
(6) aus, mit AO bezeichnen wir die Amplitude (AO = √(A1²+B1²)) in
der Frequenzdarstellung. Im folgenden bedeute der Index 1 stets
Masse, Geschwindigkeit usw. des δ-Cep-Sternes.
Der Dopplereffekt führt für die wahre Bahngeschwindigkeit v1 des
Sternes 1 (M1) relativ zum Schwerpunkt zu:
v1/c = AO∙P/sin(i), (P=1/O), i=Bahnneigung. (31)
Andererseits ist der Bahnradius a1 mit v1 verknüpft:
v1 = 2π∙a1/U, so daß (32)
a1 = AO∙P∙U∙c/(2π∙sin(i)), (33)
mit U als Umlaufszeit des Doppelsterns (gleich Grundperiode Po der
Fourierkomponente in (6)). Mit dem Schwerpunktsatz
a1∙M1 = a2∙M2 (34)
und den Definitionen
a = a1 + a2 (Abstand Stern 1 - Stern 2)
(35)
M = M1 + M2 (Gesamtmasse des Systems)
folgt aber:
a1 = a∙(1-M1/M). (36)
Es sei darauf hingewiesen, daß a1 sich nach (33) aus der
Beobachtung ergibt, es muß jedoch für sin(i) ein plausibler Wert
angenommen werden, da i in unseren Fällen nicht bestimmbar ist.
Wir wollen nun folgende Maßeinheiten verwenden:
Zeiteinheit = siderisches Jahr as = 365.256366 Tage, (37)
Längeneinheit = Astronomische Einheit AU = 1.4959787066E11 m,
Masseneinheit = Sonnenmasse MS = 1.98892E30 kg,
so daß das 3. Keplersche Gesetz in diesen Einheiten lautet:
3
2 a
U = ─── . (38)
M
Im folgenden betrachten wir die Masse M1 als gegeben. Da U und a1
als aus der Beobachtung gegebene Parameter anzusehen sind (33),
können wir aus (36) und (38) eine kubische Gleichung zur
Bestimmung der Gesamtmasse M des Doppelsternes aufstellen:
3 2
(M-M1) - Mmin ∙ M = 0 , mit (39)
3
a1 3
Mmin = ───── = (AO∙P∙c/(2π∙sin(i))) ∙ U . (40)
U²
(Nebenbemerkung: AO∙P ist dimensionslos;
c = 2.99792458E08 m/s; c/2π = 10065.30545 AU/as )
Da stets Mmin > 0 ist, hat (39) nur eine reelle Nullstelle M > 0.
Bis auf sin(i) sind alle Größen in dem Ausdruck Mmin beobachtbar
und daher bekannt. Setzen wir sin(i) = 1, so hat Mmin für
gegebene AO,P,U seinen minimalen Wert. Offenbar erhalten wir dann
aus (39) als Lösung auch das kleinstmögliche M, da für sin(i)
ungleich 1 der Wert von Mmin größer wird und damit auch ein
größeres M als Lösung resultiert.
Zusammenfassend stellen wir fest, daß mit sin(i) = 1 aus (39) eine
untere Schranke für die Gesamtmasse des Doppelsternsystems folgt.
Ist diese merklich größer als die Cepheidenmasse M1, so ist die
Masse nicht mit der Hypothese eines Doppelsternes verträglich, da
die zweite Komponente dann sichtbar wäre, es sei denn, wir hätten
es mit einem unsichtbaren Schwarzen Loch zu tun!
Übrigens können wir sogar von der (unsinnigen) Annahme ausgehen,
daß M1=0 ist. Als Lösung von (39) erhalten wir dann als
Minimalmasse gerade Mmin (40), die kleiner ist als die Lösung M
bei Annahme M1 > 0. Ist also bereits Mmin zu groß für einen
plausiblen Doppelstern, so braucht man (39) nicht mehr zu lösen.
* * *
An dieser Stelle sollen einige Abschätzungen zur Genauigkeit der
Frequenz- bzw. Periodenänderungsmessungen angegeben werden.
Fall 1:
Nach (7) betrachten wir zunächst mit NP = 1 und NF = 0 den Fall:
.
Ei = Oo∙ti + Oo∙ti²/2 + Eo + vi (41)
mit i=1,Nm Beobachtungen; vi ist der scheinbare Fehler. Mit den
Varianzen der Ei: V(Ei) = σ(Ei)², der Einheitsmatrix I und der
Annahme, daß alle Varianzen V(Ei) gleich VE sind, ist die
Kovarianzmatrix der Ei:
CE = VE∙I . (42)
Nach (41) wird der Zusammenhang zwischen dem Meßvektor E und dem .
Parametervektor A =(Oo,Oo,Eo)˜͂ durch die folgende Matrix F
beschrieben:
┌ ┐
│ t1 0.5∙t1² 1 │
│ . . . │
│ . . . │
F = │ . . . │ (43)
│ . . . │
│ . . . │
│ tNm 0.5∙tNm² 1 │ .
└ ┘
Die Kovarianzmatrix CA der Parameter A ist bei einer Lösung nach
der Methode der kleinsten Quadrate gegeben durch
CA = i(F˜͂∙PE∙F) (44)
(Vorzeichen "i" bedeutet "Inverse von"), dabei ist PE = iCE die
Gewichtsmatrix der Ei .
Wir erhalten:
┌ 2 3 ┐
│ Σti 0.5∙Σti Σti │
│ │
│ 3 4 2 │
iCA = 1/VE ∙ │ 0.5∙Σti .25∙Σti 0.5∙Σti │ (45)
│ │
│ 2 │
│ Σti 0.5∙Σti Nm │ .
└ ┘
Wählt man die ti symmetrisch zum Nullpunkt und äquidistant
verteilt, so sind die Summen über ungerade Potenzen der ti Null.
Bezeichnen wir den Beobachtungszeitraum mit T = tNm - t1, und
schreibt man genähert Σti² = Nm∙T²/12, Σ(ti²)² = Nm∙(T²)²/80, so
erhält man für CA :
┌ 2 ┐
│ 12∙T 0 0 │
│ │
4 │ 2 │
CA = VE/(Nm∙T ) ∙ │ 0 720 -30∙T │ (46)
│ │
│ 2 4 │
│ 0 -30∙T 9∙T /4 │ .
└ ┘
Da mit (8) und (9) σ(Ei) = Oo∙σ(ti) ist, erhalten wir schließlich
für die Varianz bzw. Streuung der Frequenzänderung:
.
V(Oo) = 720∙Oo²∙V(ti)/(Nm∙(T²)²) bzw. (47)
.
σ(Oo) = Oo∙σ(ti)/T²∙√(720/Nm) . (48)
. .
Da dPo = -Po²∙dOo erhalten wir auch die Streuung der Perioden-
änderung:
.
σ(Po) = (Po∙σ(ti)/T²)∙√(720/Nm) ; (49)
∙
(genähert: σ(Po) = 27∙(Po∙σ(ti)/T²)/√Nm) . (49a)
Fall 2:
Zur Anwendung der Formel (49) ist es nützlich, aus der
Meßgenauigkeit σ(fi) der Funktionswerte der periodischen Funktion
sowie aus Amplitude Af, Periode Po und der Anzahl Ne von Messungen
die Genauigkeit σ(ti) eines einzelnen Phasenpunktes abzuschätzen.
Dabei sei die Änderung der Periode über eine oder wenige Phasen
(über die gemessen wird) so klein, daß man P(t) = Po setzen darf.
Als periodische Funktion nehmen wir eine einfache Sinusfunktion,
so daß wir die folgenden Meßgleichungen haben:
fi = a1 + a2∙cos(2π∙ti/Po) + a3∙sin(2π∙ti/Po) + vi . (50)
Die Gewichtsmatrix Pf der Messungen fi ist gleich
Pf = (1/V(fi))∙I , (51)
die Bedeutung von V und I ist dieselbe wie im Fall 1 . Ganz
analog zu Fall 1 haben wir die Matrix F :
┌ ┐
│ 1 cos(2π∙t1/Po) sin(2π∙t1/Po) │
│ . . . │
│ . . . │
F = │ . . . │ (52)
│ . . . │
│ . . . │
│ 1 cos(2π∙tNe/Po) sin(2π∙tNe/Po) │ .
└ ┘
Wir nehmen an, daß die Meßorte äquidistant gerade so über den
Meßzeitraum der Länge g*Po (mit g > 0 und ganzzahlig) verteilt
sind, daß
tNe-t1 = g∙Po∙(Ne-1)/Ne mit Ne > 2g . (53)
Die Inverse der Kovarianzmatrix Ca der Parameter a1,a2,a3 lautet
dann:
┌ ┐
│ 1 0 0 │
│ │
iCa = Ne/V(fi) ∙ │ 0 1/2 0 │ (54)
│ │
│ 0 0 1/2 │ ,
└ ┘
und damit erhalten wir für Ca :
┌ ┐
│ 1 0 0 │
│ │
Ca = V(fi)/Ne ∙ │ 0 2 0 │ (55)
│ │
│ 0 0 2 │ .
└ ┘
Die Lösung a2,a3 kann in der Form
a2 = Af∙cos(2π∙te/Po) (56)
a3 = Af∙sin(2π∙te/Po)
geschrieben werden, so daß sich nach (50) auch schreiben läßt:
fi = a1 + Af∙cos(2π∙te/Po-2π∙ti/Po) + vi . (57)
te beschreibt gerade den Zeitpunkt, an dem die Funktion ein
Maximum hat. Wir können ohne Einschränkung der Allgemeinheit aus
(56) die Streuung σ(te) bestimmen, die identisch ist mit der
gesuchten Streuung σ(ti) in (49),(49a). Es folgt aus (56) :
te = Po/2π∙arctg(a3/a2) . (58)
Da a2 und a3 nicht korreliert sind, kann die Streuung σ(te) über
das einfache Fehlerfortpflanzungsgesetz berechnet werden. Man
erhält so mit den aus der Matrix (55) zu entnehmenden Streuungen
σ(a2) = σ(a3) = σ(fi)∙√(2/Ne) (59)
schließlich:
σ(te) = (σ(fi)/Af)∙(Po/π)/√(2Ne) . (60)
Ersetzen wir in (49) σ(ti) durch (60), so erhalten wir:
.
σ(Po) = (σ(fi)/Af)∙(Po/T)²∙(√360/π)/√(Nm∙Ne) ; (61)
.
(genähert: σ(Po) = 6∙(σ(fi)/Af)∙(Po/T)²/√(Nm∙Ne)). (61a)