求解一阶方程的7个方程 ode45 Matlab

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) ，否则代码会创建一个列向量。没有它，您可能会遇到问题。