《matlab在科学计算中的应用线性方程组的直接求解分析方法》由会员分享，可在线阅读，更多相关《matlab在科学计算中的应用线性方程组的直接求解分析方法（56页珍藏版）》请在金锄头文库上搜索。
1、4.2.3 线性方程组的直接求解分析方法矩阵的三角分解：LU分整理课件整理课件整理课件格式L是一个单位下三角矩阵，u是一个上三角矩阵整理课件 l,u=lu(A) lP-1 Ll = 0.5000 1.0000 0 1.0000 0 0u =例：用两种方法对A进行LU分解 A=1 2 3; 2 4 1; 4 6 7;整理课件 l,u,p=lu(A)l = 1.0000 0 0 0.5000 1.0000 0u =p = 0 0 1 0 1 0 1 0 0整理课件QR分解 将矩阵A分解成一个正交矩阵与一个上三角矩阵的乘积。格式： Q,R = qr(A) 求得正交矩阵Q和上三角阵R，Q和R满足A=Q
2、R。 整理课件例： A = 1 2 3;4 5 6; 7 8 9; 10 11 12; Q,R = qr(A)Q =R = 0 0 0整理课件Cholesky(乔里斯基 )分解 整理课件格式： D=chol(A)整理课件例：进行Cholesky分解。 A=16 4 8; 4 5 -4; 8 -4 22; D=chol(A)D = 4 1 2 0 2 -3 0 0 3整理课件（1）LU分解： A*X=b 变成 L*U*X=b所以 X=U(Lb) 这样可以大大提高运算速度。例：求方程组 的一个特解。解： A=4 2 -1;3 -1 2;11 3 -1; B=2 10 8; D=det(A)D =
3、10利用矩阵的LU、QR和cholesy分解求方程组的解整理课件L,U=lu(A)L = 0.2727 1.0000 0 1.0000 0 0U =整理课件 X=U(LB)X = A*Xans =整理课件（2）Cholesky分解 若A为对称正定矩阵，则Cholesky分解可将矩阵A分解成上三角矩阵和其转置的乘积，方程 A*X=b 变成 R*R*X=b所以 X=R(Rb) （3）QR分解 对于任何长方矩阵A，都可以进行QR分解，其中Q为正交矩阵，R为上三角矩阵的初等变换形式，即：A=QR方程 A*X=b 变形成 QRX=b所以 X=R(Qb) 这三种分解，在求解大型方程组时很有用。其优点是运算
4、速度快、可以节省磁盘空间、节省内存。整理课件三个变换 在线性方程组的迭代求解中，要用到系数矩阵A的上三角矩阵、对角阵和下三角矩阵。此三个变换在MATLAB中可由以下函数实现。上三角变换： 格式 triu(A,1)对角变换： 格式 diag(A)下三角变换： 格式 tril(A,-1)例：对此矩阵做三种变换。整理课件 A=1 2 -2;1 1 1;2 2 1; triu(A,1)ans = 0 2 -2 0 0 1 0 0 0 tril(A,-1)ans = 0 0 0 1 0 0 2 2 0 b=diag(A); bans = 1 1 1整理课件4.3 线性方程组的迭代解法迭代法的一般形式 线
5、性方程组：整理课件4.3.1 Jacobi迭代法整理课件方程组 Ax=b 可写成 x=Bx+f 由此可构造迭代法： x(k+1)=Bx(k)+f 其中：B=I-D-1A=D-1(L+U) , f=D-1b. D=diag(A), L,U分别为-A的严格下三角和严格上三角矩阵整理课件function y=jacobi(a,b,x0)D=diag(diag(a); U=-triu(a,1); L=-tril(a,-1);B=D(L+U); f=Db;y=B*x0+f;n=1; x0=y; y=B*x0+f; n=n+1;endn整理课件例：用Jacobi方法求解，设x(0)=0,精度为10-6。
6、a=10 -1 0; -1 10 -2; 0 -2 10; b=9; 7; 6; jacobi(a,b,0;0;0)n = 11ans =整理课件4.3.2 Gauss-Seidel迭代法循环计算过程中存有两个向量：整理课件整理课件由此得原方程 Axb 求解的迭代方程：Ddiag(A)L和U分别为矩阵-A的严格下三角和严格上三角矩阵其中：整理课件function y=seidel(a,b,x0)D=diag(diag(a);U=-triu(a,1);L=-tril(a,-1);G=(D-L)U ;f=(D-L)b;y=G*x0+f; n=1; x0=y; y=G*x0+f; n=n+1;end
7、n整理课件例：对上例用Gauss-Seidel迭代法求解 a=10 -1 0; -1 10 -2; 0 -2 10; b=9; 7; 6; seidel(a,b,0;0;0)n = 7ans =整理课件例：分别用Jacobi和G-S法迭代求解,看是否收敛。Jacobi迭代法： a=1 2 -2; 1 1 1; 2 2 1; b=9; 7; 6; jacobi(a,b,0;0;0)n = 4ans = -27 26 8整理课件 seidel(a,b,0;0;0)n = 1011ans = 1.0e+305 * -Inf InfSeidel迭代法：对有些问题Gauss-Seidel迭代比Jacob
8、i迭代收敛快对有些问题Gauss-Seidel迭代比Jacobi迭代收敛的慢对有些问题Jacobi迭代收敛但Gauss-Seidel迭代发散整理课件4.3.3 SOR迭代法 在很多情况下，J法和G-S法收敛较慢，所以考虑对G-S法进行改进。于是引入一种新的迭代法逐次超松弛迭代法（Succesise Over-Relaxation),记为SQR法。 迭代公式为： X(k+1)= (D-wL)-1(1-w)D+wU)x(k) + w(D-wL)-1 b 其中：w最佳值在1, 2）之间，不易计算得到，因此 w通常有经验给出。整理课件function y=sor(a,b,w,x0)D=diag(dia
9、g(a);U=-triu(a,1);L=-tril(a,-1);M=(D-w*L)(1-w)*D+w*U); f=(D-w*L)b*w;y=M*x0+f; n=1; x0=y; y=M*x0+f; n=n+1;endn整理课件例：上例中，当w=1.103时，用SOR法求解原方程。 a=10 -1 0; -1 10 -2; 0 -2 10; b=9; 7; 6; sor(a,b,1.103,0;0;0)n = 8ans =整理课件4.3.4 两步迭代法 当线性方程系数矩阵为对称正定时，可用一种特殊的迭代法来解决，其迭代公式为： （D-L）x(k+1/2) =U x(k) +b （D-U）x(k+
10、1)=Lx(k+1/2) +b= x(k+1/2) =(D-L)-1 U x(k) + (D-L)-1 b x(k+1)= (D-U)-1 Lx(k+1/2) + (D-U)-1 b整理课件function y=twostp(a,b,x0)D=diag(diag(a);U=-triu(a,1);L=-tril(a,-1);G1=(D-L)U; f1=(D-L)b;G2=(D-U)L; f1=(D-U)b;y=G1*x0+f1; y=G2*y+f2; n=1; x0=y; y=G1*x0+f1; y=G2*y+f2; n=n+1;endn整理课件例：求解方程组 a=10 -1 2 0; -1 1
11、1 -1 3; 2 -1 10 3; 0 3 -1 8; b=6; 25; -11; 15; twostp(a, b, 0; 0; 0; 0)n = 7ans =整理课件4.4 线性方程组的符号解法 在MATLAB的Symbolic Toolbox中提供了线性方程的符号求解函数，如 linsolve(A,b) 等同于 X = sym(A)sym(b). solve(eqn1,eqn2,.,eqnN,var1,var2,.,varN )整理课件例： A=sym(10,-1,0;-1,10,-2;0,-2,10)； b=(9; 7; 6)； linsolve(A,b)ans = 473/475 9
12、1/95 376/475 vpa(ans)ans =整理课件例： x,y = solve(x2 + x*y + y = 3,x2 - 4*x + 3 = 0,x,y) x = 1 3 y = 1 -3/2 整理课件4.5 稀疏矩阵技术稀疏矩阵的建立：格式 S=sparse(i,j,s,m,n)生成一mxn阶的稀疏矩阵，以向量i和j为坐标的位置上对应元素值为s。例： n=5; a1=sparse(1:n, 1:n, 4*ones(1,n), n, n)a1 = (1,1) 4 (2,2) 4 (3,3) 4 (4,4) 4 (5,5) 4整理课件例： a2=sparse(2:n, 1:n-1,o
13、nes(1,n-1),n,n)a2 = (2,1) 1 (3,2) 1 (4,3) 1 (5,4) 1 full(a2)ans = 0 0 0 0 0 1 0 0 0 0 0 1 0 0 0 0 0 1 0 0 0 0 0 1 0整理课件例：n=5,建立主对角线上元素为4，两条次对角线为1的三对角阵。 n=5; a1=sparse(1:n,1:n,4*ones(1,n),n,n); a2=sparse(2:n,1:n-1,ones(1,n-1),n,n); a=a1+a2+a2a = (1,1) 4 (2,1) 1 (1,2) 1 (2,2) 4 (3,2) 1 (2,3) 1 (3,3) 4
14、 (4,3) 1 整理课件 (3,4) 1 (4,4) 4 (5,4) 1 (4,5) 1 (5,5) 4 full(a)ans = 4 1 0 0 0 1 4 1 0 0 0 1 4 1 0 0 0 1 4 1 0 0 0 1 4整理课件格式 A=spdiags(B,d,m,n) 生成一mxn阶的稀疏矩阵，使得B的列放在由d指定的位置。例： n=5 b=spdiags(ones(n,1),4*ones(n,1),ones(n,1),-1,0,1,n,n); full(b)ans = 4 1 0 0 0 1 4 1 0 0 0 1 4 1 0 0 0 1 4 1 0 0 0 1 4整理课件格式
15、： spconvert(dd) 对于无规律的稀疏矩阵，可使用此命令由外部数据转化为稀疏矩阵。调用形式为：先用load函数加载以行表示对应位置和元素值的.dat文本文件，再用此命令转化为稀疏矩阵。例：无规律稀疏矩阵的建立。首先编制文本文件sp.dat如下：5 5 0整理课件 spconvert(sp)ans = (5,1) 5 (4,4) 2 (3,5) 8 full(ans)ans = 0 0 0 0 0 0 0 0 0 0 0 0 0 0 8 0 0 0 2 0 5 0 0 0 0整理课件稀疏矩阵的计算： 同满矩阵比较，稀疏矩阵在算法上有很大的不同。具体表现在存储空间减少，计算时间减少。例：
16、比较求解下面方程组n1000时两种方法的差别。整理课件 n=1000; a1=sparse(1:n,1:n,4*ones(1,n),n,n); a2=sparse(2:n,1:n-1,ones(1,n-1),n,n); a=a1+a2+a2; b=ones(1000,1); tic; x=ab; t1=toct1 = a=full(a); tic; x=ab; t2=toct2 =整理课件一般矩阵的特征值与特征向量格式： d=eig (A) 只求解特征值。格式： V, D=eig (A) 求解特征值和特征向量。整理课件例：直接求解： A=16 2 3 13; 5 11 10 8; 9 7 6
17、12; 4 14 15 1; eig(A)ans =整理课件精确解： eig(sym(A)ans = 0 34 4*5(1/2) -4*5(1/2)高精度数值解： vpa(ans,70)ans = 0 34.28 97083588981642084整理课件同时求出特征值与特征向量：直接求解： v, d = eig(A)v =d = 34.0000 0 0 0 0 8.9443 0 0 0 0 -8.9443 0整理课件解析解： v,d=eig(sym(A)v = -1, 1, -8*5(1/2)-17, 8*5(1/2)-17 -3, 1, 4*5(1/2)+9, -4*5(1/2)+9 3,
18、 1, 1, 1 1, 1, 4*5(1/2)+7, -4*5(1/2)+7 d = 0, 0, 0, 0 0, 34, 0, 0 0, 0, 4*5(1/2), 0 0, 0, 0, -4*5(1/2) 整理课件4.6.2 矩阵的广义特征向量问题 若B=I，则化成普通矩阵特征值问题。格式： d=eig (A，B) 求解广义特征值。格式： V, D=eig (A，B) 求解广义特征值和特征向量。整理课件例：直接求解： A=5,7,6,5; 7,10,8,7; 6,8,10,9; 5,7,9,10; B=2,6,-1,-2; 5,-1,2,3; -3,-4,1,10; 5,-2,-3,8; V,
19、 D = eig(A, B)V = 0.3697 -0.3741 + 0.6259i -0.3741 - 0.6259i 1.0000 0.9948 -0.0674 - 0.2531i -0.0674 + 0.2531i -0.6090 0.7979 0.9239 + 0.0264i 0.9239 - 0.0264i -0.2316 1.0000 -0.6599 - 0.3263i -0.6599 + 0.3263i 0.1319 整理课件D = 4.7564 0 0 0 0 0.0471 + 0.1750i 0 0 0 0 0.0471 - 0.1750i 0 检验： norm(A*V-B*V*D)ans = 符号运算工具箱中的eig( )函数不支持广义特征值的运算。整理课件整理课件
