求解一阶方程的 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

我不确定你的 inital_val 应该代表什么，但在这里你至少可以运行以下片段和 fix/change 相应地。

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

Answer 2

您的代码存在两个基本问题。第一个问题是您没有使用 ode45( ) 所需的导数函数的正确签名。第二个问题是，以这种方式盲目地合并四元数元素将导致 non-unit 个四元数。我不知道在使用 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

此外，您的 angular 速率初始值对于旋转问题似乎完全不合理。它们将被解释为 [1 -1 2] radians/second。我怀疑您为此想要 degrees/second，因此您需要缩小这些值。

第二个问题，即保持四元数元素形成单位四元数的问题，使用 ode45( ) 不容易克服。所有的小积分错误都会使四元数的大小逐渐远离 1。因此，在您的其他代码中对该四元数的任何后续使用都会因此而出现问题。我唯一的建议是完全放弃 ode45( ) 并编写您自己的自定义积分器（例如 RK4），以便在每一步都可以 re-normalize 四元数元素。

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

Solving 7 equations of first order equations ode45 Matlab

matlab

differential-equations

ode45