线性滤波发生成风速时程-正定问题
有谁看过徐赵东（土木工程常用软件分析与应用）中模拟风荷载那块%线性滤波法生成风速时程
clc;
clear;
p=4;
Np=5;
Vz=50;
T=500;
dT=0.1;
load 'coor3.txt'
XYZ=coor3;
X=coor3(:,2);
Y=coor3(:,3);
Z=coor3(:,4);
k=0.4;
z0=0.03;
Ustar=k*Vz/log(10/z0);
for n=1:Np
U(n)=1./k*Ustar*log(XYZ(n,4)/z0);
end
%由功率谱估算协方差
r=T/dT;
RR=[];
for s=0:p
for n=1:Np
for m=1:n
h=2*sqrt(16^2*(X(n)-X(m)).^2+8^2*(Y(n)-Y(m)).^2+10^2*(Z(n)-Z(m)).^2)/(U(n)+U(m));
F=@(x)200.*Ustar^2.*(Z(n)-Z(m))./U(n)./U(m)./(1.+50.*x*Z(n)./U(n)).^(5./3)./(1.+50.*x*Z(m)./U(m)).^(1./2).*exp(-1.*h*x).*cos(2.*pi*x*s*dT);
Q(m,n)=quadl(F,0,10);
Q(n,m)=Q(m,n);
end
end
RR=;
end
%求解公式系数，
C=[];
for n=1:p
CC=[];
for m=1:p
idx=abs(n-m);
if idx==0
C1=RR(:,1:Np);
else
C1=RR(:,(idx+1)*Np+1:(idx+2)*Np);
end
CC=;
end
C=;
end
%求解公式左边的系数R
D=[];
for n=1:p
DD=RR(:,Np*n+1:Np*(n+1));
D=;
end
CCC=inv(C);
phi=CCC*D;%求解自回归系数
%求解系数矩阵公式
RRO=RR(:,1:Np);
Ruu=RR(:,Np+1:Np*(p+1));
BO=chol(RRO-Ruu*phi);
运行到这的乔斯莱斯分解 要求矩阵正定，我该怎么办呢？
有懂的 大侠 希望 给解释 或分享程序
页:
[1]