10实验5线性方程组.ppt

《10实验5线性方程组.ppt》由会员分享，可在线阅读，更多相关《10实验5线性方程组.ppt（37页珍藏版）》请在得力文库 - 分享文档赚钱的网站上搜索。

1、10实验实验5线性方程组线性方程组1 1 求解线性方程组的直接法与迭代法求解线性方程组的直接法与迭代法求解线性方程组nnnnnnnnnnbxaxaxabxaxaxabxaxaxa22112222212111212111有两种方法：直接法：直接法：经过有限次算术运算求出精确解，高斯消元法高斯消元法 迭代法：迭代法：从初始解出发，根据设计好的步骤用逐次求出的近似解逼近精确解，雅可比迭代雅可比迭代和高斯高斯-赛德尔迭代赛德尔迭代法。(1)1.1 直接法直接法（高斯消元法与LU分解）设方程组为 43214321)1(44)1(43)1(42)1(41)1(34)1(33)1(32)1(31)1(24)

2、1(23)1(22)1(21)1(14)1(13)1(12)1(11,bbbbbxxxxxaaaaaaaaaaaaaaaaAbAx1000100010001)1(11)1(41)1(11)1(31)1(11)1(211aaaaaaM)2(44)2(43)2(42)2(34)2(33)2(32)2(24)2(23)2(22)2(14)2(13)2(12)2(111000aaaaaaaaaaaaaAM(2)0)1(11a设0)2(22a10001000100001)2(22)2(42)2(22)2(322aaaaM)3(44)3(43)3(34)3(33)3(24)3(23)3(22)3(14)3

3、(13)3(12)3(111200000aaaaaaaaaaaAMM)4(44)4(34)4(33)4(24)4(23)4(22)4(14)4(13)4(12)4(11123)3(33)3(433000000 ,100010000100001aaaaaaaaaaAMMMaaM0)3(33a方程组(2)就可以化为(4)(4)(4)(4)1111121314(4)(4)(4)22222324321(4)(4)333334(4)4444000000 xbaaaaxbaaaM M Mxbaaxba 其中M=M3M2M1为一个单位下三角阵。这个方程组就可以依次求出x4, x3, x2, x1,这就是高斯

4、消元法的过程。 x=U -1M b 或 x=U -1L -1b 推广到一般就是：对方程Ax=b (2)，存在一个单位下三角阵M，使得 MA x=Mb, 记U=MA，它是一个上三角阵，方程组化为 U x=Mb。我们记M的逆矩阵M -1=L，有A=LU(MA) 。可以解出 若若n阶矩阵阶矩阵A可逆可逆,当且仅当当且仅当顺序主子式顺序主子式不为零，则不为零，则A可分解可分解为一个单位为一个单位下三角阵下三角阵L 和一个和一个上三角阵上三角阵U 的积的积A=LU，分解是，分解是唯唯一一的。这就是所谓的的。这就是所谓的LU 分解分解。(数值计算方法P58杜里特尔分解惟一存在的一个充要条件)(1 )2 1

5、(1 )1 1(1 )3 1(1 )1 1(1 )4 1(1 )1 11000100010001aaaaaa(2)32(2)22(2)42(2)2210000100010001aaaa(3)43(3)331 0000 1000 0100 01aa命题命题 若若n 阶矩阵阶矩阵A可逆可逆，则存在，则存在交换阵交换阵P，使，使PA=LUL，U 分别为分别为单位下三角阵单位下三角阵和和上三角阵上三角阵。命题命题 正定对称矩阵可分解成正定对称矩阵可分解成对角元素为正对角元素为正的的下三角阵下三角阵与与它的它的转置矩阵转置矩阵之积，即之积，即A=L LT或者表为A = L D LT其中其中L是是单位下三

6、角阵单位下三角阵，D是是元素为正元素为正的的对角阵对角阵。这种分解称。这种分解称三角分解三角分解或或Cholesky分解分解。(数值计算方法P63定理3.23) 求解方程组(2)对0)(kkka的假设，等价于A的顺序主子顺序主子式式Dk0。在消元过程中，为避免因)(kkka绝对值太小而造成舍入误差的过大，在进行到第k 步时，都按列选择),()(nkiakik中最大的一个，称之为列主元列主元，将列主元所在行与第k行交换再按上面的高斯消元法进行下去，称为列主元素消元法列主元素消元法。1.2 误差分析，条件数误差分析，条件数 对于一般的方程组A x=b，如果解 x 对b 或 A 的扰动敏感，就称方程

7、组是病态病态的，也称系数矩阵 A 是病态的。向量范数向量范数和矩阵矩阵范数范数是定量估计 x 对 b 或 A 的扰动敏感程度的重要指标。 x对b 扰动的敏感程度取决于Cond(A)= |A-1|A|。我们称之为矩阵A的条件数条件数。条件数越大条件数越大，由由b的误差引起的的误差引起的x的误差可的误差可能越大能越大。 类似地，x 的（相对）误差也大致上达到 A 的(相对)误差的Cond(A)倍。关于条件数COND可用help命令等查阅。( ),xxrcond ArbAxxb1.3 MATLAB中用直接法解方程组中用直接法解方程组1. 求解求解A x=b 输入A, b之后，用 x=A b，即输出方

8、程组的解。2. LU分解，用列主元素消元法分解，用列主元素消元法L,U=lu(A) L为单位下三角阵与交换阵的积，U为上三角 阵，使A=LU3. 范数与条件数范数与条件数 L,U,p=lu(A) L为单位下三角阵，U同上，p为一交换阵， 使pA=LUU=chol(A) 对正定矩阵A的Cholesky分解，输出上三角阵 U，使A=UTUn=norm(X) 矩阵或向量X的2-范数c=cond(A) 向量或矩阵X的2-条件数求解求解A x=b可以用可以用 x = U ( L b)利用LU分解还可以快速计算矩阵的行列式与逆矩阵。1.4 迭代法迭代法对于大型稀疏矩阵大型稀疏矩阵(n 较大且零元素较多的矩

9、阵)适合用迭代法迭代法例例1 求解141035310214310321321321xxxxxxxxx将方程组改写成4 . 13 . 01 . 05 . 03 . 02 . 04 . 11 . 03 . 0213312321xxxxxxxxx这样改写的目的是：从一组初始值出发，可以进行迭代过程。我们先通过一个例子说明迭代法的基本思想。0, 0, 0)0(3)0(2)0(1xxx令这里给定了一组初值。方程组的迭代形式为：4 . 13 . 01 . 05 . 03 . 02 . 04 . 11 . 03 . 0)(2)(1)1(3)(3)(1)1(2)(3)(2)1(1kkkkkkkkkxxxxxx

10、xxx当计算到k=4时，就可得到与精确解很接近的近似解。这就是迭代法的基本思想。理想的情形是 收敛。), 2 , 1 , 0()(kxki1.5 雅可比迭代雅可比迭代 将A分解为A=DLU ，其中D = diag(a11, a22, , ann ),000 ,000 11121 121nnnnnnaaaUaaaL若对角阵D非奇异，将A x=b 改写为D x=(L+U )x+b，于是如果x(k)收敛于 x , 则 x 必为方程(27)的解，即A x=b 的解。x = D -1( L + U ) x + D -1b（27）则方程组的迭代形式为： B1 = D -1 (L + U ), f 1= D

11、 -1 b（28）x(k+1) = B1 x(k) + f1 ( k = 0, 1, 2,)（29）141035310214310321321321xxxxxxxxx例例 用雅可比迭代求解方程组A=10 3 1;2 -10 3;1 3 10;b=14 -5 14;D=diag(diag(A);L=-tril(A,-1);U=-triu(A,1);B=D(L+U);f=Db;x=0;0;0;for k=1:9 x=B*x+f; xend x1 x2 x3 1.4000 0.5000 1.4000 1.1100 1.2000 1.1100 0.9290 1.0550 0.9290 0.9906 0

12、.9645 0.9906 1.0116 0.9953 1.0116 1.0003 1.0058 1.0003 0.9982 1.0001 0.9982 1.0001 0.9991 1.0001 1.0003 1.0001 1.0003 1.0000 1.0001 1.00001.6 高斯高斯-赛德尔迭代赛德尔迭代我们先来回顾一下雅可比迭代的形式4 . 13 . 01 . 05 . 03 . 02 . 04 . 11 . 03 . 0)(2)(1)1(3)(3)(1)1(2)(3)(2)1(1kkkkkkkkkxxxxxxxxx在计算)1(2kx时，)1(1kx已经算出，而计算)1(3kx)1(

13、1kx)1(2kx时，和已经算出，4 . 13 . 01 . 0 5 . 03 . 02 . 0 4 . 11 . 03 . 0)1(2)1(1)1(3)(3)1(1)1(2)(3)(2)1(1kkkkkkkkkxxxxxxxxx因此，我们的迭代就可以改进为如下形式： 对一般的方程组 D x = (L + U ) x + b原来的迭代公式D x(k+1) = L x(k) + U x(k) + b就可以改进为D x(k+1) = L x(k+1) + U x(k) + b当D-L非奇异时，改写为x(k+1) = (D L)1Ux(k) + (D L)1b令 B2 = (D L)1U，f2 =

14、(D L)1b ，x(k+1) = B2 x(k ) + f 2 (k = 0,1,2,) 这就是高斯高斯-赛德尔迭代赛德尔迭代。于是B2 = (D L)1U，f2 = (D L)1b ，x(k+1) = B2 x(k ) + f 2 (k = 0,1,2,) 141035310214310321321321xxxxxxxxx用高斯-赛德尔迭代解方程组：A=10 3 1;2 -10 3;1 3 10;b=14 -5 14;D=diag(diag(A);L=tril(A,-1);U=triu(A,1);B=(D-L)U;f=(D-L)b;x=0;0;0;for k=1:6 x=B*x+f; xe

15、ndans = 1.4000 0.7800 1.0260ans = 1.0634 1.0205 0.9875ans = 0.9951 0.9953 1.0019ans = 1.0012 1.0008 0.9996ans = 0.9998 0.9998 1.0001ans = 1.0000 1.0000 1.00001.7 迭代法的收敛性迭代法的收敛性记一般的迭代公式为 x(k+1) = B x(k ) + f (k = 0,1,2,) （*）又设原方程组的解为x*，即 x* = B x* + f ，并与此式相减得：x(k ) x* = B k ( x(0) x* )于是，x(k ) 收敛于 x

16、* 等价于 Bk 趋于零，这又等价于B 的所有特征值的模都小于1，由此导出：1max)(iiB迭代公式迭代公式(*)收敛的充要条件收敛的充要条件是B 的的谱半径谱半径数值计算方法数值计算方法P81定理定理3.4.1若若A 是严格对角占优的，则是严格对角占优的，则雅可比雅可比和和高斯高斯-赛德尔赛德尔迭代均收敛。迭代均收敛。若若A 正定对称，则正定对称，则高斯高斯-赛德尔赛德尔迭代收敛。迭代收敛。将迭代公式递推 k 次function x,m=gs(A,b,x0,tol,n)D=diag(diag(A);L=-tril(A,-1);U=-triu(A,1);B=(D-L)U;f=(D-L)b;x

17、=x0;if max(abs(eig(B)=1 error(高斯赛得尔迭代不收敛)endfor k=1:n x=B*x+f; x if norm(A*x-b)=1 error(雅可比迭代不收敛)endfor k=1:n x=B*x+f; x if norm(A*x-b)tol break endend m=k;雅可比迭代雅可比迭代高斯赛得尔迭代高斯赛得尔迭代上机作业1、P110 12、P111 2、（1）（2）（3）超定线性方程组超定线性方程组：方程个数超过了未知数的个数.5.4 超定线性代数方程组的最小二乘法超定线性代数方程组的最小二乘法 线性最小二乘法线性最小二乘法曲线拟合问题的提法曲线拟

18、合问题的提法：已知一组（二维）数据，即平面上的 n 个点(xi , yi ), i=1,2,n, xi 互不相同，寻求一个函数（曲线）y= f(x)，使 f(x)在某种准则下与所有数据点最为接近与所有数据点最为接近，即曲线拟合得最好。1、最小二乘法概念、最小二乘法概念 设rk(x)是事先选定事先选定的一组函数，令 f(x) = a1 r1(x)+ a2 r2(x)+ am rm(x) (18)其中而 ak （k=1,m, mm), 且AA可逆,则为原方程的最小二乘解. 此时x=(AA)-1(Ab)2、x，y=lu(A) x,y,p=lu(A) u=chol(A) 3、矩阵的范数、条件数、秩、特

19、征值norm(x) norm(x,1) norm(x,inf) 范数cond(x) cond(x,1) cond(x,inf) 条件数rcond(x) 向量x的条件数倒数的估计值rank(A)e=eig(x) V,D=eig(x) 特征值 xV=VD(4) 稀疏矩阵的处理稀疏矩阵的处理建立稀疏矩阵的命令为a=sparse(r,c,v,m,n)r和c为矢量，分别是指矩阵中非零元素的行号和列号，v是一个全部为非零元素的矢量，m和n分别为稀疏矩阵的行数和列数。S=sparse(2 5 3 1 3 4,1 2 4 5 1 3,1 3 -1 2 5 4,5,5)FS = 0 0 0 0 2 1 0 0

20、0 0 5 0 0 -1 0 0 0 4 0 0 0 3 0 0 0FS=full(S)S = (2,1) 1 (3,1) 5 (5,2) 3 (4,3) 4 (3,4) -1 (1,5) 2按列的增序排列例例nbAnn21,414114114解A x = bn=1000;b=1:n;a1=sparse(1:n,1:n,4,n,n);a2=sparse(2:n,1:n-1,1,n,n);a=a1+a2+a2;tic;x=ab;toc; %计时A=full(a);tic;y=Ab;toc; %计时s1=sum(x)s2=sum(y)elapsed_time = 0.9900elapsed_tim

21、e = 3.9600s1 = 83452s2 = 83452检验两个结果是否相同二、实例二、实例nnnxqxbxxxxxpqpqpqpA00,11101221, 3 , 2, 021kqxpxxkkk 一年生植物的繁殖一年生植物的繁殖 前2.2.1（教材P26）中已推出下列方程成立bcabaqbcap)1 (,121其中xk 为第k年的植物数量,且现在问题是：在繁殖过程中，若已知某一年的植物数量，并希望若干年后，它的数量达到给定的规模。求出从第2年开始以后各年这种植物的数量。则可解稀疏矩阵方程为：2 2 投入产出与输电网络投入产出与输电网络2.1 投入产出问题投入产出问题表表1 国民经济各个部

22、门间投入产出关系表国民经济各个部门间投入产出关系表（产值 亿元） 产出投入 农 业 制造业 服务业外部需求 总产出 农 业 15 20 30 35 100 制 造 业 30 10 45 115 200 服 务 业 20 60 / 70 150 初 始 投 入 35 110 75 总 投 入 100 200 150(此表见数学实验P128） 假定每个部门的产出与投入是成正比的，由表1能够确定三个部门的投入产出表如下：表表2 投入产出表投入产出表 产出投入 农 业 制 造 业 服 务 业 农 业 0.15 0.10 0.20 制 造 业 0.30 0.05 0.30 服 务 业 0.20 0.30

23、 0(表中数字称为投入系数投入系数或消耗系数消耗系数) 1) 有n个部门，已知投入系数，给定外部需求，建立求解各部门总产出模型。 2) 设投入系数如表2，若今年对农业、制造业和服务业的外部需求分别为50、150、100亿元，问这三个部门的总产值分别是多少？ 3) 如果三个部门的外部需求分别增加一个单位，它们的总产出应分别增加多少？ 4) 对任给非负外部需求，都能得到非负的总产出，模型称为可行的，为使模型可行，投入系数应满足什么条件？问题：问题： 产出投入 农 业制造 业服务 业外部需求总产 出农 业 15 20 30 35 100 制 造 业 30 10 45 115 200 服 务 业 20

24、 60 / 70 150初始投入 35110 75 总 投 入100200150 产出投入农 业制造 业服务 业 农 业0.15 0.10 0.20 制 造 业0.30 0.05 0.30 服 务 业0.20 0.30 0记n个部门中第i 个部门总产出为xi,其中对第j个部门的投入为xij,满足的外部需求为di，则：(37) ), 2 , 1( 1nidxxnjiiji在每个部门的产出与投入成正比的假设下，如记第j个部门的单位产出需要第i个部门的投入为aij，则(38) ), 2 , 1,( njixxajijij将(38)式代入(37)即得(39) ), 2 , 1( 1nidxaxinjj

25、ijiA=( aij )nn，x=( x1, xn)T ，d = ( d1,dn)T ，(39)式可写作x = A x + d 或 ( I - A)x = d (40,41)如果( I - A)可逆，则 x =( I - A)-1d (42) a=0.15 0.1 0.2 ;0.3 0.05 0.3;0.2 0.3 0;d=50 150 100;b=eye(3)-a;x=bd,c=inv(b)x = 139.2801 267.6056 208.1377c = d= 1.3459 0.2504 0.3443 50 0.5634 1.2676 0.4930 150 0.4382 0.4304 1.

26、2167 100农业、制造业、服务业的总产出分别应为：139.2801267.6056 亿元208.1377对农业的外部需求每增加 1 个单位，农业、制造业、服务业的总产出应分别增加1.3459, 0.5634, 0.4382单位。对制造业的外部需求每增加 1 个单位，农业、制造业、服务业的总产出应分别增加0.2504, 1.2676, 0.4304单位。对服务业的外部需求每增加 1 个单位，农业、制造业、服务业的总产出应分别增加0.3443, 0.4930, 1.2167单位。x=cd根据 x =( I - A)-1d (42) 要使模型可行，即对任意的外部需求d 0，都能得到总产出x 0，

27、显然只需( I - A)-1 0。因为A 0(所有的aij0)，( I - A )( I + A + A2 + + Ak ) = I Ak + 1只要Ak 0(k )，就有 0)(11kkAAI而Ak 0等价于, 0kA而,kkAA故只要1A。我们选择1-范数，即只要(43) ), 2 , 1( 11njaniij模型就是可行的。这只要初始投入非负初始投入非负就是当然成立的了。), 2 , 1( 1njxxnijij（43）等价于2.2 输电网络输电网络 简化的大型输电网络如图，其中Ri，ri分别表示负载电阻和线路内阻，设电源电压为V。r1r2rnR1R2RnI0 = i1I1i2inI2In

28、V2) 设R1=R2 = =Rn=R ，r1=r2= =rn= r , 在R = 6，r =1, V=18, n=10的情况下求I1, I2, , In及总电流I0；3) 讨论n的情形。1) 列出各负载上电流 I1, I2, , In 的方程；r1r2rnR1R2RnI0 = i1I1i2inI2InV111122221111 nnnnnnIRIRirIRIRirVIRi r 11232121nnnnniIiiIiiIiiI 11322211nnnnnnnIiIIiIIIiIIIi 0)(0)()(11222211121111nnnnnnnIrRIRIrIrRIRVIrIrIrR(47)记nn

29、nrRRrrrRRrrrrRR1222211111100方程(47)变成R I = E (48) 回路定律节点定律00 ,21VEIIIInr=1;R=6;V=18;n=10;b1=sparse(1,1,V,n,1);b=full(b1);a1=triu(r*ones(n,n);a2=diag(R*ones(1,n);a3= tril(R*ones(n,n),-1)+. tril(R*ones(n,n),-2);a=a1+a2+a3;I=abI0=sum(I)00 ,00VbrRRrrrRRrrrrRarrrrrrrra00001000000000003RRa00000000000RRRRRR

30、RRR=I0 = 5.9970I = 2.0005 1.3344 0.8907 0.5955 0.3995 0.2702 0.1858 0.1324 0.1011 0.0867n=10n=20I0 = 6.0000 I = 2.0000 1.3333 0.8889 0.5926 0.3951 0.2634 0.1756 0.1171 0.0780 0.0520 0.0347 0.0231 0.0154 0.0103 0.0069 0.0047 0.0032 0.0023 0.0018 0.0015我们还可以令n=50或100.从运行的结果可以看出，当n(10)增大时，总电流I0几乎不变，说明总

31、电阻增加很少。Ik+1差不多总是Ik的2/3倍。下面我们从理论上讨论当n时，总电阻的变化。结果分析：R0(n)R0(n+1)Rr)()() 1(000nRRnRRrnR(49)以R0(1)=R+r 代入递推公式(49)，即可求出R0(n)，令00)(limRnRn则R0满足如下方程0 ,020000rRrRRRRRRrR即由此解得：)6, 1( 32420RrrRrrR当于是总电流6318010RViI而,32, , 263111020021021001IIRRRiRRRIiRRRiiRRRI)3, 6( 32001RRIIRRRInnn这正好证实了前面的计算结果。r1r2rnR1R2RnI0 = i1I1i2inI2InV0001IRRRRRI10020iRRRRiR