在 Matlab 中保存 ode45 中的导数值

Question

我正在模拟具有质量弹簧和双摆的（有些奇怪的）系统的运动方程，为此我有一个质量矩阵和函数 f(x)，并调用 ode45 来求解

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

我有5个状态变量，q= [ QDot, phi, phiDot, r, rDot]'; （删除 Q 因为没有任何东西依赖它，QDot 是最新的。） 现在，为了计算一些力，我还想保存 rDotDot 的计算值，ode45 会为每个积分步骤计算该值，但 ode45 不会将其返回。我搜索了一下，但我发现的唯一两个解决方案是 a) 将其转换为 3 阶问题并将 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 

输出函数只计算初始和所有成功步骤的导数，所以它基本上和 Adrian Ratnapala 建议的一样；在ode45 的每个输出处重新运行导数；我认为这会更优雅（阿德里安+1）。

输出函数方法的优点是您可以在运行集成时绘制rDotDot 值（不要忘记drawnow！），这对于长时间运行的集成非常有用。