Matlab中的摆动弹簧ODE系统 - 如何使位置矢量跟随路径？

Question

我对 Matlab 很陌生。 我有一个使用ode45 和arrow.m 的脚本来显示Matlab 上带有质量的摆动弹簧在穿过3-D 空间时的运动。该程序几乎在做我想做的事。现在，菱形的密度显示了弹簧的有效速度（ode45 采用个人最喜欢的数字样本时除外），并且速度几乎可以通过函数的步长（至少使用我的计算机运行代码的速度）。我想要做的是让我在代码中注释掉的位置矢量只显示在质量的瞬时位置，即曲线的端点，而不是钻石出现的每个点。我四处寻找帮助，但我尝试的一切似乎都会导致错误。如果有人能指出我正确的方向，将不胜感激。请尝试运行程序，你会明白我的意思，也可以玩弄函数的参数。

function elasticPendulum(totalTime)
time=1;
totalTime=20
X=1
Vx=0
Y=0
Vy=0
Z=1.1
Vz=0
w=3;
g=9;
l=1;
while (time<totalTime)
    tspan=[0,time];
    x0=[X,Vx,Y,Vy,Z,Vz];
    [t,x]=ode45(@F,tspan,x0);
    plot3(x(:,1),x(:,3),x(:,5),'dk'), axis([-2 2 -2 2 -10 2]);
    grid on;
    o=[0, 0, 0];
    Xeq=[0, 0, -(g/(w^2)+l)];
    arrow(o,Xeq,'Length',5);
    %{ 
    Xt=[x(:,1) x(:,3) x(:,5)]
    arrow(o,Xt,'Length',5);
    %}
    time=time+.1*(((x(2))^2+(x(4))^2+(x(6))^2)^(1/2))
    pause(.1);
end

    function xprime=F(t,x)
    r=sqrt(x(1)^2+x(3)^2+x(5)^2);
    xprime=[x(2);-w*(r-l)/r*x(1);x(4);-w*(r-l)/r*x(3);x(6);-w*(r-l)/r*x(5)-g];
    end

end

Answer 1

我认为您只需将Xt 设置为此即可完成您想要的操作，这样每次迭代时只绘制最后一个向量：

Xt = [x(end,1) x(end,3) x(end,5)];

此外，您提到“除了 ode45 采用个人最喜欢的数字样本时”，一切正常。您可以通过更改tspan 来告诉ode45 使用固定步长输出：

tspan = 0:dt:time;

dt = 1/time（或者你可以使用linspace）。这与使用固定步长求解器不同 - 阅读 my answer to this question 可能会对此有所了解。

我认为您没有正确更新动画。您更新时间，但不更新初始条件。因此，您在 while 循环的每次迭代中都通过相同的路径进行集成。您是否有理由在每次迭代中调整结束时间？你是想找出振荡周期还是什么？

此外，您的动画制作方法相当粗糙/效率低下。您应该阅读有关可以通过 odeset 为 Matlab 的 ODE 求解器设置的 output options 的信息。特别是'OutputFcn'。您甚至可以使用和/或查看某些内置输出函数的代码——例如，在命令窗口中输入edit odeplot 或edit odephas3。这是我为您的案例创建的一个简单的输出函数：

function elasticPendulum
totalTime = 20;
stepsPerSec = 10;
X = 1; Vx = 0;
Y = 0; Vy = 0;
Z = 1.1; Vz = 0;
w = 3; g = 9; l = 1;

opts = odeset('OutputFcn',@(t,x,flag)arrowplot(t,x,flag,w,l,g));
tspan = linspace(0,totalTime,totalTime*stepsPerSec);
x0 = [X,Vx,Y,Vy,Z,Vz];
ode45(@(t,x)F(t,x,w,l,g),tspan,x0,opts);

function xprime=F(t,x,w,l,g)    %#ok<INUSL>
r=sqrt(x(1)^2+x(3)^2+x(5)^2);
xprime=[x(2);-w*(r-l)/r*x(1);x(4);-w*(r-l)/r*x(3);x(6);-w*(r-l)/r*x(5)-g];

function status=arrowplot(t,x,flag,w,l,g)   %#ok<INUSL>
persistent h; % Handle for moving arrow
status = 0;
o = [0, 0, 0];
switch flag
    case 'init'
        axis([-2 2 -2 2 -10 2]); grid on; hold on;
        Xeq = [0, 0, -(g/(w^2)+l)];
        arrow(o,Xeq,'Length',5);
        plot3(x(1,:),x(3,:),x(5,:),'dk');  % Initial position
        Xt = [x(1,end) x(3,end) x(5,end)];
        h = arrow(o,Xt,'Length',5);        % Initial arrow, get handle
    case 'done'
        hold off; status = 1;
    otherwise
        plot3(x(1,:),x(3,:),x(5,:),'dk');        % Plot new positions
        Xt = [x(1,end) x(3,end) x(5,end)];
        arrow(h,'Start',o,'Stop',Xt,'Length',5); % Update arrow
        pause(0.1);
end

在每次迭代时调用plot3（arrowplot 中的otherwise 案例）仍然效率低下，因为它添加了新的高级绘图对象，占用了更多内存。更好/更快的方法是从第一次调用plot3 中获取句柄，并使用set 和get 添加新数据点。请参阅odeplot 和odephas3 的代码，了解如何执行此操作（有点高级）。

请注意我是如何通过匿名函数而不是通过创建子函数来传递参数的。这有点风格问题。