Not enough input arguments ode45

조회 수: 2 (최근 30일)
Aleyna
Aleyna 2023년 11월 22일
편집: Sulaymon Eshkabilov 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
dv3=[ad; 0; 0]
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
Aleyna
Aleyna 2023년 11월 22일
Question: Calculate Euler angles differential equations for 1(a)-2(b)-3(c) rotation sequence

댓글을 달려면 로그인하십시오.

답변 (2개)

Dyuman Joshi
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
dv3=[ad; 0; 0];
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

Sulaymon Eshkabilov
Sulaymon Eshkabilov 2023년 11월 22일
편집: 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
dv3=[ad; 0; 0];
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 CenterFile Exchange에서 Ordinary Differential Equations에 대해 자세히 알아보기

Community Treasure Hunt

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

Start Hunting!

Translated by