matlab平方根法和改进平方根法求解线性方程组例题与程序.doc 
 
(2)设对称正定阵系数阵线方程组123456784202400136806224109432510634XX1,2,0,TX2、数学原理1、平方根法解N阶线性方程组AXB的CHOLESKLY方法也叫做平方根法，这里对系数矩阵A是有要求的，需要A是对称正定矩阵，根据数值分析的相关理论，如果A对称正定，那么系数矩阵就可以被分解为的形式，其中L是下三角矩TAL阵，将其代入AXB中，可得TLXB进行如下分解TLXBY那么就可先计算Y,再计算X，由于L是下三角矩阵，是上三角矩阵，这TL样的计算比直接使用A计算简便，同时你应该也发现了工作量就转移到了矩阵的分解上面，那么对于对称正定矩阵A进行CHOLESKY分解，我再描述一下过程吧如果你对原理很清楚那么这一段可以直接跳过的。设，即TL1211121222121NNNNNNNAALLLLLL其中,,IJJIA第1步，由矩阵乘法，故求得211,IIALLA11,,3,IILALN一般的，设矩阵L的前K1列元素已经求出第K步，由矩阵乘法得112KKMKIIMIKALALL,于是1212,3N,1,KKMKIKIILALKLLIKN2、改进平方根法在平方根的基础上，为了避免开方运算,所以用计算；其中，TLDA；11112NNNDDDD得11212122NNNDLLLAL按行计算的元素及对元素公式LD对于I,211,21JIJIKJTATLI，/,IJJLTD，I1IIKIATL计算出的第行元素后，存放在的第行相置，LDTI,2IIJT，AI然后再计算的第行元素，存放在的第行的对角元素存放在的相应AD位置对称正定矩阵按分解和按分解计算量差不多，但分解不ATLDTLTLD需要开放计算。求解,的计算公式分别如下公式。BYYX11,IIKHL2,IN1/,NNIIKIXYDLX1,IN3、程序设计1、平方根法FUNCTIONXPFPFA,B楚列斯基分解求解正定矩阵的线性代数方程ALL’先求LYB再用L’XY即可以求出解XN,NSIZEAL1,1SQRTA1,1FORK2NLK,1AK,1/L1,1ENDFORK2N1LK,KSQRTAK,KSUMLK,1K12FORIK1NLI,KAI,KSUMLI,1K1LK,1K1/LK,KENDENDLN,NSQRTAN,NSUMLN,1N12解下三角方程组LYB相应的递推公式如下，求出Y矩阵YZEROSN,1先生成方程组的因变量的位置，给定Y的初始值FORK1NJ1K1YKBKLK,JYJ/LK,KEND解上三角方程组L’XY递推公式如下，可求出X矩阵XZEROSN,1UL 求上对角矩阵FORKN11JK1NXKYKUK,JXJ/UK,KENDA4,2,4,0,2,4,0,02,2,1,2,1,3,2,04,1,14,1,8,3,5,60,2,1,6,1,4,3,32,1,8,1,22,4,10,34,3,3,4,4,11,1,40,2,5,3,10,1,14,20,0,6,3,3,4,2,19B0620239221545XPFPFA,BX1211481140112729751560152810912026796354259201852、改进平方根法FUNCTIONXIMPROVECHOLESKYA,B,N用改进平方根法求解AXBLZEROSN,NL为NN矩阵DDIAGN,0D为NN的主对角矩阵SLDFORI1NL的主对角元素均为1LI,I1ENDFORI1NFORJ1N验证A是否为对称正定矩阵IFEIGAA4,2,4,0,2,4,0,02,2,1,2,1,3,2,04,1,14,1,8,3,5,60,2,1,6,1,4,3,32,1,8,1,22,4,10,34,3,3,4,4,11,1,40,2,5,3,10,1,14,20,0,6,3,3,4,2,19B0620239221545N8XIMPROVECHOLESKYA,B,NX1211481140112729751560152810912026796354259201854、结果分析和讨论平方根法和改进平方根法求解线性方程组的解为X(1211481，1401127，297515，601528，109120，267963，54259，20185)T。与精确解相比较也存在很大的误差，虽然系数矩阵的对角元素都大于零，原则上可以不必选择主元，但由于矩阵的数值问题较大，不选主元的结果就是产生很大的误差，所以在求解的过程中还是应该选择主元以此消除误差，提高精度。5、完成题目的体会与收获对称正定矩阵的平方根法及改进平方根法是目前解决这类问题的最有效的方法之一，合理利用的话，能够产生很好的求解效果。改进平方根法较平方根法，因为不用进行开方运算，所以具有一定的求解优势。通过求解此题，学会了平方根法和改进平方根法MATLAB编程，使我受益匪浅。

 
平方根法
 
function [x]=pingfg(A,b)
 
%喬累斯基分解
 
[n,n]=size(A);
 
L=zeros(n,n);%實際上不用為 L 申請空間，使用 A 即可
 
L(1,1)=sqrt(A(1,1));
 
for k=2:n
 
L(k,1)=A(k,1)/L(1,1);
 
end
 
for k=2:n-1
 
L(k,k)=sqrt(A(k,k)-sum(L(k,1:k-1).^2));
 
for i=k+1:n
 
L(i,k)=(A(i,k)-sum(L(i,1:k-1).*L(k,1:k-1)))/L(k,k);
 
end
 
end
 
L(n,n)=sqrt(A(n,n)-sum(L(n,1:n-1).^2));
 
%解下三角方程組 Ly=b
 
y=zeros(n,1);
 
for k=1:n
 
j=1:k-1;
 
y(k)=(b(k)-L(k,j)*y(j))/L(k,k);
 
end
 
%解上三角方程組  L'x=y
 
x=zeros(n,1);
 
U=L';
 
for k=n:-1:1
 
j=k+1:n;
 
x(k)=(y(k)-U(k,j)*x(j))/U(k,k);
 
end

满足各阶顺序主子式非零的非奇异矩阵 $A$ 存在唯一的LU分解，如果矩阵 $A$ 的结构再特殊一点，线性方程组 $Ax=b$ 的求解算法是否还可以得到优化呢。
考虑对称矩阵以及对称正定矩阵，这两者在实际科学计算中出现频率是较高的。
 


 
正文
 
对称阵.
正定阵.
平方根法.

 
对称阵.
 
设 $A$ 是 $n$ 阶对称阵，并且满足各阶顺序主子式非零，那么 $A$ 可以唯一分解为$$A=LD{L}^T$$其中 $D$ 是对角阵，$L$ 是单位下三角阵。
【简证】首先 $A$ 可以进行LU分解，即 $A=LU$，由于 $D_k={\prod }_{i=1}^k{u}_{ii}\ne 0$，因此 $u_{ii}\ne 0,i=1,2,...,n$，所以可以将 $U$ 分解为 $$U=D{U}^{\prime }$$
      
       
        
         
          D
         
         
          =
         
         
          d
         
         
          i
         
         
          a
         
         
          g
         
         
          {
         
         
          
           u
          
          
           11
          
         
         
          ,
         
         
          
           u
          
          
           22
          
         
         
          ,
         
         
          .
         
         
          .
         
         
          .
         
         
          ,
         
         
          
           u
          
          
           
            n
           
           
            n
           
          
         
         
          }
         
        
        
         D=diag\{u_{11},u_{22},...,u_{nn}\}
        
       
      D=diag{u11​,u22​,...,unn​}因此$$A=LD{U}^{\prime }$$
由于 $A$ 是对称阵，所以 $$A=A^T=(LD{U}^{\prime }{)}^T=U^{\prime T}{D}^T{L}^T=U^{\prime T}D{L}^T=LD{U}^{\prime }$$根据LU分解唯一性，必然有 $U^{\prime }=L^T$，即$$A=LD{L}^T$$
上述结论的意义在于，对于对称阵的LU分解，矩阵 $L,U$ 中的元素存在对称性，大约仅需要计算 $\frac{n(n+1)}{2}$ 个，相较于一般情况下 $n^2$ 个元素减少了一半。
 
正定阵.
 
 
 
【正定阵】狭义定义：对于 $n$ 阶实对称矩阵 $A$ 而言，当且仅当 $\forall z\in R^n,z\ne 0$，都有 $z^{T}Az>0$，那么称 $A$ 是正定矩阵。正定矩阵的各阶顺序主子式为正，其逆矩阵也是正定矩阵。
 
 
【定理】对于正定矩阵 $A$ 而言，存在唯一分解$$A=L{L}^T$$其中 $L$ 是下三角阵，且主对角线元素都为正数。
【简证】根据上部分对称阵的LU分解，此处有 $A=L_{1}D{L}_1^T$，其中 $L_1$ 是单位下三角阵，又因为 $u_{ii}=\frac{D_k}{D_{k-1}}>0$，因此令
      
       
        
         
          
           D
          
          
           ′
          
         
         
          =
         
         
          d
         
         
          i
         
         
          a
         
         
          g
         
         
          {
         
         
          
           
            u
           
           
            11
           
          
         
         
          ,
         
         
          
           
            u
           
           
            22
           
          
         
         
          ,
         
         
          .
         
         
          .
         
         
          .
         
         
          ,
         
         
          
           
            u
           
           
            
             n
            
            
             n
            
           
          
         
         
          }
         
        
        
         D'=diag\{\sqrt{u_{11}},\sqrt{u_{22}},...,\sqrt{u_{nn}}\}
        
       
      D′=diag{u11​
             ​,u22​
             ​,...,unn​
             ​}则有$$L_{1}D{L}_1^T=L_1{D}^{\prime }{D}^{\prime }{L}_1^T=(L_1{D}^{\prime })(L_1{D}^{\prime }{)}^T=L{L}^T$$
上述分解称为正定矩阵的乔列斯基分解，矩阵形式表示如下：$$A=\left[\begin{array}{cccc}{l}_{11}& 0& \cdots & 0\\ l_{21}& l_{22}& \cdots & 0\\ ⋮& ⋮& \ddots & ⋮\\ l_{n1}& l_{2n}& \cdots & l_{nn}\end{array}\right]\cdot \left[\begin{array}{cccc}{l}_{11}& l_{21}& \cdots & l_{n1}\\ 0& l_{22}& \cdots & l_{n2}\\ ⋮& ⋮& \ddots & ⋮\\ 0& 0& \cdots & l_{nn}\end{array}\right]$$
考察矩阵 $A$ 的原始元素$$\left[\begin{array}{cccc}{a}_{11}& a_{12}& \cdots & a_{1n}\\ a_{21}& a_{22}& \cdots & a_{2n}\\ ⋮& ⋮& \ddots & ⋮\\ a_{n1}& a_{n2}& \cdots & a_{nn}\end{array}\right]$$可以发现$$a_{ij}=l_{jj}{l}_{ij}+\sum _{k=1}^{j-1}{l}_{ik}{l}_{jk}$$
于是得到Cholesky分解公式：$$l_{jj}=\sqrt{a_{jj}-\sum _{k=1}^{j-1}{l}_{jk}^2},j=1,2,...,n$$$$l_{ij}=\frac{a_{ij}-{\sum }_{k=1}^{j-1}{l}_{ik}{l}_{jk}}{l_{jj}},i=j+1,j+2,...,n$$
 
平方根法.
 
根据上面得出的Cholesky分解参数 $l_{ij}$，可以解出正定线性方程组 $Ax=b$ 的解为$$y_i=\frac{b_i-{\sum }_{k=1}^{i-1}{l}_{ik}{y}_k}{l_{ii}},i=1,2,...,n$$$$x_i=\frac{y_i-{\sum }_{k=i+1}^n{l}_{ki}{x}_k}{l_{ii}},i=n,n-1,...,1$$
 
 
上述方法得名于Cholesky分解过程中计算 $l_{ii}$ 的时候需要进行开方运算，如果为了避免开方运算，也可以不对矩阵 $D$ 做分解，而是将其乘入 $L$ 中，得到如下表达形式$$A=(LD)L^T=\left[\begin{array}{cccc}{d}_1& 0& \cdots & 0\\ d_1\cdot l_{21}& d_2& \cdots & 0\\ ⋮& ⋮& \ddots & ⋮\\ d_1\cdot l_{n1}& d_2\cdot l_{2n}& \cdots & d_n\end{array}\right]\cdot \left[\begin{array}{cccc}1& l_{21}& \cdots & l_{n1}\\ 0& 1& \cdots & l_{n2}\\ ⋮& ⋮& \ddots & ⋮\\ 0& 0& \cdots & 1\end{array}\right]$$
在这种表示形式下：$$a_{ij}=d_j\cdot l_{ij}+\sum _{k=0}^{j-1}{d}_k\cdot l_{ik}\cdot l_{jk}$$$$l_{jj}=1$$
由此可以得到 $$d_j=a_{jj}-d_k\cdot \sum _{k=0}^{j-1}{l}_{jk}^2,j=1,2,...,n$$$$l_{ij}=\frac{a_{ij}-{\sum }_{k=0}^{j-1}{d}_k\cdot l_{ik}\cdot l_{jk}}{d_j}m,i=j+1,j+2,...,n$$
于是线性方程组 $Ax=b$ 转化为求解 $Ly=b,D{L}^T=y$，方程的解为$$\{\begin{array}{rl}& y_1=b_1\\ & y_i=b_i-\sum _{k=1}^{i-1}{l}_{ik}{y}_k,i=2,3,...,n\end{array}$$$$\{\begin{array}{rl}& x_n=\frac{y_n}{d_n}\\ & x_i=\frac{y_i}{d_i}-\sum _{k=i+1}^n{l}_{ki}{x}_k,i=n-1,n-2,...,1\end{array}$$
这种方法被称为改进的平方根法，避免了开方运算。

平方根法 乔累斯基分解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;
}
 
数据文件 
  
 实验结果