matrix while using ODE45

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

Answer 1

[VBBV]

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

[VBBV]

Inside the function duf you are trying to assign, different size dimensions for vectors, Since size of q is only 3 and shaft speed vector of different size,

[Adedayo Odeleye]

Thank you for the advice,

is there a way i can increase the value for Shaft_speed within the loop. So when the loop is finished calculating q at a certain speed, the final values of q are then the new initial condition for the next value of Shaft_speed

[VBBV]

Shaft_speed=0.2:1000; % specify a arbitrary small value for shaft speed, bcoz a 0 speed will produce infinite time value
for i = 1:length(Shaft_speed)
	%shaft speed in rad/s
	S_Speed = Shaft_speed(i);
	Rad_speed=Shaft_speed(i)*((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;
	Iter=10;
	for   r=1:Iter
		[t,q]=ode45(@duf,[0 max(Time)],[X Y Z]);
		X=q(1); % Final values are assigned to initial values for next iteration 
		Y=q(2);
		Z=q(3);
		Max_values(r) = max(q(:,1))
		plot(t,q(:,1),'k')
		xlabel('TIme,s')
		ylabel('Relative Displacement, Rad')
		grid
	end
end
function qdot=duf(t,q,S_Speed) % q_ddot+(Cmech/J)*q_dot+(K1/J)*q+(K3/J)*x^3=f*(w^2)*cos(w*t)
...
...
end

Pass the shaft speed as function argument

[Adedayo Odeleye]

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]

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]

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;

Answer 2

[Walter Roberson]

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.

[Adedayo Odeleye]

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