精通matlab科學(xué)計(jì)算

9通過(guò)本章的介紹，讀者既能應(yīng)用中相應(yīng)的函數(shù)求解線性方程組，又能通過(guò) 求逆在中用求逆法求解，線性方程可以直接用左除法“\”求解，也可以用求逆函數(shù)inv()求解。9-1左除法和求逆法求解線性方程組實(shí)例。用左除法和求逆法求解如下線性方

3x37x3解：用左除法求解，在命令窗口中輸入求解程序>>A=[123;-137;90>>b=[14>>>>x-用逆矩陣法求解，在命令窗口中輸入求解程序>>A=[123;-137;90>>b=[14>>>>x-由計(jì)算結(jié)果可知，方程組的解為[

9.2分解矩陣分解法是指根據(jù)一定的原理用某種算法將系數(shù)矩陣分解成若干個(gè)矩陣的代數(shù)運(yùn)LU分解、QR分解、Cholesky分解、Schur分解、Hessenberg分解和奇LULUDoolittle分解，它是將矩陣分解為一luA的LU分解，常用格式為：（，使得ALU；[L,U,P]lu(A)：產(chǎn)生一個(gè)上三角矩陣U和一個(gè)下三角矩陣L及置換矩陣P，使得AxbxU\(L\b)xU\(L\P*b)9-2】LULU1.534x4

x1

310x43x442解：用LU分解法求解，在命令窗口中輸入求解程序>>A=[1.53-0.84;20910;-74.8-0.61;1412.3-4>>b=[401->>>>>>x---由計(jì)算結(jié)果可知，方程組的解為QR

x2x3x40.37210.82911.5336,1.4547QR分解就是把矩陣分解為一個(gè)正交矩陣和一個(gè)上三角矩陣的乘積。中對(duì)矩陣A進(jìn)行QR分解的函數(shù)調(diào)用格式如下：[Q,R,E]qr(A)QRE，使得AxbxR\(Q\b)xE(R\(Q\b))9-3】QRQR0.5x20.3333x30.25x4 0.3333x20.25x30.2x4 0.25x20.2x30.1667x4 0.2x20.1667x30.1429x41解：用QR分解法求解，在命令窗口中輸入求解程序A=[0.1429b=[122Q、R--------=--0--00-000Rx=1.0e+003*由計(jì)算結(jié)果可知，方程組的解為Cholesky

x2x3x41161440,36002380RTRCholesky分解。中與Cholesky分解有關(guān)的函數(shù)如下：Rchol(A)：對(duì)A進(jìn)行Cholesky分解，使得ARTR；[Rp]chol(A)ACholeskyARTRAp0；否則，pA滿秩，Rp1階的上三RTRA(1:(p1)，1:(p1))AxbxR\(RT\b)9-4】Cholesky分解法求解線性方程組實(shí)例。用Cholesky分解法求解以下線性936x230x3 192x2180x3 180x2180x31解：用Cholesky分解法求解，在命令窗口中輸入求解程序>>A=[9-3630;-36192-180;30-180>>>>>>RR0-00x由計(jì)算結(jié)果可知，方程組的解為

在在中還提供了矩陣的奇異值分解、Hessenberg分解和Schur分解等函數(shù)，可用來(lái)求解線性方程組，如表9-1所示。9-1其他的矩陣分解函數(shù)及常用 [U,S,V]奇異值分解，使得AU*S*[P,H]Hessenberg分解,使得AP*H*PT，其中P為酉矩[U,T]Schur分解，使得AU*T*UT，其中U為正交矩陣，TSchur矩9-50.5x20.3333x30.25x40.2x5 0.3333x20.25x30.2x40.1667x5 0.25x20.2x30.1667x40.1429x5 0.2x20.1667x30.1429x40.125x50.20.1667x20.1429x30.125x40.1111x51解：用奇異值分解法求解，在命令窗口中輸入求解程序A=[b=[1010U、S、VU=- -- - -- -- S=00 00 00 000 000 V------ -- x=1.0e+004* 奇異值分解很奇異值分解很有用，將系數(shù)矩陣進(jìn)行奇異值分解以后，A的奇異值都在S的對(duì)角線上，這樣就可以粗略估計(jì)A的條件數(shù)，看看A是否是 UVS是對(duì)9-6】HessenbergHessenberg分解法求解以下線 x2x3x4 2x23x34x4 3x26x310x4 4x210x320x42解：用Hessenberg分解法求解， A=11112313614>>>>b=[147->>P、HP-0--00000H-0000-00x=

x2x3x41033,3612Hessenberg分解法分解后產(chǎn)生正交矩陣和Hessenberg矩陣，正交矩陣的逆就是其轉(zhuǎn)置；如果系數(shù)矩陣對(duì)稱的話，產(chǎn)生的Hessenberg矩陣就是三對(duì)角9-7】Schur分解法求解線性方程組實(shí)例。Schur x2x3x4 2x23x34x4 3x26x310x4 4x210x320x42解：用Schur分解法求解， >> 1112313614>>b=[147-U、T-=000000000000Tx--由計(jì)算結(jié)果可知，方程組的解為

x2x3x41033,3612產(chǎn)生的Schur矩陣就是對(duì)角陣，顯然對(duì)角陣的逆矩陣非常容易得到，因此9.3迭代105 逐次近n階線性方程組，AxbA（A為非奇異）Q為非奇異。其中下，逐次近法收斂的充分條件是迭代矩陣的譜半徑小于1。里迭代 功能：用里迭代法求線性方程組Axb的解。調(diào)用格式：[x,n]=richason(Abx0,eps,其中，A為線性方程組的系數(shù)矩陣；b為線性方程組中的常數(shù)向量；x0為迭代初始向量；eps為解的精度控制（此參數(shù)可選M為迭代步數(shù)控制（此參數(shù)可選；x為線性方程組的解；n為求出所需精度的解實(shí)際的迭代步數(shù)。 function%采用里迭代法求線性方程組Ax=b的if(nargin==3)eps1.0e- M elseif(nargin==M=I=eye(size(A));nn 【例9-8】迭代法求解線性方程組實(shí)例用迭代法求解以下線性方程組，

3331解：用迭代法求解，在命令窗口中輸入求解程序>>A=[1.017- - >>b=[10>>x0=[00>>xn[x３][0.9738,-0.0047,1.001] >>ans，經(jīng)過(guò)驗(yàn)證也確定了解的正確性一般情況下迭代法難以收斂因?yàn)橐梗珹0AJacobi 中編程實(shí)現(xiàn)的Jacobi迭代法函數(shù)為：jacobi。功能：用Jacobi迭代法求線性方程組Axb的解調(diào)用格式：[x,njacobiAbx0,eps,其中，A為線性方程組的系數(shù)矩陣；b為線性方程組中的常數(shù)向量；x0為迭代初始向量；eps為解的精度控制（此參數(shù)可選；varargin為迭代步數(shù)控制（此參數(shù)可選x為線性方程組的解；Jacobi迭代法 function%采用Jacobi迭代法求線性方程組Ax=b的ifnargin==3eps=1.0e-6;M=200;elseifnargin<3elseifnarginM= %求A的對(duì)角矩陣 %求A的下三角陣 %求A的上三角矩陣 whilenorm(x-x0)>=epsx=B*x0+f;9-9】JacobiJacobi迭代法求解以下線性方程組，0.98890.0005x20.0002x3 0.9946x20.0077x3 0.0092x20.9941x31解：用Jacobi迭代法求解，在命令窗口中輸入求解程序>>A=[ -0.0005-- - >>b=[10>>x0=>>x=-n4步迭代，Jacobi JacobiJacobix0B1Gauss-SeidelAA 中編程實(shí)現(xiàn)的Gauss-Seidel迭代法函數(shù)為：gauseidel。功能：用Gauss-Seidel迭代法求線性方程組Ax=b的解。調(diào)用格式：[x,ngauseidelAbx0,eps,M)。其中，A為線性方程組的系數(shù)矩陣；x0eps為解的精度控制（此參數(shù)可選M為迭代步數(shù)控制（此參數(shù)可選；x為線性方程組的解；Gauss-Seidel function%采用Gauss-Seidel迭代法求線性方程組Ax=b的ifnargin==3eps=1.0e-6;M=200;elseifnargin==4M=200;elseifnargin<3 %求A的對(duì)角矩 %求A %求A whilenorm(x-x0)>=eps9-10】Gauss-SeidelGauss-Seidel迭代法求解以1.44490.7948x20.8801x3 1.9568x20.1730x3 0.5226x21.9797x31解：用Gauss-Seidel迭代法求解，在命令窗口中輸入求解程序A=b=[10>>>>xn[3][0.5929,0.2443,0.3835]AD、L、UGauss-Seidel是一個(gè)事先選好的常數(shù)，稱為松弛因子，當(dāng)1時(shí)叫超松弛法，也叫SOR迭代法；當(dāng)1時(shí)叫低松弛法。其迭代公式為：若系數(shù)矩陣A對(duì)稱正定，當(dāng)02，SOR迭代法收斂。 功能：用超松弛迭代法求線性方程組Ax=b的解。調(diào)用格式：[xn]=SORA,bx0,weps,其中，A為線性方程組的系數(shù)矩陣；b為線性方程組中的常數(shù)向量；x0為迭代初始向量；weps為解的精度控制（此參數(shù)可選M為迭代步數(shù)控制（此參數(shù)可選；x為線性方程組的解；n為求出所需精度的解實(shí)際的迭代步數(shù)。 function[x,n]=SOR(A,b,x0,%采用超松弛迭代法求線性方程組Ax=b的ifnargin==4eps=1.0e-6;M=200;elseifnargin<4elseifnarginM=if(w<=0||w>=2) %求A的對(duì)角矩陣 %求A的下三角矩陣 %求A的上三角矩 whilenorm(x-x0)>=epsx=B*x0+f;9-11】超松弛迭代法求解線性方程組實(shí)例。用超松弛迭代法求解以下線性方程4x13x2 4x2x3x24x3解：用超松弛迭代法求解，在命令窗口中輸入求解程序>>A=[430;34-1;0-1>>b=[2430->>x0=[00>>xn可見，經(jīng)過(guò)14步迭代，SOR迭代法求出了方程組的解為[2,x3][3,4,5](DL)xk1/2(Uxkb)(1)Dxk(DU)xk1(Lxkb)(1)Dxk1/2在中編程實(shí)現(xiàn)的對(duì)稱逐次超松弛迭代法函數(shù)為：SSOR。功能：對(duì)稱逐次超松弛迭代法求線性方程組Axb的解。調(diào)用格式：[x,n]SSORAbx0weps,M)。其中，A為線性方程組的系數(shù)矩陣；x0weps為解的精度控制（此參數(shù)可選M為迭代步數(shù)控制（此參數(shù)可選；x為線性方程組的解；SSOR算法用實(shí)現(xiàn)如下function[x,n]=SSOR(A,b,x0,%采用對(duì)稱逐次超松弛迭代法求線性方程組Ax=b的ifnargin==4eps=1.0e-6;M=200;elseifnargin<4elseifnarginM=if(w<=0||w>=2) %求A的對(duì)角矩 %求A %求A x=B2*x12+f2; whilenorm(x-x0)>=epsx=B2*x12+f2;9-12】SSOR迭代法求解以下4x13x2 4x2x3x24x3解：用SSOR迭代法求解，在命令窗口中輸入求解程序>>A=[430;34-1;0-1>>b=[2430->>x0=[00x,n]=SSOR(A,b,x0,1.25松弛因子取xn可見，經(jīng)過(guò)34步迭代，SSOR迭代法求出方程組的解為[2,x3][3,4,5]需的步驟來(lái)看，SOR迭代法（14步迭代）比SSOR迭代法（34步迭代）要快一些。(DU)xk1Lxk1/2bD、L、UDA分解的對(duì)角矩陣，LA分解的下三角矩陣，U為矩陣A分解的上三角矩陣。在中編程實(shí)現(xiàn)的兩步迭代法函數(shù)為：twostep。功能：用兩步迭代法求線性方程組Ax=b的解。調(diào)用格式：[x,ntwostepAbx0,eps,varargin)。其中，A為線性方程組的系數(shù)矩陣；x0eps為解的精度控制（此參數(shù)可選；varargin為迭代步數(shù)控制（此參數(shù)可選x為線性方程組的解；n為求出所需精度的解實(shí)際的迭代步數(shù)。兩步迭代法用實(shí)現(xiàn)如下：function%采用兩步迭代法求線性方程組Ax=b的ifnargin==3eps=1.0e-6;M=200;elseifnargin<3elseifnarginM= %求A的對(duì)角矩陣 %求A的下三角矩陣 %求A的上三角矩x=B2*x12+f2; whilenorm(x-x0)>=epsx0=x;x=B2*x12+f2;9-13】?jī)刹降ㄇ蠼饩€性方程組實(shí)例。采用兩步迭代法求解以下線性方程組，4x13x23 33x24x324解：用兩步迭代法求解，在命令窗口中輸入求解程序>>A=[430;34-1;0-1>>b=[2430->>x0=[00>>xn30步迭代，兩步迭代法求出了方程組的解為

φ(x)xATx-在中編程實(shí)現(xiàn)的最速下降法函數(shù)為：fastdown。功能：用最速下降法求線性方程組Ax=b的解。調(diào)用格式：[x,nfastdownAbx0eps)。其中，A為線性方程組的系數(shù)矩陣；x0xn為求出所需精度解的實(shí)際的搜索步數(shù)。最速下降法用實(shí)現(xiàn)如下：function%采用最速下降法求線性方程組Ax=b的if(nargin==3)eps=1.0e-6;r=b-d=dot(r,r)/dot(A*r,r);x=x0+d*r;x0=x;r=b-d=dot(r,r)/dot(A*r,r);x=x0+d*r;n=n+】且分別取迭代初始值[1,1,1,1和[2,10,-1,2]，分析迭代初值對(duì)迭代步數(shù)的影響。0.5x20.3333x30.25x4 0.3333x20.25x30.2x4 0.25x20.2x30.1667x4 0.2x20.1667x30.1429x4解：用最速下降法求解，在命令窗口中輸入求解程序A=[>>b=[101x0=[1111]';%>>x=1.0e+003*n時(shí)，經(jīng)過(guò)解為 2,x3,x4[328.8,-3608.8,8586,-5540.4]在命令窗口中接著輸入：x0=[210-12]';%>>x=1.0e+003*n可見，當(dāng)?shù)跏贾禐閇2,10,-1,2時(shí)，經(jīng)過(guò)81365的解為 2,x3,x4[328.8,-3608.8,8586,-5540.4]在中編程實(shí)現(xiàn)的共軛梯度法函數(shù)為：conjgrad。功能：共軛梯度法求線性方程組Ax=b的解。調(diào)用格式：[x,n]=conjgrad(A,b,x0)。其中，A為線性方程組的系數(shù)矩陣；b為線性方程組中的常數(shù)向量；x0x共軛梯度法用實(shí)現(xiàn)如下function%采用共軛梯度法求線性方程組Ax=b的if(nargin==3)eps=1.0e-6;r1=b-A*x0;p1=r1;d=dot(r1,r1)/dot(p1,A*p1);x=x0+d*p1;r2=r1-f=dot(r2,r2)/dot(r1,r1);p2=r2+f*p1;n=x0=x;p1=p2;r1=d=dot(r1,r1)/dot(p1,A*p1);x=x0+d*p1;r2=r1-f=dot(r2,r2)/dot(r1,r1);p2=r2+f*p1;n=n+d=dot(r2,r2)/dot(p2,A*p2);x=x+d*p2;n=n+9-15】共軛梯度法求解線性方程組實(shí)例。采用共軛梯度法求解以下線性方程組，25x1302 400x21050x118900x2 2880x2解：用共軛梯度法求解，在命令窗口中輸入求解程序>>A= - ------b=[53-0->>x0=[0000>>x=n6為 3,x4,x5[5.7667,2.9167,1.9310,1.4333,1.1349]輸入下面的表達(dá)式以驗(yàn)證結(jié)果，在命令窗口中輸入>>ans當(dāng)Axb為方程組時(shí)，共軛梯度收斂得很慢。預(yù)處理技術(shù)是在用共軛梯度法對(duì)方程Axb用共軛梯度法求解，則FTx為原方程的解。理想情況下是取FA-11步驟2，取系數(shù)矩陣A的三條對(duì)角線構(gòu)成的矩陣，再對(duì)其進(jìn)行Cholesky分解，將F還有一些其他的取法，具體可參考數(shù)值計(jì)算的書籍。中有一個(gè)對(duì)應(yīng)的處理函數(shù)pcg()，它的用法如下：[x,flag]pcg(A,b,tol,maxit,M1,M2,x0,p1,p2,…)；[x,flag,relres]pcg(A,b,tol,maxit,M1,M2,x0,p1,p2,…)；[x,flag,relres,iter]pcg(A,b,tol,maxit,M1,M2,x0,p1,p2,…)；[x,flag,relres,iter,resvec]其中，tol為迭代的精度，maxit為迭代的最大步數(shù)，M1M2內(nèi)不收斂；2M為壞條件的預(yù)處理因子；3表示兩次連續(xù)迭代完全相同；4表relres為相對(duì)誤差：norm(b-iterx9-16】預(yù)處理的共軛梯度法求解線性方程組實(shí)例。采用預(yù)處理的共軛梯度法求25x1302 400x21050x118900x2 2880x2A=-----b=[53--0-->>x0=[0000M=pascal(5);%預(yù)處理矩>>[x,flag,re,it]=pcg(A,b,1.e-8,1000,M=111111234515 xflag0rexit=可見，當(dāng)?shù)跏贾禐閇0,0,0,0,0]10程組的解為 3,x4,x5[ 9-2其他的線性方程組迭代法求解函函數(shù)及語(yǔ) x=線性方程組的LQ解x=線性方程組的雙共軛梯度解x=線性方程組的穩(wěn)定雙共軛梯度x=線性方程組的共軛梯度的LSQR解x=線性方程組的廣義最小殘差解x=線性方程組的最小殘差解x=線性方程組的準(zhǔn)最小殘差解】3x28x3 4x27x3 6x29x3解：用最小殘差法求解，在命令窗口中輸入求解程序>>A=[038;147;-26>>b=[11>>[x,flag,re,x=flag1re=xit=由計(jì)算結(jié)果可知，經(jīng)過(guò)迭代3步，flag=1,表明默認(rèn)情況下不收斂，因此需要選擇其9.4特殊追趕法首先對(duì)系數(shù)矩陣A做如下三角分解：L為下三角的兩對(duì)角矩陣，U1。求解時(shí)先求解Lyb的解y，再求解Uxy的解x。在中編程實(shí)現(xiàn)的追趕法函數(shù)為：followup。功能：用追趕法求線性方程組Ax=b的解。調(diào)用格式：xfollowupA,b)其中，A為線性方程組的系數(shù)矩陣；b為線性方程組中的常數(shù)向量；x為線性方程組的解。追趕法用實(shí)現(xiàn)如下function%采用追趕法求線性方程組Ax=b的n=rank(A);disp('Error:對(duì)角有元素為0！');d=a=ones(n-c=ones(n-d(n,1)=%求解Ly=b的解y，解保存在b中d(i,1)=d(i,1)-(a(i-1,1)/d(i-1,1))*c(i-b(i,1)=b(i,1)-(a(i-1,1)/d(i-1,1))*b(i-%求解Ux=y的解x(n,1)=x(i,1)=(b(i,1)-9-18】2.5x1x2 1.5x2x3 0.5x3x41 0.5x4x5 1.5x5x61 x52.5x6解：用追趕法求解，在命令窗口中輸入求解程序>>A=0000000000000/