pim4.m

这是一种求解微分方程的高精度数值方法

clear all;
clc;
close all;
n=7;
G=8.08e10;
mid=7820;
l=[0.05 0.19 0.20 0.085 0.23 0.092 0.113 0.0];
D1=[0.023 0.023 0.076 0.076 0.076 0.076 0.046];
D=0.0095;
mm=[0.3 0.45 0.8 0.6 0.8 0.6 2.0];
m=zeros(1,n);
for i=1:n
    if i==2
     m(i)=6.2532e-6*(l(i)+l(i+1)/2)+mm(i)*D1(i)^2/16;   
    end
m(i)=6.2532e-6*(l(i)+l(i+1)/2)+mm(i)*D1(i)^2/8;
end
c=zeros(1,n+2);%
for i=2:8
    c(i)=0;
end
k=zeros(1,n+1);
%%%%%%%%%%%以上变量不同情况要从新计算
Pk=0.25;%电机额定功率为250，由电压调节控制输入功率
zhuansu=8000;%单位rpmin
Tin=9549.3*Pk/zhuansu;
T1=zeros(n,1);
T1(1,1)=Tin;
for j=2:n
    e(j)=l(j)/(7.9964e-10*G);
end
pp=1;
abe=zeros(n);
if pp==0
    fprintf('there is no crack');
elseif pp==1
    fprintf('there is one crack');
de=0.1;
ci=2;%出现裂纹的轴段号
xe=2;
be=bianrd(D,G,de,pp,xe);
be=l(ci+1)/be;
abe(ci,ci)=be;
else 
    fprintf('there are more crack');
    de=[0.3 0.1];
G=8.08e10;
ci=[1 2];%出现裂纹的轴段号
be=bianrd(D,G,de,pp,xe);
for j=1:pp
be(j)=l(ci(j)+1)/be(j);
end
for j=1:pp
abe(ci(j),ci(j))=be(j);
end
end

[M,C,Ee]=arryf(m,c,e,n);
%dK=[10000 0;0 0];
o=zeros(n);
o1=zeros(n,1);
%k=k*1e8;
%F=[0;10];
pz=1;
if pz==1
f0=[o1;T1];
vk=zeros(2*n,1);
elseif pz==0
    f0=[o1;o1];
    vk=[0.0;0.0;0.1;0.1];
end
f1=zeros(size(f0));
tf=5;
step=0.01; % different step size
N=20;
%%%求步长要首先确定
%figure;
%hold;
dt=step/2^N;%将h进一步细分
for iter=1:tf/step
      iter;
    t(:,iter)=step*(iter-1);
    %f10=f0*sin(zhuansu/60*t(:,iter));
    if pp==1
    BEe=Ee+abe*sin(pi*zhuansu/30*t(iter));
    end
    K=inv(BEe)
    Kx=100*eye(n);
    H=[-inv(M)*C/2 inv(M);0.25*C*inv(M)*C-K-Kx -0.5*C*inv(M)];
   f0(n+1:2*n,1)=f0(n+1:2*n,1)+Kx*vk(n+1:2*n,1);
    iH=inv(H);
    I=eye(size(H));
 Ta=H*dt+(H*dt)^2*(I+(H*dt)/3+(H*dt)^2/12+(H*dt)^3/60)/2;%泰勒展开高阶项
  for jj=1:N%求总Ta
    Ta=2*Ta+Ta*Ta;
  end
  T=I+Ta;
    v(:,iter)=vk(:);%给定初值
    vk=T*(vk+iH*(f0+iH*f1))-iH*(f0+iH*f1+f1*step);%计算精细积分点处的位移量
end
v1=v;
if pp==0
save('dgd','v');
elseif pp==1
 save('bgd','v1');   
end
  % PIM end
ng=length(v(1,:));
[pyy1,pyy2,yy1,yy2,ff]=gonglp(v,ng,step);
ff=2*pi*ff;
 figure(1);
plot(t(1:tf/step),v(1,:),'b-');
figure(2);
plot(t(1:tf/step),v(2,:),'b-');
figure(3);
plot(ff,imag(yy1));
title('虚频图');
xlabel('频率');
ylabel('虚部');
figure(4);
plot(ff,imag(yy2));
title('实频图');
xlabel('频率');
ylabel('虚部');
nfft=2^nextpow2(ng);
figure(5);
plot(ff,pyy1);
title('Frequency content of pyy1');
xlabel('frequency (Hz)');
figure(6);
plot(ff,pyy2);
title('Frequency content of pyy2');
xlabel('frequency (Hz)');
load 'dgd';
figure(5);
plot(t(1:tf/step),v(1,:)-v2(1,:));