1、實驗 5：線性代數方程組的數值解法習題3：已知方程組,其中，定義為： 試通過迭代法求解此方程組，認識迭代法收斂的含義以及迭代初值和方程組系數矩陣性質對收斂速度的影響。實驗要求：（1） 選取不同的初始向量x0和不同的方程組右端向量b，給定迭代誤差要求，用雅可比迭代法和高斯-賽德爾迭代法計算，觀測得到的迭代向量序列是否均收斂？若收斂，記錄迭代次數，分析計算結果并得出結論；（2） 取定右端向量b和初始向量x0，將A的主對角線元素成倍的增長若干次，非主對角元素不變，每次用雅可比迭代法計算，要求迭代誤差滿足,比較收斂速度，分析現象并得出結論。1、 程序設計（可直接粘貼運行）1） Jacobi迭代法fun

2、ction y=jacobi(a,b,x0,e,m) %定義jacobi函數，其中：a,b為線性方程組中的矩陣和右端向量；x0為初始值；%e和m分別為人為設定的精度和預計迭代次數；運行結果y為迭代的結果和所有中間值組成的%矩陣y=0;%對y初始化d=diag(diag(a);%按雅可比迭代標準形形式取主對角元素作為矩陣Du=-triu(a,1);%取上三角矩陣ul=-tril(a,-1);%取下三角矩陣lbj=d-1*(l+u);fj=d-1*b;x=x0,zeros(20,m-1);%初始化x,其中x1=x0,即初始值 for k=1:m%人為規定迭代次數，防止不收斂迭代導致死循環x(:,k

3、+1)=bj*x(:,k)+fj;%jacobi迭代 if norm(x(:,k+1)-x(:,k),inf)<e %判斷迭代后是否滿足迭代中止條件： y=x(:,1:k+1);%賦給y所有中間值和迭代結果 sizej=k;%若去掉;號，則輸出迭代次數 break %并結束迭代 end%若不成立，繼續迭代end %以下部分為驗證迭代公式收斂的方法，僅需運行一次即可，因為收斂性完全由A矩陣決定,而A%在本題是固定不變的；通過判斷中B的譜半徑或范數大小（B在jacobi迭代法中%為矩陣bj），即可得知收斂性:e=eig(bj)%輸出全部特征值，若，則收斂n1=norm(bj,1);%計算1-

4、范數n2=norm(bj);%計算2-范數nn=norm(bj,inf)；%計算-范數q=min(n1 n2 nn)%由于譜半徑不超過人以一種范數，所以只要范數的最小值q<1,也可判斷迭代法收斂2） Gauss迭代法：與Jacobi程序結構相同，不再注釋function y=gauss(a,b,x0,e,m)y=0;d=diag(diag(a);u=-triu(a,1);l=-tril(a,-1);x=x0,zeros(20,m-1);bgs=(d-l)-1*u;fgs=(d-l)-1*b; for k=1:mx(:,k+1)=bgs*x(:,k)+fgs; if norm(x(:,k+

5、1)-x(:,k),inf)<ey=x(:,1:k+1); sizeg=k; break endende=eig(bgs)n1=norm(bgs,1);n2=norm(bgs);nn=norm(bgs,inf);min(n1 n2 nn) 3） 操作函數：%構造矩陣An=20;a1=sparse(1:n,1:n,3,n,n);%按稀疏矩陣的輸入法構造，比較方便a2=sparse(1:n-1,2:n,-0.5,n,n);a3=sparse(1:n-2,3:n,-0.25,n,n);a=a1+a2+a3+a2'+a3'a=full(a);%還原為滿矩陣%通過給定不同的初始向量

6、x0或者右端項b，以及規定不同的誤差要求，進行jacobi和gauss%迭代，得到的結果y1、y2位兩種迭代的次數，同時輸出迭代結果，便于分析b=x0=e=m=y1=jacobi(a,b,x0,e,m);y2=gauss(a,b,x0,e,m);4） 改變矩陣A：先對jacobi函數作一定修正，方便分析，命名為jacobi2,如下:function y=jacobi2(a,b,x0,e,m) d=diag(diag(a);u=-triu(a,1);l=-tril(a,-1);bj=d-1*(l+u);fj=d-1*b;x=x0,zeros(20,m-1); n1=norm(bj,1);%計算范

7、數n2=norm(bj);nn=norm(bj,inf);q=min(n1 n2 nn);y(1)=q;%輸出結果1：范數的最小值，判斷收斂速度的方法 for k=1:mx(:,k+1)=bj*x(:,k)+fj; if norm(x(:,k+1)-x(:,k),inf)<e y(2)=k;%輸出結果2：迭代次數 break endend %改變A矩陣主對角元素的值，比較jacobi迭代的收斂速度，即迭代誤差滿足%時的迭代次數b=(1:20)'%取定bx0=20*ones(20,1);%取定x0e=10-5;%取定精度r=20; %設置主對角元素擴大的最大倍數 y=0;m=50;

8、for k=1:r; a1=k*sparse(1:n,1:n,3,n,n);%只需更改A的主對角元素 a=a1+a2+a3+a2'+a3'%重新構造A a=full(a);y(k,1:3)=k,jacobi2(a,b,x0,e,m);endy%輸出對角線擴大倍數最小范數迭代次數2、 運行結果分析1） 不同初值（x0）和右端項（b）條件下的解的情況表一：收斂性判斷1-范數2-范數-范數Jacobi0.01630.01670.01630.0167Gauss0.00080.00840.00840.0084可以看到，矩陣A無論是譜半徑或是任意范數的值都小于1，可知在A不變的情況下，Ja

9、cobi和Gauss法必然收斂。2） b取不同的值，x0=20*ones(20,1), e=10-5, m=50 條件下的情況對比迭代次數B=1:20B=10:10:200B=20:-1:1B=20*ones(20,1)B=2000*ones(20,1)Jacobi2426242330Gauss1617161520根據1）分析的結果，可以證明無論b取任何值，采用兩種方法迭代均收斂，但b的值的變化會影響迭代的次數；且Gauss迭代法總是比Jacobi迭代法收斂速度更快。下表列出的是B=1:20情況下部分結算結果，可以很明顯的看到兩種迭代法的收斂速度不同:JacobiK=1K=2K=3K=4K=5

10、K=6 K=22K=23K=24標準值X15.33332.75001.53941.08740.88580.79820.72470.72470.72470.7247X29.00004.33332.66441.92481.60821.46521.34441.34441.34441.3444X311.00005.80563.71992.78012.36262.17172.00722.00722.00722.0072GaussK=1K=2K=3K=4K=5K=6K=14K=15K=16標準值X15.33332.05401.13540.85390.76570.73770.72470.72470.7247

11、0.7247X26.55562.94321.84591.50331.39491.36051.34441.34441.34441.3444X37.53703.73862.55532.18152.06272.02492.00722.00722.00722.0072由于精度問題，在迭代的最后幾次中從顯示的數位已經不能看出標準值與計算值得差別，但是若采用long顯示設定，就可以看到更多位小數的顯示，其結果符合最初設定的精度e，數據繁瑣，略。3） x0取不同的值，b=20*ones(20,1), e=10-5, m=50 條件下的情況對比迭代次數X0=1:20X0=10:10:200X0=20:-1:1

12、X0=20*ones(20,1)X0=2000*ones(20,1)Jacobi2227222331Gauss1518151520根據1）分析的結果，可以證明無論X0取任何值，采用兩種方法迭代均收斂，但x0的值的變化會影響迭代的次數；且Gauss迭代法總是比Jacobi迭代法收斂速度更快。下表列出的是x0=1:20情況下部分結算結果，可以很明顯的看到兩種迭代法的收斂速度不同:JacobiK=1K=2K=3K=4K=5K=6 K=20K=21K=22標準值X17.25008.62509.22289.45369.55419.59749.63279.63279.63279.6327X27.66679

13、.958310.81411.18311.34111.41111.468311.468311.468311.4683X38.166710.75711.81612.28212.48712.57912.656012.656012.656012.6560GaussK=1K=2K=3K=4K=5K=6K=13K=14K=15標準值X17.25008.93529.42679.57209.61509.62769.63279.63279.63279.6327X28.708310.65411.22811.39811.44811.46311.468311.468311.468311.4683X39.805611.

14、81312.40812.58412.63612.65012.656012.656012.656012.65604） b=(1:20)'x0=20*ones(20,1);e=10-5; m=50;（固定各個參數）r=20;（改變a的主對角元素，依次擴大1倍、2倍20倍）運行結果：主對角元素擴大倍數三個范數的最小值q迭代次數1.00000.489321.00002.00000.244712.00003.00000.16319.00004.00000.12238.00005.00000.09798.00006.00000.08167.00007.00000.06997.00008.00000

15、.06127.00009.00000.05446.000010.00000.04896.000011.00000.04456.000012.00000.04086.000013.00000.03766.000014.00000.03506.000015.00000.03266.000016.00000.03066.000017.00000.02886.000018.00000.02726.000019.00000.02585.000020.00000.02455.0000結果分析：根據判斷迭代是否收斂的原理，若A是嚴格對角占優的，即，則Jacobi迭代和Gauss迭代均收斂。由此推測，A的對角

16、占優程度可能是影響收斂速度的關鍵因素。由上述試驗表面，對A的對角元素數值的擴大的確可以明顯的改善Jacobi迭代的收斂速度，與預測的情況相符。主對角占優程度可以影響收斂速度，但是這個結論因為數學水平有限沒有得到嚴格的推導。但關于A的范數和收斂速度的關系在書上有比較完整的證明，即有公式：（1）q為A矩陣的三種范數中最小的值，x*為原線性方程組的解，可知q越小，序列收斂越快。從上表中可以看出，迭代次數的減小(收斂速度的增加)，是和q的減小同步發生的，公式（1）得到驗證。 還可以看到，隨著擴大倍數的增大，迭代速度的增長率在不斷的下降。可以直觀的理解為,主對角元素的相對優勢在擴大一定倍數之后已經相當明

17、顯，即使再增大若干倍其對于非對角元素仍然是具有壓倒性優勢的，優勢地位不會因此而改變很多，根據前面的假設，收斂速度與對角占優程度有關，速度的改變也將減緩。習題8 種群的繁殖與穩定收獲：種群的數量因繁殖而增加，因自然死亡而減少，對于人工飼養的種群（比如家畜）而言，為了保證穩定的收獲，各個年齡的種群數量應維持不變，種群因雌性個體的繁殖而改變，為方便起見以下種群數量均指其中的雌性。種群年齡記作k=1，2，n，當年年齡k的種群數量記作xk，繁殖率記作bk（每個雌性個體1年繁殖的數量），自然存活率記作sk（sk=1-dk，dk為1年的死亡率），收獲量記作hk，則來年年齡k的種群數量應為，要求各個年齡的種群

18、數量每年維持不變就是要使。(1)若已知bk，sk，給定收獲量hk，建立求各個年齡的穩定種群數量xk的模型（用矩陣、向量表示）。(2)設n=5，b1=b2=b5=0，b3=5，b4=3，s1=s4=0.4，s2=s3=0.6，如果要求h1h5為500，400，200，200，100，求(3)要使h1h5均為500， 如何達到？1.模型的建立根據題目給出的模型和各個參量，建立線性方程組。（1）參數b,s,h的意義如題目所述，向量，表示年齡為k的種群在第i個統計時間段中的數量。題目中還給出了使各年齡的種群數量每年維持不變的要求，即在k取一定值時，不隨i的變化而變化。方程（1）變為如下形式：（2）得到

19、了方程（2）的形式，就可以通過直接或間接（迭代）方法求解線性方程組。2.程序設計1）構造矩陣n=5;b=zeros(1,n);b(3)=5;b(4)=3;s=0.4 0.6 0.6 0.4;h=0 500 400 200 100'ss=diag(s);m=b;ss,zeros(n-1,1)-eye(5);x1=100*ones(5,1); %初值x1，計算guass迭代使用到x=gauss1(m,h,x1); %試通過gauss迭代法得到方程組的解xx=mh; %xx為matlab通過內定方法得到的解，可以認為是精確的 2）高斯函數: 試采用高斯迭代法求解方程，內部設定精度和迭代次數，

20、輸出函迭代的結果和中間值，并輸出全部特征值和范數的最小值function x=gauss1(a,b,x0) d=diag(diag(a);u=-triu(a,1);l=-tril(a,-1);bgs=(d-l)-1*u;fgs=(d-l)-1*b;e=10.-3; %設定精度m=20; %迭代次數x=x0;for k=1:mx(:,k+1)=bgs*x(:,k)+fgs;if norm(x(:,k+1)-x(:,k),inf)<e;x(:,1:k+1);breakendend e=eig(bgs) %輸出全部特征值n1=norm(bgs,1); n2=norm(bgs);nn=norm(

21、bgs,inf);q=min(n1 n2 nn) %輸出范數的最小值3、 運行結果分析1） 問題1的方程已在“模型的建立”中給出2） 將條件n=5，b1=b2=b5=0，b3=5，b4=3，s1=s4=0.4，s2=s3=0.6，h=500,400,200 ,100,100帶入方程（2）,用左除程序和gauss迭代分別求解x1x5運行結果：e = 0 0 0 0 1.6320q = 6.4974x=1.0e+007 *colume 0 1 234 16 17 18 19 20xx = 1.0e+003 * 8.4810 2.8924 1.3354 0.6013 0.1405由特征值的和范數大于

22、1的結果可以分析出，采用gauss迭代法不能收斂，具體運算結果x發散現象很明顯，迭代四次就會出現1e5數量級的值，驗證了上面的推斷。方程的正確結果為xx，即穩定狀態下，種群處于1、2、3、4、5年齡段的個體數目分別為：8481、2892、1335、601、141，數量維持不變。3） 改變h1h5均為500，再次求解x,并討論結果由于矩陣m（2）方程中左端矩陣，沒有改變，可知方程仍然不可以用gauss或者jacobi迭代法進行求解。直接由左除運算得到：xx = 1.0e+004 * 1.1772 0.4209 0.2025 0.0715 -0.0214保持自然存活率和繁殖率不變，在對每個年齡段的

23、人工收獲量都為500的情況下，達到穩定時，種群處于1、2、3、4、5年齡段的個體數目分別為：11772、4209、2025、715、-214，數量維持不變。可以從數據中看到，年齡段5的個體數目出現了負值，這在現實中是不可能出現的。分析：問題1：由方程（2）的建立過程可以看出，影響種群某年齡段的主要因素是上一年齡段的存活率和人工收獲量，兩者對于該年齡段種群數量產生相反的作用，而特殊的，對于年齡段1，即種群中剛剛出的個體數量，其僅僅受到各個種群的繁殖率的影響。出現解為負值的情況在數學上是由方程本身的性質決定的，是正常現象，但聯系本題涉及的實際情況卻反映出了較大的問題：它說明此種人工繁殖和收獲的方法

24、是不能實現的。問題2：不收斂性，由于方程（2）在題目要求的初始條件下其范數和特征根均大于1，即采用迭代法不能收斂于方程的解。這種情況在現實中意味著：若人工給定初始種群數量時，采用不同于方程（2）的解所規定的數量，在給定的繁殖率、存活率、收獲量的情況下，這種人工繁殖的方法不可能趨向于穩定。解決1:關于問題1：解的不符合現實的情況，可以通過改變參數（s、b、h）來解決。下面針對年齡段5出現負值的問題提出一種最簡單的解決辦法：方案：在現有人工飼養條件和種群質量（對應著參數b、s）下，減少對年齡段5的個體的收集量，以保證穩定解的存在。y=0;for k=200:10:500 h(5)=k;%減小h(5

25、),即對年齡段5的人工收集量 k0=(k-190)/10;xx=mh%輸出此時方程的解y(1,k0)=xx(5);%將新的新得到h(5)的值賦予yendkk=500:-10:200;plot(kk,y(1,:),grid;%作圖運行結果： 對年齡段5的收獲量取不同值時方程的部分解： 1.0e+004 *X(1)X(2)X(3)X(4)X(5)H(5)=2001.17720.42090.20250.07150.00863001.17720.42090.20250.0715-0.00144001.17720.42090.20250.0715-0.01145001.17720.42090.20250

26、.0715-0.0214由于年齡段5的數量x（5）僅僅可能對年齡段1通過繁殖率發生影響，而b5為0，所以改變年齡段5的情況不會對其他年齡段書數量造成任何影響。僅減小小對年齡段5的收集量，就可以使方程存在合理的穩定解；由圖形還可以得到，當h(5)大約為286時，方程的解x(5)為零，此值為對年齡段5的最大收獲量。可以看出，本法可以維持其他年齡段的收集兩不變，且不需增加成本改善人工飼養條件和種群質量（改變參數b、s），比較經濟。 解決2：解的的不收斂性，實際上比起1更加嚴重，因為對方程(2)的迭代不能自動收斂，對應著實際種群數量不能自發趨向穩定，這對于人工飼養造成了非常大的麻煩，因為必須保證初始的各個年齡段個體的數量分布就滿足方程2解的形式 ，且由于現實中的情況是不可能嚴格滿足方程（2）的，總會有無法控制的個體數量波動，這又要求對種群數量的隨時監測和調整，