Not enough input arguments ode45
6 Ansichten (letzte 30 Tage)
Ältere Kommentare anzeigen
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')
Antworten (2)
Dyuman Joshi
am 22 Nov. 2023
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
0 Kommentare
Sulaymon Eshkabilov
am 22 Nov. 2023
Bearbeitet: Sulaymon Eshkabilov
am 22 Nov. 2023
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);
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
0 Kommentare
Siehe auch
Kategorien
Mehr zu Ordinary Differential Equations finden Sie in Help Center und File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!