平方根法 乔累斯基分解Cholesky_解线性方程组的直接解法



标签：计算方法实验



#include <stdio.h>
#include <math.h>

const int maxn = 15;

int main(){
    double a[maxn][maxn], b[maxn], y[maxn], x[maxn], l[maxn][maxn];
    int n, sum;

    freopen("sqrt.txt", "r", stdin);
    scanf("%d", &n);
    for(int i = 1; i <= n; i++){
        for(int j = 1; j <= n; j++)  scanf("%lf", &a[i][j]);
        scanf("%lf", &b[i]);
    }
    /*
    for(int i = 1; i <= n; i++){
        for(int j = 1; j <= n; j++)  printf("%10f", a[i][j]);
        printf("%10f\n", b[i]);
    }
    */

    for(int i = 1; i <= n; i++){  //算L: 对正定矩阵A分解，A = LL^T
        sum = 0;
        for(int k = 1; k <= i - 1; k++)  sum += (l[i][k] * l[i][k]);
        l[i][i] = sqrt(a[i][i] - sum);

        sum = 0;
        for(int j = i + 1; j <= n; j++){
            for(int k = 1; k <= i - 1; k++)  sum += (l[j][k] * l[i][k]);
            l[j][i] = (a[j][i] - sum) / l[i][i];
        }
    }

    for(int i = 1; i <= n; i++){  //求y: A = LL^T -> LL^Tx = b -> Ly = b
        sum = 0;
        for(int k = 1; k <= i - 1; k++)  sum += (l[i][k] * y[k]);
        y[i] = (b[i] - sum) / l[i][i];
    }

    for(int i = n; i >= 1; i--){  //求x: L^Tx = y
        sum = 0;
        for(int k = i + 1; k <= n; k++)  sum += (l[k][i] * x[k]);
        x[i] = (y[i] - sum) / l[i][i];
    }

    for(int i = 1; i <= n; i++)  printf("%10f\n", x[i]);

    return 0;
}


数据文件 

实验结果 

定义






例子
首先判断顺序主子式。











R矩阵每一行除以对角元


平方根分解推导
因为A的分解是唯一的，可以得到$L=R^{T}$


例题











待定系数法

%定理2.2.3：对阵正定矩阵的楚列斯基(Cholesky)分解
%设A为n阶对阵正定矩阵，则存在一个可逆的下三角矩阵G，使得
%A=GG’，当限定G的对角元为正时，这种分解是唯一的
%--------A=GG’的分解算法-------
%参考教材：《数值分析》李乃成，梅立泉，科学出版社
clear;clc;
A=[9,18,9,-27;
18,45,0,-45;
9,0,126,9;
-27,-45,9,135];
b=[1 2 16 8]';
n=length(b);%方程个数n
G=zeros(n,n);
G(1,1)=sqrt(A(1,1));
G(2:n,1)=A(2:n,1)/G(1,1);
for j=2:n-1
G(j,j)=sqrt(A(j,j)-sum(G(j,1:j-1).^2));
for i=j+1:n
G(i,j)=(A(i,j)-sum(G(i,1:j-1).*G(j,1:j-1)))/G(j,j);
end
end
G(n,n)=sqrt(A(n,n)-sum((G(n,1:n-1)).^2));
G
%--------用A=GG’分解求解方程组Ax=b(平方根法)----
%Gy=b，G'x=y;
x=zeros(n,1);%未知向量
y=zeros(n,1);%中间向量
y(1)=b(1)/G(1,1);
for i=2:n
y(i)=(b(i)-sum(G(i,1:i-1)'.*y(1:i-1)))/G(i,i);
end
y
%---------由G'x=y求出x----------
%方法类似于U由x=y解出x
G=G';
x(n)=y(n)/G(n,n);
for i=n-1:-1:1
x(i)=(y(i)-sum(G(i,i+1:n)'.*x(i+1)))/G(i,i);
end
x

转自改进的平方根法

做了一些注释和修改了函数

适用于A为对称正定矩阵：A=LDL^T(L为单位下三角矩阵，D为对角矩阵，D的对角元素di>0)
function [L,D,x,y]=improvedSqareRoot(A,b)
n=length(A(:,1)); %求A矩阵第一列的长度，即矩阵维数
%% 判断输入的矩阵是不是符合要求
for k=1:n
    if (det(A(1:k,1:k)))<=0 %由于对称正定矩阵A，必有det(A)>0
        disp('A矩阵不是对称正定矩阵，请重新运行程序')
    end
end
%分解A，使A=L*D*L'
%% 初始化L（i,i) 
for i=1:n
    L(i,i)=1;
end
%% 初始化d(i) 和L(i,k)
for k=1:n
    temp_sum1=0;
    for j=1:k-1
        temp_sum1=temp_sum1+L(k,j)^2*d(j);
    end    
    d(k)=A(k,k)-temp_sum1;

    for i=k+1:n
        temp_sum2=0;
        for j=1:k-1
            temp_sum2=temp_sum2+L(i,j)*d(j)*L(k,j);
        end
        L(i,k)=(A(i,k)-temp_sum2)/d(k);
    end
end
%% 初始化t(i,j)
for k=1:n
    for i=k+1:n
        t(i,k)=L(i,k)*d(k);
    end
end

%% 转换
for k=2:n
    temp_sum3=0;
    for j=1:k-1
        temp_sum3=temp_sum3+t(k,j)*L(k,j);
    end
    d(k)=A(k,k)-temp_sum3;
    
    for i=k:n
        temp_sum4=0;
        for j=1:k-1
            temp_sum4=temp_sum4+t(i,j)*L(k,j);
        end
        t(i,k)=A(i,k)-temp_sum4;
    end
    L(i,k)=t(i,k)/d(k);
end

%% 分解Ax=b 为 Ly=b  L'x=inv(D)*y
%% 求y
for k=1:n
    temp_sum5=0;
    for j=1:k-1
        temp_sum5=temp_sum5+L(k,j)*y(j);
    end
    y(k)=b(k)-temp_sum5;
end
%% 求x
for k=n:-1:1
    temp_sum6=0;
    for j=k+1:n
        temp_sum6=temp_sum6+L(j,k)*x(j);
    end
    x(k)=y(k)/d(k)-temp_sum6;
end
D=diag(d)
end


测试数据：

A =[    1.0000   -1.0000    1.0000

-1.0000    3.0000   -2.0000

1.0000   -2.0000    4.5000];

b =[     4    -8    12];