Matlab中的擺動彈簧ODE系統 - 如何使位置矢量跟隨路徑？

Question

我對Matlab很陌生。 我有一個腳本，使用ode45和arrow.m來顯示一個擺動的彈簧在Matlab中的質量隨着它在三維空間中移動的運動。該方案是差不多做我想做的。現在鑽石的密度顯示出有效的彈簧速度（除了當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取個人最喜歡的數字樣本。「你可以告訴ode45簡單地通過改變tspan使用固定步輸出：

tspan = 0:dt:time; 

其中dt = 1/time（或者你可以使用linspace）。這與使用固定步長求解器不同 - 閱讀my answer to this question可能對此有所瞭解。

我不認爲你正在更新你的動畫。您更新時間，但不是最初的條件。因此，您在while循環的每次迭代中集成相同的路徑。是否有理由在每次迭代中調整結束時間？你是否想找到振盪週期或什麼？

此外，您的動畫方法是相當粗糙/效率低下。您應該閱讀關於output options的信息，可以通過odeset來設置Matlab的ODE解算器。特別是'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的呼叫（將otherwise情況arrowplot）仍然是低效的，因爲它增加了新的高層次情節的對象，佔用更多的內存。更好/更快的方法是從第一次撥打電話到plot3，並使用set和get來添加新的數據點。請參閱odeplot和odephas3的代碼，瞭解如何做到這一點（這有點高級）。

請注意我是如何通過匿名函數而不是通過創建子函數來傳遞參數的。這是一個風格問題。