求解7个一阶方程组ode45Matlab

Question

我正试图求解一个有七个微分方程的系统。我很难掌握ode45解算器
这些是方程式：

ω2_dot=-0.75ω1ω3

ω1_dot=0.75ω2ω3+0.2

ω3_dot=0

q1_dot=1/2(ω1q4+ω2q3-ω3q2(

q2_dot=1/2(ω2q4+ω3q1-ω1q3(

q3_dot=1/2(ω3q4+ω1q2-ω2q2(

q4_dot=-1/2(ω1q1+ω2q2+ω3q3(

初始值按相同顺序列在inital_val中。以下是我目前所拥有的：

inital_val = [1 -1 2 0 0 0 1];
timespan = [0:0.05:3.95];
[result t] = ode45(@soe,timespan,inital_val);
function [results,t]=soe(inital_val, timespan)
omega1_dot = -0.75*omega2*omega3;
omega2_dot = 0.75*omega2*omega3+0.2;
omega3_dot = 0;
q1_dot = (1/2)*(q4*omega1+omega2*q3-omega3*q2);
q2_dot = (1/2)*(q4*omega2+omega3*q1-omega1*q3);
q3_dot = (1/2)*(q4*omega3+omega1*q2-omega2*q2);
q4_dot = (1/2)*(q1*omega1+q2*omega2+q3*omega3);
results = [omega1 omega2 omega3 q1 q2 q3 q4];
end

然而，它给了我这个错误消息，这是有道理的，但我不知道如何修复它：

Unrecognized function or variable 'omega2'.
Error in ode45_part>soe (line 10)
omega1_dot = -0.75*omega2*omega3;
Error in odearguments (line 90)
f0 = feval(ode,t0,y0,args{:});   % ODE15I sets args{1} to yp0.
Error in ode45 (line 115)
odearguments(FcnHandlesUsed, solver_name, ode, tspan, y0, options, varargin);
Error in ode45_part (line 7)
[result t] = ode45(@soe,timespan,inital_val);

如有任何帮助，将不胜感激

Answer 1

您的代码存在两个基本问题。第一个问题是，您没有用于ode45((所需的导数函数的正确签名。第二个问题是，以这种方式盲目地积分四元数元素会导致非单位四元数。在使用ode45((时，我不知道有什么简单的方法可以强制执行此限制。

让我们先来解决一个简单的问题，导数函数签名。此外，您在该函数中使用变量也是不正确的。它需要这样：

function dy = soe(t,y) % Required signature is (time,state_vector)
% Pick off the named variables from the state vector y
omega1 = y(1);
omega2 = y(2);
omega3 = y(3);
q1 = y(4);
q2 = y(5);
q3 = y(6);
q4 = y(7);
% Calculate the variable derivatives
omega1_dot = -0.75*omega2*omega3; % IS THIS A TYPO? Should it be omega1*omega3?
omega2_dot = 0.75*omega2*omega3+0.2;
omega3_dot = 0;
q1_dot = (1/2)*(q4*omega1+omega2*q3-omega3*q2);
q2_dot = (1/2)*(q4*omega2+omega3*q1-omega1*q3);
q3_dot = (1/2)*(q4*omega3+omega1*q2-omega2*q2);
q4_dot = (1/2)*(q1*omega1+q2*omega2+q3*omega3);
% Group variable derivatives into a single column vector
dy = [omega1dot; omega2dot; omega3dot; q1_dot; q2_dot; q3_dot; q4_dot];
end

此外，对于旋转问题，角速率的初始值似乎完全不合理。它们将被解释为[1-1 2]弧度/秒。我怀疑你想要度/秒，所以你需要缩小这些值。

第二个问题是保持四元数元素形成一个单位四元数，使用ode45((并不容易克服。所有的小积分误差将在驱动四元数幅值远离1时缓慢变化。任何后续在其他代码中使用此四元数都会因此而出现问题。在这方面，我唯一的建议是可能完全放弃ode45((，并编写自己的自定义积分器(例如，RK4(，以便在每一步都可以重新规范四元数元素。

Answer 2

我不确定inital_val应该代表什么，但在这里，您至少可以运行以下代码段并相应地进行修复/更改。

clear; clc;
y0 = [1 -1 2 0 0 0 1];
tspan = [0:0.05:3.95];
[t, y] = ode45(@(t, y) odefcn(t, y), tspan, y0);
function dydt = odefcn(t, y)
dydt = zeros(7, 1);
dydt(2) = -0.75 * y(1) * y(3);
dydt(1) = 0.75 * y(2) * y(3) + 0.2;
dydt(3) = 0;
dydt(4) = 1 / 2 * (y(1) * y(7) + y(2) * y(6) - y(3) * y(5));
dydt(5) = 1 / 2 * (y(2) * y(7) + y(3) * y(4) - y(1) * y(6));
dydt(6) = 1 / 2 * (y(3) * y(7) + y(1) * y(5) - y(2) * y(5));
dydt(7) = -1 / 2 * (y(1) * y(4) + y(2) * y(5) + y(3) * y(6));
end

请注意，您需要dydt = zeros(7, 1)，否则代码将创建一个列向量。没有它，你可能会遇到问题。