計算方法矩陣直接法與迭代法

1、需要兩份東西1、實驗程序：輸入、輸出、注釋2、實驗報告：問題描述、方法描述、方案設計、結 果分析、結論謝謝，麻煩寫的詳細些實驗五解線性方程組的直接方法實驗5.1 (主元的選取與算法的穩(wěn)定性)問題提出：Gauss消去法是我們在線性代數中已經熟悉的。但由于計算機的數值運算是在一個有限的浮點數集合上進行的，如何才能確保Gauss消去法作為 數值算法的穩(wěn)定性呢？ Gauss消去法從理論算法到數值算法，其關鍵是主元的選 擇。主元的選擇從數學理論上看起來平凡，它卻是數值分析中十分典型的問題。實驗內容：考慮線性方程組Ax =b, A Rn n,b Rn編制一個能自動選取主元，又能手動選取主元的求解線性方程組

2、的Gauss消去過程。實驗要求：一715,b =:,則方程有解11561.14TX =（1,11，1）T。一6 18 61(1)取矩陣A二|三三|86I8取n=10計算矩陣的條件數。讓程序自動選取主元，結果如何？(2)現選擇程序中手動選取主元的功能。每步消去過程總選取按模最小或 按模盡可能小的元素作為主元，觀察并記錄計算結果。若每步消去過程總選取按 模最大的元素作為主元，結果又如何？分析實驗的結果。(3)取矩陣階數n=20或者更大，重復上述實驗過程，觀察記錄并分析不同的問題及消去過程中選擇不同的主元時計算結果的差異，說明主元素的選取在消去過程中的作用。（ 4）選取其他你感興趣的問題或者隨機生成

3、矩陣，計算其條件數。重復上述實驗，觀察記錄并分析實驗結果。1）首先編寫gauss消元法代碼；輸入為矩陣A,向量b,精度ptol,輸出為方程 組的解x。程序如下function out = gausle（ A,b,flag,ptol ）%UNTITLED 高斯消元法可選擇的主元消去法% A n x n 矩陣% b n x 1% flag 標志主元為自動選取還是手動選取，%0，自動（對角為主元）1，選取最小模為主元，%2，選取模最大，（誤差最小）3 ，次最小的模為主元% ptol 精度% out 輸出值，為nx1 的解 Ax=b 的解%if nargin<4,ptol=50*eps;endm

4、,n=size（A）;if m=n,error（'A 不是方陣'）;endout=zeros（n,1）;% 預先設定解的維數nb=n+1;Ab=A,b;%擴維矩陣% disp('開始用擴維陣計算)；% disp(Ab);RA=rank(A);RB=rank(Ab);zhica=RB-RA;if zhica>0,disp('請注意：因為RA=RB,所以此方程組無解.') returnend% 消元過程for i=1:n-1 ;i;if flag=0;pivot=Ab(i,i);ri=i;elseif flag=1;% 按照每列模最小的選取pivot,

5、min_index=min(abs(A(i:n,i);ri=i+min_index-1;elseif flag=2;% 按照模最大的選取pivot,max_index=max(abs(A(i:n,i);ri=i+max_index-1;elseif flag=3%方程最小非0 數tA=A;pivot,min_index=min(abs(tA(i:n,i);while pivot=0;tA(min_index+i-1,i)=inf;pivot,min_index=min(abs(tA(i:n,i);endri=i+min_index-1;endif (pivot=0)|(pivot<pto

6、l) ;warning('系數矩陣奇異！ ！)；return;endif ri=i; % 交換行tmp=A(i,:);A(i,:)=A(ri,:);A(ri,:)=tmp;t=b(i);b(i)=b(ri);b(ri)=t;endfor kk=i+1:nL(kk,i)=A(kk,i)/A(i,i);A(kk,i+1:n)=A(kk,i+1:n)-L(kk,i)*A(i,i+1:n);b(kk)=b(kk)-L(kk,i)*b(i);endendif A(n,n)=0warning('系數矩陣奇異！ ！)；return;end% % 回代求解x=zeros(n,1)'fo

7、r k=n:-1:1%k%A(k,k+1:n)%x(k+1:n)if k=nx(n)=b(n)/A(n,n);elsex(k)=(b(k)-sum( A(k,k+1:n).*x(k+1:n) ) )/(A(k,k); endendout=x'end主程序為N=10時，方程解為1.00001.00001.00001.00001.00001.00001.00001.00001.00001.0000A的條件數如下：1-條件數 2.5575e+0032-條件數 1.7276e+003無窮條件數2.5575e+003代碼在 里，運行計算即可。Flag=0；,時為自動選取2) Flag=1,為選取

8、最小模的值為主元，flag=2,為選取最大模為主元，代碼是一樣的，修改flag值和維數n即可。3) N=20，只需要更改n,重復上述步驟即可。4)-nIiXon 1i=0XonnXizx1i zz0n,X2，b =n3zX2i=0anxn 1n1zx：_i=0思考題一 ：(Vadermond即陣)設1XoX2d21X1X1A = 1x2x2a鼻鼻JXnX2其中，xk =1 +0.1k,k =0,1,，n ,(1)對n=2,5,8,計算A的條件數；隨n增大，矩陣性態(tài)如何變化？(2)對n=20,解方程組Ax=b;設A的最后一個元素有擾動10-4,再求解Ax=b(3)計算(2)擾動相對誤差與解的相對

9、偏差，分析它們與條件數的關系。(4)你能由此解釋為什么不用插值函數存在定理直接求插值函數而要用拉格朗日或牛頓插值法的原因嗎？(5)嘗試做出范德蒙矩陣的預優(yōu)。答1) n=2, 5, 8 ; A的條件數為如下所示，隨n增大，條件數越大，矩陣的 病態(tài)性越嚴重。n =2,cond1(A) =16060;cond2(A)=11044; cond 二(A) =16016;n =5,cond1(A) =1.0522 107; cond2(A)=0.5594 107;cond:;( A) =1.0668 107;n =8;cond1(A) =5.0565 1010; cond2(A)=2.2739 1010;

10、cond 二(A) =4.45788 1010;代碼為 sikaoti_111.m2)無擾動時x =1，,1T,x為20父1的向量；有擾動時，x = 1.0e+006 *0.0027-0.02700.1250-0.35610.6956-0.98381.0344-0.81320.4692-0.18470.03570.0093-0.01040.0043-0.00100.00010.0000-0.00000.00000.00000.0000,代碼文件為sikaoti _112.3)由此看來，當系數矩陣的條件數越大，則病態(tài)性越嚴重，因而系數矩陣很小 相對誤差也會讓解產生很大的相對誤差。4)插值函數存在

11、定理直接求插值函數而要用拉格朗日或牛頓插值法的原因是為 了補償計算中所產生偏差，而且迭代次數越多，偏差就越大。5)預優(yōu)，相關MATLAB函數提示：zeros(m,n) 生成m行，n列的零矩陣ones(m,n) 生成m行，n列的元素全為1的矩陣eye(n)生成n階單位矩陣rand(m,n) 生成m行,n列(0,1)上均勻分布的隨機矩陣diag(x)返回由向量x的元素構成的對角矩陣tril(A)提取矩陣A的下三角部分生成下三角矩陣triu(A)提取矩陣A的上三角部分生成上三角矩陣rank(A)返回貨!陣A的秩det(A)返回方陣A的行列式inv(A)返回可逆方陣A的逆矩陣V , D=eig(A)

12、返回方陣A的特征值和特征向量norm(A,p)矩陣或向量的p范數cond(A,p)矩陣的條件數L , U ,巳=lu(A)選列主元LU分解R=chol(X)平方根分解Hi=hilb(n) 生成 n 階 Hilbert 矩陣實驗六 解線性方程組的迭代法實驗6.1 (病態(tài)的線性方程組的求解)問題提出：理論的分析表明，求解病態(tài)的線性方程組是困難的。實際情況 是否如此，會出現怎樣的現象呢？實驗內容：考慮方程組Hx=b的求解，其中系數矩陣H為Hilbert矩陣，.1. . _H =(hi,j)nn, hi,j, i, j =1,2, ,ni j -1這是一個著名的病態(tài)問題。通過首先給定解(例如取為各個分

13、量均為1)再 計算出右端b的辦法給出確定的問題。實驗要求：(1)選擇問題的維數為6,分別用Gauss消去法、列主元Gauss消去法、J 迭代法、GS迭代法和SOR迭代法求解方程組，其各自的結果如何？將計算結果 與問題的解比較，結論如何？(2)逐步增大問題的維數(至少到100),仍然用上述的方法來解它們，計算的結果如何？計算的結果說明了什么？（ 3）討論病態(tài)問題求解的算法答：1) Guass消元法的程序在gausle.m程序里，如下所示。Flag=0,是自 動消元，Guass消去法的結果為：x=1 ，1，1，1，1，1flag=2為列主元Gauss消去法法，Gauss列主元消去法的結果為：x=1

14、 ，1，1，1，1，1Jacobi迭代法程序在jacobi.m,結果為：Inf ,Inf,NaN,NaN,NaN,NaN 。 也就是迭代發(fā)散，因為迭代陣B 的譜半徑大于1,為4.3085所以此方程不適合用Jacobi方法迭GS迭代法在gauseidel.m,其結果為：0.99921.01250.95901.03061.02670.9715迭代次數為2016SOR 超級松弛法程序在SOR.m 中， w 選擇為 1.89其結果為：0.99911.01800.91361.15860.88001.0304迭代次數為517 次。主程序如下。% 實驗 6.1% 實驗 6.1(病態(tài)的線性方程組的求解)clc

15、;clear;close all;n=6; % 矩陣維數A=hilb(n);eps=1e-5;R=max(abs(eig(A) ;%disp('A 的譜半徑為',num2str(R);b=sum(A,2);x_gauss=gausle(A,b,0)x_gauss_2=gausle(A,b,2)% 編寫 jacobi.m 和 gauseidel.mx_jaco_3,n3=jacobi(A,b,zeros(size(b),eps ) %( A,b,omega,x0,eps,N)x_gs_4,n4=gauseidel(A,b,zeros(size(b),eps )% SOR 法ome

16、ga=1.89;x_SOR_5,n5=SOR(A,b,omega,zeros(size(b),eps)2) 將維數增加為100，高斯消元法，結果為1.00001.00000.99931.02280.67543.5051-10.564434.3970-60.724877.6571-77.081786.1987-72.494316.549923.9588-7.319232.2195-71.5173-9.4683102.1479-4.3990-49.2259-65.428477.29858.937810.7393-2.8585-94.704818.5097112.925253.2932-77.717

17、8-45.60030.8528-20.4096-56.1519106.6452-3.433045.7688 8.1846-62.144310.476826.8102-75.727076.4337-55.6616-22.277789.8267-63.1822 111.9219 -73.864928.5538-22.5285-44.4054-20.734745.2838-39.7352-57.6491 101.373957.6956 47.3377105.3975-47.4092-159.0873-38.1213-73.969618.9134-69.1624 118.195726.6726141.

18、973669.7617-144.1574-13.454667.8139-44.4608-126.1802-39.0300 192.4826 -67.486416.6392 82.7630-146.8443-90.7014164.5643126.838711.1138-223.7063-6.505967.2569-139.3375216.16062.0236-98.0808114.9628-3.8427-55.4881-100.470398.0132-7.2417列主元高斯消元法，結果為0,即出現計算錯誤，彈出，系數矩陣奇異，jacobi B的譜半徑為86.3375, Jacobi迭代法發(fā)散，結

19、果為NaN ,，NaNGauss-S 方法 0.99841.02480.92471.02731.05281.03211.00240.97980.96830.96620.97060.97870.98820.99761.00591.01271.01771.02111.02281.02321.02241.02071.01831.01541.01221.00891.00551.00220.99900.99610.99340.99100.98890.98720.98570.98460.98380.98330.98310.98320.98340.98390.98460.98550.98650.98760.

20、98890.99020.99160.99310.99460.99610.99760.99911.00061.00201.00341.00481.00611.00731.00841.00951.01041.01131.01211.01271.01331.01371.01411.01431.01451.01451.01441.01421.01391.01341.01291.01221.01151.01061.00961.00851.00731.00611.00471.00321.00160.99990.99820.99630.99440.99240.99030.98810.98580.98350.

21、98110.97860.97600.9734SOR 迭代法結果為0.99781.03910.84921.18200.90151.15990.85851.10600.84561.08420.8658 1.0830 0.89791.08670.92811.0874 0.9510 1.08310.96601.0743 0.97441.06270.97811.0499 0.97881.0372 0.9778 1.02540.97631.0152 0.97481.00670.97381.00010.9735 0.9953 0.9740 0.9922 0.97530.99040.9772 0.9900 0

22、.97960.99050.98260.99180.98580.9936 0.9892 0.99590.9927 0.9984 0.99611.00100.99941.00361.00251.00611.00531.00831.00781.01031.01001.01201.01171.01331.01301.01421.01381.01471.01421.01471.01421.01431.01371.01351.01271.01231.01131.01061.00941.00851.00721.00601.00451.00311.00140.99990.99800.99620.99420.9

23、9230.99010.98800.98570.98340.98100.97850.97590.9733條件數來衡量問題的病態(tài)程度。條件數越大，病態(tài)可能越嚴重實驗6.2書上P211計算實習題2,其中N=100,第二小問改為用Jacobi迭 代、G-S迭代、紅黑排序的G-S迭代求解，并比較他們之間的收斂速度；進一步， 用BSOR迭代求解，試找出最優(yōu)松弛因子。習題見圖片- x=Xi=ih,y=yj = jh,解：j ,由此，可知，x萬向上0,1上分為了 11個小格，y萬向也i, j =0,1，,N,N 1一樣。五點差分方程為11rrui.j 1 " "TTui j,j hh114

24、 一、了ui1,j 一聲ui,j 1 產ui,j = f(Xi,yj),即-121ui,j 1 121ui121ui,j 1 - 121ui 1,j484ui,j = f (xi, yj)Jacobi結果為jacobi B的譜半徑為 0.95949 迭代次數=221y_jacobi =1.01391.02781.04151.05501.06811.08041.09141.10011.10531.10491.02781.05531.08251.10931.13521.15961.18131.19891.21011.21171.04151.08251.12311.16301.20171.23841.27141.29881.31761.32411.05501.10931.16301.21601.26771.31721.36251.40131.42991.44411.06811.13521.20171.26771.33281.39581.45471.50671.54771.57271.08041.15961.23841.31721.39581.47311.54711.61471.67121.71061.09141.18131.27141.36251.45471.54711.63791.72361.79931.85771.10011.1