%单自由度系统适用，初始条件任意，比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;

杜哈梅积分，求任意激励的响应，振动力学里的，采用matlab函数！
这是我写的duhamel积分，完全是按照西安交大版《振动力学》上的公式写的，调试无误。
% duhamel        duhamel integral
% y=duhamel(x,h,fs)
% inputs:     x:analysed signal(vector)
%             h:impulse response
%             fs:frequency of sampling
% outputs:    y:duhamel integral of analysed signal x
function y=duhamel(x,h,fs)
    n=length(x);
    for i=1:n
        xx=x(1:i);
        hh=fliplr(h(1:i));
        temp=xx.*hh;
        y(i)=trapz(temp)*(1/fs);
    end
FROM CHINAVIB