# matrix while using ODE45

1 view (last 30 days)
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)
f=0.005;
Zeta=0.05;%damping ratio
Cmech=2*Zeta*(sqrt(K1/J))*J; %mechanical damping coef
Rint=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;
%freq of harmonic fluctuations
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')
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)
f=0.005;
Zeta=0.05;%damping ratio
Cmech=2*Zeta*(sqrt(K1/J))*J; %mechanical damping coef
Rint=82;
ECF=sqrt(1.9e-6/Rint); %electromagnetic coupling factor
L=160e-3;%Inductance (mH)
%shaft speed in RPM
Shaft_speed=0:1000;
%freq of harmonic fluctuations
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 = [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]);'

VBBV on 26 Nov 2020
Shaft_speed=1; % specify a value instead of vector
%freq of harmonic fluctuations
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);

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 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;

Walter Roberson on 26 Nov 2020
Shaft_speed=0:1000;
That is a vector, with first element 0.
So that is a vector, with first element 0.
%freq of harmonic fluctuations
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 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]);

R2017a

### Community Treasure Hunt

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

Start Hunting!