cholesky.m

用matlab编写的cholesky分解程序

N=200;
M=100
dt=0.1;
A=zeros(2*M,2*M);
x=zeros(1,M);
y=zeros(1,M);
for i=1:N
    for m=1:M
        x(m)=100+50*cos(2*pi*2.5*(i+m-1)*dt).*exp(-0.2*(i+m-1)*dt)+20*cos(2*pi*1*(i+m-1)*dt).*exp(0.1*(i+m-1)*dt);  
        y(m)=100+50*cos(2*pi*2.5*(i+m)*dt).*exp(-0.2*(i+m)*dt)+20*cos(2*pi*1*(i+m)*dt).*exp(0.1*(i+m)*dt);
    end
        z=[x,y];
        Z=z'*z;
        A=A+(1/N)*Z;            %构造正则矩阵对之一Rzz
end
n=length(A)
lambda=eig(A);
flagdc=0;
flagzd=0;
while(flagdc~=n*(n-1)/2)|(flagzd~=n)
    for i=1:n
        for j=i+1:n
            if(A(i,j)==A(j,i))
                flagdc=flagdc+1;
            end
        end
    end
    if(flagdc~=n*(n-1)/2)
        disp('矩阵A不是对称矩阵，请重新输入矩阵A:')
        A=input('A=\n');
        disp('矩阵的阶数为:')
        n=length(A)
        flagdc=0;
    end
    if flagdc==n*(n-1)/2
        disp('矩阵A的特征值是:')
        lambda=eig(A)
        for i=1:n
            if(lambda(i)>0)
                flagzd=flagzd+1;
            end
        end
        if flagzd~=n
            disp('矩阵A不是正定矩阵。请重新输入矩阵A:')
            A=input('A=\n');
            disp('矩阵的阶数为：')
            n=length(A)
            flagzd=0;
            flagdc=0;
        end
    end
end
l(1,1)=sqrt(A(1,1));
for j=2:n
    l(j,1)=A(j,1)/l(1,1);
end
for i=2:n
    for j=i:n
        p(j,i)=0;
        for k=1:i-1
            p(j,i)=p(j,i)+l(i,k)*l(j,k);
            t(j,i)=A(j,i)-p(j,i);
            if (j==i)
                l(j,i)=A(j,i)-p(j,i);
                if (j==i)
                    l(j,i)=sqrt(t(j,i));
                else
                    l(j,i)=t(j,i)/l(i,i);
                end
            end
        end
    end
end