Основное взаимодействие двух тел с использованием ode45 Matlab

Question

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

Вот функция, я звоню с ode45:

function dy = rigid(t,y) 
dy = zeros(4,1); 

%Position 
xx=y(1); 
yy=y(2); 

%Radius 
r=(xx.^2+yy.^2).^0.5; 

%Constants 
M=10^30; 
G=6.67*10^-11; 

%dX/dt 
dy(1)=y(3); %vx 
dy(3)=-M.*G.*xx.*r.^-3; %ax 

%dY/dt 
dy(2)=y(4); %vy 
dy(4)=-M.*G.*yy.*r.^-3; %ay 

И здесь решатель линии:

%Init 
x_0=1; 
y_0=1; 
vx_0=0; 
vy_0=5; 

%ODEs 
[T,Y] = ode45(@rigid,[0 1000],[x_0 y_0 vx_0 vy_0]); 

%Vectors 
x=Y(:,1); 
y=Y(:,2); 

%Plot 
plot(x,y) 
xlabel('x'); 
ylabel('y'); 
title('y=f(x)'); 

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

Я размышлял об этом некоторое время сейчас, и у меня очень мало идей, которые сделали несколько поисков в Интернете.

Answer 1

Реферат: Существуют проблемы с интегрированием гамильтоновых систем с нормальными численными интеграторами, и ваши особые начальные условия усугубляют это до такой степени, что численное решение не имеет сходства с правильным.

---

Нет ничего плохого в вашей реализации как таковой. Однако начальные условия, которые вы используете, не самые лучшие. Константы G и M, которые вы используете, находятся в единицах СИ, что означает, что координаты находятся в м, скорости находятся в м/с, а время - в с. Эти линии

x_0=1; 
y_0=1; 
vx_0=0; 
vy_0=5; 
[T,Y] = ode45(@rigid,[0 1000],[x_0 y_0 vx_0 vy_0]); 

поэтому означает, что вы просите на орбите с радиусом около 1,4 метра и орбитальной скоростью 5 м/с, и вы хотите эту орбиту в течение периода 17 минут. Представьте себе, что на самом деле объект находился всего в нескольких метрах от массы в 10-30 килограммов!

Так давайте попробуем попросить что-то более реалистичного, похожие на Earths' orbit, и смотреть на него в течение 1 года:

x_0=149.513e9; 
y_0=0; 
vx_0=0; 
vy_0=29.78e3; 
[T,Y] = ode45(@rigid,[0 31.536e6],[x_0 y_0 vx_0 vy_0]); 

И результат

, который выглядит, как и ожидалось ,

Но здесь есть вторая проблема. Давайте посмотрим на этой орбите в течение 10 лет ([0 315.36e6]):

Теперь мы уже не получим замкнутую орбиту, но спираль! Это связано с тем, что численное интегрирование протекает с ограниченной точностью, и для этой системы уравнений это приводит (физически) к loss of energy. Точность может быть увеличена с использованием параметров до ode45, но в конечном итоге проблема будет сохраняться.

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

Эта «орбита» представляет собой прямую линию в направлении начала координат (солнце). Что может быть в порядке, поскольку прямое колебание является возможным частным случаем эллиптической орбиты.Но, глядя на координатах с течением времени:

Мы видим, что нет никаких колебаний. «Планета» падает на солнце и остается там. То, что здесь происходит, - это тот же эффект, что и с большой орбитой: Неточная интеграция приводит к потере энергии. Более того, численное интегрирование застревает: мы просили период в 1000 с, но интеграция не выходит за пределы 1.6e-10 секунд.

---

Насколько я могу судить, ни ode45, ни других стандартных интеграторами Matlab не являются достаточными для решения этой проблемы. Для этого разработаны специальные алгоритмы цифровой интеграции, в частности для гамильтоновых систем, называемых symplectic integrators. Существует file exchange entry, который обеспечивает различные реализации. Также см. this question and answers для большего количества указателей.