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

Я пытаюсь смоделировать базовую гравитацию для объекта незначительной массы вокруг массивного тела. Я следовал примерам, представленным в документации пакета 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)');

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

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


matlab ode differential-equations numerical-integration
person Ridgy    schedule 14.02.2015    source источник
comment
Возможно, вы захотите изменить свои комментарии с французского на английский на случай, если люди здесь не понимают... merci :)   -  person Benoit_11    schedule 14.02.2015
comment
Можете ли вы показать, какой результат вы ожидаете?   -  person Benoit_11    schedule 14.02.2015
comment
Я рисую y=f(x), поэтому для этой задачи я ожидаю либо гиперболы, либо параболы, либо эллипсы.   -  person Ridgy    schedule 14.02.2015

Ответы (1)


arrow_upward
2
arrow_downward

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

В вашей реализации как таковой нет ничего плохого. Однако начальные условия, которые вы используете, не самые лучшие. Константы 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 кг!

Итак, давайте попробуем задать что-то более реалистичное, например орбиту Земли, и посмотрим при этом более 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]):

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

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

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

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

Насколько я могу судить, ни ode45, ни другие стандартные интеграторы Matlab не подходят для решения этой проблемы. Существуют специальные алгоритмы численного интегрирования, разработанные специально для гамильтоновых систем и называемые симплектическими интеграторами. Существует запись обмена файлами, обеспечивающая различные реализации. Также см. этот вопрос и ответы для получения дополнительной информации.

person A. Donda    schedule 14.02.2015
comment
Блестящий ответ! Странные начальные условия связаны с тем, что конечной целью проекта является моделирование эффекта Пойнтинга-Робертсона для пылинки, вращающейся вокруг звезды (отсюда и масса), и я действительно не играл с начальными радиусами, масштабированными до планетарные орбиты (конечно, теперь я понимаю, что это должно было быть моей отправной точкой!) Если я останусь с ode45, я смогу ограничить проблему, вычислив только один период и снизив допуск на ошибку, верно? Поищу симплектические интеграторы, тывм! - person Ridgy; 14.02.2015
comment
@Ridgy Рад, что тебе понравилось! – Да, если вы сможете настолько снизить допуски ошибок, что вы действительно получите полный период. Но, учитывая вашу цель, особенно важно, чтобы вы могли различать спиральное движение из-за эффекта от спирального движения из-за ошибок интегрирования. И в этом случае система больше не является гамильтоновой, я полагаю? - person A. Donda; 14.02.2015