我正在嘗試用七個微分方程求解一個系統。而且我很難掌握 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);
任何幫助將非常感激
uj5u.com熱心網友回復：
您的代碼存在兩個基本問題。第一個問題是您沒有使用 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），以便在每一步都可以重新規范化四元數元素。
uj5u.com熱心網友回復：
我不確定您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)，否則代碼會創建一個列向量。沒有它，您可能會遇到問題。
轉載請註明出處，本文鏈接：https://www.uj5u.com/ruanti/441046.html