在 matlab 中解决“NaN”的建议 Matlab中大数和小数的处理

如何解决在 matlab 中解决“NaN”的建议 Matlab中大数和小数的处理

我正在尝试使用 Matlab 在 3d 中绘制行星运动模型。我用牛顿定律和两个物体之间的引力得到了下面的微分方程：

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' 数组是行星质量。

然后我用数值方法Runge-Kutta-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,:);
    y((i*10)+1,:);
    for j=2:10
        k1=F(t(i),y(((i-1)*10)+j,:),j);
        k2=F(t(i)+h/2,:)+(h/2).*k1',j);
        k3=F(t(i)+h/2,:)+(h/2).*k2',j);
        k4=F(t(i)+h,:)+h.*k3',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,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,2);
    z(i,3);
end

plot3(x,z)

最后，结果根本不令人满意，我得到了很多“NAN”，然后我对 RK4 方法进行了一些调整并开始得到数字，但是当我绘制它们时，结果发现我正在绘制一条线而不是轨道。

任何帮助将不胜感激。提前致谢。

解决方法

两个错误：一个物理：公式中的alpha是代码中的j，公式中的运行索引j是公式中的循环索引i。总的来说，这会导致符号错误，将引力转化为电子之间的排斥力。因此，物理学规定，只要物体的路径不相交，物体就会几乎线性地彼此远离。

其次，您正在应用 RK4 方法，使其总体上是 1 阶方法。这些也往往会很快表现出非身体行为。您需要首先在临时 StagePos 变量中将所有位置更新到第一阶段，然后使用它来计算第二阶段的所有位置更新等。每个步骤与当前实现的差异可能很小，但是这样系统误差快速总结。