Encontrar aproximaciones enteras

Puede probar el algoritmo de reducción de base de celosía de Lenstra–Lenstra–Lovász (LLL) para calcular un vector corto en una celosía. Yo uso Pari GP (p. 382, ​​sección 3.10.62 qflll) para el cálculo.

Definimos una red con una base formada por las columnas de la Matriz$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$es un factor de escala. Una combinación lineal de estos vectores base

$$\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}$$

es de pequeña longitud, si${\lambda_1} T_d+ {\lambda_2} T_s$y${\lambda_1} T_d+ {\lambda_3} T_a$es pequeño y$\frac {\lambda_1} m, \frac {\lambda_2} m , \frac {\lambda_3} m$también son pequeños. En este caso

$${\lambda_2} T_s - {\lambda_3} T_a= ({\lambda_1} T_d+ {\lambda_2} T_s )-({\lambda_1} T_d+ {\lambda_3} T_a)$$ es pequeño, también. Tenga en cuenta que${\lambda_2}$y${\lambda_3}$tener el mismo signo.

Con Pari/GP en línea obtengo el siguiente resultado

\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]]

Los coeficientes son$242, 223, 239$, las diferencias son$0.036, 0.180, 0.216$.

Agregué la duración del mes tropical ($T_t$) y el mes sideral ($Ti$), de Wiki y obtuve, para 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]]

Invertir el problema, con$F_a=1/T_a$etc. Entonces estamos buscando números enteros$a,b,c$que satisfacen

$$\frac{F_a}{a}\approx\frac{F_b}{b}\approx\frac{F_c}{c}\approx \delta.$$ Es decir, queremos encontrar algún espacio de cuadrícula$δ$para que todos$F_{a,b,c}$están muy cerca de los puntos de la cuadrícula. La nueva interpretación de la tarea es, pues, $$\text{ minimize }\|\vec F-δ\vec a\|_2^2\text{ where }\vec a\in\Bbb Z^3,$$ con$\vec a$mínimo, es decir, coprimo como vector.

Esto se puede hacer mediante algún tipo de algoritmo euclidiano, similar a cómo se relaciona con él el cálculo de la fracción continua. En el campo de las motivaciones también se encuentra el algoritmo L3.

Iterar

• Tome los dos elementos más grandes y reduzca el más grande por el segundo más grande usando cocientes enteros.
• Realizar un seguimiento de la matriz de transformación$M$, aquí con la convención de que si$\vec F'$es el vector reducido en alguna etapa, entonces$M^T\vec F'=\vec F$es igual al vector original.
• Entonces siempre que el (nuevo) elemento maximal de$\vec F'$en la posición$k$es decididamente más grande que los elementos restantes, el vector$\vec a=M_k$, el$k$la fila de$M$es un candidato para el vector de coeficientes enteros.

Como medida de error, use el residuo de la tarea de minimización, o estrechamente relacionado, una medida de cuán paralelo$\vec F$y$\vec a$son.

En el siguiente script de python se realizan estas ideas. Se imprime un informe para el caso de que el residual mejore sustancialmente.

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]

Esto da el registro

  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°

donde efectivamente el resultado de referencia está en el 4º registro.

Para$3$periodos

Dividir$aT_1 = bT_2 = cT_3$en$2$casos más simples:$aT_1 = bT_2$y$bT_2 = cT_3$.

Resulta que

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

La continuación apropiada se puede demostrar mejor con un ejemplo.

Considere el caso provisto donde$T_1 = 27.212220815, T_2 = 29.530588861$y$T_3 = 27.554549886$.

Encontramos aproximaciones apropiadas para$\frac{a}{b}$y$\frac{b}{c}$usando el método de fracción continua descrito en la pregunta, obteniendo

$$\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}$$

Ahora podemos multiplicar$\frac{a}{b}$apropiadamente para hacer el denominador en$\frac{a}{b}$igual al numerador en$\frac{b}{c}$como sigue:

$$\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}$$

Ahora podemos igualar

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

Por eso

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

Por lo tanto,$2094586148925T_1 \approx 1930145757300T_2 \approx 2068563668668T_3$.

Obviamente, estos números son más grandes que$242, 223$y$239$pero eso es porque elegí truncar la fracción continua más allá de cierto punto. Si lo desea, también puede tomar los valores encontrados para$a, b, c$y encontrar aproximaciones enteras más pequeñas para$a : b : c$.

Para$n$periodos

El argumento anterior se puede generalizar fácilmente de la siguiente manera:

Considerar$a_1T_1 = a_2T_2 = a_3T_3 = a_4T_4 = ... = a_nT_n$.

simplemente toma$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}$, conseguir

$$\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}}$$

Luego multiplica cada una de las fracciones apropiadamente para que coincidan los numeradores y denominadores necesarios. Después de eso, encuentra el valor de$a_1$igualando los numeradores. A continuación, considere la relación$a_{1} : a_{2} : a_{3} : ... : a_{n}$, multiplica para que todo sea un número entero y listo, ¡tu respuesta está lista!