Matlab中的循环错误，计算位置

Question

目标：

在这个项目中，我们计算了太阳系中每个物体对系统中其他物体施加的引力。基于这些力和物体的初始位置/速度，可以使用牛顿第二定律预测它们的运动。我们假设给出以下信息：

• 太阳系中所有天体的质量。
• 所有行星在给定时间的位置和速度。

然后我们可以使用牛顿万有引力定律来计算万有引力。每个物体都受到来自所有其他物体的重力的作用，因此该物体上的总力是所有其他物体所感应的力的总和。一旦知道物体上的力，就可以使用牛顿第二运动定律确定每个物体的加速度。基于任何时间 t 的加速度，我们可以计算出物体在时间 t + Δt 的新速度和位置。我们从 t = 0 时刻的已知位置和速度开始，然后计算 Δt 时的位置和速度，然后计算 2Δt 时的位置和速度，依此类推。

分配：

编写一个 MATLAB 函数，输入一个结构数组（其中每个元素代表太阳系中的一个物体）、一个时间步长 Δt 和一个最终时间 T。结构数组中的每个结构都应具有以下字段：

• name：包含行星名称的字符串。
• x：行星初始位置的 x 坐标。
• y：行星初始位置的 y 坐标。
• z：行星初始位置的 z 坐标。
• vx：行星初始速度的 x 分量。
• vy：行星初始速度的 y 分量。
• vz：行星初始速度的 z 分量。

那么你的函数应该使用上面描述的引力系统模型来计算所有行星的位置和速度。您的函数应该以适当的方式返回这些值，并在 3D 图中生成所有行星的位置图。情节至少应包含标题和图例。图例应使用字段名称来命名每个行星。此外，每个行星应该以不同的颜色显示，其中颜色是随机确定的。示例图如图 1 所示。

我的错误

我的最后一个循环中有一个错误，其中包含n=1:T/t; 关于矩阵尺寸。我正在尝试使用每次迭代的更新位置来计算新的力，但是我得到了错误，我的绘图显示了一堆直线，而不是太阳系轨道。**

function [x, y, z, vx, vy, vz] = solarsystemsim(F,t,T)

m = size(F,2);
G = 6.67384e-11;
j=1;
x = 1:10;
% r = zeros(m-1,1);
% gF = zeros(m-1,1);
for(i=1:m)
    for(j=1:m)
        r(i,j) = distform2(F(i).x,F(j).x,F(i).y,F(j).y,F(i).z,F(j).z);
        gF(i,j) = (G*(F(i).mass)*(F(j).mass))/((r(i,j))^2);
        gFx(i,j) = -((gF(i,j))*(F(i).x-F(j).x))/(r(i,j));
        gFy(i,j) = -((gF(i,j))*(F(i).y-F(j).y))/(r(i,j));
        gFz(i,j) = -((gF(i,j))*(F(i).z-F(j).z))/(r(i,j));

        if(i==j)
            gF(i,j) = 0; 
            gFx(i,j) = 0;
            gFy(i,j) = 0;
            gFz(i,j) = 0;
        end
    end
end

for(i=1:m)
    gFxT(i) = sum(gFx(i,:)); 
    gFyT(i) = sum(gFy(i,:)); 
    gFzT(i) = sum(gFz(i,:)); 
end


tn = 0;
n = 1;

x = zeros(T/t,m);
y = zeros(T/t,m);
z = zeros(T/t,m);
vx = zeros(T/t,m);
vy = zeros(T/t,m);
vz = zeros(T/t,m);

for(i=1:m)
    vx(1,i) = F(i).vx;
    vy(1,i) = F(i).vy;
    vz(1,i) = F(i).vz;
    x(1,i) = F(i).x;
    y(1,i) = F(i).y;
    z(1,i) = F(i).z;
end

for(n=1:T/t)
    for(i=1:m)
        for(j=1:m)
            r(i+1,j+1) = distform2(x(i,j),x(i,j),y(i,j),y(i,j),z(i,j),z(i,j));
            gF(i+1,j+1) = (G*(F(i).mass)*(F(j).mass))/((r(i,j))^2);
            gFx(i+1,j+1) = -((gF(i,j))*(x(i,j)-x(i,j+1)))/(r(i,j));
            gFy(i+1,j+1) = -((gF(i,j))*(y(i,j)-y(i,j+1)))/(r(i,j));
            gFz(i+1,j+1) = -((gF(i,j))*(z(i,j)-z(i,j+1)))/(r(i,j));

            if(i==j)
                gF(i,j) = 0; 
                gFx(i,j) = 0;
                gFy(i,j) = 0;
                gFz(i,j) = 0;
            end
        end
    end

gFxT(i) = sum(gFx(i,:)); 
gFyT(i) = sum(gFy(i,:)); 
gFzT(i) = sum(gFz(i,:)); 

    for(i=1:T/t)
        for(j=1:m)
            vx(i,j) = vx(i,j) + t*(gFxT(j))/(F(j).mass);
            vy(i,j) = vy(i,j) + t*(gFyT(j))/(F(j).mass);
            vz(i,j) = vz(i,j) + t*(gFzT(j))/(F(j).mass);
            x(i,j) = x(i,j) + t*(vx(j));
            y(i,j) = y(i,j) + t*(vy(j));
            z(i,j) = z(i,j) + t*(vz(j));
        end
    end 

plot3(x,y,z);
end

Answer 1

F 不是矩阵，因此您不能使用括号对其进行索引。它是一个元胞数组，必须使用大括号索引 {}

已根据 lmillefiori 的评论将结构更正为单元格