實(shí)驗(yàn)一 解線性方程組的直接法

1、實(shí) 驗(yàn) 報(bào) 告課程名稱 數(shù)值分析 實(shí)驗(yàn)項(xiàng)目 解線性方程組的直接法 專業(yè)班級(jí) 姓 名 學(xué) 號(hào) 指導(dǎo)教師 成 績(jī) 日 期 月 日 一. 實(shí)驗(yàn)?zāi)康?、 掌握程序的錄入和matlab的使用和操作；2、 了解影響線性方程組解的精度的因素方法與問(wèn)題的性態(tài)。3、 學(xué)會(huì)matlab提供的“”的求解線性方程組。二. 實(shí)驗(yàn)要求 1、按照題目要求完成實(shí)驗(yàn)內(nèi)容；2、寫(xiě)出相應(yīng)的matlab 程序；3、給出實(shí)驗(yàn)結(jié)果(可以用表格展示實(shí)驗(yàn)結(jié)果)；4、分析和討論實(shí)驗(yàn)結(jié)果并提出可能的優(yōu)化實(shí)驗(yàn)。5、寫(xiě)出實(shí)驗(yàn)報(bào)告。三. 實(shí)驗(yàn)步驟1、用分解及列主元高斯消去法解線性方程組a)，輸出中系數(shù)分解的矩陣和，解向量和;用列主元法的行交換次序解向

2、量和求；比較兩種方法所得結(jié)果。2、用列主高斯消元法解線性方程組。（1）、（2）、分別輸出，解向量，（1）中的條件數(shù)。分析比較（1）、（2）的計(jì)算結(jié)果3、線性方程組的和分別為,則解. 用matlab內(nèi)部函數(shù)求和的所有特征值和. 若令,求解，輸出向量和，從理論結(jié)果和實(shí)際計(jì)算兩方面分析線性方程組解的相對(duì)誤差以及的相對(duì)誤差的關(guān)系。4、 希爾伯特矩陣，其中，（1）分別對(duì)計(jì)算，分析條件數(shù)作為的函數(shù)如何變化。（2）令，計(jì)算，然后用高斯消去法解線性方程組求出，計(jì)算剩余向量以及。分析當(dāng)增加時(shí)解分量的有效位數(shù)如何隨變化，它與條件數(shù)有何關(guān)系？當(dāng)多大時(shí)連一位有效數(shù)字也沒(méi)有了？將每種情形的兩個(gè)結(jié)果進(jìn)行表格對(duì)比，如：n=

3、6時(shí)：gauss列主消去法求得的的有效數(shù)字四、實(shí)驗(yàn)結(jié)果五、討論分析（對(duì)上述算例的計(jì)算結(jié)果進(jìn)行比較分析，主要說(shuō)清matlab的算符與消去法的適用范圍不同，自己補(bǔ)充）六、改進(jìn)實(shí)驗(yàn)建議（自己補(bǔ)充）1.列主元的高斯消去法 利用列主元的高斯消去法matlab程序源代碼：首先建立一個(gè)gaussmethod.m的文件，用來(lái)實(shí)現(xiàn)列主元的消去方法。function x=gaussmethod(a,b)%高斯列主元消去法，要求系數(shù)矩陣非奇異的， %n = size(a,1);if abs(det(a) a=10 -7 0 1;-3 2.099999 6 2;5 -1 5 -1;2 1 0 2a = 10.0000

4、 -7.0000 0 1.0000 -3.0000 2.1000 6.0000 2.0000 5.0000 -1.0000 5.0000 -1.0000 2.0000 1.0000 0 2.0000 l,u=lu(a)l = 1.0000 0 0 0 -0.3000 -0.0000 1.0000 0 0.5000 1.0000 0 0 0.2000 0.9600 -0.8000 1.0000u = 10.0000 -7.0000 0 1.0000 0 2.5000 5.0000 -1.5000 0 0 6.0000 2.3000 0 0 0 5.0800 b=8 5.900001 5 1b =

5、 8.0000 5.9000 5.00001.0000 y=lby =8.00001.00008.30005.0800 x1=uxx1 =0.0000-1.00001.00001.0000det1= det(a)det1 =-762.00012、（1）在matlab窗口：a=3.01 6.03 1.99;1.27 4.16 -1.23;0.987 -4.81 9.34a = 3.0100 6.0300 1.9900 1.2700 4.1600 -1.2300 0.9870 -4.8100 9.3400 b=1 1 1b = 1 1 1x1,det1,index=gauss (a,b)x1 =

6、1.0e+03 * 1592.599624841381 -631.9113762025488 -493.6177247593899det1 = -0.0305index = 1 (2) 在matlab窗口： a=3.00 6.03 1.99;1.27 4.16 -1.23;0.990 -4.81 9.34a = 3.0000 6.0300 1.9900 1.2700 4.1600 -1.2300 0.9900 -4.8100 9.3400 b=1 1 1b = 1 1 1x2,det2,index=gauss5555(a,b)x2 = 119.5273 -47.1426 -36.8403det

7、2 = -0.4070index = 13、在matlab窗口：a=10 7 8 7;7 5 6 5;8 6 10 9;7 5 9 10;b=32 23 33 31;x=abb1=32.1 22.9 33.1 30.9;x1=ab1a1=10 7 8.1 7.2;7.08 5.04 6 5;8 5.98 9.89 9;6.99 5 9 9.98;x2=a1bdelta_b=norm(b-b1)/norm(b)delta_a=norm(a-a1)/norm(a)delta_x1=norm(x-x1)/norm(x)delta_x2=norm(x-x2)/norm(x)cond_a=cond(a)

8、x = 1.0000 1.0000 1.0000 1.0000x1 = 9.2000 -12.6000 4.5000 -1.1000x2 = -9.5863 18.3741 -3.2258 3.5240delta_b = 0.0033delta_a = 0.0076delta_x1 = 8.1985delta_x2 = 10.4661cond_a = 2.9841e+033、在matlab窗口： a=10 7 8 7;7 5 6 5;8 6 10 9;7 5 9 10; b=32 23 33 31; x=abb1=32.1 22.9 33.1 30.9;x1=ab1a1=10 7 8.1 7.

9、2;7.08 5.04 6 5;8 5.98 9.89 9;6.99 5 9 9.98;x2=a1bdelta_b=norm(b-b1)/norm(b)delta_a=norm(a-a1)/norm(a)delta_x1=norm(x-x1)/norm(x)delta_x2=norm(x-x2)/norm(x)cond_a=cond(a)x = 1.0000 1.0000 1.0000 1.0000x1 = 9.2000 -12.6000 4.5000 -1.1000x2 = -9.5863 18.3741 -3.2258 3.5240delta_b = 0.0033delta_a = 0.0

10、076delta_x1 = 8.1985delta_x2 = 10.4661cond_a = 2.9841e+034、k=13for n=2:6 a=hilb(n); co(n)=cond(a,inf);endx=1:6;plot(x,co);b=zeros(k);x11=b;x0=b;r=b;for i=2:k x=(linspace(1,1,i); x0(1:i,(i-1)=x; h=hilb(i); b0=h*x; b(1:i,(i-1)=b0; x1=gauss2(h,b0); r(1:i,(i-1)=b0-h*x1; x11(1:i,(i-1)=x1;enddx=x11-x0;結(jié)果如

11、下：co=0,27, 748,28375,943656, 29070279可見(jiàn)，條件數(shù)隨著n的增大，急劇增加，如圖1所示。將（2）求得的結(jié)果（dx，x11）整理得n=2時(shí)：xrn有效數(shù)字14.44e-160161-6.66e-16015n=3時(shí)：xrn有效數(shù)字1-1.33e-15-2.22e-161519.55e-150141-9.99e-15014n=4時(shí)：xrn有效數(shù)字1-2.35e-1401412.56e-130131-6.11e-131.11e-161213.96e-13013n=5時(shí)：xrn有效數(shù)字1-1.21e-144.44e-161416.97e-1401311.18e-132.

12、22e-16131-6.17e-1301214.59e-13-1.11e-1613n=6時(shí)：xrn有效數(shù)字1-9.28e-1301212.67e-110111-1.82e-1001014.75e-100101-5.26e-100912.08e-10010n=7時(shí)：xrn有效數(shù)字1-9.26e-1201113.71e-100101-3.59e-092.22e-1691.000000011.40e-08080.99999997-2.57e-08081.000000022.22e-081.11e-1680.99999999-7.30e-0908n=8時(shí)：xrn有效數(shù)字1-2.82e-114.44e-

13、161111.53e-09090.99999998-2.01e-082.22e-1681.000000111.09e-07070.9999997-2.95e-07-2.22e-1671.000000424.19e-07-1.11e-1670.9999997-2.99e-071.11e-1671.000000088.44e-0807n=9時(shí)：xrn有效數(shù)字1-2.75e-104.44e-16101.000000021.90e-08080.99999968-3.20e-07071.000002282.28e-064.44e-1660.99999164-8.36e-06051.000017061.7

14、1e-05050.99998042-1.96e-05-1.11e-1651.000011821.18e-05050.99999708-2.92e-0606n=10時(shí)：xrn有效數(shù)字1-9.05e-104.44e-1691.000000087.76e-08070.99999835-1.65e-06061.000014911.49e-054.44e-1650.99992911-7.09e-05041.000194430.000194426040.99968164-0.000318365-1.11e-1641.000307150.000307149040.99983898-0.00016101904

15、1.000035373.54e-0505n=11時(shí)：xrn有效數(shù)字0.999999995-5.16e-09081.0000005395.39e-07-4.44e-1660.999986041-1.40e-05051.0001558750.00015588040.999072325-0.0009277031.0032587180.00325872030.992910161-0.0070898021.0096588070.00965881-2.22e-1620.991981908-0.0080181031.0037076540.00370765030.999267974-0.000732-2.22

16、e-163n=12時(shí)：xrn有效數(shù)字0.999999965-3.49e-088.88e-1681.0000044094.41e-064.44e-1660.999861744-0.0001383041.0018797110.00187971030.986241983-0.013758021.0603759740.060375972.22e-1620.831939173-0.16806082.22e-1611.3039733630.30397336010.643862006-0.3561381.11e-1611.2606840180.260684022.22e-1610.891666924-0.1

17、0833311.11e-1611.0195107530.0195107502n=13時(shí)：xrn有效數(shù)字1.0000000696.89e-088.88e-1670.999989224-1.08e-05-4.44e-1651.0004135380.00041354-2.22e-1640.993147495-0.00685252.22e-1631.0612462080.06124621-2.22e-1620.669261244-0.33073882.22e-1612.1490741311.1490741300-1.65417119-2.6541712005.118548084.11854808-1.

18、11e-160-3.243049939-4.24304991.11e-1603.7830048962.783004900-0.051792817-1.05179281.11e-1601.1743291170.174329121.11e-161五、分析討論： 實(shí)驗(yàn)的數(shù)學(xué)原理很容易理解，也容易上手。把運(yùn)算的結(jié)果帶入原方程組，可以發(fā)現(xiàn)符合的還是比較好。這說(shuō)明列主元消去法計(jì)算這類方程的有效性。當(dāng)a可逆時(shí)，能夠?qū)⒂?jì)算進(jìn)行到底，列主元法就能確保算法的穩(wěn)定，而且計(jì)算量不大。直接三角消去過(guò)程，實(shí)質(zhì)上是將a分解為兩個(gè)三角矩陣的乘積a=lu，并求解ly=b的過(guò)程。回帶過(guò)程就是求解上三角方程組ux=y。所以在實(shí)際的

19、運(yùn)算中，矩陣l和u可以直接計(jì)算出，而不需要任何中間步驟，從而在計(jì)算過(guò)程中將高斯消去法的步驟進(jìn)行了進(jìn)一步的簡(jiǎn)略，大大提高了運(yùn)算速度，這就是三角分解法。通過(guò)以上的計(jì)算比較，2.題方程組具有嚴(yán)重的病態(tài)性。當(dāng)系數(shù)矩陣有微小的變化時(shí)，wucha =-401.8918 159.5435 124.6330，所得的解與原方程組的解有很大的相對(duì)誤差。1題方程組中當(dāng)系數(shù)矩陣a和b有微小變化時(shí)，wucha =0 0 0 0，所得的解與方程組的解沒(méi)有相對(duì)誤差。所以1題方程組是良性的。用matlab內(nèi)部函數(shù)inv通過(guò)求逆矩陣，然后通過(guò)x=inv（a）*b也可以求出方程組的解，但是沒(méi)有列主元高斯消去法具有良好的穩(wěn)定性。det函數(shù)求方程組系數(shù)矩陣的行列式時(shí)所得結(jié)果和高斯消去法和三角法所得結(jié)果相同，具有方便快捷的優(yōu)點(diǎn)。題四可以看出，條件數(shù)越大，有效位數(shù)越少，當(dāng)n=13時(shí)，出現(xiàn)有效位數(shù)為0的情況。附：高斯列主消去法源程序代碼function x,det,index=gauss(a,b)% gauss% a - 方程組矩陣% b - 方程組右端% x - 方程組的解% det -方程組行列式% index - index=0表示求解