《數(shù)值分析》上機(jī)實(shí)驗(yàn)報(bào)告

1、完美 WORD 格式數(shù)值分析上機(jī)實(shí)驗(yàn)報(bào)告專業(yè)知識(shí)分享完美 WORD 格式數(shù)值分析上機(jī)實(shí)驗(yàn)報(bào)告1.用 Newton法求方程X7 -X4 +14=0在（ 0.1 ，1.9 ）中的近似根（初始近似值取為區(qū)間端點(diǎn)，迭代 6 次或誤差小于 0.00001 ）。1.1 理論依據(jù)：設(shè)函數(shù)在有限區(qū)間 a ，b 上二階導(dǎo)數(shù)存在，且滿足條件1. f (x) f (b)02. f . ( x)在區(qū)間 a, b上不變號(hào)3. f . ( x) 0;| f (c) |.( x) |,其中 是中使.達(dá)到的一個(gè)4.|fca,bmir (| f(a), f(b) |)ba則對(duì)任意初始近似值x0 a,b,由迭代過(guò)程N(yùn)ewtonx

2、k 1(xk )xkf ( xk ) , k0,1,2,3f '( xk )所生的迭代序列 xk平方收斂于方程 f (x)0在區(qū)間 a,b上的惟一解令f ( x)x 728 x 414, f (0.1)0, f (1.9) 0f( x)7x 6112 x 37 x3 ( x316)0f( x)42 x 5336 x 242 x 2 ( x 38) 0f (1.9)f (1.9)0故以 1.9 為起點(diǎn)xk 1xkf (xk )f (xk )x01.9如此一次一次的迭代，逼近x的真實(shí)根。當(dāng)前后兩個(gè)的差<= 時(shí)，就認(rèn)為求出了近似的根。本程序用Newton 法求代數(shù)方程 (最高次數(shù)不大于

3、 10) 在（ a,b ）區(qū)間的根。專業(yè)知識(shí)分享完美 WORD 格式1.2C 語(yǔ)言程序原代碼：#include<stdio.h>#include<math.h>main()double x2,f,f1;double x1=1.9;/ 取初值為1.9dox2=x1;f=pow(x2,7)-28*pow(x2,4)+14;f1=7*pow(x2,6)-4*28*pow(x2,3);x1=x2-f/f1;while(fabs(x1-x2)>=0.00001|x1<0.1);/ 限制循環(huán)次數(shù)printf(" 計(jì)算結(jié)果： x=%fn",x1);1.

4、3 運(yùn)行結(jié)果：1.4MATLAB 上機(jī)程序functiony=Newton(f,df,x0,eps,M)專業(yè)知識(shí)分享完美 WORD 格式d=0;for k=1:Mif feval(df,x0)=0d=2; breakelsex1=x0-feval(f,x0)/feval(df,x0);ende=abs(x1-x0);x0=x1;if e<=eps&&abs(feval(f,x1)<=epsd=1; breakendendif d=1y=x1;elseif d=0y='迭代 M 次失敗 'elsey= '奇異 'endfunctiony

5、=df(x)y=7*x6-28*4*x3;專業(yè)知識(shí)分享完美 WORD 格式Endfunctiony=f(x)y=x7-28*x4+14;End>> x0=1.9;>> eps=0.00001;>> M=100;>> x=Newton('f','df',x0,eps,M);>> vpa(x,7)1.5 問(wèn)題討論：1.使用此方法求方解，用誤差來(lái)控制循環(huán)迭代次數(shù)，可以在誤差允許的范圍內(nèi)得到比較理想的計(jì)算結(jié)果。 此程序的不足之處是， 所要求解的方程必須滿足上述定理的四個(gè)條件，但是第二和第四個(gè)條件在計(jì)算機(jī)上比較難

6、以實(shí)現(xiàn)。2.Newton迭代法是一個(gè)二階收斂迭代式，他的幾何意義Xi+1 是 Xi 的切線與 x 軸的交點(diǎn) ,故也稱為切線法。它是平方收斂的，但它是局部收斂的，即要求初始值與方程的根充分接近，所以在計(jì)算過(guò)程中需要先確定初始值。3.本題在理論依據(jù)部分，討論了區(qū)間 (0.1,1.9) 兩端點(diǎn)是否能作為 Newton 迭代的初值，結(jié)果發(fā)現(xiàn) 0.1 不滿足條件，而 1.9 滿足，能作為初值。另外，該程序簡(jiǎn)單，只有一個(gè)循環(huán)，且為順序結(jié)構(gòu)，故采用 do-while 循環(huán)。當(dāng)然也可以選擇for 和 while 循環(huán)。專業(yè)知識(shí)分享完美 WORD 格式專業(yè)知識(shí)分享完美 WORD 格式2.已知函數(shù)值如下表：x12

7、345f(x)00.69314711.0986121.3862941.6094378348x678910f(x)1.7917591.94591012.0794452.1972242.302585561f (x)f (1)=1f (10)=0.1試用三次樣條插值求f(4.563) 及 f(4.563) 的近似值。2.1 理論依據(jù)(xjx)3( x x j 1)3h2j 1xj xh2j 1x xj 1S( x) M j 16hj 1M j6hj 1( yj 1M j 1 6)(hj 1) ( yj M j6)(hj 1)這里 hj 1xjx j1，所以只要求出M j ，就能得出插值函數(shù)S（x）。

8、21121求 M j 的方法為：222N1d 06 ( y1y0y0 )h0h0M 0d0M 1d12N112M NdN這里 d j6(y j 1y jy jy j1)( j 1, 2, , N 1)h j 1h jh jh j 1d N6 yN1( yNyN 1 )h N1hN1h j11h jjh jh jjjh j 1h j1專業(yè)知識(shí)分享完美 WORD 格式最終歸結(jié)為求解一個(gè)三對(duì)角陣的解。用追趕法解三對(duì)角陣的方法如下：b1c1111a2b2c2l1122Al21LUan 1bn 1cn 1n 1n 1anbnln1nd ,即 Ld , 若記1LUx,則由 Ld 得Uxn11d111x11

9、l21，n1n1l n1ndnnxnn綜上可得求解方程 Ax=d的算法：1b1 ,1d1 ,l i 1i 1 ,i 1bi 1li 1ci ,ii 1di1l i 1ii1,2,3,n1xnn,xiici xi1 ,in1, ,2,1ni2.2C 語(yǔ)言程序代碼：#include<stdio.h>#include<math.h>void main()int i,j,m,n,k,p;專業(yè)知識(shí)分享完美 WORD 格式doubleq10,p10,s4,g4,x0,x1,g0=1,g9=0.1;double s1010;double a10,b10,c10,d10,e10,x10

10、,h9,u9,r9;double f10=0,0.69314718,1.0986123,1.3862944,1.6094378,1.7917595,1.9459101,2.079445,2.1972246,2.3025851;printf(" 請(qǐng)依次輸入 xi:n");for(i=0;i<=9;i+)scanf("%lf",&ei);/ 求 h 矩陣for(n=0;n<=8;n+)hn=en+1-en;d0=6*(f1-f0)/h0-g0)/h0;d9=6*(g9-(f9-f8)/h8)/h8;for(j=0;j<=7;j+)d

11、j+1=6*(fj+2-fj+1)/hj+1-(fj+1-fj)/hj)/(hj+hj+1);for(m=1;m<=8;m+)um=hm-1/(hm-1+hm);for(k=1;k<=8;k+)rk=hk/(hk-1+hk);for(i=0;i<=9;i+)/ 求 u 矩陣for(p=0;p<=9;p+)sip=0;if(i=p)sip=2;專業(yè)知識(shí)分享完美 WORD 格式s01=1;s98=1;for(i=1;i<=8;i+)sii-1=ui;sii+1=ri;printf(" 三對(duì)角矩陣為 :n");for(i=0;i<=9;i+)f

12、or(p=0;p<=9;p+)/ 求 r 矩陣 printf("%5.2lf",sip); if(p=9) printf("n");printf(" 根據(jù)追趕法解三對(duì)角矩陣得：n");a0=s00;b0=d0;for(i=1;i<9;i+)ci=sii-1/ai-1;/ 求 d 矩陣ai=sii-si-1i*ci;bi=di-ci*bi-1;if(i=8)p10=bi;q10=ai;專業(yè)知識(shí)分享完美 WORD 格式x9=p10/q10;printf("M10=%lfn",x9);for(i=9;i>

13、;=1;i-)xi-1=(bi-1-si-1i*xi)/ai-1;printf("M%d=%lfn",i,xi-1);printf(" 可得 s(x)在區(qū)間 4,5 上的表達(dá)式 ;n");printf(" 將 x=4.563 代入得： n");x0=5-4.563;x1=4.563-4;s4=x3*pow(x0,3)/6+x4*pow(x1,3)/6+(f3-x3/6)*(5-4.563)+(f4-x4/6)*(4.563-4);g4=-x3*pow(x0,2)/2+x4*pow(x1,2)/2-(f3-x3/6)+(f4-x4/6)

14、;printf("計(jì) 算 結(jié) 果 ： f(4.563)的 函 數(shù) 值 是 ： %lfnf(4.563)的 導(dǎo) 數(shù) 值是： %lfn",s4,g4);2.3 運(yùn)行結(jié)果：專業(yè)知識(shí)分享完美 WORD 格式2.4 問(wèn)題討論1. 三次樣條插值效果比 Lagrange 插值好，沒(méi)有 Runge 現(xiàn)象，光滑性較好。2. 本題的對(duì)任意劃分的三彎矩插值法可以解決非等距節(jié)點(diǎn)的一般性問(wèn)題。3. 編程過(guò)程中由于定義的數(shù)組比較多，需要仔細(xì)弄清楚各數(shù)組所代表的參數(shù)，要注意各下標(biāo)代表的含義，特別是在用追趕法計(jì)算的過(guò)程中。專業(yè)知識(shí)分享完美 WORD 格式3.用 Romberg算法求33x x1.4 (5x

15、 7) sin x2 dx(允許誤差0.00001) .13.1 理論依據(jù)：Romberg算法的計(jì)算步驟如下：（1）先求出按梯形公式所得的積分值T1(0)b a f ( a)f (b)2（2）把區(qū)間 2 等分，求出兩個(gè)小梯形面積之和，記為T1(1) ，即T1(1) b a f (a)f (b)2 f ( ba )421這樣由外推法可得 T (0)T1(1)()2 T1(0)(1)(0), T2(0)24T1T1。21(1)2412（3）把區(qū)間再等分（即 22等分），得復(fù)化梯形公式 T1(2)，由 T1(1) 與 T1(2) 外推可(1)4T1(2)T1(1)(0) 42 T2(1)T2(0)，

16、如此，若已算出 2k等分的復(fù)化梯形公式得 T241，T3421T1( k ) ，則由 Richardson外推法，構(gòu)造新序列( k 1) 4mTm(k )Tm(k 1)Tm 1， m=1,2, ,l, k=1,2, ,l-m+1,4m1最后求得 Tl (0)1 。（4）Tl (0)Tl (0)1 或 |Tl (0)Tl (0)1 |<就停止計(jì)算，否則回到（ 3），計(jì)算 T1( l 1) ，一般可用如下算法：T1(0)b a f (a)f (b)2T (l )1 T (l 1)ba 2l1f a(2i1) ba1212l1i12lT ( k 1)4m Tm(k)Tm(k1),m1,2,l

17、,k 1,2, , l m 1m 14m1專業(yè)知識(shí)分享完美 WORD 格式其具體流程如下，并全部存入第一列(1)T1(0)( 3T)2( 0 )( 6T)3( 0 )(2) T1(1)( 5T)2( 1 )(4) T1(2)通常計(jì)算時(shí)，用固定l=N 來(lái)計(jì)算，一般 l=4 或 5 即能達(dá)到要求。3.2C 語(yǔ)言程序代碼：#include<math.h>#include<stdio.h>double f(double x)/ 計(jì)算 f(x) 的值double z;z=pow(3,x)*pow(x,1.4)*(5*x+7)*sin(x*x);return(z);main() d

18、ouble t2020,s,e=0.00001,a=1,b=3; int i,j,l,k;t01=(b-a)*(f(b)+f(a)/2;/ 下為 romberg算法t11=(b-a)*(f(b)+2*f(b+a)/2)+f(a)/4;t02=(a*t11-t01)/(4-1);j=3;for(l=2;fabs(t0j-1-t0j-2)>=e;l+)for(k=1,s=0;k<=pow(2,l-1);k+)專業(yè)知識(shí)分享完美 WORD 格式s+=f(a+(2*k-1)*(b-a)/pow(2,l);/判斷前后兩次所得的T(0)的差是否符合要求，如果符合精度要求則停止循環(huán)tl1=(tl-

19、11+(b-a)*s/pow(2,l-1)/2;for(i=l-1,j=2;i>=0;i-,j+)tij=(pow(4,j-1)*ti+1j-1-tij-1)/(pow(4,j-1)-1);if(t01<e)printf("t=%0.6fn",t01);elseprintf("用Romberg算 法 計(jì) 算 函 數(shù) 所 得 近 似 結(jié) 果 為 ：nf(x)=%0.6fn",t0j-1);3.3 運(yùn)行結(jié)果：3.4MATLAB 上機(jī)程序functionT,n=mromb(f,a,b,eps)if nargin<4,eps=1e-6;endh

20、=b-a;R(1,1)=(h/2)*(feval(f,a)+feval(f,b);n=1;J=0;err=1;while(err>eps)專業(yè)知識(shí)分享完美 WORD 格式J=J+1;h=h/2;S=0;for i=1:nx=a+h*(2*i-1);S=S+feval(f,x);endR(J+1,1)=R(J,1)/2+h*S;for k=1:JR(J+1,k+1)=(4k*R(J+1,k)-R(J,k)/(4k-1);enderr=abs(R(J+1,J+1)-R(J+1,J);n=2*n;endR;T=R(J+1,J+1);format longf=(x)(3.x)*(x.1.4)*(

21、5*x+7)*sin(x*x);T,n=mromb(f,1,3,1.e-5)3.5 問(wèn)題討論：1.Romberge算法的優(yōu)點(diǎn)是 :把積分化為代數(shù)運(yùn)算 ,而實(shí)際上只需求 T1(i) ,以后用遞推可得 .算法簡(jiǎn)單且收斂速度快 ,一般 4 或 5 次即能達(dá)到要求。2.Romberge 算法的缺點(diǎn)是 :對(duì)函數(shù)的光滑性要求較高， 計(jì)算新分點(diǎn)的值時(shí)，這些數(shù)值的個(gè)數(shù)成倍增加。3. 該程序較為復(fù)雜，涉及函數(shù)定義，有循環(huán)，而且循環(huán)中又有判斷，編寫時(shí)需要注意該判斷條件是處于循環(huán)中， 當(dāng)達(dá)到要求時(shí)跳出循環(huán)，終止運(yùn)算。4. 函數(shù)的定義可放在主函數(shù)前也可在主程序后面。 本程序采用的專業(yè)知識(shí)分享完美 WORD 格式后置方

22、式。專業(yè)知識(shí)分享完美 WORD 格式4. 用定步長(zhǎng)四階 Runge-Kutta求解dy1 / dt1dy2 / dty3dy3 / dt1000 1000 y2100 y3y1 (0)0y2 (0)0y3 (0)0h=0.0005 ，打印 yi(0.025)， yi(0.045)， yi(0.085)， yi(0.1)，(i=1 ，2，3)4.1 理論依據(jù)：Runge_Kutta 采用高階單步法，這里不是先按 Taylor 公式展開，而是先寫成 tn 處附近的值的線性組合（有待定常數(shù)）再按 Taylor公式展開，然后確定待定常數(shù)， 這就是 Runge-Kutta法的思想方法。本題采用四階古典的

23、Runge-Kutta公式：Yn 1Yn K1 3K 23K 3K4/8K1hF ( xn ,Yn )K2hF ( xnh / 3, YnhK1 / 3)K3hF( xn2h / 3,YnhK1 / 3hK 2 )K4hF ( xnh,Yn hK1hK 2hK 3)4.2C 語(yǔ)言程序代碼：#include<stdio.h>void fun(double x4,double y4,double h)y1=1*h;y2=x3*h;y3=(1000-1000*x2-100*x2-100*x3)*h;/ 微分方程向量函數(shù) void main()專業(yè)知識(shí)分享完美 WORD 格式 double

24、 Y54,K54,m,z4,e=0.0005; double y5=0,0.025,0.045,0.085,0.1; int i,j,k;for(i=1;i<=3;i+)Y1i=0;for(i=1;i<=4;i+)for(j=1;j<=3;j+)Kij=0;for(k=1;k<=5;k+)for(m=yk-1;m<=yk;m=m+e)for(i=1;i<=3;i+)zi=Yki;fun(z,K1,e);for(i=1;i<=3;i+)zi=Yki+e*K2i/2;/ 依此求 K1,K2K3 的值fun(z,K2,e);for(i=1;i<=3;i

25、+)zi=Yki+e*K2i/2;fun(z,K3,e);for(i=1;i<=3;i+)zi=Yki+e*K3i;fun(z,K4,e);專業(yè)知識(shí)分享完美 WORD 格式for(i=1;i<=3;i+)Yki=Yki+(K1i+2*K2i+2*K3i+K4i)/6;/求YiN+1 的值 if(k!=5)for(i=1;i<=3;i+)Yk+1i=Yki;printf(" 計(jì)算結(jié)果 :n");for(i=1;i<5;i+)for(j=1;j<=3;j+)printf("y%d%4.3f=%-10.8f,",j,yi,Yij)

26、;if(j=3)printf("n");printf("n");4.3 運(yùn)行結(jié)果：專業(yè)知識(shí)分享完美 WORD 格式4.4 問(wèn)題討論：1.定步長(zhǎng)四階 Runge-kutta方法是一種高階單步法法穩(wěn)定，精度較高，誤差小且程序相對(duì)簡(jiǎn)單， 存儲(chǔ)量少。不必求出起始點(diǎn)的函數(shù)值，可根據(jù)精度的要求修改步長(zhǎng)，不會(huì)由于起始點(diǎn)的誤差造成病態(tài)。2. 本程序可以通過(guò)修改主程序所調(diào)用的函數(shù)中的表達(dá)式來(lái)實(shí)現(xiàn)對(duì)其它函數(shù)的任意初值條件求微分計(jì)算。3.程序中運(yùn)用了大量的for 循環(huán)語(yǔ)句，因?yàn)樵摴街猩婕按罅康那蠛?，且有不同的函?shù)和對(duì)不同的數(shù)值求值，編程稍顯繁瑣。所以編寫過(guò)程中一定要注意各循

27、環(huán)的次數(shù)，以免出錯(cuò)。專業(yè)知識(shí)分享完美 WORD 格式5.12.384122.115237-1.0610741.112336- 0.113584 0.7187191.7423823.067813- 2.0317432.11523719.141823- 3.125432- 1.012345 2.189736 1.563849 - 0.7841651.1123483.123124- 1.061074- 3.12543215.5679143.1238482.0314541.836742- 1.0567810.336993-1.0101031.112336- 1.0123453.12384827.108

28、4374.101011- 3.7418562.101023- 0.71828- 0.037585A - 0.113584 2.189736 2.031454 4.101011 19.897918 0.431637 - 3.111223 2.121314 1.784317 0.718719 1.563849 1.836742 - 3.741856 0.431637 9.789365 - 0.103458 - 1.103456 0.2384171.742382- 0.784165-1.0567812.101023- 3.111223- 0.103458 14.71384653.123789- 2.

29、2134743.0678131.1123480.336993- 0.718282.121314-1.1034563.12378930.7193344.446782- 2.031743 3.123124-1.010103- 0.0375851.7843170.238417- 2.2134744.44678240.00001b (2.187436933.992318- 25.1734170.846716951.784317- 86.6123431.11012304.719345 - 5.6784392 )T用列主元消去法求解Ax=b 。5.1 理論依據(jù)：列主元素消元法是在應(yīng)用 Gauss 消元法的

30、基礎(chǔ)上，憑借長(zhǎng)期經(jīng)驗(yàn)積累提出的， 是線性方程組一般解法， 目的是為避免在消元計(jì)算中使誤差的擴(kuò)大， 甚至嚴(yán)重?fù)p失了有效數(shù)字使數(shù)據(jù)失真， 而在每次初等變換前對(duì)矩陣作恰當(dāng)?shù)恼{(diào)整，以提高 Gauss 消元法的數(shù)字穩(wěn)定性，進(jìn)而提高計(jì)算所得數(shù)據(jù)的精確度。 即在每主列中取絕對(duì)值最大的元素作主元，再做對(duì)應(yīng)的行交換然后消元求解的辦法。具體做法如下：將方陣 A 和向量 b 寫成 C= （A，b ）。將 C 的第 1 列中第 1 行的元素與其下面的此列的元素逐一進(jìn)行比較， 找到最大的元素 cj1 ，將第 j 行的元素與第 1 行的元素進(jìn)行交換，然后通過(guò)行變換，將第 1 列中第 2 到第 n 個(gè)元素都消成 0。將變換

31、后的矩陣 C(1) 的第二列中第二行的元素與其下面的此列的元素逐一進(jìn)行比較，找到最大的元素ck(1)2 ，將第 k 行的元素與第 2 行的元素進(jìn)行交換，然后通過(guò)行變換，將第 2 列中第 3 到第 n 個(gè)元素都消成 0。以此方法將矩陣的左下部分全都消專業(yè)知識(shí)分享完美 WORD 格式成 0 后再求解。最終形式如下：a( n )a ( n )g1111 n0(A,b ）00ann(n ) gn5.2C 語(yǔ)言程序代碼（ 1）比較該列的元素的絕對(duì)值的大小，將絕對(duì)值最大的元素通過(guò)行變換使其位于主對(duì)角線上；（ 2）進(jìn)行高斯消去法變換，把系數(shù)矩陣化成上三角形，然后回代求#include "math.

32、h"#include "stdio.h"void Householder(double A99);void expunction(double A99,double b9,double x9);void main()double A99=12.38412,2.115237,-1.061074,1.112336,-0.113584,0.718719,1.742382,3.067813,-2.031743,2.115237,19.141823,-3.125432,-1.012345,2.189736,1.563849,-0.784165,1.112348,3.1231

33、24,-1.061074,-3.125432,15.567914,3.123848,2.031454,1.836742,-1.056781,0.336993,-1.010103,1.112336,-1.012345,3.123848,27.108437,4.101011,-3.741856,2.101023,-0.71828,-0.037585,專業(yè)知識(shí)分享完美 WORD 格式-0.113584,2.189736,2.031454,4.101011,19.897918,0.431637,-3.111223,2.121314,1.784317,0.718719,1.563849,1.836742

34、,-3.741856,0.431637,9.789365,-0.103458,-1.103456,0.238417,1.742382,-0.784165,-1.056781,2.101023,-3.111223,-0.103458,14.713847,3.123789,-2.213474,3.067813,1.112348,0.336993,-0.71828,2.121314,-1.103456,3.123789,30.719334,4.446782,-2.031743,3.123124,-1.010103,-0.037585,1.784317,0.238417,-2.213474,4.446

35、782,40.00001;double b9=2.1874369,33.992318,-25.173417,0.84671695,1.784317,-86.612343,1.1101230,4.719345,-5.6784392;double x9=0.0;int i,j;Householder(A);printf("n The Results of X are:n");expunction(A,b,x);for(i=1;i<10;i+)printf("X%1d=%fn",i,xi-1);void Householder(double A99)專業(yè)

36、知識(shí)分享完美 WORD 格式double q9,u9,y9,s,a,kr;int i,j,k;for(i=0;i<7;i+)s=0;for(j=i+1;j<9;j+)s+=Aji*Aji;s=sqrt(s);a=s*s+fabs(Ai+1i)*s;for(j=0;j<9;j+)if(j<=i) uj=0;else if(j=i+1) uj=Aji+Aji/fabs(Aji)*s;else if(j>i+1) uj=Aji;for(k=0;k<9;k+)yk=0;for(j=0;j<9;j+)yk+=Akj*uj;yk/=a;kr=0;for(k=0;k<9;k+)kr+=yk*uk;kr/=2*a;for(k=0;k<