在Matlab中將衍生值保存在ode45中

Question

我正在模擬具有質量彈簧和雙擺的（有點奇怪的）系統的運動方程，爲此我有一個質量矩陣和函數f（x），並調用ode45解決

M*x' = f(x,t); 

我有5個狀態變量，q = [QDOT，PHI，phiDot，R，rDot]'; （刪除Q因爲沒有什麼依賴於它，QDot是最新的。） 現在，爲了計算一些力量，我還想保存rDotDot的計算值，ode45計算每個積分步驟的計算值，但是ode45不會給出這個結果。我已經搜索了一下，但是我發現的唯一的兩個解決方案是 a）將它變成一個三階問題，並將phiDotDot和rDotDot添加到狀態向量中。我想盡可能地避免這種情況，因爲它已經是非線性的，這實際上使事情變得更糟，並縮短了計算時間。

b）增加狀態直接計算函數，如here所述。然而，在他說的例子中，在質量矩陣中添加了一行零。它是有道理的，因爲否則它將整合導數，而不僅僅是在一點上評估它，但另一方面它使質量矩陣單數。似乎沒有爲我工作...

這似乎是這樣一個基本的東西（想要的狀態向量的派生值），有沒有很明顯，我沒有想到？ （或者不太明顯的東西也可以......）

哦，全局變量不是很好，因爲ode45幾次調用f（）函數，同時細化它的步驟，所以全局變量的大小和返回的狀態向量q完全不匹配。

萬一有人需要它，下面的代碼：

%(Initialization of parameters are above this line) 
    options = odeset('Mass',@massMatrix); 
    [T,q] = ode45(@f, tspan,q0,options); 

function dqdt = f(t,q,p) 
    % q = [qDot phi phiDot r rDot]'; 

    dqdt = zeros(size(q)); 

    dqdt(1) = -R/L*q(1) - kb/L*q(3) +vs/L; 
    dqdt(2) = q(3); 
    dqdt(3) = kt*q(1) + mp*sin(q(2))*lp*g; 
    dqdt(4) = q(5); 
    dqdt(5) = mp*lp*cos(q(2))*q(3)^2 - ks*q(4) - (mb+mp)*g; 
end 

function M = massMatrix(~,q) 
    M = [ 
     1 0 0 0 0; 
     0 1 0 0 0; 
     0 0 mp*lp^2 0 -mp*lp*sin(q(2)); 
     0 0 0 1 0; 
     0 0 mp*lp*sin(q(2)) 0 (mb+mp) 
     ]; 
end 

Answer 1

最簡單的方法是在每個由ode45返回值的，只是重新運行功能。

困難的解決辦法是嘗試將DotDots記錄到其他位置（預分配矩陣或甚至外部文件）。問題是，如果ode45祕密地在奇怪的地方進行評估，最終可能會得到不需要的數據點。

Answer 2

由於您使用的是嵌套函數，因此您可以使用它們的範圍規則來模擬全局變量的行爲。

這是最簡單的（AB）使用output function用於此目的：

function solveODE 

    % ....   
    %(Initialization of parameters are above this line) 

    % initialize "global" variable 
    rDotDot = []; 

    % Specify output function 
    options = odeset(... 
     'Mass', @massMatrix,... 
     'OutputFcn', @outputFcn); 

    % solve ODE 
    [T,q] = ode45(@f, tspan,q0,options); 

    % show the rDotDots  
    rDotDot 



    % derivative 
    function dqdt = f(~,q) 

     % q = [qDot phi phiDot r rDot]'; 

     dqdt = [... 
      -R/L*q(1) - kb/L*q(3) + vs/L 
      q(3) 
      kt*q(1) + mp*sin(q(2))*lp*g 
      q(5) 
      mp*lp*cos(q(2))*q(3)^2 - ks*q(4) - (mb+mp)*g 
     ]; 

    end % q-dot function 

    % mass matrix 
    function M = massMatrix(~,q) 
     M = [ 
      1 0 0 0 0; 
      0 1 0 0 0; 
      0 0 mp*lp^2 0 -mp*lp*sin(q(2)); 
      0 0 0 1 0; 
      0 0 mp*lp*sin(q(2)) 0 (mb+mp) 
     ]; 
    end % mass matrix function 


    % the output function collects values for rDotDot at the initial step 
    % and each sucessful step 
    function status = outputFcn(t,q,flag) 

     status = 0; 

     % at initialization, and after each succesful step 
     if isempty(flag) || strcmp(flag, 'init') 
      deriv = f(t,q); 
      rDotDot(end+1) = deriv(end); 
     end 

    end % output function 

end 

輸出功能僅計算在最初的衍生物和所有成功的步驟，所以它基本上是做一樣的東西阿德里安Ratnapala建議;在ode45的每個輸出處重新運行導數;我認爲這會更優雅（阿德里安+1）。

輸出函數方法的優點是您可以在集成運行時繪製rDotDot值（不要忘記drawnow！），這對於長時間運行的集成非常有用。