Прогноз орбиты в Гринвиче

[GreenEkatherine]

Всем доброго вечера!

Я пытаюсь прогнозировать положение и скорость спутника примерно на 10 секунд по методу Рунге-Кутты 4 порядка. Орбиту считаю невозмущенной. Особенность заключается в том, что начальные положение и скорость заданы в Гринвиче и через 10 секунд я хочу их вычислить так же в Гринвиче. Положение спутника получается довольно правдоподобным (30 м ошибка), а вот в скорости накапливается ошибка. Пробовала крутить Землю, ошибка скорости уменьшается (но все еще накапливается), а в положении спутника вылезает лишний поворот Земли за 10 секунд.

Привожу свой код, начальные данные и те, которые должны получиться через 10 секунд. Подскажите пожалуйста, в чем косяк, а то уже всю голову сломала.

Начальные:
Point(5831740.00; 401393.00; 3631890.00)
Velocity(-3862.03; -1853.99; 6384.90)

Ожидаемые:
Point(5792760.00; 382857.00; 3695520.00)
Velocity(-3935.49; -1853.16; 6339.87)

Код: [Выделить]

std::vector<double> y(6), p(6), q(6), r(6), s(6), f(6);
double h = 0.0001;

y[0] = Point.x;
y[1] = Point.y;
y[2] = Point.z;
y[3] = Velocity.x;
y[4] = Velocity.y;
y[5] = Velocity.z;

f[0] = Velocity.x;
f[1] = Velocity.y;
f[2] = Velocity.z;
f[3] = -(double)GRAVITATION_PARAM*Point.x/(AbsVector(Point)*AbsVector(Point)*AbsVector(Point));
f[4] = -(double)GRAVITATION_PARAM*Point.y/(AbsVector(Point)*AbsVector(Point)*AbsVector(Point));
f[5] = -(double)GRAVITATION_PARAM*Point.z/(AbsVector(Point)*AbsVector(Point)*AbsVector(Point));

for (unsigned long j = 0; j < 100000; j++)
{
for (int i = 0; i < 6; i++)
{
p[i] = f[i];
q[i] = f[i] + h*p[i]/2;
r[i] = f[i] + h*q[i]/2;
s[i] = f[i] + h*q[i];
y[i] += h*(p[i] + 2*q[i] + 2*r[i] + s[i])/6;
}
f[0] = y[3];
f[1] = y[4];
f[2] = y[5];
f[3] = -(double)GRAVITATION_PARAM*y[0]/(AbsVector(CIC_Point3DD(y[0], y[1], y[2]))*AbsVector(CIC_Point3DD(y[0], y[1], y[2]))*AbsVector(CIC_Point3DD(y[0], y[1], y[2])));
f[4] = -(double)GRAVITATION_PARAM*y[1]/(AbsVector(CIC_Point3DD(y[0], y[1], y[2]))*AbsVector(CIC_Point3DD(y[0], y[1], y[2]))*AbsVector(CIC_Point3DD(y[0], y[1], y[2])));
f[5] = -(double)GRAVITATION_PARAM*y[2]/(AbsVector(CIC_Point3DD(y[0], y[1], y[2]))*AbsVector(CIC_Point3DD(y[0], y[1], y[2]))*AbsVector(CIC_Point3DD(y[0], y[1], y[2])));
}


[Geen]

Во-первых, чему равно "GRAVITATION_PARAM"?
Во-вторых, откуда взялось "ожидаемое"?
В-третьих, "производную" надо считать на каждой "подитерации".

[GreenEkatherine]

Geen,

1) GRAVITATION_PARAM = 398600441899999;
2) Ожидаемое - с датчиков GPS;
3) Вы хотите сказать, что производную нужно считать для каждого i? Но ведь итерации идут по j, а i - это всего лишь проход по элементам вектора.

[Geen]

Вы хотите сказать, что производную нужно считать для каждого i?

Нет, для каждого p,q,...
https://ru.wikipedia.org/wiki/%D0%9C%D0%B5%D1%82%D0%BE%D0%B4_%D0%A0%D1%83%D0%BD%D0%B3%D0%B5_%E2%80%94_%D0%9A%D1%83%D1%82%D1%82%D1%8B

[GreenEkatherine]

Нет, для каждого p,q,...
https://ru.wikipedia.org/wiki/%D0%9C%D0%B5%D1%82%D0%BE%D0%B4_%D0%A0%D1%83%D0%BD%D0%B3%D0%B5_%E2%80%94_%D0%9A%D1%83%D1%82%D1%82%D1%8B

Да, я видела это. Но здесь специфика в том, что нет функции, есть только данные. Я программировала формулы из "Космической фотограмметрии" Урмаева, с. 139 (в прикрепленных)

То есть, в данном случае, я могу пересчитывать производные в подитерациях только для f[3], f[4], f[5]: x(t + h/2), r(t + h/2), ..., если придумать, как их рассчитать для этого шага

[phobos24]

В правых частях дифференциальных уравнений нужно учитывать вращение Земли.
По рабоче-крестьянски:
Нужно к dVx/dvt добавить omega*omega*x+omega*Vy
К dVy/dt добавить omega*omega-omega*Vx
omega - угловая скорость вращения Земли

[GreenEkatherine]

В правых частях дифференциальных уравнений нужно учитывать вращение Земли.

От рабочих и крестьян большое спасибо! Учла, пересчитала. Допустима ли в таких масштабах погрешность для скорости (-1.33; 2.77; -0.22)? С вектором положения все стало еще лучше, там (1.84; 15.27; -2.52).

[Geen]

Но здесь специфика в том, что нет функции, есть только данные

Не бывает такой "специфики". Тем более, что вычисление \(f\) у Вас явно написано....

нужно учитывать вращение Земли

Это зачем? Что бы сто тысяч раз дополнительную ошибку накопить? Если нужно учитывать "вращение", то делать это нужно только один раз, в конце...

omega - угловая скорость вращения Земли

А куда \(V_z\) дели?

Допустима ли в таких масштабах погрешность

У Вас пока не с погрешностью проблема, а с алгоритмом... Ваш внутренний цикл эквивалентен банальному \(\Delta y = f h \left(1+\frac{h}{2}+\frac{h^2}{6}\right)\)

[phobos24]

Это зачем? Что бы сто тысяч раз дополнительную ошибку накопить? Если нужно учитывать "вращение", то делать это нужно только один раз

Затем, что она интегрируем во вращающейся системе координат, а не в инерциальной
Скорости вращения по z нет Считаем, что в пределах 10 секунд перевод вращающейся ск и инерциальной связан только с вращением вокруг оси z

По поводу скоростей, мю Земли не wgs-84, проверьте

[Geen]

Затем, что она интегрируем во вращающейся системе координат

Зачем? У нас мало источников ошибок? Нам некуда девать ресурсы процессора?

Скорости вращения по z нет

Опишите оси (систему координат), пожалуйста. (и формулу для силы Кориолиса, только в читаемом виде, пожалуйста) (не по делу)

По поводу скоростей, мю Земли не wgs-84, проверьте

Это Вы сейчас о чём? Не понял.

[phobos24]

Зачем? У нас мало источников ошибок? Нам некуда девать ресурсы процессора?

Прочитайте название темы.

[Geen]

Зачем? У нас мало источников ошибок? Нам некуда девать ресурсы процессора?

Прочитайте название темы.

В названии темы "прогноз", а не " интегрирование "

В любом случае, сначала надо правильно имплементировать RK4 (вместо кривого, на данный момент, "RK1").

[GreenEkatherine]

По поводу скоростей, мю Земли не wgs-84, проверьте

Поправила GRAVITATION_PARAM = 398600500000000, принципиальных изменений нет, ошибка Point(3.28878169; 14.72033843; 1.21065085), Velocity(-1.329678343; 2.762204044; -0.221094735). В принципе, такая точность для меня достаточна, если она будет сохраняться на других данных.

[GreenEkatherine]

Не бывает такой "специфики". Тем более, что вычисление \(f\) у Вас явно написано....

Вычисление f[0] - f[2] неявно, да и для f[3] - f[5] я честно не понимаю, как будет выглядеть f(t + h/2).

Если нужно учитывать "вращение", то делать это нужно только один раз, в конце...

Я в первом сообщении написала, что пробовала поворачивать после вычислений, в итоге ошибка вектора положения показывала как раз на лишний поворот.

У Вас пока не с погрешностью проблема, а с алгоритмом...

Может, лучше тогда сразу Адамса, говорят, он быстрее?

[ForumUser]

Кто нибудь оценивал сходимость 2суточных глонасс эфемерид?

[phobos24]

Кто нибудь оценивал сходимость 2суточных глонасс эфемерид?

Поясните вопрос? О какой сходимости речь?