杜哈美积分matlab编程,Matlab讨论区 - 声振论坛 - 振动,动力学,声学,信号处理,故障诊断 - Powered by Discuz!...

%单自由度系统适用，初始条件任意，比happy提供的适应性强。对于多自由度系统，正

%则化之后适用，加循环即可(有点修改，for t=0:dTaT end之间的不用改)。

clc; clear

n=1000; m=3; T=2; si=0.02; k=2700; w=30;

dTao=T/n; wD=w*(1-si^2)^0.5;

p2=[ones(1,333),zeros(1,668)]; % p2=[1000,zeros(1,1000)];

%%%%%%%%%%%%%%%杜哈美积分开始%%%%%%%%%%%%%%%%

x0=0; v0=0; j=1;

for t=0:dTaT

sysN=0; sysP=0;

%%%%%%%%%%%%位移\速度\加速度公用部分%%%%%%%%%%%%%%

for i=0:t/dTao

sysN=sysN+p2(i+1)*exp(si*i*dTao*w)*cos(wD*i*dTao)*dTao;

sysP=sysP+p2(i+1)*exp(si*i*dTao*w)*sin(wD*i*dTao)*dTao;

end

%%%%%%%%%%%%位移\速度\加速度公用部分%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%求位移%%%%%%%%%%%%%%%%%%%

sysA=exp(-si*w*t)*(x0*cos(wD*t)+(si*w*x0+v0)/wD*sin(wD*t));

sysB=exp(-si*w*t)/(m*wD)*(sin(wD*t)*sysN-cos(wD*t)*sysP);

x(j)=sysA+sysB;

%%%%%%%%%%%%%%%%%求位移%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%求速度%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%V=C+D+E+F%%%%%%%%%%%%%%%%%%

sysC=exp(-si*w*t)*(-x0*wD*sin(wD*t)-(si*w)^2*x0*sin(wD*t)/wD);

sysD=exp(-si*w*t)*(v0*cos(wD*t)-si*w*v0*sin(wD*t)/wD);

%%%%%%%%%%%%%%%E=G+H%%%%%%%%%%%%%%%%%%%%%%

sysG=sin(wD*t)*sysN;

sysH=-cos(wD*t)*sysP;

sysE=(sysN+sysP)*-si*w*exp(-si*w*t)/(m*wD);

%%%%%%%%%%%%%%%F=I+J+K+L%%%%%%%%%%%%%%%%%%

sysI=wD*cos(wD*t)*sysN;

sysJ=sin(wD*t)*p2(j)*exp(si*w*t)*cos(wD*t);

sysK=wD*sin(wD*t)*sysP;

sysL=-cos(wD*t)*p2(j)*exp(si*w*t)*sin(wD*t);

sysF=(sysI+sysJ+sysK+sysL)*exp(-si*w*t)/(m*wD);

v(j)=sysC+sysD+sysE+sysF;

%%%%%%%%%%%%%%%%%%%%%%%%%求速度%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%求加速度%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%sysdC%%%%%%%%%%%%%%%

sysdC1=-x0*wD^2*cos(wD*t)-(si*w)^2*x0*cos(wD*t);

sysdC=-si*w*sysC+exp(-si*w*t)*sysdC1;

%%%%%%%%%%%%%%%%%%%sysdD%%%%%%%%%%%%%%%

sysdD1=-v0*wD*sin(wD*t)-si*w*v0*cos(wD*t);

sysdD=-si*w*sysD+exp(-si*w*t)*sysdD1;

%%%%%%%%%%%%%%%%%%%sysdE%%%%%%%%%%%%%%%

sysdG=wD*cos(wD*t)*sysN+sin(wD*t)*p2(j)*exp(si*w*t)*cos(wD*t);

sysdH=wD*sin(wD*t)*sysP-cos(wD*t)*p2(j)*exp(si*w*t)*sin(wD*t);

sysdE=(si*w)^2*exp(-si*w*t)*(sysG+sysH)/(m*wD)-si*w*exp(-si*w*t)*(sysdG+sysdH)/(m*wD);

%%%%%%%%%%%%%%%%%%%sysdF%%%%%%%%%%%%%%%

sysdI=-wD^2*sin(wD*t)*sysN+wD*cos(wD*t)*p2(j)*exp(si*w*t)*cos(wD*t);

%%%%%%%%%%%%怎样考虑dP2,暂时取0%%%sysdJ%%%%%%%%%%

sysdJ1=wD*cos(wD*t)^2*p2(j)*exp(si*w*t);

sysdJ2=sin(wD*t)*0*exp(si*w*t)*cos(wD*t);

sysdJ3=sin(wD*t)*p2(j)*si*w*exp(si*w*t)*cos(wD*t);

sysdJ4=-wD*sin(wD*t)^2*p2(j)*exp(si*w*t);

sysdJ=sysdJ1+sysdJ2+sysdJ3+sysdJ4;

%%%%%%%%%%%%%%%%%%%sysdK%%%%%%%%%%%%%%%

sysdK=wD^2*cos(wD*t)*sysP+wD*sin(wD*t)*p2(j)*exp(si*w*t)*sin(wD*t);

%%%%%%%%%%%%怎样考虑dP2,暂时取0%%%sysdL%%%%%%%%%%

sysdL1=wD*sin(wD*t)^2*p2(j)*exp(si*w*t);

sysdL2=cos(wD*t)*0*exp(si*w*t)*sin(wD*t);

sysdL3=-cos(wD*t)*p2(j)*si*w*exp(si*w*t)*sin(wD*t);

sysdL4=-wD*cos(wD*t)^2*p2(j)*exp(si*w*t);

sysdL=sysdL1+sysdL2+sysdL3+sysdL4;

%%%%%%%%%%%%%%%%%sysdF%%%%%%%%%%%%%%

sysdF=-si*w*exp(-si*w*t)*(sysI+sysJ+sysK+sysL)/(wD*m)+exp(-si*w*t)*(sysdI+sysdJ+sysdK+sysdL)/(wD*m);

%%%%%%%%%%%%%%%%%acc%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

acc(j)=sysdC+sysdD+sysdE+sysdF;

%%%%%%%%%%%%%%%%%%%%%%%%%求加速度%%%%%%%%%%%%%%%%%%%%%%%%%

j=j+1;

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%杜哈美积分完毕%%%%%%%%%%%%%%%%%%%%%%%%%

figure;

t=0:dTaT;

subplot(3,1,1);plot(t,x);grid on;

subplot(3,1,2);plot(t,v);grid on;

subplot(3,1,3);plot(t,acc);grid on;