Блог

Я пытаюсь построить модель движения планет в 3D с использованием Matlab. Я использовал закон Ньютона с гравитационной силой между двумя объектами, и я получил дифференциальное уравнение ниже:

код matlab:

 function dy=F(t,y,CurrentPos,j)
m=[1.98854E 30 3.302E 23 4.8685E 24 5.97219E 24 6.4185E 23 1.89813E 27 5.68319E 26 8.68103E 25 1.0241E 26 1.307E 22];
G=6.67E-11;
dy = zeros(6,1);
dy(1) = y(4);
dy(2) = y(5);
dy(3) = y(6);
for i=1:10
    if i~=j
        deltaX=(CurrentPos(j,1)-CurrentPos(i,1));
        deltaY=(CurrentPos(j,2)-CurrentPos(i,2));
        deltaZ=(CurrentPos(j,3)-CurrentPos(i,3));
        ray=sqrt((deltaX^2) (deltaY^2) (deltaZ^2));
        dy(4) = dy(4)   G*m(i)*(deltaX/(ray^3));
        dy(5) = dy(5)   G*m(i)*(deltaY/(ray^3));
        dy(6) = dy(6)   G*m(i)*(deltaZ/(ray^3));
    end
end
 

где массив «m» — это массы планет.

затем я использовал численный метод Рунге-Кутта-4 для его решения, и вот код:

 function [y,t]=RK4(F,intPos,a,b,N)

h=(b-a)/N;
t=zeros(N,1);
y = zeros(10*N,6);
y(1,:)=intPos(1,:);
y(2,:)=intPos(2,:);
y(3,:)=intPos(3,:);
y(4,:)=intPos(4,:);
y(5,:)=intPos(5,:);
y(6,:)=intPos(6,:);
y(7,:)=intPos(7,:);
y(8,:)=intPos(8,:);
y(9,:)=intPos(9,:);
y(10,:)=intPos(10,:);

t(1)=a;

for i=1:N
    
    t(i 1)=a i*h;
    CurrentPos=y((i*10)-9:i*10,:);
%     CurrentPos(1,:)=intPos(1,:);
    y((i*10) 1,:)=intPos(1,:);
    for j=2:10
        k1=F(t(i),y(((i-1)*10) j,:),CurrentPos,j);
        k2=F(t(i) h/2,y(((i-1)*10) j,:) (h/2).*k1',CurrentPos,j);
        k3=F(t(i) h/2,y(((i-1)*10) j,:) (h/2).*k2',CurrentPos,j);
        k4=F(t(i) h,y(((i-1)*10) j,:) h.*k3',CurrentPos,j);
        y((i*10) j,:)=y(((i-1)*10) j,:) (h/6)*(k1 2*k2 2*k3 k4)';
    end
end
 

Наконец-то применил функцию для начальных состояний из системы JPL HORIZONS:

 format short

intPos=zeros(10,6);
intPos(1,:)=[1.81899E 08 9.83630E 08 -1.58778E 07 -1.12474E 01 7.54876E 00 2.68723E-01];
intPos(2,:)=[-5.67576E 10 -2.73592E 10 2.89173E 09 1.16497E 04 -4.14793E 04 -4.45952E 03];
intPos(3,:)=[4.28480E 10 1.00073E 11 -1.11872E 09 -3.22930E 04 1.36960E 04 2.05091E 03];
intPos(4,:)=[-1.43778E 11 -4.00067E 10 -1.38875E 07 7.65151E 03 -2.87514E 04 2.08354E 00];
intPos(5,:)=[-1.14746E 11 -1.96294E 11 -1.32908E 09 2.18369E 04 -1.01132E 04 -7.47957E 02];
intPos(6,:)=[-5.66899E 11 -5.77495E 11 1.50755E 10 9.16793E 03 -8.53244E 03 -1.69767E 02];
intPos(7,:)=[8.20513E 10 -1.50241E 12 2.28565E 10 9.11312E 03 4.96372E 02 -3.71643E 02];
intPos(8,:)=[2.62506E 12 1.40273E 12 -2.87982E 10 -3.25937E 03 5.68878E 03 6.32569E 01];
intPos(9,:)=[4.30300E 12 -1.24223E 12 -7.35857E 10 1.47132E 03 5.25363E 03 -1.42701E 02];
intPos(10,:)=[1.65554E 12 -4.73503E 12 2.77962E 10 5.24541E 03 6.38510E 02 -1.60709E 03];

[yy,t]=RK4(@F,intPos,0,1e8,1e3);
x=zeros(101,1);
y=zeros(101,1);
z=zeros(101,1);
for i=1:1e3
    x(i,:)=yy((i-1)*10 4,1);
    y(i,:)=yy((i-1)*10 4,2);
    z(i,:)=yy((i-1)*10 4,3);
end

plot3(x,y,z)
 

Наконец, результат совсем не удовлетворил, и я получил много «NAN», затем я немного скорректировал метод RK4 и начал получать числа, но когда я построил их, оказалось, что я строю линию вместо орбиты.

Любая помощь будет оценена. Заранее спасибо.

Две ошибки: одна физическая: альфа в формуле — это j значение в коде, текущий индекс j в формулах — это индекс цикла i в формуле. В целом это приводит к ошибке знака, преобразуя силу притяжения силы тяжести в силу отталкивания, например, между электронами. Таким образом, физика диктует, что тела удаляются друг от друга почти линейно, пока их пути не пересекаются.

Во-вторых, вы применяете метод RK4 таким образом, что в целом это метод порядка 1. Они также имеют тенденцию вести себя не физически довольно быстро. Вам нужно сначала обновить все позиции до первого этапа во временной StagePos переменной, затем использовать это для вычисления всех обновлений позиций для второго этапа и т.д. Разница с текущей реализацией может быть небольшой на каждом шаге, но такие систематические ошибки быстро суммируются.