2线性方程组的数值解法.docx
- 文档编号:3877148
- 上传时间:2023-05-06
- 格式:DOCX
- 页数:32
- 大小:395.69KB
2线性方程组的数值解法.docx
《2线性方程组的数值解法.docx》由会员分享,可在线阅读,更多相关《2线性方程组的数值解法.docx(32页珍藏版)》请在冰点文库上搜索。
2线性方程组的数值解法
第二章线性方程组的数值解法
§2.0引言
在自然科学和工程技术中很多问题的解决常常归结为解线性代数方程组,例如电学中的网络问题,船体数学放样中建立三次样条函数问题,用最小二乘法求实验数据的曲线拟合问题,解非线性方程组问题,用差分法或者有限元方法解常微分方程、偏微分方程边值问题等都导致求解线性代数方程组,而这些方程组的系数矩阵大致分为两种,一种是低阶稠密矩阵(例如,阶数大约为≤150),另一种是大型稀疏矩阵(即矩阵阶数高且零元素较多)。
设有线性方程组Ax=b,其中
为非奇异阵,
,
关于线性方程组的数值解法一般有两类:
直接法与迭代法。
1.直接法
就是经过有限步算术运算,可求得方程组精确解的方法(若计算过程中没有舍入误差)。
但实际计算中由于舍入误差的存在和影响,这种方法也只能求得线性方程组的近似解。
本章将阐述这类算法中最基本的高斯消去法及其某些变形。
这类方法是解低阶稠密矩阵方程组的有效方法,近十几年来直接法在求解具有较大型稀疏矩阵方程组方面取得了较大进展。
2.迭代法
就是用某种极限过程去逐步逼近线性方程组精确解的方法。
迭代法具有需要计算机的存贮单元较少、程序设计简单、原始系数矩阵在计算过程中始终不变等优点,但存在收敛性及收敛速度问题。
迭代法是解大型稀疏矩阵方程组(尤其是由微分方程离散后得到的大型方程组)的重要方法。
§2.1Gauss消去法
高斯(Gauss)消去法是解线性方程组最常用的方法之一,它的基本思想是通过逐步消元(行的初等变换),把方程组化为系数矩阵为三角形矩阵的同解方程组,然后用回代法解此三角形方程组(简单形式)得原方程组的解。
1.消元
例:
下面我们来讨论一般的解n阶方程组的高斯消去法。
将原方程组记为A
(1)x=b
(1),其中A
(1)=(aij
(1))=(aij),b
(1)=b
(1)第一次消元。
其中
注:
若a11
(1)=0,则在第1列中至少有一个元素不为0,可交换行后再消元
(2)第k次消元。
其中
注:
为减少计算量,令
,则
(3)当k=n–1时得
完成第n–1次消元后得到与原方程组等价的三角形方程组A(n)x=b(n)
注:
当det(A)≠0时,显然有aii(i)≠0,(i=1,…,n),称为主元素。
2.回代
求解上述三角形方程组,得到求解公式
注:
求解过程称为回代过程。
3.Gauss消去法的计算量
以乘除法的次数为主
(1)消元过程:
第k步时(n–k)+(n–k)(n–k+1)=(n–k)(n–k+2)
共有
(2)回代过程:
求xi中,乘n–i次,除1次,共n–i+1次(i=1,…,n–1)
共有
(3)总次数为
注:
当n=20时约为2670次,比克莱姆法则9.71020次大大减少。
4.说明
(1)若消元过程中出现akk(k)=0,则无法继续
(2)若akk(k)≠0,但较小,则小主元做除数将产生大误差
(3)每次消元都选择绝对值最大者作主元,称为高斯主元消去法
(4)通常第k步选择
——第k列主对角元以下元素绝对值最大者作主元(该行与第k行对调),称为列主元消去法。
例1用舍入三位有效数字求解线性方程组
(准确解是x1=10.0,x2=1.00)
解:
1)不选主元的Gauss消去法计算结果:
x1=-10.0,x2=1.01,此解无效;
2)按列选主元的Gauss消去法计算结果:
x1=10.0,x2=1.00.
(5)行标度化
当线性方程组的系数矩阵的元素大小相差很大时,则可能引进因丢失有效数字而产生的误差,并且舍入误差的影响严重,即使用Gauss主元消去法得到的解也不可靠.为避免这一问题,可将系数矩阵的行元素按比例增减以使元素的变化减小.例如,用每行元素的最大模除该行各元素,使它们的模都不大于1,这叫行标度化,标度化的目的是要找到真正的主元.消元过程仍是对原方程组进行的,只不过在Gauss列主元消去法的算法中,按列选主元ck时,应修改为
其中
5.算法设计
高斯消去法
a=[5-1-1-1-4;
-110-1-112;
-1-15-18;
-1-1-11034];
x=[0;0;0;0]
n=4
fork=1:
3
fori=k+1:
n
a(i,:
)=a(i,:
)-a(k,:
)*a(i,k)/a(k,k)
end
end
forj=n:
-1:
1
x(j)=(a(j,n+1)-a(j,j+1:
n)*x(j+1:
n))/a(j,j)
end
列主元素消去法
a=[-0.040.040.123;
0.56-1.560.321;
-0.241.24-0.280];
x=[0;0;0]
n=3;
fork=1:
n-1
[c,i]=max(abs(a(k:
n,k)));
q=i+k-1
ifq~=k
m=a(q,:
);a(q,:
)=a(k,:
);a(k,:
)=m
end
fori=k+1:
n
a(i,:
)=a(i,:
)-a(k,:
)*a(i,k)/a(k,k)
end
end
forj=n:
-1:
1
x(j)=(a(j,n+1)-a(j,j+1:
n)*x(j+1:
n))/a(j,j)
end
§2.2矩阵分解在解线性方程组中的作用
1.原理
若A=LR,则Ax=bLRx=b
若其中L、R为三角形矩阵,则方程组易解
定义:
(1)称
为单位下三角阵
(2)设L为单位下三角阵,R为上三角阵,若A=LR,则称A可三角(LR)分解,并称A=LR为A的三角分解(杜利特尔分解)
定理:
(1)A=(aij)nn有唯一LR分解A的顺序主子式k≠0,k=1,2,...,n–1
(2)若A=(aij)nn可逆,则存在置换阵P,使PA的n个顺序主子式全不等于0
注:
由Ax=bPAx=Pb知,经适当行交换后,A总可三角分解
2.直接三角分解法
设A的各阶顺序主子式不为0:
A=LR,即
(1)主对角线(含)上边:
当列标m>行标k时,lkm=0
j=k,k+1,...,n;k=1,2,...,n
(2)主对角线下边:
当行标m>列标k时,rkm=0
i=k+1,k+2,...,n;k=1,2,...,n–1
欲求lik和rkj
设k=1,a1j=l11r1jr1j=a1jj=1,2,...,n
ai1=li1r11li1=ai1/r11i=2,3,...,n
一般地,
最后:
r11
r12
...
r1n
第1步
l21
r22
...
r2n
第2步
l31
l32
...
...
ln1
ln2
rnn
第n步
3.求解方程组
下三角Ly=b:
上三角Rx=y:
紧凑格式:
a=[2100-3;
-3-4-1213;
123-4;
4149-13];
b=[10;5;-2;7]
y=[0000]'
x=[0000]'
n=4;
forr=1:
n
a(r,r:
n)=a(r,r:
n)-a(r,1:
r-1)*a(1:
r-1,r:
n)
a(r+1:
n,r)=(a(r+1:
n,r)-a(r+1:
n,1:
r-1)*a(1:
r-1,r))/a(r,r)
end
R=triu(a,0)
L=eye(n)+tril(a,-1)
forr=1:
n
y(r)=b(r)-L(r,1:
r-1)*y(1:
r-1)
end
forj=n:
-1:
1
x(j)=(y(j)-R(j,j+1:
n)*x(j+1:
n))/R(j,j)
end
紧凑格式
a=[2100-310;
-3-4-12135;
123-4-2;
4149-137];
x=[0000]'
n=4;
forr=1:
n
a(r,r:
n+1)=a(r,r:
n+1)-a(r,1:
r-1)*a(1:
r-1,r:
n+1)
a(r+1:
n,r)=(a(r+1:
n,r)-a(r+1:
n,1:
r-1)*a(1:
r-1,r))/a(r,r)
end
R=triu(a,0)
forj=n:
-1:
1
x(j)=(a(j,n+1)-R(j,j+1:
n)*x(j+1:
n))/R(j,j)
end
4.选主元的三角分解法
5.平方根法
(1)定理:
若A正定,则A可唯一分解为A=LDLT。
其中L为单位下三角,D为元素全为正数的对角阵。
(2)设A=LDLT
则有:
(3)为了避免重复计算,引进中间量tik=likdk,将上式改写为:
(4)Ax=bLDLTx=bLy=b,Dz=y,LTx=z
求解公式:
上述方法称为“改进的平方根法”
注:
不能选主元作行交换——破坏对称性,必是数值稳定的
a=[5-41;
-46-4;
1-46]
b=[2;-1;-1]
d=[000]';
t=[000;000;000];
L=[000;000;000];
n=3
fork=1:
n
d(k)=a(k,k)-t(k,1:
k-1)*L(k,1:
k-1)'
t(k+1:
n,k)=a(k+1:
n,k)-t(k+1:
n,1:
k-1)*L(k,1:
k-1)'
L(k+1:
n,k)=t(k+1:
n,k)/d(k)
end
x=[000]';
y=[000]';
z=[000]';
fori=1:
n
y(i)=b(i)-L(i,1:
i-1)*y(1:
i-1)
end
fori=1:
n
z(i)=y(i)/d(i)
end
fori=n:
-1:
1
x(i)=z(i)-L(i+1:
n,i)'*x(i+1:
n)
end
6.追赶法
在实际问题中,经常遇到以下形式的方程组
这种方程组的系数矩阵称为三对角矩阵。
以下针对这种方程组的特点提供一种简便有效的算法——追赶法。
设直接三角分解为:
A=LR
其中
,
容易看出pi=ci(1,2,…,n–1),计算L和R的其余元素的公式为
有了A的LR分解后,线性方程组Ax=d等价于两个简单的方程组:
Ly=d,Rx=y
于是,计算y的公式是
y1=d1,yi=di–liyi–1(i=2,3,…,n)
计算过程:
r1→y1→l2→r2→y2→...→ln→rn→yn(追)
计算x的公式是
xn=yn/rn,xi=(yi–cixi+1)/ri(i=n–1,n–2,…,1)(赶)
定理:
若A为对角占优/*diagonallydominant*/的三对角阵,且满足
|b1|>|c1|>0,|bn|>|an|>0,ai≠0,ci≠0,则追赶法可解以A为系数矩阵的方程组
注:
1.如果A是严格对角占优阵,则不要求三对角线上的所有元素非零。
2.根据不等式|
可知:
分解过程中,矩阵元素不会过分增大,算法保证稳定。
运算量为O(6n)。
7.QR算法
补充向量的范数与矩阵的范数
在线性代数方程组的数值解法中,经常需要分析解向量的误差,需要比较误差向量的“大小”或“长度”。
那么怎样定义向量的长度呢?
我们在初等教学里知道,定义向量的长度,实际上就是对每一个向量按一定的法则规定一个非负实数与之对应,这一思想推广到
维线性空间里,就是向量的范数或模。
用Rn表示n维实向量空间,用Cn表示n维复向量空间,首先将向量长度概念推广到Rn(或Cn)中。
1.向量的范数
范数的最简单的例子,是绝对值函数:
有三个熟知的性质:
(1)x0x>0x=0当且仅当x=0
(2)ax=axa为常数
(3)x+y≤x+y
范数的另一个简单例子是三维欧氏空间的长度
设x=(x1,x2,x3),则x的欧氏范数定义为:
欧氏范数也满足三个条件:
x,yR3,a为常数
(1)x≥0,且x=0x=0
(2)ax=ax
(3)x+y≤x+y
前两个条件显然,第三个条件在几何上解释为三角形一边的长度不大于其它两边长度之和。
因此,称为三角不等式。
向量范数的一般概念:
定义1:
设V是数域F上的向量空间,对V中任一向量α,都有唯一实数||α||与之对应,满足如下三个条件:
1)正定性:
α≥0,且α=0α=0
2)齐次性:
kα=|k|||α||,这里kF
3)三角不等式:
||α+||||α||+||||
则称α为α的范数。
Cn上的常见范数有:
1)1-范数
2)2-范数
3)-范数
不难验证,上述三种范数都满足定义的条件。
注:
上述形式的统一:
1p
定理5:
定义在Cn上的向量范数||x||是变量x分量的一致连续函数。
(f(x)=||x||)
定理6:
在Cn上定义的任何两个范数都是等价的。
即存在正数k1与k2(k1≥k2>0)对一切xCn,不等式
k1||x||b||x||ak2||x||b
成立。
对常用范数,容易验证下列不等式:
例1设A为正定矩阵,xCn,令
,
则||x||A为向量范数,称为加权范数
证明:
(1)显然xA≥0,且xA=0x=0
(2)
(3)因为A正定,所以可逆P,使A=PTP
所以||x+y||A=||P(x+y)||2=||Px+Py)||2||Px||+||Py||2=||x||A+||y||A
综上,||x+y||A为范数。
有了范数的概念,我们就可以讨论向量序列的收敛性问题。
定义2:
设给定Cn中的向量序列{xk},即
x0,x1,…,xk,…
其中
若对任何i(i=1,2,…,n)都有
则向量
称为向量序列{xk}的极限,或者说向量序列{xk}依坐标收敛于向量x*,记为
定理7:
向量序列{xk}依坐标收敛于x*的充分条件是
如果一个向量序列{xk}与向量x*满足上式,就说向量序列{xk}依范数收敛于x*,于是便得:
向量序列依范数收敛与依坐标收敛是等价的。
2.矩阵的范数
下面我们给出矩阵范数的一般定义。
为简单起见,仅讨论实数域的情形。
定义3若矩阵ARnn的某个非负的实值函数N(A)=||A||,满足条件
1)正定性:
A≥0,且A=0A=0
2)齐次性:
kA=|k|||A||kR
3)三角不等式:
||A+B||||A||+||B||
4)相容性:
||AB||||A||||B||
则称N(A)是Rnn上的一个矩阵范数(或模)。
由于在大多数与估计有关的问题中,矩阵和向量会同时参与讨论,所以希望引进一种矩阵的范数,它是和向量范数相联系而且和向量范数相容的,即
5)||Ax||||A||||x||
对任何向量xRn及矩阵ARnn都成立。
如由Rnn上2-范数可以得到Rnn中矩阵的一种范数
称为A的Frobenius范数。
||A||F显然满足正定性、齐次性及三角不等式。
为此我们再引进一种矩阵的范数。
定义4设xRn及矩阵ARnn,||x||v为一种向量范数(如v=1,2或∞),相应地定义一个矩阵的非负函数
可验证||A||v满足定义3,所以||A||v是Rnn上矩阵的一个范数,||A||v还满足相容性条件5),称为A的算子范数。
定理8:
设n阶方阵A=(aij)nn,则
与||x||1相容的矩阵范数是
与||x||2相容的矩阵范数是
其中1为矩阵ATA的最大特征值。
与||x||相容的矩阵范数是
注:
A的范数||A||与A的特征值间的关系
设为矩阵A的任一特征值,向量e为相应的特征向量,则e=Ae
因为||||e||=||Ae||||A||||e||
所以||||A||
从而得
定理9:
矩阵A的任一特征值的绝对值不超过A的范数||A||。
定义6:
矩阵A的诸特征值的最大绝对值称为A的谱半径,记为:
由定理9可知:
§2.3直接法的误差分析
1.误差向量与残差向量
设x为线性方程组Ax=b的准确解,用数值算法得到的近似解为x*,则称向量e=x*-x为误差向量。
显然,||e||可以表示近似解x*的准确程度:
||e||越小,近似解x*越准确。
但由于x不知道,所以无法计算e。
变通的办法是考虑残差向量:
r=b–Ax*的范数的大小。
但是,在某些情况下,尽管||r||很小,||e||也可能很大。
例线性方程组
的准确解是(1,1)T,若用某种方法得到计算解(2.000,0.500),则残差向量
其-范数||r||=0.0015,而误差向量
其-范数||e||=1,它是残差向量的667倍。
能否建立残差向量与误差向量之间的关系?
由于r=b–Ax*=Ax–Ax*=A(x–x*)=Ae,故有
e=A–1r,||e||||A–1||||r||
另一方面,由b=Ax得||b||||A||||x||,因此
其中,
是近似解x*的相对误差,而
表示相对残差
上式表明近似解x*的相对误差大约是相对残差的||A||||A–1||倍。
2.条件数与病态方程组
定义5:
设A为n阶非奇矩阵,称数||A||||A–1||为矩阵A的条件数,记为cond(A)。
显然,当A的条件数cond(A)较小时,只要残差向量很小,近似解x*的相对误差就很小,近似解的精度高,但若cond(A)较大时,则即使相对残差较小,近似解x*的相对误差也可能很大,近似解的精度就很低。
上例中线性方程组系数矩阵A的条件数:
||A||=3,||A-1||1000cond(A)3000很大!
这就解释了为什么残差向量小未必保证近似解的误差小。
定义6:
当cond(A)相对地大时,称方程组Ax=b为病态的,否则称为良态的。
3.病态方程组的判断
理论上虽然可以用矩阵的条件数来度量线性方程组求解问题的病态程度,但是,条件数大到怎样才算大,并没有客观的适用的标准,且当方阵的阶n较大时,A-1的计算也非易事,所以一般地说,遇到下列几种情况之一就应考虑线性方程组可能是病态的:
(1)系数矩阵的某两行(列)几乎近似相关;
(2)系数矩阵的行列式的值很小;
(3)用主元消去法解线性方程组时出现小主元;
(4)近似解x*已使残差向量r=b-Ax*的范数很小,但该近似解仍不符合问题的要求.
4.解的精度估计与改善
我们知道,影响计算精度的一个重要因素是舍人误差的传播.如前所述,按计算过程逐步地分析舍人误差的传播是极为困难又不实用的.一种常用的办法是向后误差分析方法,它将计算过程中舍入误差的影响归结为原始数据的小扰动,也就是说,将计算所得的近似解x*看成是某个带有小扰动的方程组(A+A)x=b+b的准确解.
定理:
设x*是线性方程组(A+A)x=b+b的解,x是线性方程组Ax=b的解,且A的扰动A满足||A-1||||A||<1,则有
上式表明了近似解x*的相对误差与A及b的相对误差之间的关系,其中A的条件数起着重要的作用.
估计式只是理论性的,实际应用有很大的困难.因此,—般采用下述方法估计近似解的精确程度.
对于线性方程组Ax=b,任取一个非零向量y,计算出d=Ay;然后采用与解Ax=b同一算法解线性方程组Ay=d,得到它的近似解y*,y*的相对误差为
,由于求解Ay=d与求解Ax=b的算法相同,且系数矩阵同为A,所以可以认为Ax=b的近似解x*的相对误差大约也是
。
当然,对病态方程组来说,这种估计方法是不大可靠的.
病态方程组的求解问题是计算方法研究的一个重要课题.这里介绍一下改善近似解的精度的方法:
1)采用双倍字长进行运算,特别是在做形如∑aibi的计算时,每一乘积项均取双倍字长的数,并与其他项累加后再舍入成单字长,这样能使结果的精度提高,同时也比全都用双倍字长进行运算节省内存和机时.
2)采用迭代改善的办法.设已得到Ax=b的近似解x(0),采用双倍字长计算Ax(0)后再舍入成单字长,并计算残差向量r(0)=b–Ax(0),记x=x–x(0),则x应满足方程组Ax=r(0).采用原单字长运算的算法求出x的近似解x(0)(在A=LR或A=QR分解已完成下,只需作顺代和回代求解),从而得到Ax=b的改善了的近似解x
(1)=x(0)+x(0),若||x(0)||<,则x
(1)就是满足精度()要求的解。
否则再重复上述做法。
§2.4线性方程组迭代解法
1.基本思想
若方程组Ax=b可写成等价形式
x=Bx+g,(2.4-1)
则给定一个初始向量x(0),可以得到迭代公式
x(k+1)=Bx(k)+g,k=0,1,…(2.4-2)
若上式确定的向量序列{x(k)}收敛于x,则x显然是(2.4-1)的解,从而为方程组Ax=b的解。
形如(2.4-2)的逐次逼近的方法称为简单迭代法,B称为该迭代法的迭代矩阵。
2.迭代法的收敛性
定理1:
对任意初始向量x(0)及常向量g,迭代法(2.4-2)收敛的充分必要条件是迭代矩阵B的谱半径(B)<1。
这一结论在理论上是颇为重要的,但实际用起来不甚方便。
当n较大时,谱半径的计算并非易事,但由于(B)||B||,所以只要B的某种范数小于1,迭代法(2.4-2)必收敛:
定理2:
若迭代矩阵B的某种范数||B||<1则(2.4-2)确定的迭代法对任意初值x(0)均收敛于方程组x=Bx+g的唯一解x。
定理3:
若||B||=q
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 线性方程组 数值 解法