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.

question from:https://stackoverflow.com/questions/65851290/suggestion-to-solve-nan-in-matlab-dealing-with-large-and-small-numbers-in-mat

Answer 1

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.