ode23 , matrix equation , 2nd order DE

How I will solve this matrix in Matlab ?

w=10.02;
p= 7850;     %%%Kg/m^3;
g=9.81;
G=77*10^9;        %%% N/m^2
E=206*10^9;       %%% N/m^2
L=9.6;   %%%m
D=0.15 ; %%m
m=pi*(D/2)^2*L*p;
A=pi*D^2/4 ; %%m2
J=m*D^2/8;  %%Kg*m^2
Ip=pi*D^4/32 ;   %% m^4
I=pi*D^4/64 ; 
cx=0.05*m*w;
cy=0.05*m*w;
ct=0.08*J*w;
kx=3*E*I/L;
ky=3*E*I/L;
kt=G*Ip/L;

Answer 1

[darova]

You can find roots for x'', y'' and dtheta'' every iteration solving matrix equation

function du = func(t,u)
    u1 = u(1:3)';                        % [x y theta]'
    du1 = u(4:6)';                       % [dx dy dtheta]'
    dth = u(6);                          % dtheta
    C = [cx 0 0; 0 cy 0; 0 0 cz];        % C matrix
    K = ...                              % K matrix
    F = [me(w+dth)^2+Fx; ...]            % force vector
    A = [m 0 -mesin(wt)...]              % mass matrix
    B = -C*du1 -K*u1 + F;                
    du = zeros(6,1);
    du(1:3) = u(4:6);
    du(4:6) = A\B;
end

[Gloria]

Undefined function or variable 'u'.

%function du = func(t,u)
clc;
clear all;

%%%%%%   Shaft Speed %%%%%%%%%%
rpm=100;   %rpm
w=rpm*2*pi/60;  % rad/s
%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%   coupled Vibration %%%%%%%%%%
e=0.01;
%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%  Time Span (Seconds)  %%%%%%%%
t0=0;   tf=10;  art=0.001;
t=[t0:art:tf];
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%%%%%% Forces %%%%%%%%%
Fx=@(t)  0*sin(w*t);
Fy=@(t)  0*sin(w*t); 
Mt=@(t)  184.8*sin(w*t); 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
p= 7850;     %%%Kg/m^3;
g=9.81;
G=77*10^9;        %%% N/m^2
E=206*10^9;       %%% N/m^2
L=9.6;   %%%m
D=0.15 ; %%m
m=pi*(D/2)^2*L*p;
A=pi*D^2/4 ; %%m2
J=m*D^2/8;  %%Kg*m^2
Ip=pi*D^4/32 ;   %% m^4
I=pi*D^4/64 ; 
cx=0.05*m*w;
cy=0.05*m*w;
ct=0.08*J*w;
kx=3*E*I/L;
ky=3*E*I/L;
kt=G*Ip/L;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   
    u1 = u(1:3)';                        % [x y theta]'
    du1 = u(4:6)';                       % [dx dy dtheta]'
    dth = u(6);                          % dtheta
    C = [cx ,0 ,0; 0, cy, 0; 0, 0, ct];        % C matrix
    K = [kx, 0, 0;0 ,ky, 0;0, 0 ,kt];                             % K matrix
    F = [m*e*(w+dth)^2*cos(w*t)+Fx(t); -m*g+m*e*(w+dth)^2*sin(w*t)+Fy(t);-m*e*g*cos(w*t)+Mt(t)];            % force vector
    A = [m ,0, -m*e*sin(w*t);0, m, m*e*cos(w*t);-m*e*sin(w*t), m*e*cos(w*t), J+m*e^2];              % mass matrix
    B = -C*du1 -K*u1 + F;                
    du = zeros(6,1);
    du(1:3) = u(4:6);
    du(4:6) = A\B;
    
    
figure(2), clf
subplot 311
plot(t,u(1),'b')
hold on
subplot 312
plot(t,u(2),'b')
hold on
subplot 313
plot(t,u(3),'b')
hold off
%end

[darova]

The function must be in a separate file .m (file name must be the same as function)

In main code (read about how ode45 works)

[t,u] = ode45(@func,[t0 tf],u0);
x = u(:,1);
% ...
dy = u(:,5);
dth = u(:,6);
plot(t,x)

Please use buttons for code inserting

[Jan]

@Gloria: I've edited the codes in your messages and applied a proper code formatting to improve the readability. As darova has mentioned already, you can do this by your own.

Please post the complete error message, which contains also the location of the problem.

[Gloria]

function [du] = dokuz( t,u )

clc;
clear all;

%%%%%%   Shaft Speed %%%%%%%%%%
rpm=100;   %rpm
w=rpm*2*pi/60;  % rad/s
%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%   coupled Vibration %%%%%%%%%%
e=0.01;
%%%%%%%%%%%%%%%%%%%%%%%%%%%

 %%%%%% Forces %%%%%%%%%
Fx=@(t)  0*sin(w*t);
Fy=@(t)  0*sin(w*t); 
Mt=@(t)  184.8*sin(w*t); 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
p= 7850;     %%%Kg/m^3;
g=9.81;
G=77*10^9;        %%% N/m^2
E=206*10^9;       %%% N/m^2
L=9.6;   %%%m
D=0.15 ; %%m
m=pi*(D/2)^2*L*p;
A=pi*D^2/4 ; %%m2
J=m*D^2/8;  %%Kg*m^2
Ip=pi*D^4/32 ;   %% m^4
I=pi*D^4/64 ; 
cx=0.05*m*w;
cy=0.05*m*w;
ct=0.08*J*w;
kx=3*E*I/L;
ky=3*E*I/L;
kt=G*Ip/L;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%   



    u1 = u(1:3)';                        % [x y theta]'
    du1 = u(4:6)';                       % [dx dy dtheta]'
    dth = u(6);                          % dtheta
    C = [cx 0 0; 0 cy 0; 0 0 cz];        % C matrix
    K = [kx 0 0; 0 ky 0; 0 0 kz];        % K matrix
    F = [m*e*(w+dth)^2*cos(w*t)+Fx; -m*g+m*e*(w+dth)^2*sin(w*t)+Fy;-m*e*g*cos(w*t)+Mt];            % force vector
    A = [m ,0, -m*e*sin(w*t);0, m, m*e*cos(w*t);-m*e*sin(w*t), m*e*cos(w*t), J+m*e^2];                 % mass matrix
    B = -C*du1 -K*u1 + F;                
    du = zeros(6,1);
    du(1:3) = u(4:6);
    du(4:6) = A\B;


end



%%%for m file

%%%%%%%%  Initial Conditions  %%%%%%%%%%%%
IC=[0,0,0,0,0,0];
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%  Time Span (Seconds)  %%%%%%%%
t0=0;   tf=10;  art=0.001;
tspan=[t0:art:tf];
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

[t,u] = ode45(@dokuz,tspan,IC);
x = u(:,1);
y=u(:,2);
th=u(:,3);
dx = u(:,4);
dy = u(:,5);
dth = u(:,6);
plot(t,x)


%%%%%%%%%%% I need to learn how I will use ODE45  to solve matris form equations

[darova]

Main mistake is using clear all in function. You clearing all variables including t and u (function arguments)

[Gloria]

Error using ode45

Too many input arguments.

Error in Untitled2 (line 10)

[t,u] = ode45(@dokuz,tspan,IC);

BUT WITH ODE23 IT GIVES THE ANSWERS, THANK YOU SO MUCH

[darova]

Can you accept the answer please?

Answer 2

[Torsten]

Convert the system to a first order system and write it as

M*y' = f(t,y)

Then use the mass-matrix option of the ODE solvers to supply M, define f in a function file and use ODE45, ODE15S, ... to solve.

[Gloria]

But the part I didn't understand is all the samples with ode45 ;

x''=( - c*x' -k * x) /m and for example x''= ( - c*s(2) -k * s(1)) /m

for here x'' , y'' and z'' connect eachother.

x'' =( - c*x' -k * x -m* y'')/m ( - c*s(2) -k * s(1) -m* y'')/m

how I will solve this?

[Torsten]

s1' = s4
s2' = s5
s3' = s6
m*s4' - m*e*sin(w*t)*s6' = -c_x*s4 - k_x*s1 + m*e*(w+s6)^2*cos(w*t) + F_x
m*s5' + m*e*cos(w*t)*s6' = -c_y*s5 -k_y*s2 - m*g + m*e*(w+s6)^2*sin(w*t) + F_y
-m*e*sin(w*t)*s4' + m*e*cos(w*t)*s5' + (J+m*e^2)*s6' = -c_theta*s6 -k_theta*s3 -m*e*g*cos(w*t) + M_theta

where

s1 = x, s4 = x', s2 = y, s5 = y', s3 = theta, s6 = theta'

Now write the system as

M*s' = f(s,t)

and you are nearly done.

[Gloria]

Thank you so much for your help ;

I didn't know that I can write s6',s4'.. in ODE code.

[Jan]

@Gloria: You can't. s6' is not Matlab code, but a mathematical expression.

[Gloria]

But I got the figures;

 sdot=@(t,s)...
      [s(2);
      (m*e*sin(w*t)*s(6)'-cx*s(2)-kx*s(1)+m*e*(w+s(6))^2*cos(w*t)+Fx(t))/m;
      s(4);
      (-m*e*cos(w*t)*s(6)'-cy*s(4)-ky*s(3)-m*g+m*e*(w+s(6))^2*sin(w*t)+Fy(t))/m;
      s(6);
      (m*e*sin(w*t)*s(2)'-m*e*cos(w*t)*s(4)'-ct*s(6)-kt*s(5)-m*e*g*cos(w*t)+Mt(t))/(J+m*e^2)];

So Can I write s(6)',s(2)'s(4)' or not?

[Jan]

You can write s(6)', because ' is the Matlab operator for the complex conjugate transposition:

a = [1, 1i; ...
     2, 2i]
a'

You provide real scalars. Then the ' operaotr does not change the value. So it is valid, but simply a confusing waste of time.