平方根法和改进平方根法求解线性方程组例题与程序
数值分析15-平方根法
误差分析
定 义 ：如 果 矩 阵 A 或 常 数 项 b 的 微 小 变 化 ， 引 起 方 程 组 A x b解 的 巨 大 变 化 ， 则 称 此 方 程 组 为 “ 病 态 ” 方 程 组 ， 矩 阵 A称 为 “ 病 态 ” 矩 阵 （相对于方程组而言）。 否 则 称 方 程 组 为 “ 良 态 ” 方 程 组 ， A称 为 “良态”矩阵。 矩阵的“病态”性质是矩阵本身的特性。 为了定量刻划方程组的“病态”程度，下面 对 方 程 组 A x b就 系 数 矩 阵 或 右 端 项 分 别 有 扰 动 的两种情形进行讨论。
LU 分解唯一
A = LLT
U = LT
用平方根法解线性代数方程组的算法
（1）对矩阵A进行三角分解，即A=LLT，由矩阵乘法：
对于 i = 1, 2,…, n 计算
2 lii aii lik k 1
i 1
1 2
1 i l a l l l ij ij ik jk jj 1 k
1 2 v
A v ( v 1, 2 或
A
A
2
m ax ( A T A ) . T m in ( A A ) 1 , n
当 A为 对 称 矩 阵 时 ， c o n d ( A ) 2
其 中 1， n 为 A的 绝 对 值 最 大 和 绝 对 值 最 小 的 特 征 值 。
x A 1 b A 1 b
x x
b
b
此 式 表 明 ,当 右 端 项 有 扰 动 时 ,解 的 相 对 误 差 不 超 过 右 端 项 的 相 对 误 差 的 A A 1 倍 。
误差分析
平方根法计求解线性方程组
解线性n 阶方程组直接法—Cholesky 方法解n 阶线性方程组Ax=b 的choleskly 方法也叫做平方根法，这里对系数矩阵A 是有要求的，需要A 是对称正定矩阵，根据数值分析的相关理论，如果A 对称正定，那么系数矩阵就可以被分解为的T A=L L ∙形式，其中L 是下三角矩阵，将其代入Ax=b 中，可得：T LL x=b进行如下分解：T L x L b y y ⎧=⎨=⎩那么就可先计算y,再计算x ，由于L 是下三角矩阵，是T L 上三角矩阵，这样的计算比直接使用A 计算简便，同时你应该也发现了工作量就转移到了矩阵的分解上面，那么对于对称正定矩阵A 进行Cholesky 分解，我再描述一下过程吧： 如果你对原理很清楚那么这一段可以直接跳过的。
设T A=L L ∙，即1112111112112122221222221212....................................n n n n n n nn n n nn nn a a a l l l l a a a l l l l a a a l l l l ⎡⎤⎡⎤⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦⎣⎦⎣⎦其中,,1,2,...,ij ji a a i j n ==第1步，由矩阵乘法，211111111,i i a l a l l == 故求得111111,2,3,...i i a l l i n a === 一般的，设矩阵L 的前k-1列元素已经求出第k 步，由矩阵乘法得112211k k kk km kkik im km ik kk m m a l l a l l l l --===+=+∑∑, 于是11(2,3,...,n)1(),1,2,...kk k ik ik im km m kk l k l a l l i k k n l -=⎧=⎪⎪=⎨⎪=-=++⎪⎩∑ 注意到21kkk km m a l ==∑，于是有21max km kk ii i nl a a ≤≤≤≤ 这充分说明分解过程中元素2km l 的平方不会超过系数矩阵A 的最大对角元，因而分解过程中舍入误差的放大收到了控制，用平方根法解对称正定方程组时可以不考虑选主元问题。
平方根法
j 2, l21
i2
a a21 1 2 1, j 3, l31 31 2 l11 1 l11 1
2 l22 a22 l21 2 1 1
j 3, l32 (a32 l31l21 )
i3
1 (0 1 2) 2 l22 1
2 2 l33 a33 l31 l32 11 22
i2
1 3 1 5 , j 3, l31 a31 l11 l11 3 3
2
2 l22 a22 l21 53
j 3, l32 (a32 l31l21 )
i3
1 5 1 (9 3) 2 2 l22 3 2
a11 l11 l11 a21 l21 l11
ai 1 li 1 l11
i 1, 2 , , n
L的第一列元素 li 1可以求出 假设L的第1 ~ r 1列已求出 , 考察A的第r列元素air
2 2 arr lrk lrk lrk lrr k 1 r k 1 r r 1
3．3 平方根法 系数矩阵为对称正定矩阵的方程组称为对称正定方程组。 对 称正定方程组可用高斯消去法、LU 分解法求解，但可导出计算 量更小的平方根法。 利用对称正定矩阵的三角分解（乔累斯基分解）求解对称正 定方程组的方法称为平方根法。 3.3.1 对称正定矩阵 对称矩阵 A AT 对称正定矩阵 A AT ，且对任意非零向量 x R n 有 ( Ax, x ) x T Ax 0
air lik lrk lik lrk lir lrr
k 1 k 1
r 1
i r , r 1, , n
可得L的元素的计算公式
平方根法解线性方程组
a (k 1) ij
a / a (k1) ik
(k 1) kk
*
a(k kj
1)
,
i,
j
k
1,
,n
bi(k)
b(k 1) i
a / a (k1) ik
(k 1) kk
*bk(k1) , i k 1,
, n.
(3.9)
由（3.8）逐个回代，可得（3.1）的解
b(2) i
b (1) i
b2(1)mi2 ,
i,
j
3, 4,
, n.
重复上述过程，可以得到与（3.1）等价的上三角方程组：
a11x1
a12 x2 a12 x3
a (1) 22
x2
a (1) 22
x3
a (2) 33
x3
a1n xn b1
k 2
k 1
3
(n 1) (n 2)
2
1
1 2
n(n
1); (for
right
item)
除：(n 1) (n 2)
• 回代过程：
2 1
1 2
n(n
1)
1 2
(n
1)(multiply)
n(divide)
1 2
n(n
1).
• 运算量：1 n(n2 1) 1 n(n 1) 1 n(n 1) n3 3 n2 5 n
2. 即使高斯消元法可行，如果 ak(kk很1)小，运算中用它作分 母会导致其它元素数量计的严重增长和舍入误差的扩散。
线性方程组的四种数值解法
线性方程组的四种数值解法(电子科技大学物理电子学院，四川 成都 610054)摘要：本文介绍了四种求解线性方程组的数值解法: 雅克比迭代法、高斯赛德尔迭代法、高斯消去法和改进的平方根法的基本原理和算法流程，通过求解具体方程，对四种求解方法进行了对比。
对于雅克比迭代法和高斯赛德尔迭代法，研究了两种算法对求解同一方程组的迭代效率差异，结果表明高斯赛德尔迭代法达到同样精度所需迭代次数较少。
对于高斯消去法，通过选择列主元的方法提高算法的准确度，计算结果表明高斯消去法计算精确，且运算复杂度也不是很高。
对于改进的平方根法，其运算复杂度低，但对于给定的方程组有着严苛的要求。
关键词：雅克比迭代法；高斯赛德尔迭代法；高斯消去法；改进的平方根法；线性方程组引言线性方程组的求解在日常生活和科研中有着极其重要的应用，但在实际运算中，当矩阵的维数较高时，用初等方法求解的计算复杂度随维数的增长非常快，因此，用数值方法求解线性方程组的重要性便显现出来。
经典的求解线性方程组的方法一般分为两类：直接法和迭代法。
前者例如高斯消去法,改进的平方根法等，后者的例子包括雅克比迭代法，高斯赛德尔迭代法等。
这些方法的计算复杂度在可以接受的范围内，因此被广泛采用。
一般来说，直接法对于阶数比较低的方程组比较有效；而后者对于比较大的方程组更有效。
在实际计算中，几十万甚至几百万个未知数的方程组并不少见。
在这些情况下，迭代法有无可比拟的优势。
另外，使用迭代法可以根据不同的精度要求选择终止时间，因此比较灵活。
在问题特别大的时候，计算机内存可能无法容纳被操作的矩阵，这给直接法带来很大的挑战。
而对于迭代法，则可以将矩阵的某一部分读入内存进行操作，然后再操作另外部分。
本文使用上述四种算法求解对应的方程组，验证各种算法的精确度和计算速度。
1 算法介绍1.1 雅克比迭代法 1.1.1 算法理论设线性方程组(1)b Ax的系数矩阵A 可逆且主对角元素 均不为零,令并将A 分解成 (2)从而(1)可写成令其中. (3)以B 1为迭代矩阵的迭代法(公式)(4)称为雅克比(Jacobi)迭代法(公式),用向量的分量来表示,(4)为(5)其中为初始向量.1.1.2 算法描述 1给定迭代初始向量X 0以及误差要求delta 2根据雅克比迭代公式计算出下一组向量 3判断X 是否满足误差要求，即||X k+1 – X k || < delta4若误差满足要求，则停止迭代返回结果；若否，则返回第二步进行下一轮迭代1.2 高斯赛德尔迭代法nna ,...,a ,a 2211()nna ,...,a ,a diag D 2211=()D D A A +-=()b x A D Dx +-=11f x B x +=b D f ,A D I B 1111--=-=()()111f x B x k k +=+⎩⎨⎧[],...,,k ,n ,...,i x a ba xnij j )k (j j i iii)k (i21021111==∑-=≠=+()()()()()Tn x ,...x ,x x 002010=1.2.1 算法理论由雅克比迭代公式可知,在迭代的每一步计算过程中是用的全部分量来计算的所有分量,显然在计算第i 个分量时,已经计算出的最新分量没有被利用，从直观上看,最新计算出的分量可能比旧的分量要好些.因此，对这些最新计算出来的第次近似的分量加以利用,就得到所谓解方程组的高斯—塞德尔（Gauss-Seidel ）迭代法.把矩阵A 分解成(6)其中,分别为的主对角元除外的下三角和上三角部分,于是,方程组(1)便可以写成即其中(7)以为迭代矩阵构成的迭代法(公式)(8)称为高斯—塞德尔迭代法(公式),用变量表示的形式为(9)1.2.2 算法描述 1给定迭代初始向量X 0以及误差要求delta2根据高斯赛德尔迭代公式计算出下一组向量()k x ()1+k x ()1+k ix ()()1111+-+k i k x ,...,x 1+k()1+k x()1+k jx U L D A --=()nna ,...,a ,a diag D 2211=U ,L --A ()b Ux x L D +=-22f x B x +=()()b L D f ,U L D B 1212---=-=2B ()()221f x B x k k +=+⎩⎨⎧[],...,,k ,n ,,i x a x a b a xi j n i j )k (j ij )k (j ij i ii)k (i21021111111==∑∑--=-=+=++3判断X 是否满足误差要求，即||X k+1 – X k || < delta4若误差满足要求，则停止迭代返回结果；若否，则返回第二步进行下一轮迭代1.3 高斯消去法 1.3.1 算法理论下面三种变换称为初等行变换：1.对调两行；2.以数k ≠0乘某一行中的所有元素；3.把某一行所有元素的k 倍加到另一行对应的元素上去。
改进平方根法的C++实现
实验名称：改进平方根法学院：___数学学院______________________班级姓名：学号：实验日期 2015 年 05 月 26 日自评成绩：97一、实验目的(1)熟练掌握改进平方根法和共轭梯度法的迭代过程(2)尝试使用自己熟悉的计算机语言解决数学中的问题(3)通过上机实验来巩固课本中所学的知识二、实验内容与结果题目1：改进平方根法源程序1#include<iostream>using namespace std;int main(){double a[100][100],l[100][100],u[100][100],b[10],y[10],x[10];int i,j,k,n;cout<<"请输入矩阵的行数: ";cin>>n;cout<<"请输入右端项： "<<endl;for(j=0;j<n;j++){cin>>b[j];}cout<<"请输入矩阵： "<<endl;for(i=0;i<n;i++){for(j=0;j<n;j++){cin>>a[i][j];}}for(j=0;j<n;j++){u[0][j]=a[0][j];}for(j=1;j<n;j++){l[j][0]=u[0][j]/u[0][0];}for(i=1;i<n-1;i++){double s=0;for(k=0;k<i-1;k++){s=s+l[i][k]*u[k][i];}u[i][i]=a[i][i]-s;for(j=i+1;j<n;j++){double s=0;for(k=0;k<i-1;k++){s=s+l[i][k]*u[k][j];}u[i][j]=a[i][j]-s;l[j][i]=u[i][j]/u[i][i];}}double s=0;for(k=0;k<n-1;k++){s=s+l[n-1][k]*u[k][n-1];}u[n-1][n-1]=a[n-1][n-1]-s;y[0]=b[0];for(i=1;i<n;i++){double s=0;for(k=0;k<i-1;k++){s=s+l[i][k]*y[k];}y[i]=b[i]-s;}x[n-1]=y[n-1]/u[n-1][n-1];for(i=n-2;i>=0;i--){double s=0;for(k=i+1;k<n;k++){s=s+u[i][k]*x[k];}x[i]=(y[i]-s)/u[i][i];}cout<<"输出结果："<<endl;for(i=0;i<n;i++){cout<<"x("<<i<<")"<<"="<<x[i]<<endl;}return 0;}运行结果1。
(完整word版)线性方程组的平方根解法
浅析线性方程组的平方根解法在求解线性方程组时,直接解法有顺序高斯消元法、列主元高斯消元法、全主元高斯消元法、高斯约当消元法、消元形式的追赶法、LU 分解法、矩阵形式的追赶法，当我们遇到对称正定线性方程组时，我们就要用到平方根法(对称LLT 分解法)来求解，为了熟悉和熟练运用平方根法求解线性方程组，下面对运用平方根法求解线性方程组进行解析。
一、运用平方根法求解线性方程组涉及到的定理及定义我们在运用平方根法求解线性方程组时，要判定线性方程组Ax=b 的系数矩阵A 是否是对称正定矩阵，那么我们就要了解正定矩阵的性质和如下定理及定义：1、由线性代数知，正定矩阵具有如下性质：1) 正定矩阵A 是非奇异的2) 正定矩阵A 的任一主子矩阵也必为正定矩阵 3) 正定矩阵A 的主对角元素均为正数 4) 正定矩阵 A 的特征值均大于零 5) 正定矩阵A 的行列式必为正数定义一 线性方程组Ax=b 的系数矩阵A 是对称正定矩阵，那么Ax=b 是对称正定线性方程组。
定义二 如果方阵A 满足A=AT ，那么A 是对称阵。
2.1.4 平方根法和改进的平方根法如果A 是n 阶对称矩阵，由定理2还可得如下分解定理：定理2 若A 为n 阶对称矩阵，且A 的各阶顺序主子式都不为零，则A 可惟一分解为：A ＝LDLT ，其中L 为单位下三角阵，D 为对角阵。
证明 因为A 的各阶顺序主子式都不为零，所以A 可惟一分解为：A ＝LU 因为 ，所以可将 U 分解为：⎪⎪⎪⎪⎪⎭⎫ ⎝⎛=nn u u u U 2211⎪⎪⎪⎪⎪⎪⎪⎭⎫ ⎝⎛11122211112 u u u u u u n nn n 1DU =其中 D 为对角矩阵，U1为单位上三角阵．于是：A ＝LDU1＝L(DU1)因为A 为对称矩阵，所以，A ＝AT ＝U1TDTLT ＝U1T(DLT)，由 A 的 LU 分解的惟一性即得：L ＝U1T ，即U1＝LT ，故A ＝LDLT 。
平 方 根 法
L是单位下三角阵, U是上三角阵, 将U再分解
u11
u 22
1
u
nn
u12
u11 1
u1n
u11
u n1,n
DU 0
u n1,n1
1
其中D为对角阵, U0为单位上三角阵，于是
A = L U = L D U0
又
A = AT = U0TD LT
数值计算方法
点l11是 需由a11要此进1例, 行可开l以21方看a运l1出211 算，11。平 1为,方避根免l法31 开解al方13正11 运定12算方 ，2程我组们的改缺
A LDL 用 解ll3232单成位a3三3a2l角2321 阵ll32222作1为11分24解T1阵4，13即的l把32 形对a式称32 l，正22l3其1定l21中矩 0阵11A分2 2
l11 l21 ln1
l22
l
n
2
l
nn
按矩阵乘法展开，可逐行求出分解矩阵L的元素，计
算公式是对于i=1,2,…,n
i1
1
lii (aii li2k ) 2
k 1
i 1
a ji l jk lik
l ji
k 1
lii
j=i+1, i+2,…,n
这一方法称为平方根法,又称乔累斯基(Cholesky)分
数值计算方法
平方根法 工程实际计算中,线性方程组的系数矩阵
常常具有对称正定性，其各阶顺序主子式及 全部特征值均大于0。矩阵的这一特性使它的 三角分解也有更简单的形式，从而导出一些 特殊的解法，如平方根法与改进的平方根法。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
平方根法和改进平方根法求解线性方程组例题
与程序
2、数学原理
1、平方根法解n阶线性方程组Ax=b的choleskly方法也叫做平方根法，这里对系数矩阵A是有要求的，需要A是对称正定矩阵，根据数值分析的相关理论，如果A对称正定，那么系数矩阵就可以被分解为的形式，其中L是下三角矩阵，将其代入Ax=b 中，可得：进行如下分解：那么就可先计算y,再计算x，由于L
是下三角矩阵，是上三角矩阵，这样的计算比直接使用A计算简便，同时你应该也发现了工作量就转移到了矩阵的分解上面，那么对于对称正定矩阵A进行Cholesky分解，我再描述一下过程吧：如果你对原理很清楚那么这一段可以直接跳过的。
设，即其中第1步，由矩阵乘法，故求得一般的，设矩阵L的前k-1列元素已经求出第k步，由矩阵乘法得于是
2、改进平方根法在平方根的基础上，为了避免开方运算,所以用计算；其中，；得按行计算的元素及对元素公式对于、、计算出的第行元素后，存放在的第行相置，然后再计算的第行元素，存放在的第行、的对角元素存放在的相应位置、对称正定矩阵按分解和按分解计算量差不多，但分解不需要开放计算。
求解, 的计算公式分别如下公式。
3、程序设计
1、平方根法function
[x]=pfpf(A,b)%楚列斯基分解求解正定矩阵的线性代数方程
A=LL’ 先求LY=b 再用L’X=Y 即可以求出解
X[n,n]=size(A);L(1,1)=sqrt(A(1,1));for k=2:n
L(k,1)=A(k,1)/L(1,1);endfor 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);
endendL(n,n)=sqrt(A(n,n)-sum(L(n,1:n-1)、^2));%解下三角方
程组Ly=b 相应的递推公式如下，求出y矩阵y=zeros(n,1);%先
生成方程组的因变量的位置，给定y的初始值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矩阵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=[4,2,-4,0,2,4,0,02,2,-1,-
2,1,3,2,01,14,1,-8,-3,5,6 0,-2,1,6,-1,-4,-3,32,1,-8,-
1,22,4,-10,-34,3,-3,-4,4,11,1,-4 0,2,5,-3,-10,1,14,2
0,0,6,3,-3,-4,2,19];>> b=[0;-6;20;23;9;-22;-15;45];>>
x=pfpf(A,b)x =1
21、1481
60、1528
10、91202、018
52、改进平方根法function
[x]=improvecholesky(A,b,n)
%用改进平方根法求解Ax=bL=zeros(n,n); %L为n*n矩阵
D=diag(n,0); %D为n*n的主对角矩阵S=L*D;for i=1:n %L的主对角元素均为1 L(i,i)=1;endfor i=1:n for j=1:n %验证A是否为对称正定矩阵if (eig(A)<=0)|(A(i,j)~=A(j,i))
%A的特征值小于0或A非对称时，输出wrong
disp(wrong);break;endendendD(1,1)=A(1,1); %将A分解使得
A=LDLTfor i=2:n for j=1:i-1 S(i,j)=A(i,j)-sum(S(i,1:j-
1)*L(j,1:j-1)); L(i,1:i-1)=S(i,1:i-1)/D(1:i-1,1:i-1); end D(i,i)=A(i,i)-sum(S(i,1:i-1)*L(i,1:i-
1));endy=zeros(n,1); % x，y为n*1阶矩阵x=zeros(n,1);for i=1:n y(i)=(b(i)-sum(L(i,1:i-1)*D(1:i-1,1:i-1)*y(1:i-1)))/D(i,i); %通过 LDy=b解得y的值endfor i=n:-1:1
x(i)=y(i)-sum(L(i+1:n,i)*x(i+1:n)); %通过LTx=y解得x的值end>> A=[4,2,-4,0,2,4,0,02,2,-1,-2,1,3,2,01,14,1,-8,-
3,5,6 0,-2,1,6,-1,-4,-3,32,1,-8,-1,22,4,-10,-34,3,-3,-
4,4,11,1,-4 0,2,5,-3,-10,1,14,2 0,0,6,3,-3,-4,2,19];>>
b=[0;-6;20;23;9;-22;-15;45];>> n=8;>>
x=improvecholesky(A,b,n)x =1
21、1481
60、1528
10、91202、018
54、结果分析和讨论平方根法和改进平方根法求解线性方程组的解为x=（1
21、1481，-1
40、1127，
29、7515，-
60、1528，
10、9120，-
26、7963，
5、4259，-
2、0185）T。
与精确解相比较也存在很大的误差，虽然系数矩阵的对角元素都大于零，原则上可以不必选择主元，但由于矩阵的数值问题较大，不选主元的结果就是产生很大的误差，所以在求解的过程中还是应该选择主元以此消除误差，提高精度。
5、完成题目的体会与收获对称正定矩阵的平方根法及改进平方根法是目前解决这类问题的最有效的方法之一，合理利用的话，能够产生很好的求解效果。
改进平方根法较平方根法，因为不用进行开方运算，所以具有一定的求解优势。
通过求解此题，学会了平方根法和改进平方根法matlab编程，使我受益匪浅。