數(shù)值計(jì)算方法[共35頁]

《數(shù)值計(jì)算方法[共35頁]》由會員分享，可在線閱讀，更多相關(guān)《數(shù)值計(jì)算方法[共35頁]（34頁珍藏版）》請?jiān)谘b配圖網(wǎng)上搜索。

1、1.題目造倒數(shù)表，并例求 18 的倒數(shù)。（精度為 0.0005)2.算法原理2.1 牛頓迭代法牛頓迭代法是通過非線性方程線性化得到迭代序列的一種方法。對于非線性方程 f x() = 0 ，若已知根 x* 的一個(gè)近似值 xk ，將 f (x) 在 xk 處展成一階泰勒公式后忽略高次項(xiàng)可得：f (x) f x( k ) + f (xk )(x xk )右端是直線方程，用這個(gè)直線方程來近似非線性方程 f (x) 。將非線性方程 f x( ) = 0的根 x*代入 f x( *) = 0 ，即f x( k ) + f (xk )(x* xk ) 0*xk f (xk ) 解出x f (xk )將右端取

2、為 xk+1 ，則 xk+1 是比 xk 更接近于 x* 的近似值，即33f (xk )xk+1 xk f (xk ) 這就是牛頓迭代公式，相應(yīng)的迭代函數(shù)是f (x)(x) = x f (x)2.2 牛頓迭代法的應(yīng)用11計(jì)算 是求cx =1 0的解，解出 x，即得到 。取 cc 有牛頓迭代公式cxk 11 xk+1 = xk = cc 這樣就失去了迭代的意義，達(dá)不到迭代的效果。1f (x) = cx1， f (x) = c，故重新構(gòu)造方程： cx2 x = 0 ， 也是該式的解。故取 f (x) = cx2 x ，cf (x) = 2cx 1，則有牛頓迭代公式xk+1 = xk cxk2 xk

3、 = cxk2,k = 0,1,.2cxk 12ck 111的值在 之間，取初值 x0 = 0.1。20103.流程圖0,Nx讀入1k()0?0xf=1x輸出011kkxx+()()0100fxxxfx10?xx=4.輸出結(jié)果5.結(jié)果分析當(dāng)k= 3時(shí)，得 5 位有效數(shù)字 0.05 564。此時(shí)， x3 x4 = 0.00 000 x* 時(shí)必收斂，但是當(dāng) x0 0) 時(shí)迭代結(jié)果發(fā)散，較小尚不確定。6.心得體會起初對題目的理解并不是很透徹，另外對構(gòu)建牛頓迭代公式理論依據(jù)不是特別充分，比如說為什么在原有直接得到的式子兩邊各乘一個(gè) x，只是試出來的。在編程方面不夠成熟。當(dāng)然也加深了對牛頓迭代法的理解和

4、應(yīng)用的具體實(shí)現(xiàn)。實(shí)驗(yàn)二 例 3-41.題目用列主元消去法求解方程組 12x1 3x2 + 3x3 =1518x1 3x2 x3 = 15 x1 + x2 + x3 = 6 并求出系數(shù)矩陣 A的行列式的值det A。2.算法原理2.1 順序高斯消去法順序高斯消去法是利用線性方程組初等變換中的一種變換，即用一個(gè)不為零的數(shù)乘一個(gè)方程后加至另一個(gè)方程，使方程組變成同解的上三角方程組，然后再自下而上對上三角方程組求解。這樣，順序高斯消去法可分成“消去”和“回代” 兩個(gè)過程。在用順序高斯消去法時(shí)，在消元之前檢查方程組的系數(shù)矩陣的順序主子式，當(dāng)階數(shù)較高時(shí)是很難做到的。若線性方程組的系數(shù)具有某種性質(zhì)時(shí)，如常遇

5、到的對角占優(yōu)方程組，自然能夠用高斯消去法求解。2.2 列選主元消去法線性方程組只要系數(shù)矩陣非奇異，就存在惟一解，但是按順序消元過程中可能出現(xiàn)主元素akk( )k = 0，這時(shí)盡管系數(shù)矩陣非奇異，消元過程無法再進(jìn)行，或者即使akk( )k 0，但如果其絕對值很小，用它作除數(shù)也會導(dǎo)致其他元素的數(shù)量級急劇增大和使舍入誤差擴(kuò)大，將嚴(yán)重影響計(jì)算的精度。為避免在校園過程確定乘數(shù)時(shí)的所用除數(shù)是零或絕對值小的數(shù)，即零主元或小主元，在每一次消元之前，要增加一個(gè)選主元的過程，將絕對值大的元素交換到主對角線的位置上來。列選主元是當(dāng)高斯消元到第 k 步時(shí)，從 k 列的akk 以下（包括akk ）的各元素中選出絕對值

6、大的，然后通過行交換將其交換到akk 的位置上。交換系數(shù)矩陣中的兩行（包括常數(shù)項(xiàng)），只相當(dāng)于兩個(gè)方程的位置交換了，因此，列選主元不影響求解的結(jié)果。列選主元消去法常用來求行列式。設(shè)有矩陣a11 L a1n A = MM an1 L ann 用列主元消去法將其化為上三角形矩陣，對角線上元素為a11(1) ,a22(2) ,L,a33(3),于是行列式det A = (1)ma a11(1) 22(2) Lann( )n其中 m 為所進(jìn)行的行交換次數(shù)。這是實(shí)際中求行列式值的可靠方法。3.流程圖4.輸出結(jié)果5.結(jié)果分析采用計(jì)算機(jī)運(yùn)算在計(jì)算大數(shù)據(jù)時(shí)有明顯的優(yōu)點(diǎn)，另外也需要考慮到存儲。高斯消去法的使用條件

7、是akk( )k 0,k =1,2,L,n，而列選主元法可以保證這一條件。并且可以避免在消元過程確定乘數(shù)時(shí)所用除數(shù)是絕對值小的數(shù)，相對全選主元的運(yùn)算量小，一般也可以滿足精度要求。6.心得體會此次上機(jī)不僅需要對原理了解透徹，而且要求的編程能力較強(qiáng)。在定義和思路上沒問題，只是在編程軟件的使用上遇到些不穩(wěn)定的問題，如頭文件的使用。在存儲空間上得到了新的認(rèn)識，另外發(fā)現(xiàn)了當(dāng)代碼多時(shí)流程框圖的好處。編程是一件很需要耐心的事，自己還有很大進(jìn)步空間。實(shí)驗(yàn)三 例 3-101.題目用杜里特爾分解法求矩陣 A的逆矩陣 A1。 11 1A = 12 2211 2.算法原理2.1杜里特爾分解法設(shè)線性方程組Ax b= ，

8、對系數(shù)矩陣 A進(jìn)行除不交換兩行位置得初等行變換相當(dāng)于用初等矩陣M1 左乘 A，在對方程組第一次消元后，A(1) 和b(1) 分別化為 A(2) 和b(2) ，即M1 A(1) = A(2)M1b(1) = b(2) 1m211其中M1 = m311 MOmn11第 k 次消元時(shí)， A( )k 和b( )k 分別化為 A(k+1) 和b(k+1) ，即M Ak ( )k = A(k+1) M b( )k = b(k+1)k1其中 M1 = O1mk+1,k MmnkOO1消元過程是對k=1 n1進(jìn)行的，因此有Mn1LM M A21(1) = A( )nM LM M b2(1) = b( )n n

9、 11 將上三角形矩陣 A( )n 記為 U，于是有A = M M1121LMn11U = LU其中 1m211mm1L = M M1121LMn11 = 3132 MMm43 O MMMOmn1mn2mn3 L 1為單位下三角形矩形。這樣高斯消去法的實(shí)質(zhì)是將系數(shù)矩陣 A 分解為兩個(gè)三角形矩陣 L 和 U 相乘，即A=LU在上述矩陣描述中遇到了下三角形矩陣運(yùn)算。主對角線以上元素全為零的方陣稱為下三角形矩陣。下三角形矩陣的乘積仍是下三角形矩陣。若下三角形矩陣可逆，其逆矩陣仍是下三角形矩陣，而且下三角形矩陣的乘積和逆矩陣很容易求得。把 A 分解成一個(gè)單位下三角陣和一個(gè)上三角陣 U 的乘積成為杜里特

10、爾分解。這種分解是惟一的。2.2 高斯-約當(dāng)法高斯消去法有消元和回代兩個(gè)過程，當(dāng)對消元過程稍加改變便可以使方程組化為對角形方程組Dx b= 的形式，其中矩陣D為對角形矩陣，即a11(1)(2)D = a22Oann( )n 當(dāng)高斯-約當(dāng)消去法消元的每一步都先用主元去除其所在行的各元素（包括常數(shù)項(xiàng)）時(shí)，方程組便可化成1 x1 b1( )n 1 O xM2 = bM2( )n 1 xn ( )n bn這是等號右端即為方程組的解。高斯-約當(dāng)消去法每一步都用主元去除其所在行的各元素（包括常數(shù)項(xiàng)），這個(gè)個(gè)過程成為歸一化，這時(shí)方程組的系數(shù)陣轉(zhuǎn)化為單位陣。為減小誤差，高斯-約當(dāng)消去法還常用列選主元技術(shù)。4.

11、輸出結(jié)果5.結(jié)果分析采用杜里特爾分解法求解方程組時(shí)，由于把對系數(shù)矩陣的計(jì)算和對右端項(xiàng)的計(jì)算分開，這使計(jì)算線性方程組系非常方便。只需進(jìn)行一次矩形三角分解，然后再解多個(gè)三個(gè)方程組，且多解一個(gè)方程組僅需要增加大約n2 次乘除法運(yùn)算。采用高斯約當(dāng)法僅需要進(jìn)行消元?dú)w一，而不需要回代，為編程實(shí)現(xiàn)提供了便利。6.心得體會步驟很重要，審題-確定算法-解題步驟-流程圖-程序-代入簡單值進(jìn)行驗(yàn)證。在編程時(shí)先在代碼輸入?yún)^(qū)打好框架，并且盡量在每一命令后注釋，方便檢查錯(cuò)誤和日后復(fù)習(xí)。定義和變量存儲很靈活，如我把單位向量直接賦給了 A 矩陣變量中，還有根據(jù) 終的目的直接簡化計(jì)算。另外賦值前，確定存儲空間并且要定義初值為零

12、。實(shí)驗(yàn)四 例 4-61.題目已知 f (x) 的觀測數(shù)據(jù)x1234f( )x0-5-63構(gòu)造插值多項(xiàng)式。2.算法原理首先構(gòu)造基函數(shù)l xk ( )=i=n0 xxk xxii ，可以證明基函數(shù)滿足下列條件：i k0i kl xk ( i ) = ， 1 i = k對于給定n+1個(gè)節(jié)點(diǎn)，n次拉格朗日插值多項(xiàng)式由下式給出：nL x() = l x yk ()kk=0其中l(wèi) xk ( )=i=n0 xxk xxiii k由于l xk ( ) 是一個(gè)關(guān)于 x的n次多項(xiàng)式，所以L x( ) 為關(guān)于 x的不高于n次的代數(shù)多項(xiàng)式。當(dāng) x = xi 時(shí)，L x( i ) = yi ，滿足插值條件。3.流程圖1

13、t00yk,n0,1,iii=L輸入（x,y）0,k1,k1,njkjxxttxxj+=yktyy+=4.輸出結(jié)果5.結(jié)果分析由于所知的拉格朗日計(jì)算機(jī)算法只能實(shí)現(xiàn)計(jì)算某一特定值的近似函數(shù)值，而不知如何導(dǎo)出表達(dá)式，故例求 x=2.5 處的函數(shù)值以說明表達(dá)式以得出，只是在計(jì)算機(jī)程序中。并且也能達(dá)到拉格朗日插值法使用的目的。6.心得體會編程不夠細(xì)心，程序沒問題，卻因?yàn)椴恢朗禽斎胛募e(cuò)了而檢查了好長時(shí)間。但同時(shí)也加深了對拉格朗日基函數(shù)性質(zhì)的認(rèn)識和理解。實(shí)驗(yàn)五 習(xí)題 5-21.題目給出平面函數(shù) z x y( ,) = ax +by +c的數(shù)據(jù)i12345xi0.10.20.40.60.9yi0.20.

14、30.50.70.8zi0.580.630.730.830.92按 小二乘原理確定a b c， 。2.算法原理2.1 小二乘原理設(shè)已知某物理過程 y = f x( )的一組觀測數(shù)據(jù)(xi , f x( i ),i =1,2,L,m要求在某特定函數(shù)類( )x 中尋找一個(gè)函數(shù)( )x 作為 y = f x( )的近似函數(shù)，使得二者在 xi 上的誤差或稱殘差i =(xi ) f x( i ) , i =1,2,L,m按某種度量標(biāo)準(zhǔn)為 小，這就是擬合問題。1)ii要求殘差i 按某種度量標(biāo)準(zhǔn)為 小，即要求由殘差i 構(gòu)成的殘差向量 =0, 1,L, m T 的某種范數(shù) 為 小。例如，要求 ，或 即mm 1

15、 = i = (xi ) f x( i=0i=0 = max i = max ( )xi f x( )ii為 小，這本來都是很自然的，可是計(jì)算不太方便。所以通常要求：11 2 =m i2 2 =m ( )xi f x( )i 2 2 i=0 i=0或者mm 22= i2 =(xi ) f x( i ) 2i=0i=0為 小。這種要求誤差平方和 小的擬合稱為曲線擬合的 小二乘法。2.2多變量數(shù)據(jù)擬合對于給定的一組數(shù)據(jù) (xi , yi ) ，i=1，2，m ，尋求做n次多項(xiàng)式ny =a xkkk=0使性能指標(biāo)J a a( 0,1,L,an ) = (yi a xkik )2 為 小。i=1k=0

16、由于性能指標(biāo)J可以被看做關(guān)于ak ，k=0，1，n的多元函數(shù)，故上述擬合多項(xiàng)式的構(gòu)造問題可轉(zhuǎn)化為多元函數(shù)的極值問題。令J = 0從而有正則方程組 m xixi23 xixi2xi MMMLLLxin a0 yi n+1 a1 xi yi xi M = MM ak xinxin+1xin+2 L xi2n anxin yi 對多變量（或稱多元）線性模型y* = a0 +a x1 1 +a x22 +L+a xnn進(jìn)行了m次觀測y1* = a0 +a x1 11 +a x221 +L+a xn1ny2* = a0 +a x1 12 +a x222 +L+a xn2n Myn* = a0 +a x1

17、 1n +a x22n +L+a xnnn這個(gè)稱為回歸方程組，寫成矩陣形式y(tǒng)=Ay1* 1x11 * 其中y= 2 , A=1x12yM MM * 1x1mym x21 L xn1 a0 x22 L xn2 , =a1 。M M M M x2m L xnm an 當(dāng)m n 時(shí)，要確定一組a ii , = 0,1,L,n, 使之精確地滿足 m 個(gè)方程，這是超定方程組的問題，只能在 小平方誤差的基礎(chǔ)上確定i 。定義殘差向量= 1,2,L,m T , 則= y-A其中 y = y y1,2,L, ym T 代替輸出向量。取性能指標(biāo)J = T使之 小，以此確定出。由J = T = (y A) (T y

18、 A) =y yT TAy yTA+TA AT 利用向量和矩陣的運(yùn)算公式，有A AT =AT y此即為正則方程組，當(dāng)A AT 非奇異時(shí)，可求得 = (A AT)1AT y3.流程圖nizyxiii,2,1,L=輸入n5niATyATxATiiiii,.,2,1;1321=k1kniikikjnijikiYzajXaa=113,2,13?=kkk+1=輸出參數(shù)4.輸出結(jié)果5.結(jié)果分析曲線擬合的 小二乘法是反映所給數(shù)據(jù)點(diǎn)的總的趨勢，并不是嚴(yán)格的通過每個(gè)數(shù)據(jù)點(diǎn)，這樣就避免了大量數(shù)據(jù)插值時(shí)需要高次多項(xiàng)式，同時(shí)又去掉了數(shù)據(jù)所含的測量誤差。6.心得體會整理思路越來越熟練了，所以執(zhí)行各個(gè)步驟也相對簡單了很多

19、。另外對原理也加深了認(rèn)識。附錄：1.造倒數(shù)表1源程序：#includestdio.h#includemath.h #define N 30 void main()int i; float xN,c; FILE *fp1,*fp2; fp1=fopen(input1.txt,r); fp2=fopen(output1.txt,w); fscanf(fp1,%f,&c); fscanf(fp1,%f,&x0);/初值 fprintf(fp2, *倒數(shù)表*n); for(i=0;iN;i+)/牛頓迭代法 xi+1=xi*xi*c/(2*c*xi-1); fprintf(fp2,k=%dtx(%d)=

20、%.5fn,i,i,xi); if(fabs(xi+1-xi)輸入文件：input1.txt180.13輸出文件：“output1.txt”*倒數(shù)表*k=0x(0)=0.10000 k=1x(1)=0.06923 k=2x(2)=0.05781 k=3x(3)=0.05564 k=4x(4)=0.05564計(jì)算結(jié)果：1/18.000000=0.0562.例 3-41源程序：#includeiostream #includecmath using namespace std; #define N 10 void main()int i,j,k,l,n; float bN,aNN,t,d,det=

21、1.0;/*數(shù)據(jù)輸入*/FILE *fp1,*fp2; fp1=fopen(input2.txt,r); fp2=fopen(output2.txt,w); fscanf(fp1,%d,&n); for(i=0;in;i+) for(j=0;jn;j+) fscanf(fp1,%f,&aij); for(i=0;in;i+) fscanf(fp1,%f,&bi);/*數(shù)據(jù)輸入*/*高斯消去*/*消元*/ /*列選主元函數(shù)*/ for(k=0;kn-1;k+)/從第一次消元到第 N-1 次消元 d=akk;l=k;for(i=k+1;ifabs(d)d=aik;i=l;if(i=n)/判斷是否奇

22、異，不奇異進(jìn)行行交換 if(d=0)fprintf(fp2,奇異);/如果所有行的“首列”都為 0，為奇異elseif(l!=k)/如果第 k 行的“首列”并不是 大det=det*(-1);for(j=k;j=n;j+)/交換系數(shù)矩陣中的兩行t=alj;alj=akj;akj=t;t=bl;bl=bk;bk=t;/交換右端常向量中的兩行/*列選主元函數(shù)*/for(i=k+1;in;i+)/第（k+1）次消元要得到 N-(k+1)個(gè)乘數(shù)aik=aik/akk;for(j=k+1;j=0;i-)/從倒數(shù)第二項(xiàng)開始依次回代 N-1 次 t=0; for(j=i+1;jn;j+) t=t+aij*b

23、j; bi=(bi-t)/aii;/* *高斯消去*/*數(shù)據(jù)輸出*/ for(i=0;in;i+)/輸出方程組的解fprintf(fp2,x(%d)=%.4fn,i+1,bi); for(i=0;i輸入文件：“input2.txt”312 -3 3-18 3 -11 1 115 -15 63輸出文件：“output2.txt”x(1)=1.0000 x(2)=2.0000 x(3)=3.0000 detA=-66.00003. 例 3-101源程序：#includeiostream #includecmath using namespace std; #define N 30 void mai

24、n()int i,j,r,k,n;float aNN=0,s;/*數(shù)據(jù)輸入*/FILE *fp1,*fp2; fp1=fopen(input3.txt,r); fp2=fopen(output3.txt,w); fscanf(fp1,%d,&n); for(i=0;in;i+) for(j=0;jn;j+) fscanf(fp1,%f,&aij);for(i=0;in;i+)/輸入單位陣 j=i+n; aij=1;/*數(shù)據(jù)輸入*/*LU 分解*/ for(i=1;in;i+)/r=0 時(shí)：1區(qū)間不變化，2區(qū)間變化。ai0=ai0/a00;for(i=1;in;i+)/第二行列至末行列進(jìn)行變化時(shí)

25、 for(j=i;j2*n;j+)/第2（r+1）-1區(qū)間的變化行 s=0.0;for(k=0;ki;k+)/對應(yīng)行列元素的乘積進(jìn)行求和s+=aik*akj; aij=aij-s;for(j=i+1;jn;j+)/2（r+1）區(qū)間的變化列 s=0.0;for(k=0;kLUY/*LU 分解*/*高斯約當(dāng)法解 Ux=Y*/*提取 UY 減少計(jì)算*/for(i=1;in;i+) for(j=0;ji;j+) aij=0;/*提取 UY 減少計(jì)算*/*消元*/ for(j=0;j2*n;j+)/首行歸一化a0j=a0j/a00;a00=1;/第一列其余行已為零for(i=1;in;i+)/(n-1)

26、次歸一消元 for(j=i+1;j2*n;j+)/第 i+1 行的各列進(jìn)行歸一化aij=aij/aii; aii=1;for(r=0;ri;r+)/對第一行至第 i-1 行的 i 列進(jìn)行置零for(j=r+2;j2*n;j+)/r 行的各列與第 r+1 行的對應(yīng)列進(jìn)行減法運(yùn)算 arj=arj-aij*arr+1;/第 r+1 行歸一后，乘數(shù)即為要置零處的值 arr+1=0;/乘數(shù)用完之后置零即可/*消元*/*高斯約當(dāng)法解 Ux=Y*/*數(shù)據(jù)輸出*/ fprintf(fp2,A 的逆矩陣為nn);for(i=0;in;i+) fprintf(fp2,|); for(j=0;j輸入文件：input

27、3.txt31 1 -11 2 -2-2 1 13輸出文件：output3.txtA 的逆矩陣為|2.0000-1.00000.0000|1.5000-0.50000.5000|2.5000-1.50000.5000|4.例 4-61源程序：#includestdio.h#includemath.h#define N 50int n,i,j,k; float xx=0.0,yy=0.0,t,xN,yN,cN,AN;main() FILE *fp1,*fp2; fp1=fopen(input4.txt,r); fp2=fopen(output4.txt,w); fscanf(fp1,%d,&n)

28、;/n=4 for(i=0;in;i+) fscanf(fp1,%f,%f,&xi,&yi);/1,0,2,-5,3,-6,4,3fscanf(fp1,%f,&xx);/要計(jì)算的值 for(k=0;kn;k+)/拉格朗日插值t=1.0; for(j=0;jk;j+) t=(xx-xj)*t/(xk-xj); for(j=k+1;j輸入文件：input4.txt41,02,53,64,32.53輸出文件：output4.txtx=2.5000000 處的函數(shù)值為:y=-6.37500005.習(xí)題 5-21源程序：#include stdio.h#include math.h #define N

29、30 void main()int i,n,k,j,l;float xN,yN,zN;/定義輸入變量 float AT3N;/定義 A 的轉(zhuǎn)置 float X33,Y3;/定義中間變量 ATA 和 ATyfloat s,t,d;FILE *fp1,*fp2;/*數(shù)據(jù)輸入*/fp1=fopen(input5.txt,r); fp2=fopen(output5.txt,w); fscanf(fp1,%d,&n); for(i=0;in;i+) fscanf(fp1,%f,%f,%f,&xi,&yi,&zi);/*數(shù)據(jù)輸入*/for(i=0;in;i+)/為 AT 賦值；AT0i=1;AT1i=xi

30、;AT2i=yi; for(k=0;k3;k+)/中間變量 ATAfor(j=0;j3;j+) s=0.0; for(i=0;in;i+) s+=ATki*ATji;Xkj=s; for(k=0;k3;k+)/中間變量 ATy s=0.0; for(i=0;in;i+)s+=ATki*zi;Yk=s;/*計(jì)算/*高斯消去*/*消元*/ /*列選主元函數(shù)*/n=3;for(k=0;kn-1;k+)/從第一次消元到第 N-1 次消元 d=Xkk;l=k;for(i=k+1;ifabs(d)d=Xik;i=l;if(i=n)/判斷是否奇異，不奇異進(jìn)行行交換 if(d=0)fprintf(fp2,奇異

31、);/如果所有行的“首列”都為 0，為奇異elsefor(j=k;j=n;j+)/交換系數(shù)矩陣中的兩行t=Xlj;Xlj=Xkj;Xkj=t;t=Yl;Yl=Yk;Yk=t;/交換右端常向量中的兩行/*列選主元函數(shù)*/for(i=k+1;in;i+)/第（k+1）次消元要得到 N-(k+1)個(gè)乘數(shù)Xik=Xik/Xkk;for(j=k+1;j=0;i-)/從倒數(shù)第二項(xiàng)開始依次回代 N-1 次 t=0; for(j=i+1;j輸入文件：input5.txt50.1,0.2,0.580.2,0.3,0.630.4,0.5,0.730.6,0.7,0.830.9,0.8,0.923輸出文件：output5.txt所求參數(shù)為：a=0.2000 b=0.3000 c=0.5000