實(shí)驗線性代數(shù)方程組的數(shù)值解法

1、實(shí)驗 5 線性代數(shù)方程組的數(shù)值解法精 94 2009010586實(shí)驗?zāi)康?.學(xué)會用軟件數(shù)值求解線性代數(shù)方程組，對迭代法的收斂性和解的穩(wěn)定性作初步分析；2.通過實(shí)例學(xué)習(xí)線性代數(shù)方程組解決簡化的實(shí)際問題。實(shí)驗內(nèi)容3.已知方程組 ，其中 A ，定義見書。試通過迭代法求解此方程組，認(rèn)識迭代法收斂的含義以及迭代初值和方程組系數(shù)矩陣性質(zhì)對收斂速度的影響。實(shí)驗要求：(1). 選取不同的初始向量 和不同的方程組右端項向量 b，給定迭代誤差-賽要求，用雅可比迭代法和迭代法計算，觀測得到的迭代向量序列是否均收斂？若收斂，迭代次數(shù)，分析計算結(jié)果并得出你的結(jié)論；(2). 取定右端向量 b 和初始向量 ，將 A 的主對

2、角元素成倍增加若干次，非主對角線元素不變，每次用雅可比迭代法計算，要求迭代誤差滿足 ，比較收斂速度，分析現(xiàn)象并得出你的結(jié)論。 分析設(shè)計分析題目，主要有如下關(guān)鍵點(diǎn)：雅可比迭代法和-賽迭代法的編寫；選取不同數(shù)據(jù)進(jìn)行比較分析。采用函數(shù)的方式編寫雅可比和告訴迭代法（參見 參考程序），說明如下：function x y=jacobi(A,b,x0,e,limit)參數(shù)解釋：x返回值，為解的最終結(jié)果；1y返回值，迭代的過程，y(:,1)存放初始向量初值，y(:,2)存放右端項，后依次儲存每次迭代結(jié)果；A，b，x0依次為矩陣，右端項和初始向量；e迭代誤差要求，這里全部采用前后迭代向量之差的 -范數(shù)定義；li

3、mit迭代次數(shù)上限，超過則推出迭代過程； r屏顯，特征值的最大值，即譜半徑sum迭代次數(shù)（上限為 limit）編程思路：A 分解為 D、L、U迭代計算判斷迭代誤差條件反復(fù)以上兩步直至滿足誤差條件或者超出迭代次數(shù)界限function x y=gaussseideil(A,b,x0,e,limit)參數(shù)設(shè)置，思路均同 jacobi，不再贅述。結(jié)果 1：Jacobi：r=0.4893 gaussseideil：r=0.2523分析 1：對于改變 和 b 的情況，因為 A 不變，所以譜半徑 r 均不變，如上所示。結(jié)果 2：(1). 如下表2x0b0.0000 1.00000.33330.33330.3

4、4490.34650.34740.34770.0000 0.00000.00000.05560.06020.06440.06550.06600.0000 0.00000.00000.02780.03700.04050.04190.04250.0000 0.00000.00000.00000.00930.01160.01300.01350.0000 0.00000.00000.00000.00230.00460.00550.00600.0000 0.00000.00000.00000.00000.00120.00180.00210.0000 0.00000.00000.00000.00000.0

5、0020.00060.00080.0000 0.00000.00000.00000.00000.00000.00010.0002(2). 如下表3X0b01 0.3333 0.3333 0.3681 0.3927 0.4026 0.4078 0.4100 0.4110 0.4115 0.411700 0.0000 0.0556 0.1528 0.1840 0.2034 0.2111 0.2148 0.2165 0.2173 0.217600 0.0000 0.3056 0.4074 0.4630 0.4862 0.4973 0.5022 0.5045 0.5055 0.506000 0.000

6、0 0.5556 0.6574 0.7292 0.7550 0.7680 0.7737 0.7763 0.7775 0.77800.0000 0.00000.00000.00000.00000.00000.00000.00010.0000 0.00000.00000.00000.00000.00000.00000.00000.0000 0.00000.00000.00000.00000.00000.00000.00000.0000 0.00000.00000.00000.00000.00000.00000.00000.0000 0.00000.00000.00000.00000.00000.0

7、0000.00000.0000 0.00000.00000.00000.00000.00000.00000.00000.0000 0.00000.00000.00000.00000.00000.00000.00000.0000 0.00000.00000.00000.00000.00000.00000.00000.0000 0.00000.00000.00000.00000.00000.00000.00000.0000 0.00000.00000.00000.00000.00000.00000.00000.0000 0.00000.00000.00000.00000.00000.00000.0

8、0000.0000 0.00000.00000.00000.00000.00000.00000.00000.0000 1.00000.33330.34570.34760.34790.0000 0.00000.05560.06470.06600.06620.0000 0.00000.03700.04180.04270.04290.0000 0.00000.01080.01330.01380.01390.0000 0.00000.00490.00600.00630.00640.0000 0.00000.00170.00230.00240.00240.0000 0.00000.00070.00090

9、.00100.00100.0000 0.00000.00030.00040.00040.00040.0000 0.00000.00010.00010.00020.00020.0000 0.00000.00000.00010.00010.00010.0000 0.00000.00000.00000.00000.00000.0000 0.00000.00000.00000.00000.00000.0000 0.00000.00000.00000.00000.00000.0000 0.00000.00000.00000.00000.00000.0000 0.00000.00000.00000.000

10、00.00000.0000 0.00000.00000.00000.00000.00000.0000 0.00000.00000.00000.00000.00000.0000 0.00000.00000.00000.00000.00000.0000 0.00000.00000.00000.00000.00000.0000 0.00000.00000.00000.00000.0000同理對于不同的 X0，B，能得到迭代結(jié)果(3)-(9)，在中不再贅述，詳見附件包中數(shù)據(jù)。40 10 3.3333 3.3333 3.5671 3.6157 3.6485 3.6607 3.6669 3.6696 3.

11、6709 3.671500 0.0000 0.5556 0.6481 0.7188 0.7438 0.7567 0.7624 0.7651 0.7663 0.766900 0.0000 0.2778 0.3704 0.4226 0.4449 0.4559 0.4610 0.4635 0.4646 0.465200 0.0000 0.0000 0.0926 0.1196 0.1390 0.1472 0.1515 0.1536 0.1546 0.155000 0.0000 0.0000 0.0231 0.0463 0.0572 0.0637 0.0669 0.0685 0.0693 0.06970

12、0 0.0000 0.0000 0.0000 0.0116 0.0180 0.0222 0.0244 0.0256 0.0263 0.026600 0.0000 0.0000 0.0000 0.0019 0.0058 0.0080 0.0095 0.0103 0.0107 0.011000 0.0000 0.0000 0.0000 0.0000 0.0013 0.0025 0.0033 0.0038 0.0041 0.004300 0.0000 0.0000 0.0000 0.0000 0.0002 0.0007 0.0011 0.0014 0.0016 0.001700 0.0000 0.0

13、000 0.0000 0.0000 0.0000 0.0001 0.0003 0.0005 0.0006 0.000600 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0001 0.0001 0.0002 0.000200 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0001 0.000100 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.000000 0.0000 0.0000 0.0000 0

14、.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.000000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.000000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.000001 0.3333 0.3457 0.3707 0.4008 0.4094 0.4113 0.411700 0.0556 0.0647 0.1702 0.2076 0.2156 0.2174 0.217800 0.03

15、70 0.3196 0.4690 0.4979 0.5045 0.5060 0.506300 0.0108 0.6615 0.7518 0.7727 0.7772 0.7782 0.77840 10 3.3382 3.5940 3.6554 3.6682 3.6711 3.6718 3.672000 0.5573 0.7250 0.7576 0.7652 0.7669 0.7673 0.767400 0.3711 0.4425 0.4602 0.4644 0.4654 0.4656 0.465600 0.1083 0.1438 0.1526 0.1548 0.1553 0.1555 0.155

16、500 0.0490 0.0643 0.0686 0.0697 0.0700 0.0701 0.070100 0.0172 0.0241 0.0261 0.0267 0.0268 0.0269 0.026900 0.0069 0.0099 0.0108 0.0111 0.0112 0.0112 0.011200 0.0026 0.0039 0.0043 0.0044 0.0045 0.0045 0.004500 0.0010 0.0015 0.0017 0.0018 0.0018 0.0018 0.001800 0.0004 0.0006 0.0007 0.0007 0.0007 0.0007

17、 0.000700 0.0001 0.0002 0.0003 0.0003 0.0003 0.0003 0.000300 0.0001 0.0001 0.0001 0.0001 0.0001 0.0001 0.000100 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.000000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.000000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.000000 0.0000 0.0000 0.0000 0.0000 0.0000

18、0.0000 0.0000分析 2：采用不同的初值向量和右端向量，發(fā)現(xiàn)兩種迭代方法均收斂，并看到G-S 迭代法較 JACOBI 迭代法收斂速度快，這與兩者的譜半徑大小有關(guān)。譜半徑小于 1 時，迭代法能夠收斂，且對同組數(shù)據(jù)，譜半徑越小收斂速度越快。結(jié)論：兩種方法均收斂收斂速度 G-S 法更快譜半徑小于 1 時，迭代法收斂對于同組 X0 和 B，譜半徑越小，收斂速度越快（猜測）對于同一個矩陣，G-S 法的譜半徑一定小于 JACOBI 法的譜半徑結(jié)果 3：分析 3：對對角線元素依次乘以 1-10，得到上表，可以發(fā)現(xiàn)譜半徑越來越小，迭代次數(shù)越來越少。結(jié)論：譜半徑小于 1 時，迭代法收斂對于同組 X0

19、和 B，譜半徑越小，收斂速度越快此處是嚴(yán)格對角占優(yōu)情況：對角占優(yōu)趨勢越嚴(yán)重，則譜半徑越小，收斂速度越快53*12345678910R* 100048932447816699612544489sum201298877776組數(shù)123456789迭代次數(shù)/jacobi61013131313141314迭代次數(shù)/g-s479999999參考程序%jacobifunction x y=jacobi(A,b,x0,e,limit) D=diag(diag(A);L=-tril(A,-1);U=-triu(A,1); B=D(L+U);%output radius r=max(eig(B); r%f=Db

20、; n=length(b);temp1=x0;key=false; if abs(r)1sum=0; elsesum=limit; end y(:,1)=temp1;y(:,2)=b;while (not(key)&(sumlimit) temp2=B*temp1+f;if norm(temp1-temp2,inf)e key=true;elsekey=false; end temp1=temp2; sum=sum+1;y(:,sum+2)=temp1; endsum x=temp1;%gaussseideilfunction x y=gaussseideil(A,b,x0,e,limit)

21、D=diag(diag(A);L=-tril(A,-1);U=-triu(A,1);B=(D-L)(U);%output radius r=max(eig(B); r%6f=(D-L)b;n=length(b); temp1=x0; key=false;if abs(r)1 sum=0;elsesum=limit; end y(:,1)=temp1;y(:,2)=b;while (not(key)&(sumlimit) temp2=B*temp1+f;if norm(temp1-temp2,inf)e key=true;elsekey=false; end temp1=temp2; sum=s

22、um+1;y(:,sum+2)=temp1; endsum x=temp1;clear allclc%input matrix A n=input(input the dim : ); A1=sparse(1:n,1:n,3,n,n);A2=sparse(1:n-1,2:n,-1/2,n,n);A3=sparse(1:n-2,3:n,-1/4,n,n); A=A1+A2+A3+A2+A3;AA=full(A);%input result value & original value b=eye(1,n); b(5)=10;b(6:15)=15;b(10:20)=5;b(4*(1:5)=20;b

23、=b x0=zeros(1,n);x0(5)=100;x0(6:15)=15;x0(10:20)=5;x0(4*(1:5)=20; x0=x0;%prosx y1=jacobi(AA,b,x0,0.001,200);y1 pausex y2=gaussseideil(AA,b,x0,0.001,200);y27clear allclcfor j=1:10%input matrix A% n=input(input the dim : ); n=20;A1=sparse(1:n,1:n,3*j,n,n);A2=sparse(1:n-1,2:n,-1/2,n,n);A3=sparse(1:n-2,

24、3:n,-1/4,n,n); A=A1+A2+A3+A2+A3;AA=full(A);%input result value & original value b=eye(1,n); b(5)=10;b(6:15)=15;b(10:20)=5;b(4*(1:5)=20;b=b; x0=zeros(1,n);x0(5)=100;x0(6:15)=15;x0(10:20)=5;x0(4*(1:5)=20; x0=x0;%proe=1e-5;sx y1=jacobi(AA,b,x0,e,200);y1; pauseend5.設(shè)國民經(jīng)濟(jì)由農(nóng)業(yè)、制造業(yè)和服務(wù)業(yè)三個部門，已知某年它們之間的投入產(chǎn)出關(guān)系、外

25、部需求、初始投入等如書上表所示。根據(jù)表回答下列問題：1)如果今年對農(nóng)業(yè)、制造業(yè)和服務(wù)業(yè)的外部需求分別為 50,150,100 億元，問這三個部門的總產(chǎn)出分別應(yīng)為多少？2)如果三個部門的外部需求分別增加 1 個多少？，問它們的總產(chǎn)出應(yīng)分別增加分析設(shè)計模型分析由于投入產(chǎn)出表的編制是受時間限制的，不同時期表中數(shù)字可能不同，因此需要依靠一些相對穩(wěn)定的來分析其中的內(nèi)在規(guī)律。投入產(chǎn)出法中相對穩(wěn)定的最主要是直接消耗系數(shù)，記作 ，表示第 j 部門 1 個的產(chǎn)出對第8i 部門的直接消耗量（均已產(chǎn)值計算），由投入產(chǎn)出表直接計算得到表 A。則建立模型，記 x 為個部門在一定時期內(nèi)總產(chǎn)出，d 為外部需求，則有（I-A

26、）x=d由上式可以得到 x=（I-A）-1*d，表明總產(chǎn)出 x 對外部需求 d 是線性的，所以當(dāng) d 增加 1 個，記作 時，x 的增量為 。若農(nóng)業(yè)外部 增加一個，即 ， 為 的第一列，依此類推即可。編程思路輸入投入產(chǎn)出表，生成直接消耗量表 A利用直接法計算總產(chǎn)出求 A 逆矩陣結(jié)果 1：x0 =(139,268,208)分析 1：所以對應(yīng)外部需求為 50,150,100 億元的情況，農(nóng)業(yè)、制造業(yè)和服務(wù)業(yè)三部門對應(yīng)總產(chǎn)出分別為 139,268,208 億元。結(jié)果 2：dx =1.34590.56340.43820.25041.26760.43040.34430.49301.2167分析 2：農(nóng)業(yè)

27、需求增加 1 個，則三個部門總產(chǎn)出分別增加 1.3459,0.5634和 0.4382 億元。制造業(yè)需求增加 1 個，則三個部門總產(chǎn)出分別增加0.2504,1.2676 和 0.4304 億元。服務(wù)業(yè)需求增加 1 個，則三個部門總產(chǎn)出分別增加0.3443,0.4930 和 1.2167 億元。參考程序clear all;clc;9%input matrix AA=15 20 30; 30 10 45; 20 60 0;s=100 200 150;for i=1:3 A(:,i)=A(:,i)./s(i);end A%b=50 150 100;% for i=1:3%b(i)=input(input requirement: );% end b=bx=round(eye(3)-A)b) pauseclear allclc%input matrix AA=15 20 30; 30 10 45; 20 60 0;s=100 200 150;for i=1:3 A(:,i)=A(:,i)./s(i);end A%prosdx=inv(eye(3)-A);dx7.輸電網(wǎng)絡(luò)：一種大型輸電網(wǎng)絡(luò)可簡化為如書