%微分求积法求解直管的流致振动
clc;
clear;
format long;
a=0.002;
tao=0;
c=100;
kt=10;
u=4;
N=10;%截点数目
beta=0.2;
zx=ones(1,N);
% 定义delat，截点坐标
delta = 0.0001;
coord_x(1,1) = 0;
coord_x(1,2) = delta;
coord_x(1,N-1) = 1-delta;
coord_x(1,N) = 1;
%中间截点
for i=3:(N-2)
coord_x(1,i)=(1/2)*(1-cos((i-2)/(N-3)*pi));
end
%计算一阶权系数
for i=1:N
for k=1:N
if i~=k
zx(1,i)=zx(1,i)*(coord_x(1,i)-coord_x(1,k));
end
end
end
for i=1:N
for j=1:N
if i~=j
wei_x_1(i,j)=zx(1,i)/zx(1,j)/(coord_x(1,i)-coord_x(1,j));
end
end
end
for i=1:N
for k=1:N
if i~=k
wei_x_1(i,i)=wei_x_1(i,i)-wei_x_1(i,k);
end
end
end
%计算二阶权系数
for i=1:N
for j=1:N
if i~=j
wei_x_2(i,j)=2*(wei_x_1(i,i)*wei_x_1(i,j)-wei_x_1(i,j)/(coord_x(1,i)-coord_x(1,j)));
end
end
end
for i=1:N
for k=1:N
if i~=k
wei_x_2(i,i)=wei_x_2(i,i)-wei_x_2(i,k);
end
end
end
%计算三阶权系数
for i=1:N
for j=1:N
if i~=j
wei_x_3(i,j)=3*(wei_x_2(i,i)*wei_x_1(i,j)-wei_x_2(i,j)/(coord_x(1,i)-coord_x(1,j)));
end
end
end
for i=1:N
for k=1:N
if i~=k
wei_x_3(i,i)=wei_x_3(i,i)-wei_x_3(i,k);
end
end
end
%计算四阶权系数
for i=1:N
for j=1:N
if i~=j
wei_x_4(i,j)=4*(wei_x_3(i,i)*wei_x_1(i,j)-wei_x_3(i,j)/(coord_x(1,i)-coord_x(1,j)));
end
end
end
for i=1:N
for k=1:N
if i~=k
wei_x_4(i,i)=wei_x_4(i,i)-wei_x_4(i,k);
end
end
end
%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for ii=1:100
c=ii
rm=-1
im=0
cc=0
for cc=1:100
u=0.1*cc
K=zeros(N-4,N-4);
G=zeros(N-4,N-4);
M=zeros(N-4,N-4);
%kbb
kbb=[1,0,0,0;wei_x_1(2,1),wei_x_1(2,2),wei_x_1(2,N-1),wei_x_1(2,N);wei_x_2(N-1,1)+kt*wei_x_1(N-1,1),wei_x_2(N-1,2)+kt*wei_x_1(N-1,2),wei_x_2(N-1,N-1)+kt*wei_x_1(N-1,N-1),wei_x_2(N-1,N)+kt*wei_x_1(N-1,N);wei_x_3(N,1),wei_x_3(N,2),wei_x_3(N,N-1),(wei_x_3(N,N)-c)];
%kbb=[wei_x_3(1,1),wei_x_3(1,2),wei_x_3(1,N-1),(wei_x_3(1,N)-r);wei_x_2(2,1),wei_x_2(2,2),wei_x_2(2,N-1),wei_x_2(2,N);wei_x_2(N-1,1),wei_x_2(N-1,2),wei_x_2(N-1,N-1),wei_x_2(N-1,N);wei_x_3(N,1),wei_x_3(N,2),wei_x_3(N,N-1),(wei_x_3(N,N)-r)];
%kbd
for j=1:N-4
kbd(1,j)=0;
end
for j=1:N-4
kbd(2,j)=wei_x_1(2,j+2);
end
for j=1:N-4
kbd(3,j)=wei_x_2(N-1,j+2)+kt*wei_x_1(N-1,j+2);
end
for j=1:N-4
kbd(4,j)=wei_x_3(N,j+2);
end
%kdb
for i=1:N-4
kdb(i,1)=wei_x_4(i+2,1)+u^2*wei_x_2(i+2,1);
end
for i=1:N-4
kdb(i,2)=wei_x_4(i+2,2)+u^2*wei_x_2(i+2,2);
end
for i=1:N-4
kdb(i,3)=wei_x_4(i+2,N-1)+u^2*wei_x_2(i+2,N-1);
end
for i=1:N-4
kdb(i,4)=wei_x_4(i+2,N)+u^2*wei_x_2(i+2,N);
end
%kdd
for i=1:N-4
for j=1:N-4
kdd(i,j)=wei_x_4(i+2,j+2)+u^2*wei_x_2(i+2,j+2);
end
end
%Gdb
for i=1:N-4
Gdb(i,1)=2*u*sqrt(beta)*wei_x_1(i+2,1)+a*wei_x_4(i+2,1);
end
for i=1:N-4
Gdb(i,2)=2*u*sqrt(beta)*wei_x_1(i+2,2)+a*wei_x_4(i+2,2);
end
for i=1:N-4
Gdb(i,3)=2*u*sqrt(beta)*wei_x_1(i+2,N-1)+a*wei_x_4(i+2,N-1);
end
for i=1:N-4
Gdb(i,4)=2*u*sqrt(beta)*wei_x_1(i+2,N)+a*wei_x_4(i+2,N);
end
%组装GDD矩阵
for i=1:N-4
for j=1:N-4
GDD(i,j)=2*u*sqrt(beta)*wei_x_1(i+2,j+2)+a*wei_x_4(i+2,j+2);
end
end
%G
G=-Gdb*inv(kbb)*kbd+GDD;
%K
K=-kdb*(inv(kbb))*kbd+kdd;
%组装M矩阵
for i=1:N-4
M(i,i)=1;
end
MI=inv(M);
MIK=MI*K;
MIG=MI*G;
R=zeros(2*(N-4),2*(N-4));
%组 装R
for i=1:N-4
R(i+N-4,i)=1;
end
for i=1:N-4
for j=1:N-4
R(i,j)=(-1)*MIG(i,j);
R(i,N-4+j)=(-1)*MIK(i,j);
end
end
eigr=eig(R);
[PIM,sor]=sort(real(eigr));
for j=1:2*(N-4)
numb=sor(j);
mmr(j,1)=eigr(numb,1);
end
mmr
rm=real(mmr(2*(N-4),1))
im=imag(mmr(2*(N-4),1))
cc=cc+1
if rm>0&&im~=0
plot(c,u,'or')
hold on
elseif rm>0&&im==0
plot(c,u,'*b')
hold on
else
plot(c,u,'+k')
end
end
end
