matrix while using ODE45
4 views (last 30 days)
Show older comments
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
Accepted Answer
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
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
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;
More Answers (1)
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.
See Also
Categories
Find more on Ordinary Differential Equations in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!