Doubt in Coupled Ode

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

[Torsten]

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

Answer 1

[Torsten]

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

[Susmita Panda]

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

Answer 2

[Sam Chak]

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

[Sam Chak]

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]

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.