MATLAB Answers

matrix while using ODE45

1 view (last 30 days)
Adedayo Odeleye
Adedayo Odeleye on 26 Nov 2020
Commented: VBBV on 26 Nov 2020
I am trying to wrtie Shaft_speed as a function of time , so the shaft speed should increase to from 1 to 4000 as time goes on.
I also need the initial condition values X,Y and Z of the last speed to be the new intial condition of the new shaft speed. so far i have:
clc,clear,close all
J=3.37e-5; %inertia of rotor (kg/m2)
K1=0.7243; %Linear stiffness coef (N.m/rad)
K3=1409.7; %Nonlinear stiffness coef (N.m/rad3)
f=0.005;
Zeta=0.05;%damping ratio
Cmech=2*Zeta*(sqrt(K1/J))*J; %mechanical damping coef
Rint=82;
Rload=82;
ECF=sqrt(1.9e-6/Rint); %electromagnetic coupling factor
L=160e-3;%Inductance (mH)
%initial conditions
X = 0; %initial displacement
Y = 0; %initial Velocity
Z = 0; %initial .....
%shaft speed in RPM
Shaft_speed=0:1000;
%shaft speed in rad/s
Rad_speed=Shaft_speed.*((2*pi)/60);
%freq of harmonic fluctuations
w=2.*Rad_speed;
w_Hz=w./2*pi; %convert to Hertz
nt=1./w_Hz;
Time=0:nt/100:100*nt;
Iter=10;
for r=1:Iter
[t,q]=ode45(@duf,Time,[X Y Z]);
X=q(end,1);
Y=q(end,2);
Z=q(end,3);
Max_values(r) = max(q(:,1))
plot(t,q(:,1),'k')
xlabel('TIme,s')
ylabel('Relative Displacement, Rad')
grid
end
max=max(Max_values)
function qdot=duf(t,q) % q_ddot+(Cmech/J)*q_dot+(K1/J)*q+(K3/J)*x^3=f*(w^2)*cos(w*t)
J=3.37e-5; %inertia of rotor (kg/m2)
K1=0.7243; %Linear stiffness coef (N.m/rad)
K3=1409.7; %Nonlinear stiffness coef (N.m/rad3)
f=0.005;
Zeta=0.05;%damping ratio
Cmech=2*Zeta*(sqrt(K1/J))*J; %mechanical damping coef
Rint=82;
Rload=82;
ECF=sqrt(1.9e-6/Rint); %electromagnetic coupling factor
L=160e-3;%Inductance (mH)
%shaft speed in RPM
Shaft_speed=0:1000;
%shaft speed in rad/s
Rad_speed=Shaft_speed.*((2*pi)/60);
%freq of harmonic fluctuations
w=2.*Rad_speed;
w_Hz=w./2*pi; %convert to Hertz
qdot(1)=q(2);
qdot(2)=f.*w.^2.*cos(w.*t)-(Cmech/J).*q(2)-(K1/J).*q(1)-(K3/J).*q(1).^3- ECF.*q(3);
qdot(3)=(-ECF.*q(2)-(Rint+Rload).*q(3))./L;
qdot = [qdot(1);qdot(2);qdot(3)];
end
but i geet keep getting '-Error using odearguments (line 21)
When the first argument to ode45 is a function handle, the tspan argument must have at least two elements.
Error in ode45 (line 115)
odearguments(FcnHandlesUsed, solver_name, ode, tspan, y0, options, varargin);
Error in mainduffing (line 33)
[t,q]=ode45(@duf,Time,[X Y Z]);'
Thanks in advance

  0 Comments

Sign in to comment.

Accepted Answer

VBBV
VBBV on 26 Nov 2020
Shaft_speed=1; % specify a value instead of vector
%shaft speed in rad/s
Rad_speed=Shaft_speed*((2*pi)/60);
%freq of harmonic fluctuations
w=2*Rad_speed;
w_Hz=w/2*pi; %convert to Hertz
nt=1/w_Hz;
Time=0:nt:100*nt; % use coarse step length to reduce time of run
inside the for loop specify time as initial and final states
[t,q]=ode45(@duf,[0 max(Time)],[X Y Z]);
omit the element wise operation inside the function duf
qdot(2)=f*w^2*cos(w*t)-(Cmech/J)*q(2)-(K1/J).*q(1)-(K3/J)*q(1)^3- ECF*q(3);

  6 Comments

Show 3 older comments
Adedayo Odeleye
Adedayo Odeleye on 26 Nov 2020
Thank you very much! Really apericiate your help. The script is running fine.
Just one last question, what is the implication of having a different value for Shaft_speed in the function because am not able to mirror 'Shaft_speed=0.2:1000' into the function part, I can only write a single value unless I get error:
'Error using ^
One argument must be a square matrix and the other must be a scalar. Use POWER (.^) for elementwise power.
Error in Model>duf (line 65)
qdot(2)=f*w^2*cos(w*t)-(Cmech/J)*q(2)-(K1/J)*q(1)-(K3/J)*q(1)^3- ECF/J*q(3);
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 Model (line 32)
[t,q]=ode45(@duf,[0 max(Time)],[X Y Z]);'
VBBV
VBBV on 26 Nov 2020
Since you are using Shaft_speed as vector inside duf function comment that line. And use S_Speed in place of Shaft_speed.
If you use it as vector inside the function the matrix dimensions do not match with output vector q size I.e.3
VBBV
VBBV on 26 Nov 2020
Instead pass shaft spped as function argument I.e. S_Speed
Also define a coarse step size for Shaft_speed in the beginning to reduce computation time e.g.
%if true
Shaft_speed = 0.2:200:1000;

Sign in to comment.

More Answers (1)

Walter Roberson
Walter Roberson on 26 Nov 2020
Shaft_speed=0:1000;
That is a vector, with first element 0.
%shaft speed in rad/s
Rad_speed=Shaft_speed.*((2*pi)/60);
So that is a vector, with first element 0.
%freq of harmonic fluctuations
w=2.*Rad_speed;
So w is a vector, with first element 0.
w_Hz=w./2*pi; %convert to Hertz
so w_Hz is a vector, with first element 0.
nt=1./w_Hz;
so nt is a vector, with first element 1/0 -> inf
Time=0:nt/100:100*nt;
and the right hand side is then 0:nt(1):100*nt(1) -> 0:inf/100:100*inf -> 0:inf:inf which is just the element 0 .
When you use a vector in a colon operation, the value used is the first one, which happens to be inf in this case.

  1 Comment

Adedayo Odeleye
Adedayo Odeleye on 26 Nov 2020
Thank you for the help, i have replaced Shaft_speed with 1:1000 now i am getting :
In an assignment A(:) = B, the number of elements in A and B must be the same.
Error in mainduffing>duf (line 72)
qdot(2)=f.*w.^2.*cos(w.*t)-(Cmech/J).*q(2)-(K1/J).*q(1)-(K3/J).*q(1).^3- (ECF/J).*q(3);
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 mainduffing (line 32)
[t,q]=ode45(@duf,[Time],[X Y Z]);

Sign in to comment.

Products


Release

R2017a

Community Treasure Hunt

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

Start Hunting!