虚拟激励法程序求助

hustdd

学习林家浩教授和张亚辉教授的专著《随机振动的虚拟激励法》时，编制的程序和书上的结果对不上，请大家帮忙看看。
clear;
clc;
n=4;
II=sqrt(-1);
%主结构质量、阻尼、刚度矩阵

M=eye(n)*1.0e+4;
K=eye(n)*1.6*1.0e+7;
%主结构刚度矩阵聚合
zk=zeros(n);

for j=1:(n-1)
zk(j,j)=K(j,j)+K(j+1,j+1);
zk(j,j+1)=-K(j+1,j+1);
zk(j+1,j)=-K(j+1,j+1);
end
zk(n,n)=K(n,n);
k=zk;
m=M;
%求解各阶模态频率
[tzxl,tzz]=eig(k,m);
d=diag(sqrt(tzz));
%振型规一化
for i=1:n
   [tzz1(i),j]=min(d);
   tzxl1(:,i)=tzxl(:,j);
   d(j)=max(d)+1;
end
%振型归一化取第一层振型为1
for j=1:n
   tzxl1(:,j)=tzxl1(:,j)/tzxl1(1,j);
end

w0=tzz1;
w=tzz1/(2*pi);
zhx=tzxl1;
广义阻尼矩阵   
各阶模态阻尼比都取
%阻尼比
ks0=0.05;
ks=ones(n,1)*ks0;
第n阶广义质量：   
%求广义质量
Mn=zhx'*m*zhx;
阻尼矩阵为：
%求阻尼矩阵
C=zeros(n);
for i=1:n
   C(i,i)=2*ks(i)*w0(i)*Mn(i,i);
end
c=(zhx')\C/zhx;
参数
eg即
%过滤白噪声参数
eg=0.6;
wg=15.708;
S0=0.001574;
%按照书上的要求，取频率和时间的最大值和步长
%频率间隔
dw=0.3;
%最大频率范围
maxw=45;
%最大时间值
maxt=40;
%时间间隔
dt=0.2;

%各层各时间点频率点的功率谱密度,循环变量：层数，时间点，频率点
Pwt=zeros(n,maxt/dt,maxw/dw);
%频率点数循环变量wn
wn=1;

%对频率进行循环，求解各频率点的时间历程
for w=0:dw:maxw
x1=1+4*eg^2*(w/wg)^2;
x2=(1-(w/wg)^2)^2+4*eg^2*(w/wg)^2;
Sgw=x1*S0/x2;
s=sqrt(Sgw);
%采用精细积分法进行求解时间历程，得到位移和速度时程
[disp,velp]=JINGXI67(M,zk,c,dt,maxt,w,s,n);
Ywt=disp;
for kkk=1:maxt/dt
%求确定频率下各时间点的功率谱
Yw=Ywt(:,kkk);
  每一时刻和频率点的位移向量，对其进行求共轭和装置得到协方差矩阵，对角上的元素即是每一时刻的各层的功率谱
y1=conj(Yw);
y2=transpose(Yw);

%确定时间点确定频率下的功率谱Yw,取对角线元素
Syyw=y1*y2;
   for kk=1:n
   Pwt(kk,kkk,wn)=Syyw(kk,kk);
   end
end
wn=wn+1;
end






%求解完成后实际上wn为maxw/dw+2，实际频率点个数为maxw/dw+1
%各层的时变方差，循环变量为：层数，时间点
Fangcha=zeros(n,maxt/dt);
for tn=1:maxt/dt
%求解各层的时变方差
for kk=1:n
xx1=zeros(wn-1,1);
%每一个时刻的方差对各频率点进行积分，频率点数取maxw/dw+1，即wn-1
   for wn0=1:wn-1
   xx1(wn0)=Pwt(kk,tn,wn0);
   end
%采用复合梯形求积公式对功率谱进行积分得到方差
Fangcha(kk,tn)=(xx1(1)+xx1(wn-1)+2*sum(xx1(2:wn-1-1)))*dw;
end
end

%画图
c1=(1:maxt/dt)*dt;
d1=Fangcha(1,:)/S0;
d2=Fangcha(2,:)/S0;
d3=Fangcha(3,:)/S0;
d4=Fangcha(4,:)/S0;
figure(3)
plot(c1,d1,'k',c1,d2,'r',c1,d3,'m',c1,d4,'r-')



精细积分的程序
function [disp,velp]=JINGXI67(m,k,c,dt,maxt,w,s,n)
%虚数单位
II=sqrt(-1);
%  中的
IIW=II*w;
I=eye(n);
Z=zeros(n);
离散化n自由度结构受均匀调制演变随机激励 时的运动微分方程可表示为：

其中 为平稳高斯白噪声随机过程向量， 为调制函数。该方程可表示为

其中I为单位矩阵，记
， ，

程序中的A即上面的矩阵A，AF即G
A=[Z,I;-m\k,-m\c];
AF=-1*[zeros(n,1);ones(n,1)];
对线性体系来说，方程为线性向量常微分方程其解为

对上式进行离散化 为时间步长

矩阵T 可用数值方法精细算得 采用函数EXPHDT计算，H即输入的矩阵A，dt即为时间间隔
%采用科茨公式求解非齐次项
%求解积分点的矩阵指数
EXPDT=EXPHDT(A,dt);
EXP0_75DT=EXPHDT(A,0.75*dt);
EXP0_50DT=EXPHDT(A,0.50*dt);
EXP0_25DT=EXPHDT(A,0.25*dt);
第二项 等价于单个函数的积分，采用科茨公式进行积分，具有5阶精度


zeroq=zeros(n,maxt/dt);
disp=zeroq;  
velp=zeroq;  

zerop=zeros(n,1);
dispt=zerop;
velpt=zerop;
UK=[dispt;velpt];

for i=1:maxt/dt
%给出各积分点的X坐标，即时间
t= ；t0_25= ；t0_50= ；t0_75= ；t1=
t=(i-1)*dt;
t0_25=t+0.25*dt;
t0_50=t+0.50*dt;
t0_75=t+0.75*dt;
t1=i*dt;
%调制函数
求解gt= ；gt0_25= ；gt0_50= ；gt0_75= ；gt1= 。
gt=4*(exp(-0.0995*t)-exp(-0.199*t));
gt0_25=4*(exp(-0.0995*t0_25)-exp(-0.199*t0_25));
gt0_50=4*(exp(-0.0995*t0_50)-exp(-0.199*t0_50));
gt0_75=4*(exp(-0.0995*t0_75)-exp(-0.199*t0_75));
gt1=4*(exp(-0.0995*t1)-exp(-0.199*t1));
即AF*exp(IIW*t)； 即AF*exp(IIW*t0_25)
FTK=AF*exp(IIW*t)*gt*s;
FT0_25K=AF*exp(IIW*t0_25)*gt0_25*s;
FT0_50K=AF*exp(IIW*t0_50)*gt0_50*s;
FT0_75K=AF*exp(IIW*t0_75)*gt0_75*s;
FTK1=AF*exp(IIW*t1)*gt1*s;

% 即TUK= EXPDT*UK   EXPDT即
TUK=EXPDT*UK;
%采用科茨公式求解非齐次项

X1=7*EXPDT*FTK;
X2=32*EXP0_75DT*FT0_25K;
X3=12*EXP0_50DT*FT0_50K;
X4=32*EXP0_25DT*FT0_75K;
X5=7*FTK1;
TKTK1=1/90*dt*(X1+X2+X3+X4+X5);
即TKTK1

UK1=TUK+TKTK1;
%保存各时刻的位移和速度向量
disp(:,i)=UK1(1:n,:);
velp(:,i)=UK1(n+1:2*n,:);
%将UK1赋给UK进入下一轮循环
UK=UK1;
end




求解
function [T]=EXPHDT(A,dt)
%计算矩阵指数
N=20;
m=2^N;
t=dt/m;
I=eye(size(A));
Ta=A*t+(A*t)^2*(I+(A*t)/3+(A*t)^2/12)/2;
%计算指数矩阵T
for i=1:20
  Ta=2*Ta+Ta^2;
end
T=I+Ta;

[ 本帖最后由 hustdd 于 2008-12-14 11:48 编辑 ]

星矢

楼主说“程序中的A即上面的矩阵A，AF即G”，但实际上程序中写的G为[0;I]，未包含g(t)，但楼主在程序中积分时是包含了g(t)，不知道有没有影响。

wuzhenyu522

楼主和楼上的能不能留个联系方式，大家交流一下
我的QQ353999102
E-mail：wzy5221027@163.com

Tony1027

我们学部的林家浩教授，呵呵。。。

sunjianpeng2001

好帖子，学习了！

星矢

好久没上，我的邮箱：iamcsw@126.com，讨论讨论

liqing324

顶一个，有用

chnzhn

学习了，还是蛮有用的

yxliuliu

这本书我也有看过！很麻烦

zhangyi1169

欧阳李

内容不错，谢谢分享

mzbys

需要逆虚拟激励法 这个怎么编呢

mzbys

需要逆虚拟激励法 这个怎么编呢？

zizhubingxian

我怎么弄不明白呢