克劳特分解.txt

数值计算方法的一些题

%输入线性方程组的增广矩阵A
function X=myCrout2(A)
[n,m]=size(A);
if n~=m-1
    'error'
    break
end                                   %判定线性方程组的系数矩阵是否为方阵
[l,m]=max(abs(A(:,1)));
if A(m,1)==0
    'error'
    break
end                                   %判定线性方程组的系数矩阵是否奇异
if m~=1
    S=A(1,:);
    A(1,:)=A(m,:);
    A(m,:)=S;
end                                                            
A(1,2:n)=A(1,2:n)/A(1,1);
for k=2:n-1
    for i=k:n
        A(i,k)=A(i,k)-sum(A(i,1:k-1)*A(1:k-1,k));
    end
    [l,m]=max(abs(A(k:n,k)));
    if A(k+m-1,k)==0
       'error'
       break
    end
    if m~=1
      S=A(k,:);
      A(k,:)=A(k+m-1,:);
      A(k+m-1,:)=S;
    end
    for j=k+1:n
        A(k,j)=(A(k,j)-sum(A(k,1:k-1)*A(1:k-1,j)))/A(k,k);
    end
end
    A(n,n)=A(n,n)-sum(A(n,1:n-1)*A(1:n-1,n));
    if A(n,n)==0
        'error'
        break
    end                                %判定线性方程组的系数矩阵是否奇异
    A(1,n+1)=A(1,n+1)/A(1,1);
for k=2:n
    A(k,n+1)=(A(k,n+1)-sum(A(k,1:k-1)*A(1:k-1,n+1)))/A(k,k);
end
X(n)=A(n,n+1);
for k=n-1:-1:1
    X(k)=A(k,n+1)-sum(A(k,k+1:n).*X(k+1:n));
end