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

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

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

2、0 * xk ? f (xk ) 解出 x ≈ f (xk ) 將右端取為 xk+1 ，則 xk+1 是比 xk 更接近于 x* 的近似值，即 33 f (xk ) xk+1 ≈ xk ? f (xk ) 這就是牛頓迭代公式，相應(yīng)的迭代函數(shù)是 f (x) ?(x) = x ? f (x) 2.2 牛頓迭代法的應(yīng)用 1 1 計算 是求cx? =1 0的解，解出 x，即得到 。取 c c 有牛頓迭代公式 cxk ?1 1 xk+1 = xk ? = c c 這樣就失去了迭代的意義，達不到迭代的效果。 1 f (x) = cx?1， f (x) = c，

3、 故重新構(gòu)造方程： cx2 ? x = 0 ， 也是該式的解。故取 f (x) = cx2 ? x ， c f (x) = 2cx ?1，則有牛頓迭代公式 xk+1 = xk ? cxk2 ? xk = cxk2 , k = 0,1,... 2cxk ?1 2ck ?1 1 1 的值在 ~ 之間，取初值 x0 = 0.1。 20 10 3.流程圖 0 ,,N x ε 讀入 1 k ? ( ) 0 ?0 x f ′ = 1 x 輸出 0 1 1 k k x x ? + ? ( ) ( ) 0 1

4、0 0 f x x x f x ? ? ′ 1 0 ? x x ε ? < ≠ = < = ≥ ≠ 4.輸出結(jié)果 5.結(jié)果分析 當k= 3時，得 5 位有效數(shù)字 0.05 564。此時， x3 ? x4 = 0.00 000 < 0.0 005， 故取 x* = x3 = 0.05 564 ≈ 0.056 。 此種迭代格式仍存在一定的缺陷，經(jīng)實驗后發(fā)現(xiàn)當初值 x0 > x* 時必收斂，但是當 x0 < x* ???( > 0) 時迭代結(jié)果發(fā)散，?較小尚不確定。 6.心得體會 起初對題目的理解并不是很透徹，另外對構(gòu)建牛頓迭代公

5、式理論依據(jù)不是特別充分，比如說為什么在原有直接得到的式子兩邊各乘一個 x，只是試出來的。 在編程方面不夠成熟。當然也加深了對牛頓迭代法的理解和應(yīng)用的具體實現(xiàn)。 實驗二 例 3-4 1.題目 用列主元消去法求解方程組 ? 12x1 ? 3x2 + 3x3 = 15? ? ? ??18x1 ? 3x2 ? x3 = ?15? ?? x1 + x2 + x3 = 6?? 并求出系數(shù)矩陣 A的行列式的值det A。 2.算法原理 2.1 順序高斯消去法 順序高斯消去法是利用線性方程組初等變換中的一種變換，即用一個不為零的數(shù)乘一個方程后加至另一個方程，使方程組變成同解的上三角方

6、程組，然后再自下而上對上三角方程組求解。這樣，順序高斯消去法可分成“消去”和“回代” 兩個過程。 在用順序高斯消去法時，在消元之前檢查方程組的系數(shù)矩陣的順序主子式，當階數(shù)較高時是很難做到的。若線性方程組的系數(shù)具有某種性質(zhì)時，如常遇到的對角占優(yōu)方程組，自然能夠用高斯消去法求解。 2.2 列選主元消去法 線性方程組只要系數(shù)矩陣非奇異，就存在惟一解，但是按順序消元過程中可能出現(xiàn)主元素akk( )k = 0，這時盡管系數(shù)矩陣非奇異，消元過程無法再進行，或者 即使akk( )k ≠0，但如果其絕對值很小，用它作除數(shù)也會導(dǎo)致其他元素的數(shù)量級急劇 增大和使舍入誤差擴大，將嚴重影響計算的精度。 為

7、避免在校園過程確定乘數(shù)時的所用除數(shù)是零或絕對值小的數(shù)，即零主元或小主元，在每一次消元之前，要增加一個選主元的過程，將絕對值大的元素交換到主對角線的位置上來。 列選主元是當高斯消元到第 k 步時，從 k 列的akk 以下（包括akk ）的各元素中選出絕對值 大的，然后通過行交換將其交換到akk 的位置上。交換系數(shù)矩陣中的兩行（包括常數(shù)項），只相當于兩個方程的位置交換了，因此，列選主元不影響求解的結(jié)果。 列選主元消去法常用來求行列式。設(shè)有矩陣 ?a11 L a1n ? A =?? M M ?? ??an1 L ann ?? 用列主元消去法將其化為上三角形矩陣，對角線上元素為a11(1

8、) ,a22(2) ,L,a33(3),于是行列式 det A = ?( 1)ma a11(1) 22(2) Lann( )n 其中 m 為所進行的行交換次數(shù)。這是實際中求行列式值的可靠方法。 3.流程圖 4.輸出結(jié)果 5.結(jié)果分析 采用計算機運算在計算大數(shù)據(jù)時有明顯的優(yōu)點，另外也需要考慮到存儲。 高斯消去法的使用條件是akk( )k ≠ 0, k =1,2,L,n，而列選主元法可以保證這一 條件。并且可以避免在消元過程確定乘數(shù)時所用除數(shù)是絕對值小的數(shù)，相對全選主元的運算量小，一般也可以滿足精度要求。 6.心得體會 此次上機不僅需要對原理了解透徹，而且要求的編程能

9、力較強。在定義和思路上沒問題，只是在編程軟件的使用上遇到些不穩(wěn)定的問題，如頭文件的使用。在存儲空間上得到了新的認識，另外發(fā)現(xiàn)了當代碼多時流程框圖的好處。編程是一件很需要耐心的事，自己還有很大進步空間。 實驗三 例 3-10 1.題目 用杜里特爾分解法求矩陣 A的逆矩陣 A?1。 ? 1 1 ?1? A =?? 1 2 ?2?? ???2 1 1 ?? 2.算法原理 2.1杜里特爾分解法設(shè)線性方程組Ax b= ，對系數(shù)矩陣 A進行除不交換兩行位置得初等行變換相當于用初等矩陣M1 左乘 A，在對方程組第一次消元后，A(1) 和b(1) 分別化為 A(2) 和 b(2) ，即

10、 ??M1 A(1) = A(2) ???M1b(1) = b(2) ? 1 ? ???m21 1 ? 其中 M1 = ???m31 1 ??? ? M O ? ? ? ???mn1 1?? 第 k 次消元時， A( )k 和b( )k 分別化為 A(k+1) 和b(k+1) ，即 ??M Ak ( )k = A(k+1) ??M b( )k = b(k+1) ? k ?1 ? 其中 ? ? ? M1 = ?? ? ? ? ?? O 1 ?mk+1,k M ?mnk ? ? ?? ? O ? O ? ? 1?? 消元

11、過程是對k=1~ n?1進行的，因此有 ??Mn?1LM M A2 1 (1) = A( )n ???M ? LM M b2 (1) = b( )n n 1 1 將上三角形矩陣 A( )n 記為 U，于是有 A = M M1?1 2?1LMn??11U = LU 其中 ? 1 ?m21 1 ? ?m m 1 L = M M1?1 2?1LMn??11 = ? 31 32 ? M M m43 O ? M M M O ??mn1 mn2 mn3 L ? ? ? ? ? ? ? ? ? 1?? 為單位下三角形矩形。 這樣高斯消去法的實質(zhì)是將系

12、數(shù)矩陣 A 分解為兩個三角形矩陣 L 和 U 相乘，即 A=LU 在上述矩陣描述中遇到了下三角形矩陣運算。主對角線以上元素全為零的方陣稱為下三角形矩陣。下三角形矩陣的乘積仍是下三角形矩陣。若下三角形矩陣可逆，其逆矩陣仍是下三角形矩陣，而且下三角形矩陣的乘積和逆矩陣很容易求得。 把 A 分解成一個單位下三角陣和一個上三角陣 U 的乘積成為杜里特爾分解。 這種分解是惟一的。 2.2 高斯-約當法 高斯消去法有消元和回代兩個過程，當對消元過程稍加改變便可以使方程組化為對角形方程組 Dx b= 的形式，其中矩陣D為對角形矩陣，即 ?a11(1) ? ? (2) ? D = ?

13、a22 ? ? O ? ?? ann( )n ??? ? 當高斯-約當消去法消元的每一步都先用主元去除其所在行的各元素（包括常數(shù)項）時，方程組便可化成 ?1 ?? x1 ? ?b1( )n ? ??? 1 O ?????? xM2 ??? = ???bM2( )n ??? ?? 1???? xn ?? ??? ( )n ?? bn ? 這是等號右端即為方程組的解。高斯-約當消去法每一步都用主元去除其所在行的各元素（包括常數(shù)項），這個個過程成為歸一化，這時方程組的系數(shù)陣轉(zhuǎn)化為單位陣。為減小誤差，高斯-約當消去法還常用列選主元技術(shù)。 4.輸出結(jié)果 5.結(jié)果分析

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

15、為零。 實驗四 例 4-6 1.題目已知 f (x) 的觀測數(shù)據(jù) x 1 2 3 4 f( )x 0 -5 -6 3 構(gòu)造插值多項式。 2.算法原理 首先構(gòu)造基函數(shù)l xk ( )=∏i=n0 xxk ??xxii ，可以證明基函數(shù)滿足下列條件： i k≠ ?0 i ≠ k l xk ( i ) = ? ， ?1 i = k 對于給定n+1個節(jié)點，n次拉格朗日插值多項式由下式給出： n L x( ) = ∑l x yk ( ) k k=0 其中 l xk ( )=∏i=n0 xxk ??xxii i k≠ 由于l xk ( ) 是一

16、個關(guān)于 x的n次多項式，所以L x( ) 為關(guān)于 x的不高于n次的代數(shù)多項式。當 x = xi 時，L x( i ) = yi ，滿足插值條件。 3.流程圖 1 t ? 0 0 y k ? ? ,n 0 ,1, i i i = L 輸入（x,y） 0 , ,k1,k1, ,n j k j x x t t x x j ? ? ? ??? ? + ??? = y k ty y + ? ≠ = 4.輸出結(jié)果 5.結(jié)果分析 由于所知的拉格朗日計算機算法只能實現(xiàn)計算某一特定值的近似函數(shù)值，而不

17、知如何導(dǎo)出表達式，故例求 x=2.5 處的函數(shù)值以說明表達式以得出，只是在計算機程序中。并且也能達到拉格朗日插值法使用的目的。 6.心得體會 編程不夠細心，程序沒問題，卻因為不知道是輸入文件錯了而檢查了好長時間。但同時也加深了對拉格朗日基函數(shù)性質(zhì)的認識和理解。 實驗五 習題 5-2 1.題目 給出平面函數(shù) z x y( , ) = ax +by +c的數(shù)據(jù) i 1 2 3 4 5 xi 0.1 0.2 0.4 0.6 0.9 yi 0.2 0.3 0.5 0.7 0.8 zi 0.58 0.63 0.73 0.83 0.92 按 小二

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

19、δ ，或 δ ∞即 m m δ 1 =∑ ∑δi = ?(xi ) ? f x( i=0 i=0 δ ∞ = max δi = max ?( )xi ? f x( ) i i 為 小，這本來都是很自然的，可是計算不太方便。所以通常要求： 1 1 δ 2 =??∑m δi2 ??2 =??∑m [?( )xi ? f x( )i ]2 ??2 ? i=0 ? ? i=0 ? 或者 m m δ 22=∑ ∑δi2 = [?(xi ) ? f x( i ) ]2 i=0 i=0 為 小。這種要求誤差平方和 小的擬合稱為曲線擬合的 小二乘法。 2.2多

20、變量數(shù)據(jù)擬合 對于給定的一組數(shù)據(jù) (xi , yi ) ，i=1，2…，m ，尋求做n次多項式 n y =∑a xk k k=0 使性能指標 J a a( 0, 1,L,an ) = (yi ? a xk ik )2 為 小。 i=1 k=0 由于性能指標J可以被看做關(guān)于ak ，k=0，1，…，n的多元函數(shù)，故上述擬合多項式的構(gòu)造問題可轉(zhuǎn)化為多元函數(shù)的極值問題。令 ?J = 0 從而有正則方程組 ? m ∑ ∑xi xi2 ? 3 ?∑ ∑ ∑xi xi2 xi ? ? M M M L L L ∑xin ?? ?a0 ?∑ yi

21、? n+1?? ?a1 ??∑xi yi ?? ∑xi ?? ?? ?M =? M ?? ? M ?? ? ? ?ak ?? xin xin+1 xin+2 L ∑xi2n ??? ?an ??∑xin yi ?? ∑ ∑ ∑ 對多變量（或稱多元）線性模型 y* = a0 +a x1 1 +a x2 2 +L+a xn n 進行了m次觀測 ?y1* = a0 +a x1 11 +a x2 21 +L+a xn1 n ??y2* = a0 +a x1 12 +a x2 22 +L+a xn2 n ? ? M ??yn* = a0 +a x1 1n +a x2 2

22、n +L+a xnn n 這個稱為回歸方程組，寫成矩陣形式 y=Aα ?y1* ? ?1 x11 ? * ? ? 其中 y=? 2 ? , A=?1 x12 y ?M ? ?M M ? * ? ?1 x1m ??ym ?? ? x21 L xn1 ? ?a0 ? x22 L xn2 ??, α=??a1 ??。 M M M ? ?M ? x2m L xnm ?? ??an ?? 當m > n 時，要確定一組a ii , = 0,1,L,n, 使之精確地滿足 m 個方程，這是超定方程組的問題，只能在 小平方誤差的基礎(chǔ)上確定αi 。 定義殘差向量δ= [δ

23、δ δ1, 2,L, m ]T , 則 δ= y-Aα 其中 y = [y y1, 2,L, ym ]T 代替輸出向量。取性能指標 J =δ δT 使之 小，以此確定出α。由 J =δ δT = (y ? Aα) (T y ? Aα) =y yT ?αTAy ? yTAα+αTA AT α 利用向量和矩陣的運算公式，有 A AT α=AT y 此即為正則方程組，當A AT 非奇異時，可求得 α = (A AT )?1AT y 3.流程圖 n i z y x i i i , , 2 , 1 , , L = 輸入 n ? 5 n

24、i AT y AT x AT i i i i i ,..., 2 , 1 ; ; 1 3 2 1 = ? ? ? k ? 1 k n i i ki kj n i ji ki Y z a j X a a ? = ? ∑ ∑ = = 1 1 3 , 2 , 1 3 ? = k k k ? + 1 ≠ = 輸出參數(shù) 4.輸出結(jié)果 5.結(jié)果分析 曲線擬合的 小二乘法是反映所給數(shù)據(jù)點的總的趨勢，并不是嚴格的通過每個數(shù)據(jù)點，這樣就避免了大量數(shù)據(jù)插值時需要

25、高次多項式，同時又去掉了數(shù)據(jù)所含的測量誤差。 6.心得體會 整理思路越來越熟練了，所以執(zhí)行各個步驟也相對簡單了很多。另外對原理也加深了認識。 附錄： 1.造倒數(shù)表 1>源程序： #include"stdio.h" #include"math.h" #define N 30 void main() { int i; float x[N],c; FILE *fp1,*fp2; fp1=fopen("input1.txt","r"); fp2=fopen("output1.txt","w"); fscanf(fp1,"%f",&c); fscanf(fp1,"%f",&x[0])

26、;//初值 fprintf(fp2," ****倒數(shù)表****\n"); for(i=0;i

27、fp1); fclose(fp2); } 2>輸入文件："input1.txt" 18 0.1 3>輸出文件：“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 計算結(jié)果： 1/18.000000=0.056 2.例 3-4 1>源程序： #include"iostream" #include"cmath" using namespace std; #define N 10 void main() { int

28、i,j,k,l,n; float b[N],a[N][N],t,d,det=1.0; //***數(shù)據(jù)輸入*/ FILE *fp1,*fp2; fp1=fopen("input2.txt","r"); fp2=fopen("output2.txt","w"); fscanf(fp1,"%d",&n); for(i=0;i

29、************************************************ *************高斯消去*/ //*****************************************消元*/ //****************************列選主元函數(shù)*/ for(k=0;kfabs(d)) {d=a[i][k];i=

30、l;} } if(i==n)//判斷是否奇異，不奇異進行行交換 { if(d==0) fprintf(fp2,"奇異");//如果所有行的“首列”都為 0，為奇異 else { if(l!=k)//如果第 k 行的“首列”并不是 大 { det=det*(-1); for(j=k;j<=n;j++)//交換系數(shù)矩陣中的兩行 {t=a[l][j];a[l][j]=a[k][j];a[k][j]=t;} t=b[l];b[l]=b[k];b[k]=t;//交換右端常向量中的兩行 } } } //****************************列選主元函數(shù)*/

31、 for(i=k+1;i

32、 for(i=n-2;i>=0;i--)//從倒數(shù)第二項開始依次回代 N-1 次 { t=0; for(j=i+1;j

33、"x(%d)=%.4f\n",i+1,b[i]); for(i=0;i輸入文件：“input2.txt” 3 12 -3 3 -18 3 -1 1 1 1 15 -15 6 3>輸出文件：“output2.txt” x(1)=1.0000 x(2)=2.0000 x(3)=3.0000

34、 detA=-66.0000 3. 例 3-10 1>源程序： #include"iostream" #include"cmath" using namespace std; #define N 30 void main() { int i,j,r,k,n; float a[N][N]={0},s; //***數(shù)據(jù)輸入*/ FILE *fp1,*fp2; fp1=fopen("input3.txt","r"); fp2=fopen("output3.txt","w"); fscanf(fp1,"%d",&n); for(i=0;i

35、 fscanf(fp1,"%f",&a[i][j]); for(i=0;i

36、(j=i;j<2*n;j++)//第"2（r+1）-1"區(qū)間的變化行 { s=0.0; for(k=0;kLUY //**************************

37、*************************************************** *************LU 分解*/ //***************************************************************************** *********高斯約當法解 Ux=Y*/ //*****************************************************提取 UY 減少計算*/ for(i=1;i

38、**************************************************提取 UY 減少計算*/ //*****************************************消元*/ for(j=0;j<2*n;j++)//首行歸一化 { a[0][j]=a[0][j]/a[0][0]; } a[0][0]=1;//第一列其余行已為零 for(i=1;i

39、 for(r=0;r

40、*********** *********高斯約當法解 Ux=Y*/ //************************************數(shù)據(jù)輸出*/ fprintf(fp2,"A 的逆矩陣為\n\n"); for(i=0;i

41、輸入文件："input3.txt" 3 1 1 -1 1 2 -2 -2 1 1 3>輸出文件："output3.txt" A 的逆矩陣為 | 2.0000 -1.0000 0.0000 | | 1.5000 -0.5000 0.5000 | | 2.5000 -1.5000 0.5000 | 4.例 4-6 1>源程序： #include"stdio.h" #include"math.h" #define N 50 int n,i,j,k; float xx=0.0,yy=0.0,t,x[N],y[N],c[N],A[N]; main() { FILE *f

42、p1,*fp2; fp1=fopen("input4.txt","r"); fp2=fopen("output4.txt","w"); fscanf(fp1,"%d",&n);//n=4 for(i=0;i

43、j])*t/(x[k]-x[j]); yy=yy+t*y[k]; } fprintf(fp2,"\nx=%.7f 處的函數(shù)值為:y=%.7f",xx,yy); fclose(fp1); fclose(fp2); } 2>輸入文件："input4.txt" 4 1,0 2, 5 3, 6 4,3 2.5 3>輸出文件："output4.txt" x=2.5000000 處的函數(shù)值為:y=-6.3750000 5.習題 5-2 1>源程序： #include "stdio.h" #include "math.h" #define N 30 void main()

44、 { int i,n,k,j,l; float x[N],y[N],z[N];//定義輸入變量 float AT[3][N];//定義 A 的轉(zhuǎn)置 float X[3][3],Y[3];//定義中間變量 ATA 和 ATy float 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;i

45、 } //***數(shù)據(jù)輸入*/ // for(i=0;i

46、************************計算 //***************************************************************************** *************高斯消去*/ //*****************************************消元*/ //****************************列選主元函數(shù)*/ n=3; for(k=0;k

47、找出絕對值 大的 X[i][k]和 i 行 { if(fabs(X[i][k])>fabs(d)) {d=X[i][k];i=l;} } if(i==n)//判斷是否奇異，不奇異進行行交換 { if(d==0) fprintf(fp2,"奇異");//如果所有行的“首列”都為 0，為奇異 else { for(j=k;j<=n;j++)//交換系數(shù)矩陣中的兩行 {t=X[l][j];X[l][j]=X[k][j];X[k][j]=t;} t=Y[l];Y[l]=Y[k];Y[k]=t;//交換右端常向量中的兩行 } } //********************

48、********列選主元函數(shù)*/ for(i=k+1;i

49、];//計算 x(N)的解 for(i=n-2;i>=0;i--)//從倒數(shù)第二項開始依次回代 N-1 次 { t=0; for(j=i+1;j輸入文件："input5.txt" 5 0.1,0.2,0.58 0.2,0.3,0.63 0.4,0.5,0.73 0.6,0.7,0.83 0.9,0.8,0.92 3>輸出文件："output5.txt" 所求參數(shù)為： a=0.2000 b=0.3000 c=0.5000