哈工大-數(shù)值分析上機(jī)實(shí)驗(yàn)報(bào)告.doc

實(shí)驗(yàn)報(bào)告 16SD33223 陳凱實(shí)驗(yàn)報(bào)告一題目：非線(xiàn)性方程求解摘要：非線(xiàn)性方程的解析解通常很難給出，因此線(xiàn)性方程的數(shù)值解法就尤為重要。本實(shí)驗(yàn)采用兩種常見(jiàn)的求解方法二分法和Newton法及改進(jìn)的Newton法。前言：（目的和意義）掌握二分法與Newton法的基本原理和應(yīng)用。數(shù)學(xué)原理：對(duì)于一個(gè)非線(xiàn)性方程的數(shù)值解法很多。在此介紹兩種最常見(jiàn)的方法：二分法和Newton法。對(duì)于二分法，其數(shù)學(xué)實(shí)質(zhì)就是說(shuō)對(duì)于給定的待求解的方程f(x)，其在a,b上連續(xù)，f(a)f(b)0，且f(x)在a,b內(nèi)僅有一個(gè)實(shí)根x*，取區(qū)間中點(diǎn)c，若，則c恰為其根，否則根據(jù)f(a)f(c)5e-6) ; c=(a+b)/2; if f12(a)*f12(c)0; a=c; else b=c; end R=b-a;%求出誤差k=k+1;endx=c%給出解Newton法及改進(jìn)的Newton法源程序：clear% 輸入函數(shù)f=input(請(qǐng)輸入需要求解函數(shù),s)%求解f(x)的導(dǎo)數(shù)df=diff(f);%改進(jìn)常數(shù)或重根數(shù)miu=2;%初始值x0x0=input(input initial value x0);k=0;%迭代次數(shù)max=100;%最大迭代次數(shù)R=eval(subs(f,x0,x);%求解f(x0)，以確定初值x0時(shí)否就是解while (abs(R)1e-8) x1=x0-miu*eval(subs(f,x0,x)/eval(subs(df,x0,x); R=x1-x0; x0=x1; k=k+1;if (eval(subs(f,x0,x)max;%如果迭代次數(shù)大于給定值，認(rèn)為迭代不收斂，重新輸入初值 ss=input(maybe result is error,choose a new x0,y/n?,s); if strcmp(ss,y) x0=input(input initial value x0); k=0; else break end endendk;%給出迭代次數(shù)x=x0；%給出解結(jié)果分析和討論：1. 用二分法計(jì)算方程在1，2內(nèi)的根。(,下同)計(jì)算結(jié)果為x= 1.40441513061523；f(x)= -3.797205105904311e-007；k=18；由f(x)知結(jié)果滿(mǎn)足要求，但迭代次數(shù)比較多，方法收斂速度比較慢。2. 用二分法計(jì)算方程在1，1.5內(nèi)的根。計(jì)算結(jié)果為x= 1.32471847534180；f(x)= 2.209494846194815e-006；k=17；由f(x)知結(jié)果滿(mǎn)足要求，但迭代次數(shù)還是比較多。3. 用Newton法求解下列方程a) x0=0.5；計(jì)算結(jié)果為x= 0.56714329040978；f(x)= 2.220446049250313e-016；k=4；由f(x)知結(jié)果滿(mǎn)足要求，而且又迭代次數(shù)只有4次看出收斂速度很快。b) x0=1；c) x0=0.45, x0=0.65； 當(dāng)x0=0.45時(shí)，計(jì)算結(jié)果為x= 0.49999999999983；f(x)= -8.362754932994584e-014；k=4；由f(x)知結(jié)果滿(mǎn)足要求，而且又迭代次數(shù)只有4次看出收斂速度很快，實(shí)際上該方程確實(shí)有真解x=0.5。當(dāng)x0=0.65時(shí)，計(jì)算結(jié)果為x= 0.50000000000000；f(x)=0；k=9；由f(x)知結(jié)果滿(mǎn)足要求，實(shí)際上該方程確實(shí)有真解x=0.5，但迭代次數(shù)增多，實(shí)際上當(dāng)取x00.68時(shí)，x1，就變成了方程的另一個(gè)解，這說(shuō)明Newton法收斂與初值很有關(guān)系，有的時(shí)候甚至可能不收斂。4. 用改進(jìn)的Newton法求解，有2重根，取 x0=0.55；并與3.中的c)比較結(jié)果。當(dāng)x0=0.55時(shí)，程序死循環(huán)，無(wú)法計(jì)算，也就是說(shuō)不收斂。改時(shí)，結(jié)果收斂為x=0.50000087704286；f(x)=4.385198907621127e-007；k=16；顯然這個(gè)結(jié)果不是很好，而且也不是收斂至方程的2重根上。當(dāng)x0=0.85時(shí)，結(jié)果收斂為x= 1.00000000000489；f(x)= 2.394337647718737e-023；k=4；這次達(dá)到了預(yù)期的結(jié)果，這說(shuō)明初值的選取很重要，直接關(guān)系到方法的收斂性，實(shí)際上直接用Newton法，在給定同樣的條件和精度要求下，可得其迭代次數(shù)k=15，這說(shuō)明改進(jìn)后的Newton法法速度確實(shí)比較快。結(jié)論： 對(duì)于二分法，只要能夠保證在給定的區(qū)間內(nèi)有根，使能夠收斂的，當(dāng)時(shí)收斂的速度和給定的區(qū)間有關(guān)，二且總體上來(lái)說(shuō)速度比較慢。Newton法，收斂速度要比二分法快，但是最終其收斂的結(jié)果與初值的選取有關(guān)，初值不同，收斂的結(jié)果也可能不一樣，也就是結(jié)果可能不時(shí)預(yù)期需要得結(jié)果。改進(jìn)的Newton法求解重根問(wèn)題時(shí)，如果初值不當(dāng)，可能會(huì)不收斂，這一點(diǎn)非常重要，當(dāng)然初值合適，相同情況下其速度要比Newton法快得多。實(shí)驗(yàn)報(bào)告二題目： Gauss列主元消去法摘要：求解線(xiàn)性方程組的方法很多，主要分為直接法和間接法。本實(shí)驗(yàn)運(yùn)用直接法的Guass消去法,并采用選主元的方法對(duì)方程組進(jìn)行求解。前言：（目的和意義）1. 學(xué)習(xí)Gauss消去法的原理。2. 了解列主元的意義。3. 確定什么時(shí)候系數(shù)陣要選主元數(shù)學(xué)原理：由于一般線(xiàn)性方程在使用Gauss消去法求解時(shí)，從求解的過(guò)程中可以看到，若=0，則必須進(jìn)行行交換，才能使消去過(guò)程進(jìn)行下去。有的時(shí)候即使0，但是其絕對(duì)值非常小，由于機(jī)器舍入誤差的影響，消去過(guò)程也會(huì)出現(xiàn)不穩(wěn)定得現(xiàn)象，導(dǎo)致結(jié)果不正確。因此有必要進(jìn)行列主元技術(shù)，以最大可能的消除這種現(xiàn)象。這一技術(shù)要尋找行r，使得并將第r行和第k行的元素進(jìn)行交換，以使得當(dāng)前的的數(shù)值比0要大的多。這種列主元的消去法的主要步驟如下：1. 消元過(guò)程對(duì)k=1,2,n-1,進(jìn)行如下步驟。1) 選主元，記若很小，這說(shuō)明方程的系數(shù)矩陣嚴(yán)重病態(tài)，給出警告，提示結(jié)果可能不對(duì)。2) 交換增廣陣A的r，k兩行的元素。 (j=k,n+1)3) 計(jì)算消元 (i=k+1,n; j=k+1,n+1)2. 回代過(guò)程對(duì)k= n, n-1,1,進(jìn)行如下計(jì)算至此，完成了整個(gè)方程組的求解。程序設(shè)計(jì)：本實(shí)驗(yàn)采用Matlab的M文件編寫(xiě)。 Gauss消去法源程序：cleara=input(輸入系數(shù)陣：n)b=input(輸入列陣b：n)n=length(b);A=a bx=zeros(n,1);%函數(shù)主體for k=1:n-1;%是否進(jìn)行主元選取if abs(A(k,k)abs(t) p=r; else p=k; end end %交換元素 if p=k; for q=k:n+1; s=A(k,q); A(k,q)=A(p,q); A(p,q)=s; end end end %判斷系數(shù)矩陣是否奇異或病態(tài)非常嚴(yán)重if abs(A(k,k) yipusilongdisp(矩陣奇異，解可能不正確)end %計(jì)算消元,得三角陣 for r=k+1:n; m=A(r,k)/A(k,k); for q=k:n+1; A(r,q)=A(r,q)-A(k,q)*m; end endend %求解x x(n)=A(n,n+1)/A(n,n); for k=n-1:-1:1; s=0; for r=k+1:n; s=s+A(k,r)*x(r); end t=(A(k,n+1)-s) x(k)=(A(k,n+1)-s)/A(k,k)end結(jié)果分析和討論：例：求解方程。其中為一小數(shù)，當(dāng)時(shí)，分別采用列主元和不列主元的Gauss消去法求解，并比較結(jié)果。記Emax為求出的解代入方程后的最大誤差，按要求，計(jì)算結(jié)果如下：當(dāng)時(shí)，不選主元和選主元的計(jì)算結(jié)果如下，其中前一列為不選主元結(jié)果，后一列為選主元結(jié)果，下同。 0.99999934768391 0.99999934782651 2.00000217421972 2.00000217391163 2.99999760859451 2.99999760869721Emax= 9.301857062382624e-010，0此時(shí)，由于不是很小，機(jī)器誤差就不是很大，由Emax可以看出不選主元的計(jì)算結(jié)果精度還可以，因此此時(shí)可以考慮不選主元以減少計(jì)算量。當(dāng)時(shí)，不選主元和選主元的計(jì)算結(jié)果如下 1.00001784630877 0.99999999999348 1.99998009720807 2.00000000002174 3.00000663424731 2.99999999997609Emax= 2.036758973744668e-005，0此時(shí)由Emax可以看出不選主元的計(jì)算精度就不好了，誤差開(kāi)始增大。當(dāng)時(shí)，不選主元和選主元的計(jì)算結(jié)果如下 1.42108547152020 1.00000000000000 1.66666666666666 2.00000000000000 3.11111111111111 300000000000000Emax= 0.70770085900503，0此時(shí)由Emax可以看出，不選主元的結(jié)果應(yīng)該可以說(shuō)是不正確了，這是由機(jī)器誤差引起的。當(dāng)時(shí)，不選主元和選主元的計(jì)算結(jié)果如下NaN 1NaN 2NaN 3Emax=NaN, 0不選主元時(shí)，程序報(bào)錯(cuò)：Warning: Divide by zero.。這是因?yàn)闄C(jī)器計(jì)算的最小精度為10-15，所以此時(shí)的就認(rèn)為是0，故出現(xiàn)了錯(cuò)誤現(xiàn)象。而選主元時(shí)則沒(méi)有這種現(xiàn)象，而且由Emax可以看出選主元時(shí)的結(jié)果應(yīng)該是精確解。結(jié)論：采用Gauss消去法時(shí)，如果在消元時(shí)對(duì)角線(xiàn)上的元素始終較大（假如大于10-5），那么本方法不需要進(jìn)行列主元計(jì)算，計(jì)算結(jié)果一般就可以達(dá)到要求，否則必須進(jìn)行列主元這一步，以減少機(jī)器誤差帶來(lái)的影響，使方法得出的結(jié)果正確。實(shí)驗(yàn)報(bào)告三題目： Rung現(xiàn)象產(chǎn)生和克服摘要：由于高次多項(xiàng)式插值不收斂，會(huì)產(chǎn)生Runge現(xiàn)象，本實(shí)驗(yàn)在給出具體的實(shí)例后，采用分段線(xiàn)性插值和三次樣條插值的方法有效的克服了這一現(xiàn)象，而且還取的很好的插值效果。前言：（目的和意義）1. 深刻認(rèn)識(shí)多項(xiàng)式插值的缺點(diǎn)。2. 明確插值的不收斂性怎樣克服。3. 明確精度與節(jié)點(diǎn)和插值方法的關(guān)系。數(shù)學(xué)原理：在給定n+1個(gè)節(jié)點(diǎn)和相應(yīng)的函數(shù)值以后構(gòu)造n次的Lagrange插值多項(xiàng)式，實(shí)驗(yàn)結(jié)果表明（見(jiàn)后面的圖）這種多項(xiàng)式并不是隨著次數(shù)的升高對(duì)函數(shù)的逼近越來(lái)越好，這種現(xiàn)象就是Rung現(xiàn)象。解決Rung現(xiàn)象的方法通常有分段線(xiàn)性插值、三次樣條插值等方法。分段線(xiàn)性插值：設(shè)在區(qū)間a, b上，給定n+1個(gè)插值節(jié)點(diǎn)a=x0x1xn=b和相應(yīng)的函數(shù)值y0，y1，yn，求作一個(gè)插值函數(shù)，具有如下性質(zhì)：1) ，j=0，1，n。2) 在每個(gè)區(qū)間xi, xj上是線(xiàn)性連續(xù)函數(shù)。則插值函數(shù)稱(chēng)為區(qū)間a, b上對(duì)應(yīng)n個(gè)數(shù)據(jù)點(diǎn)的分段線(xiàn)性插值函數(shù)。三次樣條插值：給定區(qū)間a, b一個(gè)分劃 ：a=x0x1xN=b 若函數(shù)S(x)滿(mǎn)足下列條件：1) S(x)在每個(gè)區(qū)間xi, xj上是不高于3次的多項(xiàng)式。2) S(x)及其2階導(dǎo)數(shù)在a, b上連續(xù)。則稱(chēng)S(x)使關(guān)于分劃的三次樣條函數(shù)。程序設(shè)計(jì)：本實(shí)驗(yàn)采用Matlab的M文件編寫(xiě)。其中待插值的方程寫(xiě)成function的方式，如下function y=f(x);y=1/（1+25*x*x）;寫(xiě)成如上形式即可，下面給出主程序 Lagrange插值源程序：n=input(將區(qū)間分為的等份數(shù)輸入：n);s=-1+2/n*0:n;%給定的定點(diǎn)，Rf為給定的函數(shù)x=-1:0.01:1;f=0;for q=1:n+1; l=1;%求插值基函數(shù) for k=1:n+1; if k=q; l=l.*(x-s(k)./(s(q)-s(k); else l=l; end end f=f+Rf(s(q)*l;%求插值函數(shù)endplot(x,f,r)%作出插值函數(shù)曲線(xiàn)grid on hold on分段線(xiàn)性插值源程序clearn=input(將區(qū)間分為的等份數(shù)輸入：n);s=-1+2/n*0:n;%給定的定點(diǎn)，Rf為給定的函數(shù)m=0;hh=0.001;for x=-1:hh:1; ff=0; for k=1:n+1;%求插值基函數(shù) switch k case 1 if xs(n); l=(x-s(n)./(s(n+1)-s(n); else l=0; end otherwise if x=s(k-1)&x=s(k)&xR);%精度控制 j=j+1; s=0; for p=1:2(j-2); s=s+f(a+(2*p-1)*h/(2(j-1); end T(1,j)=T(1,j-1)/2+h*s/(2(j-1); %梯形公式應(yīng)用 for m=2:j; k=(j-m+1); T(m,k)=(4(m-1)*T(m-1,k+1)-T(m-1,k)/(4(m-1)-1); endend%給出 Romberg積分法的函數(shù)表I=T(m,1)結(jié)果分析和討論： 進(jìn)行具體的積分時(shí)，精度取R=1e-8。1. 求積分。精確解I= 24999676。運(yùn)行程序得Romberg積分法的函數(shù)表為1.0e+007 * 4.70101520000000 3.05022950000000 2.63753307500000 2.49996760000000 2.49996760000000 0 2.49996760000000 0 0由函數(shù)表知Romberg積分給出的結(jié)果為2.4999676*107，與精確沒(méi)有誤差，精度很高。2. 求積分。精確解I=ln3= 1.09861228866811。運(yùn)行程序得Romberg積分法的函數(shù)表為1.33333333333333 1.16666666666667 1.11666666666667 1.10321067821068 1.09976770156303 1.09890151516846 1.098684618785591.11111111111111 1.10000000000000 1.09872534872535 1.09862004268048 1.09861278637027 1.09861231999130 01.09925925925926 1.09864037197371 1.09861302227749 1.09861230261625 1.09861228889937 0 01.09863054836600 1.09861258815533 1.09861229119306 1.09861228868164 0 0 01.09861251772313 1.09861229002850 1.09861228867179 0 0 0 01.09861228980593 1.09861228867046 0 0 0 0 01.09861228867019 0 0 0 0 0 0從積分表中可以看出程序運(yùn)行的結(jié)果為1.09861228867019，取8位有效數(shù)字，滿(mǎn)足要求。3. 求積分。直接