Finden ganzzahliger Approximationen

Sie können den Lenstra-Lenstra-Lovász (LLL)-Gitterbasisreduktionsalgorithmus ausprobieren , um einen kurzen Vektor in einem Gitter zu berechnen. Zur Berechnung verwende ich Pari GP (S. 382, ​​Abschnitt 3.10.62 qflll).

Wir definieren ein Gitter mit einer Basis, die aus den Spalten der Matrix besteht$M$

$$M:=\begin{pmatrix} \frac 1 m& 0& 0 \\0&\frac 1 m & 0 \\0 & 0& \frac 1 m \\T_d &T_s & 0\\ T_d&0&T_a\end{pmatrix}$$ $m$ist ein Skalierungsfaktor. Eine Linearkombination dieser Basisvektoren

$$\begin{pmatrix} \frac {\lambda_1} m \\ \frac {\lambda_2} m \\ \frac {\lambda_3} m \\{\lambda_1} T_d+ {\lambda_2} T_s \\{\lambda_1} T_d+ {\lambda_3} T_a \end{pmatrix}$$

ist von geringer Länge, wenn${\lambda_1} T_d+ {\lambda_2} T_s$Und${\lambda_1} T_d+ {\lambda_3} T_a$ist klein u$\frac {\lambda_1} m, \frac {\lambda_2} m , \frac {\lambda_3} m$sind auch klein. In diesem Fall

$${\lambda_2} T_s - {\lambda_3} T_a= ({\lambda_1} T_d+ {\lambda_2} T_s )-({\lambda_1} T_d+ {\lambda_3} T_a)$$ ist auch klein. Beachten Sie, dass${\lambda_2}$Und${\lambda_3}$haben das gleiche Vorzeichen.

Bei Pari/GP online bekomme ich folgende Ausgabe

\p20
Td= 27.212220815;
Ts = 29.530588861;
Ta = 27.554549886;
m=1000
M=[1/m,0,0;0,1/m,0;0,0,1/m;Td,Ts,0;Td,0,Ta]
s=qflll(M)
[s[,1],[s[1,1]*Td+s[2,1]*Ts, s[1,1]*Td+s[3,1]*Ta, s[2,1]*Ts-s[3,1]*Ta]]

%2 = [[242, -223, -239]~, [0.036121227000000000000, -0.17998552400000000000, 0.21610675100000000000]]

Die Koeffizienten sind$242, 223, 239$, die Unterschiede sind$0.036, 0.180, 0.216$.

Ich habe die Länge des tropischen Monats hinzugefügt ($T_t$) und der Sternmonat ($Ti$), von Wiki und bekam, für m=100000

? \p20
Td= 27.212220815;
Ts = 29.530588861;
Ta = 27.554549886;
Tt=27.321582252
Ti=27.321661554
m=100000
M=[1/m,0,0,0,0;0,1/m,0,0,0;0,0,1/m,0,0;0,0,0,1/m,0;0,0,0,0,1/m;Td,Ts,0,0,0;Td,0,Ta,0,0;Td,0,0,Tt,0;Td,0,0,0,Ti]
s=qflll(M)
[s[,1],[s[1,1]*Td+s[2,1]*Ts, s[1,1]*Td+s[3,1]*Ta, s[1,1]*Td+s[4,1]*Tt,  s[1,1]*Td+s[5,1]*Ti]]

%27 = [[4993, -4601, -4931, -4973, -4973]~, [0.37917983400000000000, -0.86695857100000000000, 0.38999009900000000000, -0.0043787470000000000000]]

Kehren Sie das Problem um, mit$F_a=1/T_a$usw. Dann suchen wir nach ganzen Zahlen$a,b,c$das befriedigt

$$\frac{F_a}{a}\approx\frac{F_b}{b}\approx\frac{F_c}{c}\approx \delta.$$ Das heißt, wir wollen einen Gitterabstand finden$δ$damit alle$F_{a,b,c}$sind sehr nah an Gitterpunkten. Die Neuinterpretation der Aufgabe ist somit zu $$\text{ minimize }\|\vec F-δ\vec a\|_2^2\text{ where }\vec a\in\Bbb Z^3,$$ mit$\vec a$minimal, d. h. teilerfremd als Vektor.

Dies kann durch eine Art euklidischen Algorithmus erfolgen, ähnlich wie die Kettenbruchberechnung damit zusammenhängt. Im Bereich der Motivationen findet sich auch der L3-Algorithmus.

Iterieren

• Nehmen Sie die beiden größten Elemente und reduzieren Sie das größte durch das zweitgrößte, indem Sie ganzzahlige Quotienten verwenden.
• Behalten Sie die Transformationsmatrix im Auge$M$, hier mit der Konvention, dass wenn$\vec F'$ist dann irgendwann der reduzierte Vektor$M^T\vec F'=\vec F$gleich dem ursprünglichen Vektor ist.
• Dann, wenn das (neue) maximale Element von$\vec F'$an Stelle$k$deutlich größer ist als die übrigen Elemente, der Vektor$\vec a=M_k$, Die$k$Reihe von$M$ein Kandidat für den ganzzahligen Koeffizientenvektor ist.

Verwenden Sie als Fehlermaß das Residuum der Minimierungsaufgabe oder ein eng verwandtes Maß dafür, wie parallel$\vec F$Und$\vec a$Sind.

Im folgenden Python-Skript werden diese Ideen umgesetzt. Für den Fall, dass sich das Residuum wesentlich verbessert, wird ein Bericht gedruckt.

T = [ 27.212220815, 29.530588861, 27.554549886 ]

def largest(F):
    k0=0; k1=1
    if F[0]<F[1]: k0=1; k1=0
    for k in range(2,len(F)):    
        if F[k]>F[k0]: k1=k0; k0=k; continue
        if F[k]>F[k1]: k1=k    
    return k0,k1

def defect(a):
    n2_F=n2_a=0
    s_Fa=0
    for ak, Tk in zip(a,T): n2_F+=1/Tk**2; n2_a+=ak**2; s_Fa+=ak/Tk
    delta = s_Fa/n2_a
    return (sum((1/Tk-delta*ak)**2 for ak, Tk in zip(a,T))/n2_F)**0.5

F = [ 1/Ta for Ta in T ]
tol=1e-1
max_F = max(F)
M = np.eye(len(F), dtype=int)
while True:
    k0,k1 = largest(F)
    err = defect(M[k0])
    if err < tol:
        tol=0.2*err
        print("  a: ",list(M[k0]))
        print("a*T: ",[a*Ta for a,Ta in zip(M[k0],T)])
        print(" ^F: ",[Fk/F[k0] for Fk in F])
        print(f"angle (deg): {np.rad2deg(np.arcsin(err)):.5g}°\n")
    if err < 5e-10: break
    q = int(round(F[k0]/F[k1]))
    F[k0] -= q*F[k1]
    M[k1] += q*M[k0]
    if F[k0]<0: F[k0]=-F[k0]; M[k0]=-M[k0]

Das ergibt das Protokoll

  a:  [1, 1, 1]
a*T:  [27.212220815, 29.530588861, 27.554549886]
 ^F:  [0.013482132197497906, 1.0, 0.07171370910340992]
angle (deg): 2.035°

  a:  [15, 14, 15]
a*T:  [408.18331222499995, 413.42824405399995, 413.31824829]
 ^F:  [0.18799937091605323, 0.055664774527038254, 1.0]
angle (deg): 0.34535°

  a:  [76, 70, 75]
a*T:  [2068.12878194, 2067.14122027, 2066.59124145]
 ^F:  [1.0, 0.29609021698213056, 0.31916673511916593]
angle (deg): 0.018052°

  a:  [242, 223, 239]
a*T:  [6585.35743723, 6585.321316003, 6585.537422754]
 ^F:  [0.14353663918949092, 1.0, 0.0779374555912061]
angle (deg): 0.00082299°

  a:  [3783, 3486, 3736]
a*T:  [102943.831343145, 102943.63276944599, 102943.798374096]
 ^F:  [0.15830991529461017, 0.06102937657324801, 1.0]
angle (deg): 4.7152e-05°

  a:  [20928, 19285, 20668]
a*T:  [569497.35721632, 569497.406184385, 569497.437043848]
 ^F:  [1.0, 0.3855057117532664, 0.3167237386175655]
angle (deg): 3.3867e-06°

  a:  [66325, 61118, 65501]
a*T:  [1804850.545554875, 1804850.530006598, 1804850.572082886]
 ^F:  [0.49417557377594934, 0.21716709153510322, 1.0]
angle (deg): 5.4577e-07°

  a:  [285986, 263534, 282433]
a*T:  [7782314.18199859, 7782314.204894774, 7782314.187952638]
 ^F:  [0.2755545984556782, 1.0, 0.05364004447339751]
angle (deg): 6.9819e-08°

  a:  [1255666, 1157087, 1240066]
a*T:  [34169460.46188779, 34169460.4734079, 34169460.458932474]
 ^F:  [1.0, 0.3709551370058305, 0.19466212784696205]
angle (deg): 1.0192e-08°

wobei das Referenzergebnis tatsächlich im 4. Datensatz steht.

Für$3$Perioden

Teilen$aT_1 = bT_2 = cT_3$hinein$2$einfachere Fälle:$aT_1 = bT_2$Und$bT_2 = cT_3$.

Es folgt dem

$$\frac{a}{b} = \frac{T_2}{T_1} \; \text{and} \; \frac{b}{c} = \frac{T_3}{T_2}$$

Die sachgerechte Weiterführung lässt sich am besten an einem Beispiel demonstrieren.

Betrachten Sie den bereitgestellten Fall wo$T_1 = 27.212220815, T_2 = 29.530588861$Und$T_3 = 27.554549886$.

Wir finden geeignete Näherungen für$\frac{a}{b}$Und$\frac{b}{c}$Verwenden Sie die in der Frage beschriebene Kettenbruchmethode und erhalten Sie

$$\frac{a}{b} \approx [1; 11, 1, 2, 1, 4, 3, 5, 1, 31, 1, 1, 6] = \frac{2381497}{2194532}$$ $$\frac{b}{c} \approx [0; 1, 13, 1, 16, 1, 27, 3, 5, 1, 1] = \frac{879525}{942599}$$

Wir können jetzt multiplizieren$\frac{a}{b}$den Nenner passend zu machen$\frac{a}{b}$gleich dem Zähler in$\frac{b}{c}$folgendermaßen:

$$\frac{a}{b} \times \frac{2194532}{879525} \approx \frac{2381497}{2194532} \times \frac{2194532}{879525} = \frac{2381497}{879525}$$ $$\implies \frac{\frac{2194532}{879525}a}{b} \approx \frac{2381497}{879525}$$

Wir können jetzt gleichsetzen

$$\frac{2194532}{879525}a = 2381497$$ $$\implies a = \frac{2094586148925}{2194532}$$

Somit

$$a : b : c : : \frac{2094586148925}{2194532} : 879525 : 942599$$ $$\implies a : b : c : : 2094586148925 : 1930145757300 : 2068563668668$$

Deshalb,$2094586148925T_1 \approx 1930145757300T_2 \approx 2068563668668T_3$.

Offensichtlich sind diese Zahlen größer als$242, 223$Und$239$aber das liegt daran, dass ich mich dafür entschieden habe, den fortgesetzten Bruch über einen bestimmten Punkt hinaus abzuschneiden. Wenn Sie möchten, können Sie auch die gefundenen Werte für übernehmen$a, b, c$und finden Sie kleinere ganzzahlige Annäherungen für$a : b : c$.

Für$n$Perioden

Das obige Argument lässt sich leicht wie folgt verallgemeinern:

In Betracht ziehen$a_1T_1 = a_2T_2 = a_3T_3 = a_4T_4 = ... = a_nT_n$.

Einfach nehmen$a_1T_1 = a_2T_2, \, a_2T_2 = a_3T_3, \, a_3T_3 = a_4T_4, ... , a_{n-1}T_{n-1} = a_{n}T_{n}$, bekommen

$$\frac{a_1}{a_2} = \frac{T_2}{T_1}, \, \frac{a_2}{a_3} = \frac{T_3}{T_2}, ..., \frac{a_{n-1}}{a_n} = \frac{T_n}{T_{n-1}}$$

Multiplizieren Sie dann jeden der Brüche entsprechend, sodass die erforderlichen Zähler und Nenner übereinstimmen. Suchen Sie danach den Wert für$a_1$durch Gleichsetzen der Zähler. Betrachten Sie als Nächstes das Verhältnis$a_{1} : a_{2} : a_{3} : ... : a_{n}$, multipliziere, um alles zu einer Ganzzahl zu machen und voila, deine Antwort ist fertig!