计计算算方方法法第七章第七章 求解方程组的数值解法求解方程组的数值解法计计算算方方法法计计算算方方法法对线性方程组：或者：我们有Cramer法则：当且仅当时，有唯一的解，而且解为：计计算算方方法法 运算量运算量 （Amount of Computation）用克莱姆用克莱姆（CramerCramer）法则求解法则求解n n阶线性方程组阶线性方程组 每个行列式由每个行列式由n!n!项相加，而每项包含了项相加，而每项包含了n n个因子个因子相乘，乘法运算次数为相乘，乘法运算次数为( (n-1)n !n-1)n !次次. .仅考虑乘仅考虑乘( (除除) )法运算法运算, ,计算解向量包括计算计算解向量包括计算n+1n+1个行列式和个行列式和n n次除法运算次除法运算, ,乘乘( (除除) )法运算次数法运算次数N=(n+1)(n-1)n!+n. .20阶的方程组用克莱姆法则进行运算阶的方程组用克莱姆法则进行运算,如用每如用每秒秒1亿次乘法运算的计算机要亿次乘法运算的计算机要30万年。万年。 计计算算方方法法怎么办呢?解线性方程组的两类方法:直接法: 经过有限次运算后可求得方程组精确解的方法(不计舍入误差!)迭代法：从解的某个近似值出发，通过构造一个无穷序列去逼近精确解的方法(一般有限步内得不到精确解)计计算算方方法法7.1 解线性方程组的直接法解线性方程组的直接法 高斯消元法高斯约当法计计算算方方法法7.1 解线性方程组的直接法解线性方程组的直接法 7.1.1 回代法回代法(Backsubstitution )计计算算方方法法且aii0,i=1,2,n计计算算方方法法计计算算方方法法计计算算方方法法7.1 解线性方程组的直接法解线性方程组的直接法 一、高斯消去法一、高斯消去法(Gaussian Elimination )思思路路首先将方程组首先将方程组Ax=b 化为上三角方程组化为上三角方程组，此过程称为此过程称为消去过程消去过程，再求解上三角方程，再求解上三角方程组，此过程称为组，此过程称为回代过程回代过程. .7.1.2 高斯消去法和选主元高斯消去法高斯消去法和选主元高斯消去法计计算算方方法法计计算算方方法法其中其中将增广矩阵将增广矩阵的的第第 i 行行 + li1 第第1 1行行，得到：，得到：消去过程：消去过程：第一步第一步：设设 ，计算因子，计算因子计计算算方方法法第第k步：步：设设 ，计算因子，计算因子且计算且计算计计算算方方法法第第k步：步：设设 ，计算因子，计算因子且计算且计算共进行共进行 n 1步，得到步，得到计计算算方方法法-=-=-=+9262214282332321xxxxxx计计算算方方法法消元过程中出现了除法很小时计计算算方方法法二、二、 选主元消去法选主元消去法为避免这种情况的发生, 可通过交换方程的次序，选取绝对值大的元素作主元. 基于这种思想导出了主元素法在高斯消去法消去过程中可能出现 的情况，这时高斯消去法将无法进行；即使主因素 但很小，其作除数 ，也会导致其它元素数量级的严重增长和舍误差的扩散计计算算方方法法v 列主元消去法列主元消去法在第在第k k步消元前，在系数矩阵第步消元前，在系数矩阵第k k列的对角线以下列的对角线以下的元素中找出绝对值最大的元。的元素中找出绝对值最大的元。 若若pkpk，交换第交换第k k个与第个与第p p个方程后，再继续消去计算个方程后，再继续消去计算. . 这种方法称为这种方法称为列主元列主元GaussGauss消去法。消去法。 列主元列主元GaussGauss消去法保证了消去法保证了L Likik1 1 (i=k+1,k+2,(i=k+1,k+2,，n).n).计计算算方方法法计计算算方方法法例例 用列主元消去法解方程组用列主元消去法解方程组解解 第一次消元对第一次消元对 因列主元素为因列主元素为a a3131(1)(1), ,故先作行交换故先作行交换E E1 1 E E3 3, ,然后进行然后进行消元计算可得消元计算可得-0.002x1+2x2+2x3 =0.4 x1+0.78125x2 =1.38163.996x1+5.5625x2+4x3=7.4178 -0.002 2 2 0.4 A(1) |b(1) = 1 0.78125 0 1.3816 3.996 5.5625 4 7.4178 3.996 5.5625 4 7.4178 A(2) |b(2) = 0 -0.61077 -1.0010 -0.47471 0 2.0029 2.0020 0.40371计计算算方方法法 由此回代由此回代, ,得得x x=(1.9272,-0.69841,0.90038)=(1.9272,-0.69841,0.90038)T T与精确解与精确解 x x=(1.9273,-0.69850,0.90042)=(1.9273,-0.69850,0.90042)T T相相比较是比较准确的比较是比较准确的. . 3.996 5.5625 4 7.4178A(3) |b(3) = 0 2.0029 2.0020 0.40371 0 0 -0.39050 -0.35160 第二次消元对第二次消元对 A A(2)(2) | |b b(2)(2) , ,因列主元素为因列主元素为a a3232(2)(2) , ,故先故先作行交换作行交换E E2 2 E E3 3, ,然后进行消元计算可得然后进行消元计算可得计计算算方方法法 if nargini t=a(i,:);a(i,:)=a(r,:);a(r,:)=t; end %消元计计算算方方法法 a(i+1):n,(i+1):(n+1)= a(i+1):n,(i+1):(n+1)- a(i+1):n,i)/a(i,i)* a(i,(i+1):(n+1); a(i+1):n,i)=zeros(n-i,1); if flag=0,a,end end %回代 x=zeros(n,1); x(n)=a(n,n+1)/a(n,n); for i=n-1:-1:1 x(i)= (a(i,n+1)-a(i,(i+1):n)*x(i+1):n)/a(i,i); end计计算算方方法法v 全主元消去法全主元消去法在第在第k k步消去前，步消去前， 在系数矩阵右下角的在系数矩阵右下角的n-k+1阶主子阵中，阶主子阵中，选绝对值最大的元素作为主元素选绝对值最大的元素作为主元素。(1) If p k then 交换第交换第 k 行与第行与第p行行; If q k then 交换第交换第 k 列与第列与第 q 列列;(2) 消元消元注注注注：列交换改变了列交换改变了 x xi i 的顺序，须记录的顺序，须记录交换次序交换次序，解完后再换回来。解完后再换回来。计计算算方方法法 例： 用全主元素消去法求解方程组 解：这里主元为-18,交换方程与方程得 计计算算方方法法 再全选主元,主元为2.333,交换x2和x3所在的两列,同时改变两未知数的排列号得 -0.944/2.333得 计计算算方方法法已经化为三角方程组,回代求解 x1=1.000,x2=3.000,x3=2.000这里未知数x2与x3已对调,所以应恢复解的顺序,方程组的实际精确解为 x1=1.000,x2=2.000,x3=3.000 计计算算方方法法计计算算方方法法7.1.37.1.3 Gauss-Jordan消元法消元法计计算算方方法法计计算算方方法法计计算算方方法法算法:将在Gauss消元第k步，变为计计算算方方法法约当法求逆阵计计算算方方法法计计算算方方法法计计算算方方法法迭代终止条件？7.2向量范数计计算算方方法法设n维向量xRn,记对应向量x的一个实数为x,若x满足下面三个性质: (1)非负性,即x0,当且仅当x=0时,x=0 (2)齐次性,x=x,为任意实数 (3)三角不等式性,x+yx+y,yRn。 则称该实数x为向量向量x的范数的范数。计计算算方方法法 具有向量范数的三条基本性质,称之为向量x中的p范数,记为xp在Rn中,常用的几种向量范数有: 计计算算方方法法距离量度欧几里德距离（欧几里德距离（Eulidean）:街区距离（街区距离（city-block）:棋盘距离（棋盘距离（chess board）:计计算算方方法法7.3解线性方程组的迭代法 直接法得到的解是理论上准确的，但是我们可以看得出，它们的计算量都是n3数量级，存储量为n2量级，这在n比较小的时候还比较合适（n400），但是对于现在的很多实际问题，往往要我们求解很大的n的矩阵，而且这些矩阵往往是系数矩阵就是这些矩阵含有大量的0元素。对于这类的矩阵，在用直接法时就会耗费大量的时间和存储单元。因此我们有必要引入一类新的方法：迭代法。 迭代法具有的特点是速度快。计计算算方方法法格式很简单：1 雅克比（雅克比（Jacobi）法）法计计算算方方法法 迭代法解方程组 雅克比迭代公式计计算算方方法法2 GaussSeidel迭代迭代在Jacobi迭代中，使用最新计算出的分量值计计算算方方法法 迭代法解方程组 赛德尔迭代公式计计算算方方法法Jacobi迭代算法function X=jacobi(A,B,P,delta, max1) % Input - A is an N x N nonsingular matrix% - B is an N x 1 matrix% - P is an N x 1 matrix; the initial guess% - delta is the tolerance for P% - max1 is the maximum number of iterations% Output - X is an N x 1 matrix: the jacobi approximation to% the solution of AX = B N = length(B); for k=1:max1 for i=1:N X(i)=(B(i)-A(i,1:i-1,i+1:N)*P(1:i-1,i+1:N)/A(i,i); end err=abs(norm(X-P); relerr=err/(norm(X)+eps); P=X; if (errdelta)|(relerrdelta) break endend X=X;计计算算方方法法function X=gseid(A,B,P,delta, max1)N = length(B);for k=1:max1 for j=1:N if j=1 X(1)=(B(1)-A(1,2:N)*P(2:N)/A(1,1); elseif j=N X(N)=(B(N)-A(N,1:N-1)*(X(1:N-1)/A(N,N); else %X contains the kth approximations and P the (k-1)st X(j)=(B(j)-A(j,1:j-1)*X(1:j-1)-A(j,j+1:N)*P(j+1:N)/A(j,j); end end err=abs(norm(X-P); relerr=err/(norm(X)+eps); P=X; if (errdelta)|(relerrdelta) break endendX=X;高斯赛德尔迭代算法计计算算方方法法作业：H61，H62，H67，H65改为向量x的无穷范数、2范数和1范数