Doubt in Coupled Ode

Question

Susmita Panda 2023-11-15

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2047717-doubt-in-coupled-ode

评论： Torsten 2023-11-16

I have written a code for coupled ode using ode45; however my results seems errornous when compared with analytical results:

%% Analytical

L=24;

A=9;

Iz=2.43;

J=21.18;

E=33.2e9;

EI=E*Iz;

rho=2.4e4;

volume=A*L;

m=(rho*A);

Fv=29.9e4;

theta=30*(pi/180);

R=L/theta;

v=40;

mu=0.2;

G=E/(2*(1+mu));

GJ=G*J;

g=10;

a1=(1/m)*(pi/L)^2*(EI*((pi/L)^2)+(GJ/(R^2)));

a2=(1/(m*R))*((pi/L)^2)*(EI+GJ);

b1=-(1/(rho*J))*((EI/R^2)+(GJ*((pi/L)^2)));

b2=-(1/(rho*J))*(1/R)*((pi/L)^2)*(EI+GJ);

wv=32.10; %% book formula provided in snip shot with constants

% factor_1=sqrt(((a1-b1)^2)+(4*a2*b2));

% wv=sqrt((a1+b1+factor_1)/2) ;% however I am getting different with my constants

g=10;

t=0.01:0.01:(L/v);

Sv=(pi*v)/(L*wv);

beta=(b1-(pi*v/L)^2)/(b1+a1-wv^2-(pi*v/L)^2);

xi=sin(pi*v*t/L)-(Sv*sin(wv*t));

u_deflect=-((2*Fv*g)/(m*L))*(1/wv^2)*(1/(1-Sv^2))*beta*xi;

figure();plot(t,u_deflect,'o-');title('Validation');

My code is:

%% Analysis using ode45

tspan_1=[0:0.001:0.6];%time range

y0_1=[0;0;0;0];%initial conditions f

[t1,y1]=ode45(@diffeqn11,tspan_1,y0_1);

%plot displacement

figure(1)

plot(t1, y1(:,1), 'r', 'LineWidth',2);title('Displacement of the beam');

hold on

plot(t1, y1(:,3), 'b', 'LineWidth',2);

%plot velocity_forced+free

figure(2)

plot(t1, y1(:,2), 'g', 'LineWidth',2);title('Velocity of the beam');

figure(3)

plot(t1, y1(:,4), 'r', 'LineWidth',2);

function f= diffeqn11(t,y)

L=24;

A=9;

Iz=2.43;

J=21.18;

E=33.2e9;

EI=E*Iz;

rho=2.4e4;

volume=A*L;

m=(rho*A);

Fv=29.9e4;

theta=30*(pi/180);

R=L/theta;

v=40;

mu=0.2;

G=E/(2*(1+mu));

GJ=G*J;

g=10;

a1=(1/m)*(pi/L)^2*(EI*((pi/L)^2)+(GJ/(R^2)));

a2=(1/(m*R))*((pi/L)^2)*(EI+GJ);

b1=-(1/(rho*J))*((EI/R^2)+(GJ*((pi/L)^2)));

b2=-(1/(rho*J))*(1/R)*((pi/L)^2)*(EI+GJ);

f=zeros(4,1);

f(1)=y(2);

f(2)=((2*Fv)/(m*L))*sin((pi*v*t)/L)-(a1*y(1))-(a2*y(3));

f(3)=y(4);

f(4)=-(b1*y(3))-(b2*y(1));

end

1 个评论
显示 -1更早的评论隐藏 -1更早的评论

Torsten 2023-11-16

Either the solution you took from the book is wrong or your differential equations are wrong.

请先登录，再进行评论。

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

Answer 1

Torsten 2023-11-15

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2047717-doubt-in-coupled-ode#answer_1353227

编辑：Torsten 2023-11-15

在 MATLAB Online 中打开

Your code is correct, but it seems there is a singularity near 0.32....

syms t x(t) y(t) a1 b1 a2 b2 Constant

eqn1 = diff(x,t,2)+a1*x+a2*y==Constant*sin(0.5*t);

eqn2 = diff(y,t,2)+b1*y+b2*x==0;

Dx = diff(x,t);

Dy = diff(y,t);

sol = dsolve([eqn1,eqn2],[x(0)==0,y(0)==0,Dx(0)==0,Dy(0)==0]);

sol.x

ans =

sol.y

ans =

figure(1)

fplot(subs(sol.x,[a1 a2 b1 b2 Constant],[120.7217,1.3587e+06,-4.4643e+06,-1.2327e+05,0.1154]),[0 0.6])

figure(2)

fplot(subs(sol.y,[a1 a2 b1 b2 Constant],[120.7217,1.3587e+06,-4.4643e+06,-1.2327e+05,0.1154]),[0 0.6])

Note that the eigenvalues of the characteristic polynomial for the system are tremendous:

a1 = 120.7217;
a2 = 1.3587e+06;
b1 = -4.4643e+06;
b2 = -1.2327e+05;
M = [0 0 1 0;0 0 0 1;-a1 -a2 0 0;-b2 -b1 0 0];
eig(M)
ans = 4×1
1.0e+03 *

   -2.1039
    2.1039
    0.1942
   -0.1942

1 个评论
显示 -1更早的评论隐藏 -1更早的评论

Susmita Panda 2023-11-16

You are correct sir. I must check the values of constants.

请先登录，再进行评论。

Answer 2

Sam Chak 2023-11-16

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2047717-doubt-in-coupled-ode#answer_1353392

在 MATLAB Online 中打开

Hi @Susmita Panda

Your 4th-order coupled system has two eigenvalues with positive real parts. Therefore, in theory, the state responses will diverge when the system is forced.

% parameters
a1 =  120.7217;
a2 =  1.3587e+06;
b1 = -4.4643e+06;
b2 = -1.2327e+05;
% state matrix
A  = [ 0 1   0 0;
     -a1 0 -a2 0;
       0 0   0 1;
     -b2 0 -b1 0];
% check eigenvalues
eig(A)
ans = 4×1
1.0e+03 *

    2.1039
   -2.1039
   -0.1942
    0.1942

2 个评论
显示无隐藏无

Sam Chak 2023-11-16

在 MATLAB Online 中打开

Hi @Susmita Panda

Please note that in your beam-coupled system, there is typically no damping considered. However, in the real physical world, damping components exist. As a structural engineer, you can either incorporate other beams with stabilizing properties, or measure the deflections of existing beams and then counteract destabilizing vibration modes.

%% Analysis using ode45

tspan_1 = [0:0.001:4*pi]; % time range

y0_1 = [0; 0; 0; 0]; % initial conditions f

options = odeset('RelTol', 1e-8, 'AbsTol', 1e-8);

[t1, y1] = ode45(@diffeqn11, tspan_1, y0_1, options);

% plot displacement

figure(1)

plot(t1/pi, y1(:,1:2:3)); grid on,

title('Displacements of the beams'), xlabel('t/\pi')

legend('y_{1}', 'y_{3}', 'location', 'SE')

% A pair of Beam couples

function f = diffeqn11(t, y)

% parameters

a1 = 120.7217;

a2 = 1.3587e+06;

b1 = -4.4643e+06;

b2 = -1.2327e+05;

Constant= 0.1154;

% measure the deflections and velocities, and feed them back into the

% system in the same channel, where the sinusoidal force is injected.

k1 = 10562734.7949266;

k2 = 4596.24524333348;

k3 = 382536034.238245;

k4 = 181690.721406603;

K = [k1 k2 k3 k4];

f = zeros(4, 1);

u = Constant*sin(0.5*t) - K*y;

f(1) = 0*y(1) + 1*y(2) + 0*y(3) + 0*y(4);

f(2) = - a1*y(1) + 0*y(2) - a2*y(3) + 0*y(4) + u;

f(3) = 0*y(1) + 0*y(2) + 0*y(3) + 1*y(4);

f(4) = - b2*y(1) + 0*y(2) - b1*y(3) + 0*y(4);

end

Susmita Panda 2023-11-16

编辑：Susmita Panda 2023-11-16

Thank you for the idea sir; however the problem in my case I am getting solution using analytical solution and not getting using ode45. I am also supposing that my constants are errorneous. I should provide the actual code and actual values of constant instead of demo code. Sir do check my revised question.

请先登录，再进行评论。