Описание проблемы

Проблемой Плеяд является астрономическая проблема механики, описывающая гравитационные взаимодействия семи звезд [1]. Этот кластер звезд также упоминается как Эти Семь Сестер, и он отображается к человеческому глазу в ночном небе из-за его близости, чтобы Заземлить [2]. Система уравнений, описывающая движение звезд в кластере, состоит из 14 нежестких дифференциальных уравнений второго порядка, которые создают систему 28 уравнений, когда переписано в форме первого порядка.

Второй закон ньютона движения имеет отношение, сила применялась к телу своей скорости изменения в импульсе в зависимости от времени,

${\mathit{F}}_{\mathit{i}}=\frac{\mathrm{d}}{\mathrm{d}\mathit{t}}{\mathit{p}}_{\mathit{i}}$.

Импульс (${\mathit{p}}_{\mathit{i}}={\mathit{m}}_{\mathit{i}}{\mathit{v}}_{\mathit{i}}$) имеет отдельный x-и y-компоненты. В то же время универсальный закон гравитации описывает силу, работающую над телом i от тела j как

${\mathit{F}}_{\mathrm{ij}}=\mathit{g}\frac{{\mathit{m}}_{\mathit{i}}{\mathit{m}}_{\mathit{j}}}{{\Vert {\mathit{p}}_{\mathit{i}}-{\mathit{p}}_{\mathit{j}}\Vert }^2}{\mathit{d}}_{\mathrm{ij}}$.

Термин ${\mathit{d}}_{\mathrm{ij}}=\frac{{\mathit{p}}_{\mathit{j}}-{\mathit{p}}_{\mathit{i}}}{\Vert {\mathit{p}}_{\mathit{j}}-{\mathit{p}}_{\mathit{i}}\Vert }$ дает направление расстояния между телами, и массы тел ${\mathit{m}}_{\mathit{i}}=\mathit{i}$ для $\mathit{i}=1,2,...,7$. Для системы со многими телами сила применилась к любому отдельному телу, сумма его взаимодействий со всеми другими, таким образом,

Установка ускорения свободного падения $\mathit{g}$ равняйтесь 1 и решающие выражения система уравнений второго порядка, описывающих эволюцию x-и y-компонентов в зависимости от времени.

где ${\mathit{r}}_{\mathrm{ij}}=\sqrt{{\left({\mathit{x}}_{\mathit{i}}-{\mathit{x}}_{\mathit{j}}\right)}^2+{\left({\mathit{y}}_{\mathit{i}}-{\mathit{y}}_{\mathit{j}}\right)}^2}$. Поскольку эти два уравнения запрашивают каждую из этих семи звезд в системе, 14 дифференциальных уравнений второго порядка ($\mathit{i}=1,2,...,7)$ опишите целую систему.

Начальные условия для системы, как дали в [1]:

Чтобы решить эту систему ОДУ в MATLAB®, необходимо закодировать уравнения в функцию прежде, чем вызвать решатели ode78 и ode89. Можно или включать необходимые функции как локальные функции в конце файла (как сделано здесь) или сохранить их как отдельные, именованные файлы в директории на пути MATLAB.

Уравнения кода

Решатели ОДУ в MATLAB требуют, чтобы уравнения были написаны в форме первого порядка, ${\mathit{q}}^{\prime }=\mathit{u}\left(\mathit{t},\mathit{q}\right)$. Для этой проблемы система уравнений может быть переписана в форме первого порядка с помощью замен $\mathit{w}={\mathit{x}}^{\prime }$ и $\mathit{z}={\mathit{y}}^{\prime }$. С этими заменами система содержит 28 уравнений первого порядка, организованных в четыре группы из семи уравнений:

$$\left[\begin{array}{c}{x}_i^{\prime }\\ y_i^{\prime }\\ w_i^{\prime }\\ z_i^{\prime }\end{array}\right]=\left[\begin{array}{c}{w}_i\\ z_i\\ f_i(x,y)\\ h_i(x,y)\end{array}\right]$$.

Вектор решения, данный путем решения системы, имеет форму

Поэтому при записи исходных уравнений в терминах вектора решения $\mathit{q}$ выражения

где$\mathit{x}=\left({\mathit{q}}_1,{\mathit{q}}_2,...,{\mathit{q}}_7\right)$ и $\mathit{y}=\left({\mathit{q}}_8,{\mathit{q}}_9,...,{\mathit{q}}_{14}\right)$. Уравнениями, написанными в форме первого порядка ${\mathit{q}}^{\prime }=\mathit{u}\left(\mathit{t},\mathit{q}\right)$, можно теперь записать функцию, которая вычисляет компоненты во время каждого временного шага процесса решения. Функция, которая кодирует эту систему уравнений:

function dqdt = pleiades(t,q)
  x = q(1:7);
  y = q(8:14);
  xDist = (x - x.');
  yDist = (y - y.');
  r = (xDist.^2+yDist.^2).^(3/2);
  m = (1:7)';
  dqdt = [q(15:28);
          sum(xDist.*m./r,1,'omitnan').';
          sum(yDist.*m./r,1,'omitnan').'];
end

В этой функции, $\mathit{x}$ и $\mathit{y}$ значения извлечены непосредственно из вектора решения $\mathit{q}$, как первые 14 компонентов выхода. Затем функция использует значения положения, чтобы вычислить расстояния между всеми семью звездами, и эти расстояния используются в коде для ${\mathit{f}}_{\mathit{i}}\left(\mathit{x},\mathit{y}\right)$ и ${\mathit{h}}_{\mathit{i}}\left(\mathit{x},\mathit{y}\right)$.

Установите дополнительные параметры

Используйте odeset функционируйте, чтобы установить значение нескольких дополнительных параметров:

opts = odeset("RelTol",1e-13,"AbsTol",1e-15,"Stats","on");

Решите систему уравнений

Создайте вектор-столбец с начальными условиями и временной вектор с расположенными с равными интервалами точками в области значений $\left[0,15\right]$. Когда вы задаете временной вектор больше чем с двумя элементами, решатель возвращает решения в то время точки, которые вы задаете.

init = [3 3 -1 -3 2 -2 2 ...
        3 -3 2 0 0 -4 4 ...
        0 0 0 0 0 1.75 -1.5 ...
        0 0 0 -1.25 1 0 0]';
tspan = linspace(1,15,200);

Теперь решите уравнения с помощью обоих ode78 и ode89 путем определения уравнений, отрезка времени, начальных значений и опций. Используйте tic и toc ко времени каждый решатель для сравнения (отмечают, что синхронизации могут отличаться в зависимости от вашего компьютера).

tic
[t78,q78] = ode78(@pleiades,tspan,init,opts);

14963 successful steps
549 failed attempts
201899 function evaluations

toc

Elapsed time is 6.994958 seconds.

tic
[t89,q89] = ode89(@pleiades,tspan,init,opts);

4181 successful steps
56 failed attempts
68726 function evaluations

toc

Elapsed time is 2.169416 seconds.

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

Постройте кривые решения

Первые 14 компонентов q89 содержите положения X и Y для каждой из этих семи звезд, как получено ode89. Постройте эти компоненты решения, чтобы видеть траектории всех звезд в зависимости от времени.

plot(q89(:,1),q89(:,8),'--',...
     q89(:,2),q89(:,9),'--',...
     q89(:,3),q89(:,10),'--',...
     q89(:,4),q89(:,11),'--',...
     q89(:,5),q89(:,12),'--',...
     q89(:,6),q89(:,13),'--',...
     q89(:,7),q89(:,14),'--')
title('Position of Pleiades Stars, Solved by ODE89')
xlabel('X Position')
ylabel('Y Position')

Создайте анимацию результатов

Поскольку траектории звезд значительно перекрываются, лучший способ визуализировать результаты состоит в том, чтобы создать анимацию, показывающую перемещение звезд в зависимости от времени. Функциональный AnimateOrbits, включенный как локальная функция в конце этого примера, принимает выход от решателя для этой проблемы и создает анимированный файл GIF в текущей папке, которая показывает, что звезды проходят свои траектории.

Например, можно сгенерировать анимацию с выходом от ode89 решатель с помощью команды

AnimateOrbits(t89,q89);

Следующим является демонстрационный выход GIF.

Ссылки

[1] Hairer, E., и др. Решение Обыкновенных дифференциальных уравнений I: Нежесткие проблемы. 2-е исправленное издание, Спрингер, 2009.

[2] “Плеяды”. Википедия, 21 июня 2021. Википедия, https://en.wikipedia.org/wiki/Pleiades.

Локальные функции

Этот раздел включает локальные функции, вызванные решателем ОДУ, чтобы вычислить решение. В качестве альтернативы можно сохранить эти функции как их собственные файлы в директории на пути MATLAB.

function dqdt = pleiades(t,q)
  x = q(1:7);
  y = q(8:14);
  xDist = (x - x.');
  yDist = (y - y.');
  r = (xDist.^2+yDist.^2).^(3/2);
  m = (1:7)';
  dqdt = [q(15:28);
          sum(xDist.*m./r,1,'omitnan').';
          sum(yDist.*m./r,1,'omitnan').'];
end
%-----------------------------------------------------------------
function AnimateOrbits(t,q)
for k = 1:length(t)
    plot(q(:,1),q(:,8),'--',q(:,2),q(:,9),'--',...
         q(:,3),q(:,10),'--',q(:,4),q(:,11),'--',...
         q(:,5),q(:,12),'--',q(:,6),q(:,13),'--',...
         q(:,7),q(:,14),'--')
    hold on
    xlim([-20 20])
    ylim([-10 10])
    sz = 15;
    plot(q(k,1),q(k,8),'o','MarkerSize',sz,'MarkerFaceColor','r')
    plot(q(k,2),q(k,9),'o','MarkerSize',sz,'MarkerFaceColor','k')
    plot(q(k,3),q(k,10),'o','MarkerSize',sz,'MarkerFaceColor','b')
    plot(q(k,4),q(k,11),'o','MarkerSize',sz,'MarkerFaceColor','m')
    plot(q(k,5),q(k,12),'o','MarkerSize',sz,'MarkerFaceColor','c')
    plot(q(k,6),q(k,13),'o','MarkerSize',sz,'MarkerFaceColor','y')
    plot(q(k,7),q(k,14),'o','MarkerSize',sz,'MarkerFaceColor',[210 120 0]./255)
    hold off
    drawnow
    M(k) = getframe(gca);
    im{k} = frame2im(M(k));
end

filename = "orbits.gif";
for idx = 1:length(im)
    [A,map] = rgb2ind(im{idx},256);
    if idx == 1
        imwrite(A,map,filename,'gif','LoopCount',Inf,'DelayTime',0);
    else
        imwrite(A,map,filename,'gif','WriteMode','append','DelayTime',0);
    end
end
close all
end