Getting NaN in Rung-Kutta 4, Simulation of 10 planets
2 Ansichten (letzte 30 Tage)
Ältere Kommentare anzeigen
I am trying to make a model of planets' movement plot it in 3d using Matlab. I used Newton's law with the gravitational force between two objects and I got the differential equation which can be found in the images below.
then I used the numerical method Runge-Kutta-4 to solve it and here's the code:
function [y,t]=RK4(pos,a,b,N)
h=(b-a)/N;
y = zeros(N*10,6);
t = zeros(N,1);
y(1,:)=pos(1,:);
y(2,:)=pos(2,:);
y(3,:)=pos(3,:);
y(4,:)=pos(4,:);
y(5,:)=pos(5,:);
y(6,:)=pos(6,:);
y(7,:)=pos(7,:);
y(8,:)=pos(8,:);
y(9,:)=pos(9,:);
y(10,:)=pos(10,:);
t(1)=a;
for i=1:N
t(i+1)=a+i*h;
CurrentPos=y((10*i)-9:i*10,:);
for j=1:10
k1=F(t(i),y(i+j-1,:),CurrentPos,j);
k2=F(t(i)+h/2,y(i+j-1,:)+(h/2).*k1',CurrentPos,j);
k3=F(t(i)+h/2,y(i+j-1,:)+(h/2).*k2',CurrentPos,j);
k4=F(t(i)+h,y(i+j-1,:)+h.*k3',CurrentPos,j);
y(i+9+j,:)=y(i+j-1,:)+(h/6)*(k1+2*k2+2*k3+k4)';
end
end
Where F is the (dY/dt):
function dy=F(t,y,CurrentPos,j)
m=[1.98854E-11 3.302E-18 4.8685E-17 5.97219E-17 6.4185E-18 1.89813E-14 5.68319E-15 8.68103E-16 1.0241E-15 1.307E-19];
G=6.67;
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(i,1)-CurrentPos(j,1));
deltaY=(CurrentPos(i,2)-CurrentPos(j,2));
deltaZ=(CurrentPos(i,3)-CurrentPos(j,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
Finally applied the function for the Initial States from JPL HORIZONS System:
pos=zeros(10,6);
pos(1,:)=[1.81899E-2 9.83630E-2 -1.58778E-3 -1.12474E-11 7.54876E-10 2.68723E-11];
pos(2,:)=[-5.67576 -2.73592 2.89173E-1 1.16497E-6 -4.14793E-6 -4.45952E-7];
pos(3,:)=[4.28480 1.00073E+1 -1.11872E-10 -3.22930E-6 1.36960E-6 2.05091E-7];
pos(4,:)=[-1.43778E+1 -4.00067 -1.38875E-3 7.65151E-7 -2.87514E-6 2.08354E-10];
pos(5,:)=[-1.14746E+1 -1.96294E+1 -1.32908E-1 2.18369E-6 -1.01132E-6 -7.47957E-8];
pos(6,:)=[-5.66899E+1 -5.77495E+1 1.50755 9.16793E-7 -8.53244E-7 -1.69767E-8];
pos(7,:)=[8.20513 -1.50241E+2 2.28565 9.11312E-8 4.96372E-8 -3.71643E-8];
pos(8,:)=[2.62506E+2 1.40273E+2 -2.87982 -3.25937E-7 5.68878E-7 6.32569E-9];
pos(9,:)=[4.30300E+2 -1.24223E+2 -7.35857 1.47132E-7 5.25363E-7 -1.42701E-8];
pos(10,:)=[1.65554E+2 -4.73503E+2 2.77962 5.24541E-7 6.38510E-8 -1.60709E-7];
[yy,t]=RK4(pos,0,100,5);
2 Kommentare
James Tursa
am 22 Jan. 2021
Could you please use comments to indicate the units used for all of your numbers? E.g., why is G listed as 6.67 instead of 6.67e-11? What are the units of the masses and initial positions and velocities?
Antworten (1)
James Tursa
am 22 Jan. 2021
Bearbeitet: James Tursa
am 22 Jan. 2021
OK, I will take a look at things. But your attempt at "simplifying" your numbers by using a custom distance unit and apparently also custom mass unit only adds confusion IMO. I would advise abandoning this right away and go back to using standard units. In fact, that is the first thing I would do to check this is to convert everything to standard units and run your code as is to see if it is a units problem or a code logic problem.
2 Kommentare
Siehe auch
Produkte
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!