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

Posted 2023-02-17

【中文标题】在 matlab 中解决“NaN”的建议。在 Matlab 中处理大数和小数

【英文标题】：Suggestion to solve 'NaN' in matlab. Dealing with large and small numbers in Matlab

【发布时间】：2021-01-22 19:05:41

【问题描述】：

我正在尝试使用 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,:)=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

最终应用了 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,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)

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

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

【问题讨论】：

请将代码（和命令窗口结果）发布为文本，而不是图像。否则想尝试的人将不得不手动输入，因此您不太可能获得帮助

复制粘贴文本肯定比制作和上传屏幕截图要容易得多。我无法阅读这些图像。请复制粘贴代码并输出为文本。

代码已更新并作为文本上传。抱歉上传截图。

【参考方案1】：

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

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

【讨论】：

关于第二点，我正在传递所有行星的 currentPos，以便系统以正确的方式工作。或者我没有理解你的意思。关于第一点，符号问题将与 deltaX^2 一起出现。或者我也没有明白你的想法。

如果您计算第二阶段，所有演变的数据必须在时间t+h/2，在第四阶段t+h，如果使用插值或外推，则错误顺序为4。否则，您会得到一个在一段时间内看起来正确但从长远来看会迅速偏离、崩溃或爆炸的 order 1 方法。 CurrentPos 在时间 t+h 用于较小的索引，在时间 t 用于较大的索引。请注意，deltaX 也以第一次幂的形式进入力计算，其中符号很重要。