Not enough input arguments ode45

2 ビュー (過去 30 日間)
Aleyna 2023 年 11 月 22 日

function dy = rigid(t,y)
% a column vector
dy = zeros(3,1);
w1=sin(t)*exp(-t);
w2=sin(2*t)*exp(-2*t);
w3=sin(3*t)*exp(-3*t);
dy(1) =(1/cos(y(2)))* (cos(y(3))*w1-sin(y(3))*w2);
dy(2) =(sin(y(3))*w1+cos(y(3))*w2);
dy(3) =w3+(1/cos(y(2)))*(-w1*cos(y(3))*sin(y(2))+w2*sin(y(2))*sin(y(3)));
syms a b c ad bd cd om1 om2 om3 t real
%rotations for 1(a)-2(b)-3(c)
X=[1 0 0;0 cos(a) sin(a);0 -sin(a) cos(a)];
Y=[cos(b) 0 -sin(b);0 1 0;sin(b) 0 cos(b)];
Z=[cos(c) sin(c) 0;-sin(c) cos(c) 0;0 0 1];
%Last rotation 3(c)
L=Z
%Last 3(c) and middle rotation 2(b)
LM=Z*Y
%Last rotation angle derivative
dv1=[0; 0; cd]
%Middle rotation angle derivative
dv2=[0; bd; 0]
%First rotation angle derivative
w11=dv1;
w22=L*dv2;
w33=LM*dv3;
%summation
w=w11+w22+w33
[S] = equationsToMatrix([w], [ad, bd, cd])
%Convert to differential equation form
%kinematic equation
IS=simplify (inv(S))
%Calculate differential equations
%angular velocity vectors
om1=sin(t)*exp(-t);
om2=sin(2*t)*exp(-2*t);
om3=sin(3*t)*exp(-3*t);
%Differential equation
angd=IS*[om1 om2 om3]'
options = odeset('RelTol',1e-4,'AbsTol',[1e-4 1e-4 1e-5]);
[T,Y] = ode45(@rigid,[0 2],[0 0 0],options);
plot(T,Y(:,1),'*',T,Y(:,2),'.',T,Y(:,3),'+')
legend('a','b','c')
1 件のコメント-1 件の古いコメントを表示-1 件の古いコメントを非表示
Aleyna 2023 年 11 月 22 日
Question: Calculate Euler angles differential equations for 1(a)-2(b)-3(c) rotation sequence

サインインしてコメントする。

回答 (2 件)

Dyuman Joshi 2023 年 11 月 22 日
I am not sure how the different symbolic variables and expressions come into play here for solving the ODE.
Nonetheless, Any functions defined in a script must be at the end of the script.
Also, the call to the ODE solver (ode45 in this case) should be outside the ODE funciton.
syms a b c ad bd cd om1 om2 om3 t real
%rotations for 1(a)-2(b)-3(c)
X=[1 0 0;0 cos(a) sin(a);0 -sin(a) cos(a)];
Y=[cos(b) 0 -sin(b);0 1 0;sin(b) 0 cos(b)];
Z=[cos(c) sin(c) 0;-sin(c) cos(c) 0;0 0 1];
%Last rotation 3(c)
L=Z;
%Last 3(c) and middle rotation 2(b)
LM=Z*Y;
%Last rotation angle derivative
dv1=[0; 0; cd];
%Middle rotation angle derivative
dv2=[0; bd; 0];
%First rotation angle derivative
w11=dv1;
w22=L*dv2;
w33=LM*dv3;
%summation
w=w11+w22+w33;
[S] = equationsToMatrix([w], [ad, bd, cd]);
%Convert to differential equation form
%kinematic equation
IS=simplify (inv(S));
%Calculate differential equations
%angular velocity vectors
om1=sin(t)*exp(-t);
om2=sin(2*t)*exp(-2*t);
om3=sin(3*t)*exp(-3*t);
%Differential equation
angd=IS*[om1 om2 om3]';
%% Call to ODE
options = odeset('RelTol',1e-4,'AbsTol',[1e-4 1e-4 1e-5]);
[T,Y] = ode45(@rigid,[0 2],[0 0 0],options);
%% Plotting the ODE
plot(T,Y(:,1),'*',T,Y(:,2),'.',T,Y(:,3),'+')
legend('a','b','c')
%% Move the ODE function to the bottom/end of the script
function dy = rigid(t,y)
% a column vector
dy = zeros(3,1);
w1=sin(t)*exp(-t);
w2=sin(2*t)*exp(-2*t);
w3=sin(3*t)*exp(-3*t);
dy(1) =(1/cos(y(2)))* (cos(y(3))*w1-sin(y(3))*w2);
dy(2) =(sin(y(3))*w1+cos(y(3))*w2);
dy(3) =w3+(1/cos(y(2)))*(-w1*cos(y(3))*sin(y(2))+w2*sin(y(2))*sin(y(3)));
end
0 件のコメント-2 件の古いコメントを表示-2 件の古いコメントを非表示

サインインしてコメントする。

Sulaymon Eshkabilov 2023 年 11 月 22 日

Don't use the MATLAB's built in function names as a variable name - cd. Why to display all iterations that can be also skipped with ";". Here is the corrected code:
options = odeset('RelTol',1e-4,'AbsTol',[1e-4 1e-4 1e-5]);
[T,Y] = ode45(@rigid,[0 2],[0 0 0],options);
angd =
angd =
angd =
angd =
angd =
angd =
angd =
angd =
angd =
angd =
angd =
angd =
angd =
angd =
angd =
angd =
angd =
angd =
angd =
angd =
angd =
angd =
angd =
angd =
angd =
angd =
angd =
angd =
angd =
angd =
angd =
angd =
angd =
angd =
angd =
angd =
angd =
angd =
angd =
angd =
angd =
angd =
angd =
angd =
angd =
angd =
angd =
angd =
angd =
angd =
angd =
angd =
angd =
angd =
angd =
angd =
angd =
angd =
angd =
angd =
angd =
plot(T,Y(:,1),'*',T,Y(:,2),'.',T,Y(:,3),'+')
legend('a','b','c')
function dy = rigid(t,y)
% a column vector
dy = zeros(3,1);
w1=sin(t)*exp(-t);
w2=sin(2*t)*exp(-2*t);
w3=sin(3*t)*exp(-3*t);
dy(1) =(1/cos(y(2)))* (cos(y(3))*w1-sin(y(3))*w2);
dy(2) =(sin(y(3))*w1+cos(y(3))*w2);
dy(3) =w3+(1/cos(y(2)))*(-w1*cos(y(3))*sin(y(2))+w2*sin(y(2))*sin(y(3)));
syms a b c ad bd cdd t real % om1 om2 om3 should be skipped
%rotations for 1(a)-2(b)-3(c)
X=[1 0 0;0 cos(a) sin(a);0 -sin(a) cos(a)];
Y=[cos(b) 0 -sin(b);0 1 0;sin(b) 0 cos(b)];
Z=[cos(c) sin(c) 0;-sin(c) cos(c) 0;0 0 1];
%Last rotation 3(c)
L=Z;
%Last 3(c) and middle rotation 2(b)
LM=Z*Y;
%Last rotation angle derivative
dv1=[0; 0; cdd];
%Middle rotation angle derivative
dv2=[0; bd; 0];
%First rotation angle derivative
w11=dv1;
w22=L*dv2;
w33=LM*dv3;
%summation
w=w11+w22+w33;
[S] = equationsToMatrix([w], [ad, bd, cdd]);
%Convert to differential equation form
%kinematic equation
IS=simplify (inv(S));
%Calculate differential equations
%angular velocity vectors
om1=sin(t)*exp(-t);
om2=sin(2*t)*exp(-2*t);
om3=sin(3*t)*exp(-3*t);
%Differential equation
angd=IS*[om1 om2 om3]'
end

サインインしてコメントする。

カテゴリ

Help Center および File ExchangeOrdinary Differential Equations についてさらに検索

R2022b

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by