线性方程组的求解
目录
摘要 ................................................................................................ 1 一.用列主元消去法解方程组 ................................................................... 2 1.问题的提出 ......................................................................................................................... 2 2.问题的分析 ......................................................................................................................... 2 3.问题的解决 ......................................................................................................................... 3 二.编写一个列主元消去法求逆矩阵的程序 ................................................... 4 1.问题的提出 ......................................................................................................................... 4 2.问题的分析 ......................................................................................................................... 4 3.问题的解决 ......................................................................................................................... 5 三.用LU分解法解方程组Axb ............................................................. 5 1.问题的提出 ......................................................................................................................... 5 2.问题的分析 ......................................................................................................................... 5 3.问题的解决 ......................................................................................................................... 6 四.用改进平方根法解方程组 ................................................................... 7 1.问题的提出 ......................................................................................................................... 7 2.问题的分析 ......................................................................................................................... 7 3.问题的解决 ......................................................................................................................... 8 五.用追赶法解方程组 ........................................................................... 9 1.问题的提出 ......................................................................................................................... 9 2.问题的分析 ......................................................................................................................... 9 3.问题的解决 ....................................................................................................................... 10 六.分别用雅可比迭代法与高斯-赛德尔迭代法解方程组 ................................... 11 1.问题的提出 ....................................................................................................................... 11 2.问题的分析 ....................................................................................................................... 11 3.问题的解决 ....................................................................................................................... 12 参考文献 ......................................................................................... 14 个人体会 ......................................................................................... 15 附录:程序代码 ................................................................................. 16
摘要
在科技研究和工程技术所提出的计算问题中,经常会遇到线性方程组的求解问题,这里主要是有关线性方程组的直接解法。解线性方程组的直接法是用有限次运算求出线性方程组 Ax=b 的解的方法。线性方程组的直接法主要有Gauss消元法及其变形、LU(如Doolittle、Crout方法等)分解法和一些求解特殊线性方程组的方法(如追赶法、LDLT法等)。这里主要有列主元消元法,LU分解法,改进的平方根法,追赶法和雅可比迭代,高斯—塞德尔迭代 的构造过程及相应的程序。线性方程的解法在数值计算中占有极重要的地位,因此,线性方程组的求解是数值分析课程中最基本的内容之一。
关键词:列主元消元法;LU分解 ;改进平方根法;追赶法;雅可比迭代;高斯—塞德尔迭代
一.用列主元消去法解方程组:
1.问题的提出:
x1x22x3x48x1x23x44
2x2x3x3x202xxxx1
11234234
(1) (2)
3x1x2x32x43x1x2x32x12x23x3x44x1x24x33x44
2.问题的分析:
列主元消去算法主要分为两个过程:消去过程和回代过程
1. 消去过程
对k1,2,,n1
(1)选主元 找ikk,,n,使a
(k)ikk
(k)
aik kin
(k)
(2)若aikk0则停止计算(detA=0)
(3)若ikk 则换行EkEik (4)消元
对i=k1,,n1
lik
(k)
aik(k)akk
对jk1,,n1 aij
(k1)
(k)(k)
aijlikakj
2.回代过程
(1)若ann0则停止计算(detA=0) (2) xn
(n)
an,n1(n)ann
(n)
(3)对in1,,1
)ai(,nn1
xi
ji1
(n)aii
aij(n)xj
n
3.问题的解决:
11034(1)解:对于A(1)|b
=A|b=21111
31123
1234
第1步选列主元为a(1)
313,i13,作变换E1E3,然后计算
l2120.667, l3110.333,l411
0.333
再作变换
E2l21E1E2,E3l31E1E3,E4l41E1E4,
31123A
(2)
|b(2)
=01.6670.3330.333301.3330.3332.3335 01.6672.6670.3333
第2步,对A(2)
选列主元为a(2)
22
31.667,i22,计算
l4
32
5
0.8, l421, 再做变换(E3l32E2)E3,(E4l42E2)E4,得到
31123A
(3)
b(3)
05/31/31/33009/1539/1513/5 00300
消去过程结束,回代计算得到解
x11x22x30
x41
所以原方程组的解为X(1,2,0,1)T
。
11
2
18
(2)解:对于
A
(1)
|b
A|b
22
3320
11102
11434
得到
第1步选列主元为a21
(1)
, 2,i12,作变换E1E2,然后计算l21l31l411再作变换E2l21E1E2,E3l31E1E3,E4l41E1E4,得到
A
(2)
|b(2)
223320
1/21/2200
021/23/28005/29/214
(2)
第2步,对A
(2) 选列主元为a322,i23,作变换(E2)(E3)得到
A
(3)
b(3)
223320021/23/28 001/21/22005/29/214
(3)
第3步,对A
(3)
选列主元为a435/2,i34,作换行(E3)(E4),计算l431/5,作
变换(E4l43E3)E4,则(A
(4)
223320
8021/23/2
|b(4))
005/29/2140002/54/5
消去过程结束,回代计算得到x17所以原方程组的解为X(7,3,2,2)
T
x23x32x42
二.编写一个列主元消去法求逆矩阵的程序,并用以求如下矩阵的逆矩阵。
1.问题的提出:
152
A283
136
2.问题的分析:
设A是n阶非奇异矩阵,则A1 存在,令 A1X(X1,X2,,Xn)
则
AXE
其中, E ( 1 , E 2 , , E n )是单位矩阵 E
从而 ( AX 1 , AX 2 , , AX n ) ( E 1 , E 2, n ) ,E即 AX k k , k 1 , 2 , , n E
求解上述n个方程组时,将它们的增广矩阵 ( A k )合并后有 ,E
a11a1ra1n100
(A,E)ar1arrarn010
a aa001nrnnn1
按列选主元无回代高斯消元法得 100b11b1r
(E,B)(E1,E2,,En,B1,B2,,Bn)010br1brr
001bb
n1nr
显然列向量 B k E k , k 1 ,2 , , n 的解。 AXk是方程 因此 A1B
b1n
brn
bnn
3.问题的解决:
1521
解:由上述过程可得283B1=0,
1360
1520
283B30,解得BB11361
1520
283B21, 1360
0.2030.0320.05
B30.0320.110.06
0.050.060.128
B2
0.2030.0320.05
1
0.110.06 所以A0.032
0.050.060.128
三.用LU分解法解方程组Axb
1.问题的提出:
48240124242412124
A b2 06202
626216
2.问题的分析:
求解线性方程组的LU分解法:将系数矩阵A分解成两个三角矩阵的乘积,
即: A=LU的形式
其中,L为下三解矩阵,U为上三解矩阵,则线性方程组:Ax=b 可改写为 LUx=b 令 Ux=y 得 Ly=b
然后,用前代方法求Ly=b得列矩阵y,再用回代方法求Ux=b,得到的列矩阵x即为所求的线性方程组的解.
3.问题的解决:
1l21
解:设L
l31l41u11ulLU1121
ul1131ul1141
1l32l42
1l43u12
u12l21u22u12l31u22l32u12l41u22l42
u11
U
1
u13
u13l21u23
u13l31u23l32u33u13l41u23l42u33l43
u12u22
u13u23u33
u14u24
得到 u34u44
u14
u14l21u24
u14l31u24l32u34
u14l41u24l42u34l43u44
4824012
24241212
06202
66216
解得:
010.5100L00.5100.1250.250.0711
subma(tMri0x3811)
0048240
012126U
0014100012.929
12
令Uxy,AxLUxb,Lyb,
1y144
y0.51462
= 解得 y0y0.5125
30.1250.250.0711y23.357412x14482400.521
x1212661.0062
解得x-0.376 141x35
12.929x43.357-0.26
所以原方程组的解为X0.5211.0060.3760.26。
T
四.用改进平方根法解方程组:
1.问题的提出:
用改进平方根法解下列方程组: (1)Axb
0.500000.51
0.51.50.50.250.2500000.51.50.250.250
b A
0.250.251.50.5000
000.250.250.51.50.5000000.50.5
(2)Axb
0.352.002.000.350.351.002.001.352.001.352.001.000.350.352.002.000.352.00
2.0012.002.001.350.351.002.000.350.35
2.002.000.350.352.001.000.351.3512.00
Ab2.000.350.351.002.001.350.352.002.00
0.350.352.002.000.351.352.001.002.00
2.002.001.350.351.002.000.350.352.00
2.002.002.000.350.352.001.000.351.35
2.问题的分析:
平方根法主要用于求解对称正定矩阵方程:首先要提到的是关于正定矩阵的定理,说的是若A为n阶地对称正定矩阵,则存在一个实的非奇异下三角矩阵L,使A=LL’(L’表示L的对称矩阵)。根据这个前提,在结合前面的LU分解算法,便有了平方根算法。平方根算法的计算量约为普通三角分解算法的一半,但由于这里要用到开平方,效率不是很高,所以,便有了为效率而存在的改进的平方根算法。 对称矩阵分解法(改进平方根法)如下:
设对称矩阵,存在分解式ALDU,因A对称矩阵,即 AAT得
AATLDUUTDLT,由分解的唯一性得ULT,LUT,ALDLT
T
于是可得对称矩阵A的LR分解式
因此完成三角分解ALR,只需计算出R的元素,然后由
由紧凑格式计算公式得
,求出L的元素。
(4.21)
(4.22)
(4.23)
公式(4.21)(4.22)(4.23)为改进平方根法的计算公式。
3.问题的解决:
(1)解:由(4.21)可计算出
00000
由(4.22)可计算出 由(4.23)可计算出
10
11
00.5L
00.2500
00.25
0010.10
00010
00001
0.10.2982.5961
110.5y
0.200.14
2.311.310.47x
0.1580.3640.14
T
所以原方程组的解为X2.31,1.31,0.47,0.158,0.364,0.14 (2)解:由(4.21)可计算出
00000001
0.51000000
1000000.1750
0.1750.3330.96910000L
100.4850.3271000
10.4240.9380.4161001
0.17500.4550.9810.1230.316100.6750.3330.8270.240.7220.6070.0291
由(4.22)可计算出 y2312.351.3165.5614.6958.8531.763
T
由(4.23)可计算出 x0.75所以原方程组的解为X0.75
-0.5-0.5
3.253.25
-0.5-0.5
0.750.75
-0.5-0.5
3.253.25
T
-0.5 T
-0.5
五.用追赶法解方程组:
1.问题的提出:
用追赶法解方程组:
2x1x25
x2xx12123
(1)
x22x3x411x32x41(2)Axb
04127
1411515141
b A1010
0151411415
2.问题的分析:
矩阵Doolittle分解形式b1a2
追赶法构造过程:追赶法仍然保持LU分解特性,它是一种特殊的LU分解。追赶法充分
c1b2
c2an1bn1
an
1p2cn1
bn
1
p3
1pn
q11q22
n1
1qn
利用了系数矩阵的三对角特点,而且使之分解更简单,得到对三对角线性方程组的快速解法。
计算公式为:
y1d1
由矩阵乘法及相等定义,有:ykdkpkyk1k2,3,,n
qnyb1x1nqnpqkca,q)qb12kck1(k2,3,,n)x(ynn,11k,,1,1kkkkkxk1kpkkk于是得计算L的元素pi及U的qi和i的计算公式,为:综合以上,求解出Doolittle分解的计算公式:
可由Lyd及Uxy解出。
追赶法较简单,计算量、乘除法次数仅有5n4。追赶法的特殊求解过程,节省了计算时间和存贮空间,但是因为追赶法来源于Gauss消元法,因此也存在Gauss消元法的缺点,即当qk0时不能进行。
qqbby1d11111,pqkpaak1kkk1
bkcpk2,3,,nq1k2,3,,n)kk1k1kck(ydykqkkbkkpk11ynqn,xk(ykckxk1)qk,kn1,,1nx若记d(dd2,,dn)T表示Axd,当ALU时,用这组公式解线性方程组的方法亦称为追赶法。
3.问题的解决:
21005121012
(1)解:A 由LU=A得 b012111
00121
1pLU2
1p3
q11
r1q2
r2
q3
0210
1210A r30121q40012
1p4
1
10.5
L0.6671
0.751
由Ly=b得y(5
21
1.51
U
1.3331
1.25
T
9.54.6672.5)T;由Ux=y得x1352
所以原方程组的解为X135(2)解:
由LU=A得
2
T
1
00.251L
0.2671
0.2681
00.2681
4
103.751
U
3.7331
3.732
103.732
由Lyb得
y2721.7520.820.57120.51220.49620.49220.49120.49由Ux=y得
x8.7067.8237.5867.5227.5037.4917.4627.3566.9625.49T
所以原方程组的解为
X8.7067.8237.5867.5227.5037.4917.4627.3566.9625.49T
。
六.分别用雅可比迭代法与高斯-赛德尔迭代法解方程组:
1.问题的提出:
2x1x210x311
3x2x38x4
1110x1x要求精度为0.5e-4 22x36x111x2x33x425
2.问题的分析:
(1) .雅可比迭代法:
设线性方程组Axb (1) 的系数矩阵A可逆且主对角元素a11,a22,,ann均不为零,令 Ddiag(a11,a22,,ann)
并将A分解成
AADD (2) 从而(1)可写成
DxDAxb 令xB1xf1
20.49
T
其中B1ID1A,f1D1b. (3)
以B1为迭代矩阵的迭代法(公式)
x(k1)B1x(k)f1 (4)
称为雅可比(Jacobi)迭代法(公式),用向量的分量来表示,(4)为
n(k1)(k)aiibiaijxjxi
j1 (5)
i1,2,,n,k0,1,2,
其中x
(0)
(0)(0)
x1(0),x2,,xn
T
为初始向量.
(k1)
由此看出,雅可比迭代法公式简单,每迭代一次只需计算一次矩阵和向量的乘法.在电算时需要两组存储单元,以存放x及x
(2). 高斯—塞德尔迭代法: 的所有分量,显然在计算第i个分量xi的第k1次近似x的分量(Gauss-Seidel)迭代法.
把矩阵A分解成
k1
(k)
.
(k)
由雅可比迭代公式可知,在迭代的每一步计算过程中是用x
k1
的全部分量来计算x
k1
(k1)
没有
被利用,从直观上看,最新计算出的分量可能比旧的分量要好些.因此,对这些最新计算出来
时,已经计算出的最新分量x1
,...,xi1
k1
xj
k1
加以利用,就得到所谓解方程组的高斯—塞德
ADLU (6)
其中Ddiaga11,a22,...,ann,L,U分别为A的主对角元除外的下三角和上三
角部分,于是,方程组(1)便可以写成
即 其中
DLxUxb
xB2xf2
B2DLU,
1
f2DLb (7)
1
以B2为迭代矩阵构成的迭代法(公式)
xk1B2xkf2 (8)
称为高斯—塞德尔迭代法(公式),用 量表示的形式为
(9)
由此看出,高斯—塞德尔迭代法的一个明显的优点是,在电算时,只需一组存储单元(计算出xi
k1
x
(k1)i
i1n1(k1)
biaijxjaijx(jk)
j1ji1aii
i1,2,n,k0,1,2,...
k
后xi
k
不再使用,所以用xi
k1
冲掉xi
,以便存放近似解.
3.问题的解决:
(1).用雅可比迭代法求解下列方程组
解:将方程组按雅可比方法写成
x10.5x25x35.5
x0.33x2.67x3.67234
x35x10.5x23x40.33x13.67x20.33x38.33
取初始值x
(0)
T(0)(0)(0)
x1(0),x2,x3,x4=0,0,0,0按迭代公式
T
x1(k1)0.5x2(k)5x3(k)5.5
(k1)
(k)(k)
0.33x32.67x43.67x2
(k1)(k)(k)
5x10.5x23x3
(k1)(k)(k)(k)
0.33x13.67x20.33x38.33x4
解:取初始值x
(0)
(0)(0)(0)
x1(0),x2,x3,
x4,按迭代公式
T
x1(k1)0.5x2(k)5x3(k)5.5(k1)
(k)(k)
0.33x32.67x43.67x2
(k1)(k1)(k1)
5x10.5x23x3
(k1)(k1)(k1)(k1)
0.33x13.67x20.33x38.33x4
数少),但这个结论,在一定条件下才是对的,甚至有的方程组,雅可比方法收敛,而高斯—塞德尔迭代法却是发散的.
参考文献:
[1] 杨大地,涂光裕.数值分析.重庆大学出版社.1998年1月 [2] 韩国强.数值分析.华南理工大学出版社.2005年3月 [3] 李有法.数值计算方法.高等教育出版社.1996年12月
[4] 关治,陆金甫.数值分析基础.高等教育出版社.1998年5月
个人体会:
通过这次数值计算综合训练,我从中学到了很多的东西,最重要的一点就是,我对数值计算的思想有了更深层次的认识,同时,也对C及MATLAB在编程方面更熟悉了一些。我做的主要是有关线性方程组求解问题,解线性方程组的直接法是用有限次运算求出线性方程组 Ax=b 的解的方法,主要有Gauss消元法及其变形、LU(如Doolittle、Crout方法等)分解法和一些求解特殊线性方程组的方法(如追赶法、LDLT法等)。我感觉到数值计算是一门应用性很强的基础课程。同时,在编程的过程中我遇到了很多的的问题,经过这次的训练,我解决了一些主要的问题,还有一些问题急待解决。我想在以后的学习中,我一定会将它们很好的解决的。
附录:程序代码:
1.程序代码:
function x=gauss2(A,b) n=length(b);A=[A,b]; for k=1:(n-1)
[Ap,p]=max(abs(A(k:n,k)));p=p+k-1; if p>k,
t=A(k,:);A(k,:)=A(p,:);A(p,:)=t; end
A((k+1):n,(k+1):(n+1))=A((k+1):n,(k+1):(n+1))... -A((k+1):n,k)/A(k,k)*A(k,(k+1):(n+1)); A((k+1):n,k)=zeros(n-k,1); end
x=zeros(n,1);
x(n)=A(n,n+1)/A(n,n); for k=n-1:-1:1
x(k,:)=(A(k,n+1)-A(k,(k+1):n)*x((k+1):n))/A(k,k); end
2.程序代码: # include "stdio.h" # define M 3
void main ( ) {
float MAT[M][2*M]; float MAT1[M][M]; float t; int i,j,k,l;
/***********************************************/ /*对矩阵进行初始化*/ for(i=0;i
for(j=0;j
/*对MAT1矩阵赋初值 */
for(i=0;i
scanf("%f",&MAT1[i][j]); /*打印目标矩阵?*/ printf("原矩阵为:\n"); for (i=0;i
for (j=0;j
/********************************************/
/*对MAT1矩阵进行扩展,MAT1矩阵添加单位阵,由M*M变成2M*2M矩阵 for(i=0;i
if (j
/*对M矩阵进行变换,使得前半部分矩阵成为单位阵,则 */ /*后半部分矩阵即为所求矩阵逆阵 */ for(i=0;i
/*对第i行进行归一化 */ for (j=0;j
MAT[i][j]=MAT[i][j]+MAT[k][j]; t=MAT[i][i];
for(j=i;j
/*对矩阵进行行变换,使得第i 列只有一个元素不为零,且为1*/ for(k=0;k
t=MAT[i][k];
for (l=i;l
MAT[k][l]=MAT[k][l]-MAT[i][l]*t; } }
/*将后半部分矩阵即所求矩阵逆阵存入MAT2矩阵。*/ for(i=0;i
for(j=0;j
MAT1[i][j]=MAT[i][j+M];
*/
printf("\n"); }
/*********************************************/ /*输出所求的逆阵*/ printf("逆阵为:\n"); for(i=0;i
for(j=0;j
printf("%5.2f",MAT1[i][j]); printf("\n"); } }
3.程序代码: #include #include #include
int LUTriangle(double** a, double* b, int n); void Solve(double** a, double* b, int n);
void main() {
double **a, *b; int i, n;
n=3;
a = malloc(n*sizeof(double *)); for (i=0;i
a[i]=malloc(n*sizeof(double)); b = malloc(n*sizeof(double));
a[0][0]=2; a[0][1]=0; a[0][2]=0; a[1][0]=3; a[1][1]=1; a[1][2]=0; a[2][0]=0; a[2][1]=4; a[2][2]=1; b[0]=2; b[1]=4; b[2]=5;
if (LUTriangle(a, b, n)==1) { Solve(a, b, n); for (i=0; i
{ printf("%f ",b[i]); }
printf("\n"); } else
{ printf("This is a singular matrix!\n"); }
for (i=0;i
getchar(); }
int LUTriangle(double** a, double* b, int n)
{ //LU Factorisation, coded by www.dayi.net btef (please let this line remain) int maxI, ii, i, j;
double maxV, d1, *temp;
//Find the row with the max element in column ii for (ii=0; ii
{ maxV = fabs(a[ii][ii]); maxI = ii;
for (i=ii+1; i
{ if (fabs(a[i][ii])>maxV) { maxV = fabs(a[i][ii]); maxI = i; } }
//If the max element = 0, this matrix is singular, return if (maxV==0.0) { return(0); }
//Possibly, exchange rows so that row ii got max element if (ii!=maxI)
{ temp = a[maxI]; a[maxI] = a[ii]; a[ii] = temp; d1 = b[ii];
b[ii] = b[maxI]; b[maxI] = d1; }
//Gauss elemination for (i=ii+1; i
a[i][ii] = d1; //update element in Lower triangle for (j=ii+1; j
}
}
return(1);
}
void Solve(double** a, double* b, int n)
{
int i, j;
double d1;
//Solve Lower triangle
for (i=1; i
{ d1 = 0.0;
for (j=0; j
{ d1 += a[i][j]*b[j];
}
b[i] -= d1;
}
//Solve Upper triangle
b[n-1] /= a[n-1][n-1];
for (i=n-2; i>=0; i--)
{ d1 = 0.0;
for (j=i+1; j
{ d1 += a[i][j]*b[j];
}
b[i] = (b[i]-d1)/a[i][i];
}
}
5.程序代码:
Clear[a,b,c,a,x];
n=Input["线性方程组阶数n="];
a=Input["三对角系数向量a="];
b=Input["三对角系数向量b="];
c=Input["三对角系数向量c="];
d=Input["三对角常数项向量d="];
eps1=0.00001;
t=1;
Do[If[Abs[b[[k-1]]]
a[[k]]=a[[k]]/b[[k-1]];
b[[k]]=b[[k]]-a[[k]]*c[[k-1]];
d[[k]]=d[[k]]-a[[k]]*d[[k-1]],
{k,2,n}];
x=Table[0,{n}];
If[t==1,
x[[n]]=d[[n]]/b[[n]];
Do[x[[k]]=(d[[k]]-c[[k]]*x[[k+1]])/b[[k]], {k,n-1,1,-1}] ];
Print["Ax=d的解为 " , x]
6.程序代码:
jacobi程序:
function x=jacobi(A,b,x0,e,N)
D=diag(diag(A));G=-inv(D)*(A-D);d=inv(D)*b;
k=0;x=x0;x0=x+2*e;
while norm(x0-x,inf)>e&k
k=k+1;
x0=x;x=G*x0+d;
end
if k==N, warning('already reach maximum number of iteration.');end
gauss-seidel程序:
function x=gauss-seidel(A,b,x0,e,N)
AL=tril(A);G=-inv(AL)*(A-AL);d=inv(AL)*b;
k=0;x=x0;x0=x+2*e;
while norm(x0-x,inf)>e&k
k=k+1;
x0=x;x=G*x0+d;
end
if k==N, warning('already reach maximum number of iteration.');end