线性滤波发生成风速时程-正定问题

magic3794

有谁看过徐赵东（土木工程常用软件分析与应用）中模拟风荷载那块
%线性滤波法生成风速时程
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=[RR,Q];
   
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=[CC;C1];
    end
    C=[CC,C];
end
%求解公式左边的系数R
D=[];
for n=1:p
    DD=RR(:,Np*n+1:Np*(n+1));
    D=[D;DD];
end
CCC=inv(C);
phi=CCC*D;  %求解自回归系数
%求解系数矩阵公式
RRO=RR(:,1:Np);
Ruu=RR(:,Np+1:Np*(p+1));
BO=chol(RRO-Ruu*phi);

运行到这的乔斯莱斯分解 要求矩阵正定，我该怎么办呢？
有懂的 大侠 希望 给解释 或分享程序