Nové metody řešení Keplerovy rovnice

[vaclavII]

Keplerova rovnice ve tvaru E = M + e*sin(E) kde e - excentricita dráhy v intervalu 0-0,99999999.
M - střední anomálie, počítáme excentrickou anomálii E.
Pro hyperbolickou dráhu e je větší než 1.
e*sinh(E) - E = M , sinh - hyperbolický sinus.

Zajímavá metoda zde:
https://www.researchgate.net/publicatio ... s_Equation
a
zde metoda CORDIC
https://arxiv.org/pdf/1808.07062

[palkog]

Ja som davnejsie porovnaval rozne metody na riesenie Keplerovej rovnice. A ako najrychlejsia sa mi zdala klasicka Newton-Raphsonova metoda, pricom som ako startovaciu hodnotu pouzil vztah S9 z tejto prace https://articles.adsabs.harvard.edu/pdf ... ..38..307O - vacsinou uz po 3-4 iteraciach bola presnost na urovni 1e-10. A pre excentricitu vacsiu ako 0.8-0.9 som pouzival tuto aproximaciu https://ntrs.nasa.gov/api/citations/199 ... 021346.pdf.

[Psion]

Ta rychlost výpočtu Keplerovy bude asi dobrá při výpočtu stovek tisíc poloh těles, pro pár výpočtů to asi nic viditelně nezrychlí, můj názor. Každopádně zajímavé metody.

[vaclavII]

Tyto metody mám v jazyku C pod Visual Studio 2022 , rychlost je dobrá a přesnost e-14 při řešení metodou půlení intervalu nepotřebuji 1.derivaci Keplerovy rovnice a v 50 krocích je hotovo.
zde bisekce https://cs.wikipedia.org/wiki/P%C5%AFle ... rval%C5%AF
Zkoušel jsem řešení rovnice Halleyovou metodou, zde používám vypočtenou 1 a 2. derivaci Keplerovy rovnice. Počet kroků 4-6 pro dráhu téměř parabolickou e = 0.99999999.
Halley metoda zde https://en.wikipedia.org/wiki/Halley%27s_method
Newtonova metoda zde https://cs.wikipedia.org/wiki/Metoda_te%C4%8Den

Kus kódu v C
public static double FNKepler(double e, double m, double exc)
{
double f = 0.0;

try // Funkční hodnota Keplerovy rovnice.
{
f = exc switch // Excentricita drahy.
{
<= 1.0 => e - exc * Math.Sin(e) - m,
_ => exc * Math.Sinh(e) - e - m,
};
}
catch (Exception ex)
{
Console.WriteLine(ex.ToString());
}

return f;
}

public static double FND1Kepler(double e, double exc)
{
double f = 0.0;

try // První derivace Keplerovy rovnice.
{
f = exc switch // Excentricita drahy.
{
<= 1.0 => 1.0 - exc * Math.Cos(e),
_ => exc * Math.Cosh(e) - 1.0,
};
}
catch (Exception ex)
{
Console.WriteLine(ex.ToString());
}

return f;
}

public static double FND2Kepler(double e, double exc)
{
double f = 0.0;

try // Druha derivace Keplerovy rovnice.
{
f = exc switch // Excentricita drahy.
{
<= 1.0 => exc * Math.Sin(e),
_ => exc * Math.Sinh(e),
};
}
catch (Exception ex)
{
Console.WriteLine(ex.ToString());
}

return f;
}

Mimo jiné metodou CORDIC lze přesně spočítat goniometrické a bolometrické funkce, dokončuji na toto malý prográmek v C.

[vaclavII]

Pokud je M blízko 0, např. 0.000001 a excentricita dráhy 0.999999 je první a druhá derivace Keplerovy rovnice velikosti řádu epsilon = e-10 a metoda Newtonova i Halleyova selhávají. Metody Cordic a Fukushima, případně půlení intervalu - bisekce neselhávají, sice potřebují cca 50 výpočetních cyklů.