Getting NaN in Rung-Kutta 4, Simulation of 10 planets

2 Ansichten (letzte 30 Tage)
Bashar Al-Mohammad
Bashar Al-Mohammad am 22 Jan. 2021
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
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?
Bashar Al-Mohammad
Bashar Al-Mohammad am 22 Jan. 2021
Thanks for being here,
I use A=10^10m, for the
When using A for distance G became in ordre E-41, Since G and the mass of planet is multiplied by each other I simplified the numbers. Ex: Mass of the Earth: 5.97219E+24 after using A for distance and I multiplied by E-41 then the mass of the Earth: 5.97219E-17.
Thanks again

Melden Sie sich an, um zu kommentieren.

Antworten (1)

James Tursa
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
Bashar Al-Mohammad
Bashar Al-Mohammad am 23 Jan. 2021
The new and improve RK4:
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
The new and improved (dY/dt):
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
The new and a bit satisfying results:
>> PlanetaryTest
>> yy
yy =
1.0e+12 *
0.0002 0.0010 -0.0000 -0.0000 0.0000 0.0000
-0.0568 -0.0274 0.0029 0.0000 -0.0000 -0.0000
0.0428 0.1001 -0.0011 -0.0000 0.0000 0.0000
-0.1438 -0.0400 -0.0000 0.0000 -0.0000 0.0000
-0.1147 -0.1963 -0.0013 0.0000 -0.0000 -0.0000
-0.5669 -0.5775 0.0151 0.0000 -0.0000 -0.0000
0.0821 -1.5024 0.0229 0.0000 0.0000 -0.0000
2.6251 1.4027 -0.0288 -0.0000 0.0000 0.0000
4.3030 -1.2422 -0.0736 0.0000 0.0000 -0.0000
1.6555 -4.7350 0.0278 0.0000 0.0000 -0.0000
0.0002 0.0010 -0.0000 -0.0000 0.0000 0.0000
-0.0568 -0.0274 0.0029 0.0000 -0.0000 -0.0000
0.0428 0.1001 -0.0011 -0.0000 0.0000 0.0000
-0.1438 -0.0400 -0.0000 0.0000 -0.0000 0.0000
-0.1147 -0.1963 -0.0013 0.0000 -0.0000 -0.0000
-0.5669 -0.5775 0.0151 0.0000 -0.0000 -0.0000
0.0821 -1.5024 0.0229 0.0000 0.0000 -0.0000
2.6251 1.4027 -0.0288 -0.0000 0.0000 0.0000
4.3030 -1.2422 -0.0736 0.0000 0.0000 -0.0000
1.6555 -4.7350 0.0278 0.0000 0.0000 -0.0000
0.0002 0.0010 -0.0000 -0.0000 0.0000 0.0000
-0.0568 -0.0274 0.0029 0.0000 -0.0000 -0.0000
0.0428 0.1001 -0.0011 -0.0000 0.0000 0.0000
-0.1438 -0.0400 -0.0000 0.0000 -0.0000 0.0000
-0.1147 -0.1963 -0.0013 0.0000 -0.0000 -0.0000
-0.5669 -0.5775 0.0151 0.0000 -0.0000 -0.0000
0.0821 -1.5024 0.0229 0.0000 0.0000 -0.0000
2.6251 1.4027 -0.0288 -0.0000 0.0000 0.0000
4.3030 -1.2422 -0.0736 0.0000 0.0000 -0.0000
1.6555 -4.7350 0.0278 0.0000 0.0000 -0.0000
0.0002 0.0010 -0.0000 -0.0000 0.0000 0.0000
-0.0568 -0.0274 0.0029 0.0000 -0.0000 -0.0000
0.0428 0.1001 -0.0011 -0.0000 0.0000 0.0000
-0.1438 -0.0400 -0.0000 0.0000 -0.0000 0.0000
-0.1147 -0.1963 -0.0013 0.0000 -0.0000 -0.0000
-0.5669 -0.5775 0.0151 0.0000 -0.0000 -0.0000
0.0821 -1.5024 0.0229 0.0000 0.0000 -0.0000
2.6251 1.4027 -0.0288 -0.0000 0.0000 0.0000
4.3030 -1.2422 -0.0736 0.0000 0.0000 -0.0000
1.6555 -4.7350 0.0278 0.0000 0.0000 -0.0000
In this time I used (m),(s),(kg). and I got some numbers. I am too tired to plot these but as soon as I do, I will get backk to here and tell. I apperciate your time Mr. Tursa. I aslo appreciate any improvement suggestions on the RK4 function. Thanks again.
Bashar Al-Mohammad
Bashar Al-Mohammad am 23 Jan. 2021
Update plotting went wrong. I extracted the x y z coordinates to plot the earth and I got a straight line. 'h' step size used is 1 day variation.

Melden Sie sich an, um zu kommentieren.

Kategorien

Mehr zu Programming finden Sie in Help Center und File Exchange

Tags

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by