數(shù)值計(jì)算或計(jì)算方法試驗(yàn)教學(xué)講義

第一章實(shí)驗(yàn)的目的和要求1.1實(shí)驗(yàn)?zāi)康臑榱苏莆沼?jì)算方法的基本思想、原理和方法，要注意計(jì)算方法的處理技巧與計(jì)算機(jī)實(shí)現(xiàn)的結(jié)合，需將各種數(shù)值方法設(shè)計(jì)成算法，并編制好程序，拿到計(jì)算機(jī)上實(shí)現(xiàn)，最后得到可行性的驗(yàn)證。1.2實(shí)驗(yàn)要求⑴用C或C++、Java、FORTRAN、Matlab等計(jì)算機(jī)程序設(shè)計(jì)語(yǔ)言編寫(xiě)程序。⑵上機(jī)前充分準(zhǔn)備，復(fù)習(xí)相關(guān)知識(shí)，選用合適的數(shù)據(jù)結(jié)構(gòu)并詳細(xì)設(shè)計(jì)算法，盡量寫(xiě)出具有通用性的程序，反復(fù)檢查程序。⑶上機(jī)時(shí)快速輸入程序；首先排除語(yǔ)法錯(cuò)誤；然后采用多組數(shù)據(jù)，詳細(xì)測(cè)試，排除邏輯錯(cuò)誤；最后將程序調(diào)試成功，運(yùn)行程序得到準(zhǔn)確結(jié)果。⑷完成計(jì)算后，反復(fù)體會(huì)和分析，試著改善計(jì)算復(fù)雜性，使程序或算法更加完美。1.3實(shí)驗(yàn)環(huán)境1.3.1硬件環(huán)境CPU:Pentium4以上內(nèi)存：256MB以上1.3.2軟件環(huán)境操作系統(tǒng)：MicrosoftWindowsXP和2000編譯器：C或C++、Java、FORTRAN、Matlab1.4本實(shí)驗(yàn)課程與其它課程的關(guān)系本課程的前導(dǎo)課程有高等數(shù)學(xué)、線性代數(shù)(或高等代數(shù))、C語(yǔ)言或FORTRAN語(yǔ)言等，最好事先開(kāi)設(shè)數(shù)據(jù)結(jié)構(gòu)；后續(xù)課程有計(jì)算機(jī)圖形學(xué)、圖像處理、模式識(shí)別等。第二章實(shí)驗(yàn)的計(jì)劃和內(nèi)容2.1實(shí)驗(yàn)計(jì)劃計(jì)算方法實(shí)驗(yàn)課共安排20學(xué)時(shí)。計(jì)算方法實(shí)驗(yàn)計(jì)劃如下⑴Lagrange插值多項(xiàng)式⑵Newton插值多項(xiàng)式⑶Hermite插值多項(xiàng)式⑷最小二乘法⑸復(fù)化求積公式⑹Romberg求積公式⑺數(shù)值微分的外推算法⑻Gauss消元法⑼直接三角分解法⑽解方程組的迭代法2.2實(shí)驗(yàn)內(nèi)容共十個(gè)實(shí)驗(yàn)題目，每次課(2學(xué)時(shí))一個(gè)題目。2.2.1實(shí)驗(yàn)一實(shí)驗(yàn)題目：Lagrange插值多項(xiàng)式相關(guān)知識(shí)：通過(guò)n+1個(gè)節(jié)點(diǎn)的次數(shù)不超過(guò)n的Lagrange插值多項(xiàng)式為：L(x)=€yl(x)nkkk=0Hx?x匚，k=0,1,…，n。k x?xj=0kjj豐k另外，補(bǔ)充C語(yǔ)言繪制圖形方面的內(nèi)容如下1.屏幕坐標(biāo)系坐標(biāo)原點(diǎn)在屏幕的左上角，x軸水平向右，y軸垂直向下。2．常用的繪圖函數(shù)(繪圖庫(kù)函數(shù)所在的頭文件graphics.h)初始化圖形系統(tǒng)的函數(shù)voidinitgraph(int*graphdriver,int*graphmode,char*pathtodriver);畫(huà)點(diǎn)函數(shù)voidputpixel(intx,inty,intpixelcolor);移“畫(huà)筆”函數(shù)voidmoveto(intx,inty);畫(huà)直線函數(shù)voidline(intx1,inty1,intx2,inty2);voidlineto(intx,inty);設(shè)置前景顏色函數(shù)voidsetcolor(intcolor);設(shè)置背景顏色函數(shù)voidsetbkcolor(intcolor);設(shè)置畫(huà)線寬度和類型函數(shù)voidsetlinestyle(intlinestyle,unsignedupattern,intthickness);關(guān)閉圖形系統(tǒng)函數(shù)voidclosegraph(void);3.繪圖程序的設(shè)計(jì)模式#include""graphics.h"main(){intgraphdriver二DETECT,graphmode;initgraph(&graphdriver，&graphmode,""");調(diào)用繪圖函數(shù)進(jìn)行繪圖closegraph();}數(shù)據(jù)結(jié)構(gòu)：兩個(gè)一維數(shù)組或一個(gè)二維數(shù)組算法設(shè)計(jì)：(略)編寫(xiě)代碼：(略)實(shí)驗(yàn)用例：已知函數(shù)y=f(x)的一張表：x0102030405060708090100110120y517.534.5&815.56.5-5-10-24.57試驗(yàn)要求：利用Lagrange插值多項(xiàng)式L(x)求被插值函數(shù)f(x)在點(diǎn)x=65處的近似值。建n議：畫(huà)出Lagrange插值多項(xiàng)式L(x)的曲線。n實(shí)驗(yàn)題目：Newton插值多項(xiàng)式相關(guān)知識(shí)：通過(guò)n+1個(gè)節(jié)點(diǎn)的次數(shù)不超過(guò)n的Newton插值多項(xiàng)式為:N(x)，f(x)+f[x,x](x一x)+f[x,x,x](x一x)(x一x)+—n 0 01 0 012 0 1+f[x,x，…,x](x一x)(x一x).…(x一x)0 1 n 0 1 n—1數(shù)據(jù)結(jié)構(gòu)：兩個(gè)一維數(shù)組或一個(gè)二維數(shù)組算法設(shè)計(jì)：(略)編寫(xiě)代碼：(略)實(shí)驗(yàn)用例：已知函數(shù)y=f(x)的一張表(同上一個(gè)試驗(yàn))試驗(yàn)要求：利用Newton插值多項(xiàng)式N(x)求被插值函數(shù)f(x)在點(diǎn)x=65處的近似值。建議:n畫(huà)出Newton插值多項(xiàng)式N(x)的曲線。n2.2.3實(shí)驗(yàn)三實(shí)驗(yàn)題目：Hermite插值多項(xiàng)式相關(guān)知識(shí)：通過(guò)n+1個(gè)節(jié)點(diǎn)的次數(shù)不超過(guò)2n+1的Hermite插值多項(xiàng)式為:H (x)，工[y…(x)+m?(x)]2n+1 jjjjj，o其中，Hermite插值基函數(shù)…(x)，[1€2(x一x)1'(x)]l2(x)jjjjj?(x)，(x一x)12(x)jjj.—…(JT.sH)CjT.—JT—_)…h(huán)-R1(x)，工―1■jjx—xk，0jkk土jj，0,1,?…，n數(shù)據(jù)結(jié)構(gòu)：三個(gè)一維數(shù)組或一個(gè)二維數(shù)組算法設(shè)計(jì)：(略)編寫(xiě)代碼：(略)實(shí)驗(yàn)用例：已知函數(shù)y=f(x)的一張表(其中m，f'(x))：x0.100.200.300.400.50y0.9048370.8187310.7408180.6703200.606531m-0.904837-0.818731-0.740818-0.670320-0.606531x0.600.700.800.901.00y0.5488120.4965850.4493290.4065700.367879m-0.548812-0.496585-0.449329-0.406570-0.367879實(shí)驗(yàn)用例：利用Hermite插值多項(xiàng)式H十‘x)求被插值函數(shù)f(x)在點(diǎn)x=0.55處的近似值。建議：畫(huà)出Hermite插值多項(xiàng)式H (x)的曲線。2n+12.2.4實(shí)驗(yàn)四實(shí)驗(yàn)題目：曲線擬合的最小二乘法相關(guān)知識(shí)：已知C[a,b]中函數(shù)f(x)的一組實(shí)驗(yàn)數(shù)據(jù)(x,y)(i=O,l,…,m),其中y=f(x)。ii i i設(shè)9(x)(j，0,1,?…,n;n<m)是C[a,b]上線性無(wú)關(guān)函數(shù)族。在j=Span{90(X)，91(X)，…，9n(X)}中找函數(shù)f(x)曲線擬合的最小二乘解S*(x)二…a9(x)，其法方程(組)為:jjj，0…(9,9)a，d(k，0,1，?…,n)kjjkj，0其中，(9,9)，…①(x)9(x)9(x)jk ijikii，0k=0,1，…,n…①(x)f(x)9 (x)三dk=0,1，…,nTOC\o"1-5"\h\zi iki k特別是，求函數(shù)f(x)曲線擬合的線性最小二乘解S*(x)，ax+b的計(jì)算公式為:(…x2)(…y)-(…x)(…xy)i i i iib= z j40 (m+1)…x2-(…x)2

iii，0 i，0

(m(m+1)€xyii-(€x)(€y)ii—7x)2i i-—7x)2i(m+1)€x2一Ci數(shù)據(jù)結(jié)構(gòu)：兩個(gè)一維數(shù)組或一個(gè)二維數(shù)組算法設(shè)計(jì)：(略)編寫(xiě)代碼：(略)實(shí)驗(yàn)用例：已知函數(shù)y=f(x)的一張表:x0102030405060708090y6867.166.465.664.661.861.060.860.460試驗(yàn)要求：利用曲線擬合的線性最小二乘法求被逼近函數(shù)f(x)在點(diǎn)x=55處的近似值，并畫(huà)出實(shí)驗(yàn)數(shù)據(jù)和直線。2.2.5實(shí)驗(yàn)五實(shí)驗(yàn)題目：復(fù)化求積公式,b一a相關(guān)知識(shí)：將積分區(qū)間[a,b]n等分，節(jié)點(diǎn)x=a+kh,k=O,l,…,n,步長(zhǎng)h= 。復(fù)化k n梯形公式為：T=h[f(a)+2€f(a+kh)+f(b)]n2k=i再將每個(gè)小區(qū)間二等分，即整個(gè)積分區(qū)間[a,b]2n等分，此時(shí)復(fù)化梯形公式為T(mén)。復(fù)2n化梯形公式的遞推關(guān)系為T(mén)=-(T+H)2n2nn€ (2i—1)h b—a其中，H=h€f(a+ )，h=n 2 ni=1數(shù)據(jù)結(jié)構(gòu)：(略)算法設(shè)計(jì)：復(fù)化梯形公式的算法如下第一步：n=1,h=b-a；h第二步：T=-[f（a）€f（b）]；n2第三步：計(jì)算H；n第四步：計(jì)算T；2n第五步：n=2n,h=h/2；第六步：若T-T：>,（事先給定的誤差精度），則轉(zhuǎn)第三步；n n2第七步：輸出Tn和等分?jǐn)?shù)n/2，結(jié)束算法。編寫(xiě)代碼：（略）實(shí)驗(yàn)用例：I=卜excosxdxo試驗(yàn)要求：利用復(fù)化梯形求積公式求I二卜excosxdx的近似值（積分的精確值0I=T2.0703463164,兀二3.14159265358979323846...）,誤差精度，=10-6。2.2.6實(shí)驗(yàn)六實(shí)驗(yàn)題目：Romberg求積公式相關(guān)知識(shí)：用兩個(gè)相鄰的近似公式（其中后一個(gè)公式是由前一個(gè)公式的分半得到的）的線性組合而得到更好的近似公式的方法，就是近代電子計(jì)算機(jī)上常用的Romberg求積方法，也叫逐次分半加速（收斂）法。設(shè)以T（k）表示二分k次后求得的梯形值，且以T（k-j）表示序列｛T（k）｝的）次加速值。o j oRomberg求積公式的T表如下khT(k)0T(k)1T(k)2T(k)3T(k)4???0b-aT⑼0???1b一a2T⑴0T(0)1???2b一a4T⑵0T⑴1T(0)2???3b一a8T⑶0T⑵1T⑴2T⑼3???4b一a16T⑷0T⑶1T(2)2T⑴3T⑼4???■■■

Romberg求積公式(逐次分半加速公式)如下b—aT0(o)=丁[f(a)+f(b)]f b—a萌"/ .八b—a、、，<T(k)=—[T(k—i)+ 乙f(a+(2i—1) )],k=1,2,?…0 2 0 2k-1 2ki=14jT(k—j+l)—T(k—j)T(k—j)= j—l— 舊——,j=l,2,€,k,j 4j—l數(shù)據(jù)結(jié)構(gòu)：一個(gè)二維數(shù)組算法設(shè)計(jì):Romberg求積公式的算法如下h第一步：取k=0,h=b—a,求T(0)= [f(a)數(shù)據(jù)結(jié)構(gòu)：一個(gè)二維數(shù)組算法設(shè)計(jì):Romberg求積公式的算法如下h第一步：取k=0,h=b—a,求T(0)= [f(a)+f(b)]02第二步：令1—k(k記區(qū)間［a,b］的二分次數(shù))⑴求梯形值T(k)'按梯形的遞推公式;⑵求加速值，按公式逐個(gè)求出T表的第k行其余各元素T(k-j)(j=l,2,…,k)；jTk(0)…化(預(yù)先給定的誤差精度)，k+l_k，則轉(zhuǎn)⑴；第三步：輸出Tk(0)和等分?jǐn)?shù)2k(或二分次數(shù)k)，結(jié)束算法。編寫(xiě)代碼:實(shí)驗(yàn)用例:(略)I=flcos(x)dx02試驗(yàn)要求:利用Romberg求積公式求上述定積分(I=<2沁0.636619772),誤差精度<8=10-6。2.2.7實(shí)驗(yàn)七實(shí)驗(yàn)題目:數(shù)值微分相關(guān)知識(shí):數(shù)值微分的中點(diǎn)公式為相關(guān)知識(shí):2hf'(x)沁G(h)=f(X+h)—f(X—h)2h應(yīng)用理查森(Richardson)外推對(duì)h逐次分半，計(jì)算過(guò)程如下表(G0(hO二G(h))G0(h)???Gi(h)???

G(f)022G1(|)G2(h)???g(<)023g_(22)g2(2)G3(h)???????計(jì)算公式為,k€,k€0,1,2,…2k-1h h4jG —)-G()G(――)€——j-12k-尸1 j-12k…丿,j€1,2,，,k,j2k-j 4j—1數(shù)據(jù)結(jié)構(gòu)：一個(gè)二維數(shù)組算法設(shè)計(jì)：（略）編寫(xiě)代碼：（略）實(shí)驗(yàn)用例：f（x）€ex試驗(yàn)要求：利用數(shù)值微分的外推算法求f'（1）的近似值2.2.8實(shí)驗(yàn)八實(shí)驗(yàn)題目：用Gauss消元法求解線性代數(shù)方程組相關(guān)知識(shí)：在做除法運(yùn)算時(shí)，分母的絕對(duì)值越小，舍入誤差就越大。因此，消元的每一步都先選取絕對(duì)值比較大的元素（稱作主元），用它作分母再消元。這就是主元素消去法的基本思想。數(shù)據(jù)結(jié)構(gòu)：用一個(gè)二維數(shù)組存儲(chǔ)線性代數(shù)方程組的增廣矩陣；線性代數(shù)方程組的解最后存儲(chǔ)在增廣矩陣的最后一列上。算法設(shè)計(jì)：用列主元Gauss消元法求解線性代數(shù)方程組（同時(shí)求出系數(shù)行列式的值det）的算法為第一步：detT;第二步：對(duì)于k=l,2,…，n-1⑴按列選取主元，找ik⑴按列選取主元，找ik丘{k，k+1,，，n}ikk=max|aik⑵如果a€0（最好是a 為預(yù)先給定的一個(gè)非常小的數(shù)），則det=0,ik ikkk⑶如果i€k(或i豐k),則換行kka?a,j=k,k,1,…,n,1kj ikjdet--det⑷消元計(jì)算，對(duì)于i=k+l,k+2,…，na①a…一^ikakk②對(duì)于j=k+1,k+2,…，n+1a…a一axaij jikkj(5)det…detxakk第二步:如果a—0（最好是l第二步:如果a—0（最好是laInn nnanan,n+1⑴an,n+1 ann⑵對(duì)于i=n-l,n-2,⑵對(duì)于i=n-l,n-2,…，1對(duì)于j=i+l,i+2,…，nai:n+1aii第四步：det…detxann編寫(xiě)代碼:（略）編寫(xiě)代碼:（略）實(shí)驗(yàn)用例：實(shí)驗(yàn)用例：「10一701_x1-8_-32.099999621x5.9000012=5一15一1x2521023x41線性代數(shù)方程組為試驗(yàn)要求:利用列主元的Gauss消元法求解上述線性代數(shù)方程組（精確解為（0-111》），并同時(shí)求出系數(shù)行列式的值試驗(yàn)要求:2.2.9實(shí)驗(yàn)九實(shí)驗(yàn)題目：直接三角分解法相關(guān)知識(shí)：有矩陣A的三角LU分解，則求解線性代數(shù)方程組Ax=b的問(wèn)題就等價(jià)于求解兩個(gè)三角方程組Ly=b和Ux=y。而利用矩陣相等則對(duì)應(yīng)元素相等的事實(shí)，可逐一求出系數(shù)矩陣A的三角分解中L和U的各元素。a11a21a12a22?…a1n€?…a2n€=…1」u：?…u1n€?…u2n€/211aa?…a//…1〔?uLn1n2nnLn1n2nn分解過(guò)程的計(jì)算公式如下u=a一藝lu,i=r,r+1,，,n;r=1,2,，,nri ri rkkik=1l=（a一藝lu）.u,i=r,r+1,，,n;r=1,2,，,n一1ir ir ikkrFrrk=1 :數(shù)據(jù)結(jié)構(gòu)：一個(gè)一維數(shù)組先后存儲(chǔ)b、y和x；一個(gè)二維數(shù)組先存儲(chǔ)A,后被L和U覆蓋（者中的零元素不用存，L的對(duì)角元1亦不用存）算法設(shè)計(jì)：利用直接三角分解法求解線性代數(shù)方程組的算法第一步：分解對(duì)于r=1,2,…，n-1⑴求L的r列對(duì)于i=r+1,r+2,…，naJ（a-藝aa）/air ir ikkr rrk=1⑵求U的r+1行①l~r+l②