數(shù)值分析實(shí)驗(yàn)報(bào)告

1、湖北民族學(xué)院理學(xué)院數(shù)值分析課程實(shí)驗(yàn)報(bào)告（一）課程名稱數(shù)值分析班級(jí)0212404實(shí)驗(yàn)日期2014/12/15姓名向義俊學(xué)號(hào)021240418實(shí)驗(yàn)成績實(shí)驗(yàn)名稱實(shí)驗(yàn)一算法的穩(wěn)定性分析與誤差分析（驗(yàn)證性實(shí)驗(yàn)）實(shí)驗(yàn)?zāi)康募耙?、 理解算法的數(shù)值穩(wěn)定性概念2、 會(huì)進(jìn)行算法的舍入誤差和截?cái)嗾`差分析實(shí)驗(yàn)環(huán)境matlab軟件實(shí)驗(yàn)內(nèi)容設(shè)計(jì)兩個(gè)算法計(jì)算下面的定積分并估計(jì)誤差。要求兩個(gè)算法理論正確，并且一個(gè)算法是不穩(wěn)定的，另外一個(gè)算法是數(shù)值穩(wěn)定的。算法描述及實(shí)驗(yàn)步驟算法描述：由于初始數(shù)據(jù)誤差在計(jì)算中傳播使計(jì)算結(jié)果誤差增長很快，就數(shù)值不穩(wěn)定。實(shí)驗(yàn)步驟：由分部積分可得計(jì)算的遞推公式 （1）若計(jì)算出，代入（1）式，可逐次求

2、出 的值。要算出就要先算出，若用泰勒多項(xiàng)式展開部分和 并取k=19,用4位小數(shù)計(jì)算，則得，截?cái)嗾`差.計(jì)算過程中小數(shù)點(diǎn)后第5位的數(shù)字按四舍五入原則舍入，由此產(chǎn)生的舍入誤差這里先不討論。當(dāng)初值取為時(shí)，用（1）式遞推的計(jì)算公式為，n=1,2，。計(jì)算結(jié)果見表1的列。用近似產(chǎn)生的誤差就是初值誤差，它對(duì)后面計(jì)算結(jié)果是有影響的.從表1中看到出現(xiàn)負(fù)值，這與一切相矛盾。實(shí)際上，由積分估值得 （2）因此，當(dāng)n較大時(shí)，用近似顯然是不正確的。這里計(jì)算公式與每步計(jì)算都是正確的，那么是什么原因合計(jì)算結(jié)果出現(xiàn)錯(cuò)誤呢？主要就是初值有誤差，由此引起以后各步計(jì)算的誤差滿足關(guān)系由此容易推得，這說明有誤差，則就是的n!倍誤差。例如，

3、n=19，若，則。這就說明完全不能近似了。它表明計(jì)算公式（a）是數(shù)值不穩(wěn)定的。我們現(xiàn)在換一種計(jì)算方案。由（2）式取n=19,得，我們粗略取，然后將公式（1）倒過來算，即由算出，公式為計(jì)算結(jié)果見表1的列。我們發(fā)現(xiàn)與的誤差不超過。記，則，比縮小了n!倍，因此，盡管較大，但由于誤差逐步縮小，故可用近似。反之，當(dāng)用方案（a）計(jì)算時(shí)，盡管初值相當(dāng)準(zhǔn)確，由于誤差傳播是逐步擴(kuò)大的，因而計(jì)算結(jié)果不可靠。此例說明，數(shù)值不穩(wěn)定的算法是不能使用的。 實(shí)驗(yàn)結(jié)果和分析總結(jié)實(shí)驗(yàn)一：實(shí)驗(yàn)二：方案一盡管相當(dāng)準(zhǔn)確，由于誤差傳播使逐步擴(kuò)大的，因而計(jì)算結(jié)果不可靠。方案二比縮小了倍，盡管較大，但誤差是逐步縮小的，故可靠。附錄程序算法

4、一i(1)=0.6321;for i=1:9 i(i+1)=1-i*i(i);endi算法二i(9)=0.0684;for i=9:-1:2 i(i-1)=1/i*(1-i(i);endi湖北民族學(xué)院理學(xué)院數(shù)值分析課程實(shí)驗(yàn)報(bào)告（二）課程名稱數(shù)值分析班級(jí)0212404實(shí)驗(yàn)日期2014/12/15姓名向義俊學(xué)號(hào)021240418實(shí)驗(yàn)成績實(shí)驗(yàn)名稱實(shí)驗(yàn)二 多項(xiàng)式插值（驗(yàn)證性實(shí)驗(yàn)）實(shí)驗(yàn)?zāi)康募耙?、 掌握拉格郎日插值多項(xiàng)式的用法，適用范圍及精確度。2、 掌握牛頓插值多項(xiàng)式的用法，適用范圍及精確度。3、掌握分段線性插值與三次樣條插值的基本原理，會(huì)利用matlab庫函數(shù)進(jìn)行分段線性插值與三次樣條插值。實(shí)驗(yàn)環(huán)

5、境matlab軟件實(shí)驗(yàn)內(nèi)容已知函數(shù)在下列個(gè)點(diǎn)的值為1進(jìn)行l(wèi)agrange插值，newton插值，畫出插值多項(xiàng)式的圖形。0.20.40.60.81.00.980.920.810.640.382. 進(jìn)行分段線性插值與三次樣條插值，畫出插值多項(xiàng)式的圖形。3. 采用四種插值方法計(jì)算處的近似值。算法描述及實(shí)驗(yàn)步驟算法分析：lagrange插值計(jì)算：lagrange插值函數(shù)newton插值公式：分段插值：簡(jiǎn)單地說，將每兩個(gè)相鄰的節(jié)點(diǎn)用直線連起來，如此形成的一條折線就是分段線性插值函數(shù)，記作它滿足 ，且在每個(gè)小區(qū)間上是線性函數(shù) ?？梢员硎緸?實(shí)驗(yàn)步驟：由于中含有3+n個(gè)待定系數(shù)，故應(yīng)需要3+n個(gè)插值條件,已

6、知插值節(jié)點(diǎn)和相應(yīng)的函數(shù)值,這里提供了1+1個(gè)條件，還需要2個(gè)邊界條件。實(shí)驗(yàn)結(jié)果和分析總結(jié)分析總結(jié)：拉格朗日插值多項(xiàng)式，公式結(jié)構(gòu)緊湊，相對(duì)于理論分析比較重要。當(dāng)插值節(jié)點(diǎn)增加或減少時(shí)，計(jì)算要全部重新開始，計(jì)算量大；其插值曲線光滑，誤差估有表達(dá)式。但當(dāng)節(jié)點(diǎn)再次增加時(shí)，它就會(huì)表現(xiàn)出runge現(xiàn)象。牛頓插值多項(xiàng)式，具有一定的繼承性，計(jì)算方便。分段低次插值分段hermit插值三次樣條插值都避免了runge現(xiàn)象的產(chǎn)生。但分段低次插值函數(shù)的導(dǎo)數(shù)是間斷的；分段hermit插值比分段低次插值效果較好。但它要給出節(jié)點(diǎn)上的導(dǎo)數(shù)值，要的信息太多，并且也不太光滑，即，光滑度不高。所以三樣條插值就避免了以上的缺點(diǎn)。附錄程序

7、1.1、lagrange插值function y=lagrange(x0,y0,x)n=length(x0);m=length(x);for i=1:m z=x(i); s=0.0; for k=1:n p=1.0; for j=1:n if j=k p=p*(z-x0(j)/(x0(k)-x0(j); end end s=p*y0(k)+s; end y(i)=s;end在運(yùn)行窗口輸入x0=0.2:0.2:1;y0=0.98 0.92 0.81 0.64 0.38;x=0.2:0.05:1;y=lagrange(x0,y0,x)plot(x0,y0,r,x,y)1.2newton插值func

8、tion f = newton(x,y,x0)syms t; if(length(x) = length(y) n = length(x); c(1:n) = 0.0;else return;end f = y(1);y1 = 0;l = 1; for(i=1:n-1) for(j=i+1:n) y1(j) = (y(j)-y(i)/(x(j)-x(i); end c(i) = y1(i+1); l = l*(t-x(i); f = f + c(i)*l; simplify(f); y = y1; if(i=n-1) if(nargin = 3) f = subs(f,t,x0); else

9、f = collect(f); f = vpa(f, 6); end endend在運(yùn)行窗口輸入x0=0.2:0.2:1; y0=0.98 0.92 0.81 0.64 0.38;x=0.2:0.05:1;y=newton(x0,y0,x)plot(x0,y0,r,x,y)1.3分段線性插值 在運(yùn)行窗口輸入x0=0.2:0.2:1;y0=0.98 0.92 0.81 0.64 0.38;x=0.2:0.05:1;y=interp1(x0,y0,x)plot(x0,y0,r,x,y)title(分段線性插值);1.4三次樣條插值在運(yùn)行窗口輸入x0=0.2:0.2:1;y0=0.98 0.92 0

10、.81 0.64 0.38;x=0.2:0.05:1;y=spline(x0,y0,x)plot(x0,y0,r,x,y)title(三次樣條插值); 湖北民族學(xué)院理學(xué)院數(shù)值分析課程實(shí)驗(yàn)報(bào)告（五）課程名稱數(shù)值分析班級(jí)0212404實(shí)驗(yàn)日期2012/12/23姓名向義俊學(xué)號(hào)021240418實(shí)驗(yàn)成績實(shí)驗(yàn)名稱實(shí)驗(yàn)四 線性方程組的直接法（驗(yàn)證性實(shí)驗(yàn)）實(shí)驗(yàn)?zāi)康募耙?、掌握高斯消去法求解線性方程組2、掌握lu分解法求解線性方程組實(shí)驗(yàn)環(huán)境matlab軟件實(shí)驗(yàn)內(nèi)容編制程序，用掌握列主元高斯消去法和lu分解法求解線性方程組，其中 ， 算法描述及實(shí)驗(yàn)步驟（1）列主元消去法： 先編寫m文件liezhu.m，然后

11、在命令窗口中輸入a=10 -7 0 1;-3 2.099999 6 2;5 -1 5 -1;2 10 0 2; b=8 5.900001 5 1;ra,rb,n,x=liezhu(a,b) 運(yùn)行程序，即可得結(jié)果（2）lu分解法： 先將a做lu分解，用l,u=lu(a)得出l和u。由于a=l*u,令ux=y,則有y=lb,x=uy,即可求解。實(shí)驗(yàn)結(jié)果和分析總結(jié)（1） 列住元高斯消去法： x = 0.4436 -0.2923 0.80151.5178（2）lu分解法：l = 1.0000 0 0 0 -0.3000 -0.0000 1.0000 0 0.5000 0.2193 0.8333 1.0

12、000 0.2000 1.0000 0 0u = 10.0000 -7.0000 0 1.0000 0 11.4000 0 1.8000 0 0 6.0000 2.3000 0 0 0 -3.8114x = 0.4436 -0.2923 0.8015 1.5178附錄程序 （1）列主元高斯消去法:function ra,rb,n,x=liezhu(a,b)b=a b; n=length(b); ra=rank(a); rb=rank(b);zhica=rb-ra;if zhica0,disp(請(qǐng)注意：因?yàn)閞a=rb，所以此方程組無解.)returnendif ra=rb if ra=ndisp

13、(請(qǐng)注意：因?yàn)閞a=rb=n，所以此方程組有唯一解.) x=zeros(n,1); c=zeros(1,n+1); for p= 1:n-1y,j=max(abs(b(p:n,p); c=b(p,:);b(p,:)= b(j+p-1,:); b(j+p-1,:)=c;for k=p+1:n m= b(k,p)/ b(p,p); b(k,p:n+1)= b(k,p:n+1)-m* b(p,p:n+1);endend湖北民族學(xué)院理學(xué)院數(shù)值分析課程實(shí)驗(yàn)報(bào)告（六）課程名稱數(shù)值分析班級(jí)0212404實(shí)驗(yàn)日期2014/12/23姓名向義俊學(xué)號(hào)021240418實(shí)驗(yàn)成績實(shí)驗(yàn)名稱實(shí)驗(yàn)四 線性方程組的迭代法（驗(yàn)

14、證性實(shí)驗(yàn)）實(shí)驗(yàn)?zāi)康募耙?、掌握雅克比迭代算法程序設(shè)計(jì)2、掌握高斯賽德爾迭代算法程序設(shè)計(jì)實(shí)驗(yàn)環(huán)境matlab軟件實(shí)驗(yàn)內(nèi)容編制程序，用雅克比迭代算法和高斯賽德爾迭代算法求解線性方程組，其中 ， 算法描述及實(shí)驗(yàn)步驟首先把所給內(nèi)容化成：x=bx+f的形式，然后給出初始向量：，并按迭代公式 進(jìn)行計(jì)算，如果按上述迭代公式所得到的向量序列收斂于某一向量，則就是方程組的解，并稱此迭代法收斂；否則稱為不收斂或發(fā)散。稱b為迭代矩陣。jacobi迭代設(shè)有線性方程組即式中，非奇異，且由式（*）得到 其相應(yīng)的迭代公式為 （*）稱迭代公式（*）為jacobi迭代。1、jacobi迭代法的程序算法：（1）取初始點(diǎn)，置k等

15、于0，精度要求和最大迭代次數(shù)n（2）用式（*）計(jì)算（3）若,則停止計(jì)算（作為線性方程組的解）。（4）若k=n，則停止計(jì)算（輸出某些信息），否則置k=k+1,轉(zhuǎn)（2）2、gauss-seidel迭代法其中。其迭代法優(yōu)點(diǎn)為只需一組存儲(chǔ)單元。實(shí)驗(yàn)結(jié)果和分析總結(jié)（1）在matlab命令窗口輸入如下命令 a = 5 2 1;-1 4 2;2 -3 10; b = -12 20 3; ep = 1e-5; x0 = 0 0 0; it_max = 90; x,k,index = jacobi(a,b,x0,ep,it_max)x = -3.9997 2.9998 1.9998k = 14index = 1

16、由index = 1可知，方程迭代14次后得到滿足要求的近似解。（1）在matlab命令窗口輸入如下命令 a = 5 2 1;-1 4 2;2 -3 10; b = -12 20 3; ep = 1e-5; x0 = 0 0 0; it_max = 90; x,k,index = jacobi(a,b,x0,ep,it_max)x = -4 3.0001 2k = 7index = 1由index = 1可知，方程迭代7次后得到滿足要求的近似解。附錄程序1、jacobi迭代法：function x,k,index = jacobi(a,b,x0,ep,it_max)if nargin 4 it_max = 100;endif nargin 3 ep = 1e-5;endn = length(a);k=0;x = zeros(n,1);y = zeros(n,1);index = 1;while 1 for i=1:n y(i)=b(i); for j=1:n if j=i y(i)=y(i)-a(i,j)*x(j); end end if abs(a(i,i