Астрономическая матчасть

[Сергей Шмальц]

(Для чтения этой темы потребуется "Включить рендеринг формул": Профиль --> Изменить профиль --> Внешний вид форума --> Включить рендеринг формул.)

Поскольку вопросы по астрономической матчасти (как минимум у меня) буду в обозримом будущем возникать всё чаще, то чтобы не плодить множество тем, я решил завести одну общую на все случаи жизни.

В этом сообщении речь пойдёт о методе наименьших квадратов, используемом для нахождения эллипса орбиты по опытным точкам в плоскости координат. Я не математик, но вроде бы разобрался с этим методом, здесь же призываю людей с хорошими знаниями математики просто проверить нижеприведённое решение на логические и прочие ошибки.

Итак, имеется множество точек проиндексованных через i=1...n, все они заведомо лежат вблизи линии такого эллипса, что невязки для них всех наименьши. Для нахождения эллипса я использую каноническое уравнение эллипса вида:

\[ \frac{x^2}{a^2}+\frac{y^2}{b^2}=1 \]

Я применил метод наименьших квадратов так, как он описан на примере линейной функции здесь. Сначала я создал два так называемых нормальных уравнения:

\[
\frac{\sum_{i=1}^n x_i^4}{a^2}+\frac{\sum_{i=1}^n x_i^2y_i^2}{b^2}=\sum_{i=1}^n x_i^2   (1)
 \]

\[
\frac{\sum_{i=1}^n x_i^2y_i^2}{a^2}+\frac{\sum_{i=1}^n y_i^4}{b^2}=\sum_{i=1}^n y_i^2   (2)
 \]

Далее для упрощения писанины заменяю все суммы временными значениями от k₁ до k₅:

\[
\sum_{i=1}^n x_i^4=k_1,  \sum_{i=1}^n x_i^2y_i^2=k_2,  \sum_{i=1}^n x_i^2=k_3,  \sum_{i=1}^n y_i^4=k_4,  \sum_{i=1}^n y_i^2=k_5
 \]

Теперь решаю систему двух уравнений с двумя неизвестными a и b:

\[
\left|
\begin{matrix}
\frac{k_1}{a^2}+\frac{k_2}{b^2}=k_3\\
\frac{k_2}{a^2}+\frac{k_4}{b^2}=k_5
\end{matrix}
\right.
 \]

Из второго уравнения вывожу a²:

\[
a^2=\frac{b^2k_2}{b^2k_5-k_4}   (3)
 \]

Подставляю значение a² в первое уравнение и получаю:

\[
\frac{b^2k_5k_1-k_4k_1}{b^2k_2}+\frac{k_2}{b_2}=k_3
 \]

Нахожу из него b²:

\[
b^2=\frac{k_4k_1-k_2^2}{k_5k_1-k_2k_3}
 \]

Так как отрицательное значение b в данном случае нас не интересует, то можно оформить итог в абсолютном значении:

\[
b=\left| \sqrt{\frac{k_4k_1-k_2^2}{k_5k_1-k_2k_3}} \right|
 \]

Подставив значение b² в уравнение (3) нахожу а:

\[
a=\left| \sqrt{\frac{k_1k_4-k_2^2}{k_3k_4-k_2k_5}} \right|
 \]

Далее для наглядности можно k₁...k₅ снова записать суммами, но здесь опущу это "косметическое" преобразование, лень писать на LaTeX кучу букав.

Вот, собственно, и всё. Найдены обе полуоси эллипса, аппроксимирующего все опытные точки.

Нет ошибок?

А нужно мне это всё как один из вычислительных шагов для нахождения орбитального периода околосолнечных комет группы Майера – задачи, о которой я писал в теме про наблюдения на SOHO.

Заранее благодарю за комментарии к вышеприведённому.

[Сергей Шмальц]

Вопрос снимается с повестки дня, мне только что подтвердили в одном математическом интернет-форуме, что всё сделано верно. Спасибо за внимание.

[xd]

Верно-то оно верно, но тут нет учёта поворотов и сдвигов эллипса. Вы предполагаете что эллипс имеет центр в начале координат, и большой полуосью ориентирован по одной из осей (заметьте, формула инвариантна относительно разворота на pi/2 и отражения относительно осей и начала координат)

[Виктор Фёдорович]

Уважаемый Сергей, это какую «матчасть» Вы имеете ввиду? Судя по приведённым символам (а это маетматические условности, и совсем не астрономические, что совсем не одно и то же), Вы решили абстрактно «поразмышлять», оторвавшись от реалий.

Во-первых, какой «эллипс» какой «орбиты» Вы в этих манипуляциях с символами выкрутили? Никаких «опытных точек» в реальном космическом пространстве нет, они только могут быть в Вашем воображении. Поэтому никакой логичностью здесь и не пахнет.

Во-вторых, реальная, относительно стабильная орбита космического объекта существенно отличается от той, которую Вы создали в своём сознании (и в приведённой формуле).

В-третьих, «множество точек», которое «заведомо лежат вблизи линии эллипса», это  заведомый постулат, в который Вы почему-то поверили без проверки на реальность и предлагаете это сделать окружающим! Вы полагаете, что от этого (от Вашей веры) можно получить реальную орбиту космического объекта? Никогда!
Более того, в своих «расчётах» Вы используете ещё один постулат: каноническое уравнение эллипса – математический вымысел, не применимый к фактической форме движения реального космического объекта.
Но, как известно, применение в серьёзной научной работе хоть одного  постулата всегда приводило к полностью ложному математическому и физическому результату и его интерпретации. А у Вас их уже, как минимум, два! Осознаёте ли Вы это?

Реальная орбита реального космического объекта имеет очень сложную форму, которая не поддаётся выражению математической формулой, потому что определяется массой факторов, ни один из которых не учтён в Ваших формулках. Да Вы их просто и не знаете. Следовательно, к сожалению, Ваше «убеждённое» математическое решение, подтверждённое на математическом форуме, полностью неверно.

[xd]

Это вообще к чему было-то?

[Сергей Шмальц]

Верно-то оно верно, но тут нет учёта поворотов и сдвигов эллипса. Вы предполагаете что эллипс имеет центр в начале координат, и большой полуосью ориентирован по одной из осей (заметьте, формула инвариантна относительно разворота на pi/2 и отражения относительно осей и начала координат)

Абсолютно верно, каноническое уравнение предполагает именно такую ориентацию в системе координат, о которой Вы говорите. Мне было важно узнать верно ли я вообще применил метод наименьших квадратов, т.к. никогда раньше не имел с ним дело. Теперь его следует применить к уравнению эллипса как кривой второго порядка, там будет уже система из шести уравнений, что не сильно изменит ход решения. Кстати, спасибо за комментарий, а то я хоть и осознавал это, но как-то не был на все 100 уверен, однако теперь сомнений не осталось.

[Сергей Шмальц]

Уважаемый Сергей, это какую «матчасть» Вы имеете ввиду?

Имею в виду ту же самую матчасть, которой пользовались и пользуются астрономы, такие как например Брайан Марсден – можете ознакомиться с его работой The sungrazing comet group (1967), в которой он тоже прибегает к использованию метода наименьших квадратов для решения системы 5 уравнений для пяти комет группы Крейца. Я прекрасно знаю и понимаю, что орбитальное движение – вещь сложная и трудноописуемая формулами, но я и не пытаюсь описать орбитальное движение группы Майера во всех его тонкостях, а лишь пытаюсь в первом приближении определить орбитальный период. Для этой цели есть не вымышленные мной опытные точки, а реальные точки (положение кометы) в плоскости снимка. В плоскости снимка методом наименьших квадратов определяется эллипс наилучшим образом аппроксимирующий их. Далее производится двухшаговая трансформация эллипса (точнее, в моём случае, только длины большой полуоси) сначала из плоскости снимка в другую плоскость с учётом орбитального наклона к эклиптике, и далее с учётом угла между осью икс плоскости снимка и линии узлов орбиты (вычисляемого из значения угла восходящего узла и солнечной долготы SOHO), при этом учитывается аргумент перигелия. После подобной трансформации из перигелийного растояния в единицах измерения и перигелийного расстояния в а.е вычисляется длина в а.е на 1 ед. измерения, и получив значение масштаба вычисляется окончательная длина большой полуоси в а.е., что тут же даёт и орбитальный период в годах. Для первого приближения этого более чем достаточно.

[xd]

Теперь его следует применить к уравнению эллипса как кривой второго порядка, там будет уже система из шести уравнений, что не сильно изменит ход решения.

Там уравнения будут отнюдь не линейными. Возможно придётся решать задачу оптимизации численными методами, тут могут быть проблемы, когда решение идёт вразнос.

[Сергей Шмальц]

Там уравнения будут отнюдь не линейными.

Да вроде бы линейные у меня получаются, только загвоздка в другом. В общем, вот что у меня выходит:

1) исходное уравнение кривой второго порядка:

\[ a_{11}x^2 + a_{22}y^2 + 2a_{12}xy + 2a_{13}x + 2a_{23}y + a_{33}=0 \]

2) после применения метода наименьших квадратов выходит система из шести уравнений с шестью неизвестными:

\[
a_{11}\sum_{i=1}^n x_i^4 + a_{22}\sum_{i=1}^n x_i^2y_i^2 + 2a_{12}\sum_{i=1}^n x_i^3y_i + 2a_{13}\sum_{i=1}^n x_i^3 + 2a_{23}\sum_{i=1}^n x_i^2y_i + a_{33}\sum_{i=1}^n x_i^2   (1)\\

a_{11}\sum_{i=1}^n x_i^2y_i^2 + a_{22}\sum_{i=1}^n y_i^4 + 2a_{12}\sum_{i=1}^n x_iy_i^3 + 2a_{13}\sum_{i=1}^n x_iy_i^2 + 2a_{23}\sum_{i=1}^n 2y_i^3 + a_{33}\sum_{i=1}^n y_i^2   (2)\\

a_{11}\sum_{i=1}^n x_i^3y_i + a_{22}\sum_{i=1}^n x_iy_i^3 + 2a_{12}\sum_{i=1}^n x_i^2y_i^2 + 2a_{13}\sum_{i=1}^n x_i^2y_i + 2a_{23}\sum_{i=1}^n x_iy_i^2 + a_{33}\sum_{i=1}^n x_iy_i   (3)\\

a_{11}\sum_{i=1}^n x_i^3 + a_{22}\sum_{i=1}^n x_iy_i^3 + 2a_{12}\sum_{i=1}^n x_i^2y_i + 2a_{13}\sum_{i=1}^n x_i^2 + 2a_{23}\sum_{i=1}^n x_iy_i + a_{33}\sum_{i=1}^n x_i   (4)\\

a_{11}\sum_{i=1}^n x_i^2y_i + a_{22}\sum_{i=1}^n y_i^3 + 2a_{12}\sum_{i=1}^n x_iy_i^2 + 2a_{13}\sum_{i=1}^n x_iy_i + 2a_{23}\sum_{i=1}^n y_i^2 + a_{33}\sum_{i=1}^n y_i   (5)\\

a_{11}\sum_{i=1}^n x_i^2 + a_{22}\sum_{i=1}^n y_i^2 + 2a_{12}\sum_{i=1}^n x_iy_i + 2a_{13}\sum_{i=1}^n x_i + 2a_{23}\sum_{i=1}^n y_i + a_{33}n   (6)
 \]

Но при решении системы классическим методом Гаусса выходит, что все неизвестные от a₁₁ до a₃₃ равны 0 – явно не искомое решение. Тут я и забуксовал...

[xd]

Совершенно верно. Проблема в том, что Вы неверно задаёте метрику для невязок (кстати в обоих случаях).
Поясню. Классический МНК определяет уравнения так, чтобы
\[ \sum_i y_i^2 \to min \]
Но для двумерной постановки задачи это неверно. Вам надо определять расстояние не до ближайшей точки в сечении по Y (рисунок, чёрные линии), а до ближайшей точки на кривой. Или расстояние до точки, лежащей на том же радиус-векторе от центра эллипса до текущей точки (не уверен что это равнозначно, надо подумать). На рисунке эти линии показаны красным, нас интересуют отрезок от точки до окружности эллипса. Я сходу затрудняюсь написать, как будут выглядеть уравнения. Вполне возможно, что вид уравнений будет проще, если ввести смещение и разворот эллипса, но в любом случае уравнения будут нелинейными или трансцендентными.
То решение, которое получилось, с нулевыми коэффициентами при членах 2 порядка, совершенно логично следуют из определения метрики, которую Вы задаёте. Вот

[Сергей Шмальц]

Я тут тоже время даром не терял и искал в интернете информацию.

1) http://dxdy.ru/topic38557.html
Алексей К. как раз затрагивает мою проблему: "Вот и с эллипсом, чтобы не получить a=b=c=d=...=0, мы должны придумать разумный constraint на коэффициенты. Что ребята и пытаются сделать. [...]"

2) http://forum.sources.ru/index.php?showtopic=287351&view=findpost&p=2460534
Tur решает задачу методом абсолютно идентичным моим представлениям – я тоже в итоге собираюсь найти инварианты D и I, из которых тут же можно посчитать длину большой полуоси (все остальные расчёты, например, центра эллипса, мне не нужны). Интересно, что он для а₃₃ произвольно берёт значение 1, я пока не совсем понял почему.

3) Несколько статей по теме (PDF):

http://arrow.dit.ie/cgi/viewcontent.cgi?article=1003&context=biodevcon
http://www.emis.de/journals/BBMS/Bulletin/sup962/gander.pdf
http://research.microsoft.com/en-us/um/people/awf/ellipse/ellipse-pami.pdf
http://autotrace.sourceforge.net/WSCG98.pdf

4) При поиске наткнулся на одну задачку (причём строго астрономического характера – речь идёт о малом теле в Солнечной системе) для студентов, где, как и в моём случае, даются n точек с заданными координатами в плоскости, такое же уравнение кривой 2-го порядка, а вот найти коэффициенты требуют не через МНК, а неким методом Хаусхолдера. Интересно, чем он лучше/хуже МНК?

В общем, задача, как я понял, вполне себе решаема, надо читать и разбираться. До связи.

[xd]

по пункту 2)
Коэффициенты эти судя по всему не являются независимыми, один лишний. Один (любой) фиксируем - получаем однозначную систему уравнений.