ode23 , matrix equation , 2nd order DE

Question

Gloria 2019-6-14

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/467010-ode23-matrix-equation-2nd-order-de

评论： darova 2019-6-19

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;

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

darova 2019-6-14

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/467010-ode23-matrix-equation-2nd-order-de#answer_379194

在 MATLAB Online 中打开

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

7 个评论
显示 5更早的评论隐藏 5更早的评论

Gloria 2019-6-16

编辑：Jan 2019-6-16

在 MATLAB Online 中打开

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

Gloria 2019-6-19

在 MATLAB Online 中打开

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

Gloria 2019-6-19

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 2019-6-19

Can you accept the answer please?

请先登录，再进行评论。

Answer 2

Torsten 2019-6-14

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/467010-ode23-matrix-equation-2nd-order-de#answer_379144

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.

6 个评论
显示 4更早的评论隐藏 4更早的评论

Gloria 2019-6-16

在 MATLAB Online 中打开

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