Suggestion to solve 'NaN' in matlab. Dealing with large and small numbers in Matlab

Question

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 below:

matlab code:

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

where the 'm' array is the planet masses.

then I used the numerical method Runge-Kutta-4 to solve it, and here's the code:

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

Finally applied the function for the Initial States from JPL HORIZONS System:

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)

Finally, the result wasn't satisfying at all and I got many 'NAN', then I did some adjustment on the RK4 method and started to get numbers, but when I plotted them it turned out I'm plotting a line instead of an orbit.

Any help would be appreciated. Thanks in advance.

Please post code (and command window results) as text, not as images. Otherwise people who want to try it will have to type it by hand, so it's unlikely you will get help — Luis Mendo, Jan 22 '21 at 20:18
Copy-pasting text must be much easier than making and uploading screenshots. I cannot read these images. Please copy-paste the code and output as text. — Cris Luengo, Jan 23 '21 at 01:52
code updated and uploaded as text. Sorry for uploading screenshots. — Bashar El-Mouhammad, Jan 23 '21 at 08:40

score 0 · Answer 1 · answered Jan 25 '21 at 09:27

0

Two errors: One physical: The alpha in the formula is the j in the code, the running index j in the formulas is the loop index i in the formula. In total this makes a sign error, transforming the attracting gravity force into a repelling force like between electrons. Thus the physics dictates that the bodies move away from each other almost linearly, as long as their paths don't cross.

Second, you are applying the RK4 method in such a way that in total it is an order 1 method. These also tend to behave un-physically rather quickly. You need to update first all positions to the first stage in a temporary StagePos variable, then use that to compute all position updates for the second stage etc. The difference to the current implementation may be small in each step, but such systematic errors quickly sum up.

answered Jan 25 '21 at 09:27

Lutz Lehmann

25,219
2
22
51

Concerning the second point, I am passing the currentPos for all planets so the system work the right way. Or maybe I didnt understand your point. Concerning the first point, the sign problem will go with the deltaX^2. Or also I didnt get your idea. – Bashar El-Mouhammad Jan 26 '21 at 13:02
If you compute the second stage, all evolving data has to be for the time `t+h/2`, in the 4th stage for `t+h`, with error order 4 if you use interpolation or extrapolation. Else you get an order 1 method that looks right for some time but in the longer term deviates rapidly, collapsing or exploding. `CurrentPos` is at time `t+h` for smaller indices and at time `t` for larger indices. Note that `deltaX` enters the force computation also in the first power, where the sign is significant. – Lutz Lehmann Jan 26 '21 at 13:08

Suggestion to solve 'NaN' in matlab. Dealing with large and small numbers in Matlab

1 Answers1