我正在模拟一个带有质量弹簧和双摆的（有点奇怪）系统的运动方程，我有一个质量矩阵和函数 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 http://www.mathworks.co.uk/matlabcentral/newsreader/view_thread/246590。然而，在这个例子中，他说要在质量矩阵中添加一行零。这是有道理的，因为否则它将对导数进行积分，而不仅仅是在一个点上对其求值，但另一方面它会使质量矩阵变得奇异。似乎对我不起作用...
这看起来像是一个基本的事情（想要状态向量的导数值），是否有一些我没有想到的非常明显的事情？ （或者一些不那么明显的东西也可以......）
哦，全局变量并不是那么好，因为 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
