貴大數(shù)值分析上機(jī)實(shí)習(xí)報(bào)告

1、數(shù)值分析上機(jī)實(shí)習(xí)報(bào)告學(xué)院：土建學(xué)院 專業(yè):學(xué)號(hào): 姓名:上機(jī)實(shí)習(xí)題一、題目：已知a與b(12. 38412, 2. 115237,-1. 061074, 1. 112336, -0. 113584, 0. 718719, 1. 742382, 3. 067813,-2. 0317432. 115237, 19. 141823,-3. 125432,-1. 012345, 2. 189736, 1. 563849, -0. 784165, 1. 11234& 3. 123124 -1. 061074, -3. 125432, 15. 567914, 3. 12384& 2. 03

2、1454, 1. 836742,-1. 056781, 0. 336993,-1. 0101031. 112336,-1.012345, 3. 123848,27. 108437, 4. 101011,-3. 741856, 2. 101023, -0.71828, -0. 037585a彳 -0. 113584, 2. 189736, 2. 031454, 4. 101011, 19. 89791 & 0. 431637,-3. 111223, 2. 121314, 1. 7841370. 718719, 1. 563849, 1. 836742, -3. 741856, 0. 43

3、1637, 9. 789365, -0. 103458,-1. 103456, 0. 2384171. 742382, -0. 784165,-1. 056781, 2. 101023, -3. 111223, -0. 103458, 14. 7138465, 3. 123789,-2. 2134743. 067813, 1. 11234& 0. 336993, -0. 7182& 2. 121314,-1. 103456, 3. 123789, 30. 719334, 4. 446782(-2. 031743, 3. 123124,-1. 010103,-0. 037585,

4、 1. 784317, 0. 238417,-2. 213474, 4. 446782, 40. 00001b=2. 1874369, 33. 992318,-25. 173417, 0. 84671695, 1. 784317,-86. 612343, 1. 1101230, 4. 719345,-5. 67843921. 用household變換，把a(bǔ)化為三對(duì)角陣b (并打印b)。2. 用超松弛法求解bx=b (取松弛因子3二1.4,対°)二0,迭代9次)。3. 用列主元素消去法求解bx=bo二、解題方法的理論依據(jù):1、用householder變換的理論依據(jù)(1)令 a0=a,

5、a(i j) l=a(i j),已知 ar_l 即 ar_l=a(i j)r(2 ) sr=sqrt(pow(a, 2)(3 ) a (r) =sr*sr+abs (a (r+1, r) *sr(4 ) y(r)=a(r_l)*u®/a®(5)kr二(/2)*ur 的轉(zhuǎn)置*yr/a®(6 ) qr=yr-kr*ur(7 ) ar=a(r-l)-(qr*ur 的轉(zhuǎn)置+ur*qr 的轉(zhuǎn)置)r=l, 2,，n22、用超松弛法求解其基本思想：在高斯方法己求出x<m),的基礎(chǔ)上，組合新的序列，從而加快收斂速度。其算式：bii -1 x z -1 + bii -xi +

6、 bii +1 x / +1 = bi卻“匹心勿一1 一処也x0l +如 biibiibii<xi二 w文i + x0ix0i = xi其中3是超松弛因子，當(dāng)3>1時(shí)，可以加快收斂速度3、用消去法求解用追趕消去法求bx二b的方法：2h + 1tb5s c迂 + 二 hb5h +二qlovo " u1ovo -pen c1s/&1e+ alsdlsh i 丄2:.&msn (d 一s 巳sms)/(£s+ a 一s s 二)9 i n 一 2：9x9vul9§ h gls + 二 + h 8z:jr4犧諭曲“藝nclude = maul.

7、 it-include =stdio.h=define ge xvoid main。im sign(doubox)；double a=9n(一 2.3x4122 二 523711061074" 1 1123361p二 3500407100719 j .74230023067001312 031743 =(2二 5237j 9一±8233 125432二-.0123452189736l5638496.784165r二 23483123124 廠 .1.061074312543215.5679143123848203145418367421.056700103369931.01

8、003廠 p112336工01234531238482708437400lr37418562120236.718286037585l 0二 35oc421oc9736 2 0314544ss二9oc9791oc.0.431637:3一二 2232 12131417*4317 廠 (0.718719 一 .563849" 1oo36742j3.7418560.43163797893656. 10345811 1034560.238417 (1.7423826.784165 工.056781210323:3 三 22310. - 03458" 14.71300463123789

9、:2.213474 -" -3o67813a112348o.336993o.718282121314l1o345631237893o.719334 4446782l (203 一 74331231241 -.0-0-0310.037585.1.7843170.23841722 一 34744446782.40.0000 二xdoublerhsw-inp.jsmg;double 三9>二 9ws互 9虧二 9=10一>9一；double bs=(2 一 x74369.33.9923 一 孑 25 173417ooc4671695.7x43176c6.612343j 二

10、63;23p4.7 - 9345 二 5.67oc4392 - for<h)；a7+j) household 煨滯芒sh0ofor(ij!.+l;a9;+i)shs+m=jpa 三三；shsqn(s)八haau+二 uv05s 詫 s+s*au±=j>?(s*ss 詫呂+二三)八for(gupga9+g)if (g<=j) ug=o;else if (g=j+l)ug=a|j+lj+s*sign(a|j+lj);else ug=agj;for(m=0;m<9;+m)ym=o;for(n=0;n<9;+n)ym=ym4-amfnl*un;ym=yml/h;

11、k=0;for(i=0;i<9;+i)k=k+ufi*yi;k=0.5*k/h;for(i=0;i<9;+i)qi=yi-k*ui;for(n=0;n<9;+n)for(m=0;m<9;4-+m) amn=amn-(qm*un+um*qn);printf(,'household:n,');for(i=0;i<9;+i)for(j=0;j<9 ;+j)if(j%9=0)pnntf(unh);printf(n%-9.5r,aij);printf(nnh);w=1.4;/*超松弛法*/for(i=0;i<9;i+)xli=o;for(i=0;

12、i<9;i+)for(j=0;j<9;j+)if(i=j)bliul=o； elseblij=-aiu/aii;for(i=0;i<9;i+)bli9=bi/aii;for(n=0;n<9;n+) for(i=0;i<9;i+)s=0;for(j=0;j<9;j+)s=s+blij*xirj;s=s+bli9;x 1 i=x 1 i *( l-w)+w*s; for(i=0;i<9;i+)if(i=5) printf(”n“)； printf(hx%d=%-10.6f',i,x 1 i);嚴(yán)以下是消去法*/ printfcxn'1);u

13、o=aoo;yo=bo; for(i=l;i<9;i+)qi=airi-l/ui-l; ui=aii-qi*ai-li; yi=bi-qi*yi-l;xge=yge/uge;fbr(i=ge-1 ;i>=o;i-)xi=(yi-aii+l*xi+l)/ui;for(i=0;i<9;i+)if(i=5)printf(hnn);printfc' x%d=%-10.6f*,i,xi);)int sign(double x)int z;z=(x>=(le-6)?l:-l);return(z);2.打印結(jié)果：household:12.38412-4.89308 0.000

14、00 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000-4.89308 25.39842 6.494100.00000 0.000000.000000.000000.000000.000000.000006.4941020.61150 8.24393 0.000000.000000.000000.000000.000000.000000.00000&2439323.42284-13.880070.000000.000000.00000-0.000000.000000.000000.00000-13.8800729.69828 4.53450

15、0.000000.000000.000000.000000.000000.000000.00000 4.5345016.00612 4.881440.000000.000000.000000.000000.00000().00000 0.000004.8814426.01332 -4.50363-0.000000.000000.000000.000000.00000 0.000000.00000450363 21.25406 4.504500.000000.000000.00000-0.00000 0.000000.00000-0.00000 4.5045014.53412x0= 1.0734

16、09 x 1=2.272579x2= -2.856601x3 二2.292514 x4=2.112165x5= -6.422586x6= 1.357802 x7=0.634259x8= -0.587042x0=l075799x3=2.293099x6= 1.357923x 1=2.275744x4=2.112634x7=0.634244x2= 2855515x5= -6.423838x8= -0.587266四.問題討論:此程序具有很好的通用性。在gs方法的基礎(chǔ)上，已經(jīng)求出x的第m解，第m-1基礎(chǔ)上，經(jīng)過重新組合 得新的序列，而在此新序列收斂速度加快。上機(jī)實(shí)習(xí)題二用對(duì)分法求b（即a）的小于23

17、且最接近23的特征值的近似值（誤差小于0.1）在用反幕法求b的該特征值的更精 確近似值及相應(yīng)的特征向量.二. 程序算法1. 對(duì)分法:對(duì)于對(duì)稱的三對(duì)角方陣，其特征值必落在卜ii b ii , ii b ii沖，其中ii b ii是其任意范數(shù).通過同號(hào)數(shù)的計(jì) 算,可知在0, ii b ii沖特征根的個(gè)數(shù)，在不斷對(duì)分，就可得特征根所在的區(qū)間，如在ik,uk中有一根入k,且lk,uk很小, 就可取lk+uk/2作為ak的近似值.2. 反幕法:基于乘幕法屮求不太容易,因此不直接用a-1xk=xk+1,而使用axk+1=xk.用求解線性方程組求出 xk+1,再標(biāo)準(zhǔn)化，就可得到絕對(duì)值最小的特征值與相應(yīng)的特征

18、向量,這就是反代法（也稱反幕法）.三. 程序清單1. 對(duì)分法：#include<stdio.h>#include<math.h>#define n 9int f（float x）int n=0,i;double sn,w;floatbnn=12.384120,-4.893077,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000, -4.893077,25.398416,6.494097,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0

19、.000000,6.494097,20.611499,8.243925,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,8.243925,23.422838, 13.880071,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,/3.880071,29.698278,4.534502,0.000000,0.000000,0.000000, 0.000000,0.000000,0.000000,0.000000,4.534502,16.006

20、1 仃,4.881435,0.000000,0.000000, 0.000000,0.000000,0.000000,0.000000,0.000000,4.881435,26.013315,4503635,0.000000, 0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,-4.503635,21.254061,4.504498, 0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,4.504498,14.534122;s0=b00-x;s1=(b11-x)-

21、b01*b01/(b00-x);for(i=2;i<n;i+)w=si-1*si-2;if(w!=0)si=bii-x-bi-1i*bi-1i/si-1;if(s1=0)si=bii-x;if(si-1=0)si=-1;for(i=0;i<n;i+)if(si>=0)n+;retum(n);void main()int f(float x);int m=0;double x 1=0.0,x2=23.0,t1;while(fabs(x2-x1)>0.1)t1= (x1+x2)/2.0;if(f(t1)-f(x2)>=1)x 仁i;if(f(t1)-f(x2)=0)x

22、2=t1;if(f(t1)-f(x2)>=0)t1= (x1+x2)/2.0;m+;printfc't 1=%fnh,t1);printf(nm=%dn",m);2. 反幕法：#include<stdio.h>#include<math.h>#define n 9main()float bnn=12.384120,-4.893077,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000, -4.893077,25.398416,6.494097,0.000000,0.0000

23、00,0.000000,0.000000,0.000000,0.000000, 0.000000,6.494097,20.611499,8.243925,0.000000,0.000000,0.000000,0.000000,0.000000, 0.000000,0.000000,8.243925,23.422838,3.880071,0.000000,0.000000,0.000000,0.000000, 0.000000,0.000000,0.000000,3.880071,29.698278,4.534502,0.000000,0.000000,0.000000, 0.000000,0.

24、000000,0.000000,0.000000,4.534502,16.006117,4.881435,0.000000,0.000000, 0.000000,0.000000,0.000000,0.000000,0.000000,4.881435,26.013315,-4.503635,0.000000, 0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,4503635,21.254061,4.504498, 0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000

25、,4.504498,44.534122;int i,j,k;float max,t1 ,t2,qn,un,yn,xn;float zn=1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0;t仁21.876953;/*對(duì)分法所求得的特征值*/for(i=0;i<n;i+)bii=bii-t1;for(k=1;k<=n+1;k+)q0=-b01/b00； u0=z0/b00;for(i=1;i<n-1;i+)qi=-bii+1/(bii+bii-1*qi-1);for(i=1;i<n;i+) ui=(zibi/ru/)/(bii+bi/rq/);yn-1

26、=un-1;for(i=7;i>=0;i-)yi=qi*yi+1+ui;max=0;for(i=0;i<n;i+)if(fabs(max)<=fabs(yi)max=yi;for(i=0;i<n;i+) zi=yi/max;t2=t1+1.0/max;printf(,'t2=%fnh,t2);for(i=0;i<n;i+)xi=zi;for(i=0;i<n;i+)printf(,x%d=%9.6fn,i,zi);四. 運(yùn)行結(jié)果1. 對(duì)分法求強(qiáng)特征值為：t1=21.9183652. 用反幕法求特征值為：t2=21.918293特征向量為：x 0=0.1

27、57189x1=-0.306283x2=0.282570x3=0.286065x4=0.198838x 5=0.534490x6=0.462646x7=1.000000刈 8=0.610016五. 問題討論1對(duì)分法簡(jiǎn)單可靠，數(shù)值穩(wěn)定性較高，對(duì)于求少量幾個(gè)特征值特別適宜.但收斂速度較慢.2.反幕法是結(jié)合對(duì)分法使用的，所以近似值較恰當(dāng)，其收斂速度是驚人的只耍迭代兩次就可得到較 滿意的結(jié)果但在運(yùn)用中需把一般矩陣化為hessebberg陣才可計(jì)算.上機(jī)實(shí)習(xí)題三一、題目：己知函數(shù)值如下表:x19345f(x)00.69314781.09861231.38629441.6094378x678910f(x)

28、1.7917951.94591012.0794452.19722462.3025851f'(x)f(l)=lf(0)=0.1試用三次樣條插值求f(4.563)及f (4.563)的近似值。二、解題方法的理論依據(jù)：任意劃分的三彎矩插值法以及方程組解法中的三對(duì)角陣追趕算法。應(yīng)用三次樣條插值法能夠?qū)瘮?shù)產(chǎn)生很好的逼近效果。而追趕算法又具有計(jì)算量少、方法簡(jiǎn)單、算法穩(wěn)定的特 點(diǎn)。方法應(yīng)用條件：適用于求復(fù)雜函數(shù)在給定區(qū)間內(nèi)某一點(diǎn)的函數(shù)值，給出函數(shù)f(x)在區(qū)間a,b中的n個(gè)插值點(diǎn)，并 且給出函數(shù)在區(qū)間端點(diǎn)處的值。三* 1計(jì)算程序:#include nstdio.hn#include nmath.h

29、m#define n 11#define ge 11void main()int i,m;double e,s,e,q 12,u 12,y 12,c 12,w 12;double b12=2,0,4.15888308,6.5916738,8.3177664,9.6566268,10.750557,11.6754606,12.47667,13.1833476,13.8155106,14.0155106;double 班1212=2,4,0,0,0,0,0,0,0,0,0,0,0;ann-l=4;ann=2;for(i=l;i<ll;i+)aii-l=l；ariiril=4;aii+l=l;

30、u0=a00;/*消去法求 ci*/y0=b0;for(i=l;i<12;i+)qi=aii-l/ui-l;ui=aii-qi*ai-li;yi=bi-qi*yi-l;cge=yge/uge;for(i=ge-1 ;i>=0;i)ci=(y i-ai i+1 *ci+1 )/ui;printf(”請(qǐng)輸入要插的值:”);scanf(”f,&e);for(i=0;i<12;i+)e=fabs(e-i);if(e>=2)wil=0;else if(e<=l)wi=0.5*fabs(e*e*e)-e*e+2.0/3.0;elsewlij=(-1.0/6.0)*fa

31、bs(e*e*e)+e*e-2*fabs(e)+4.0/3.0;s=0.0;for(i=0;i<12;i+)s=s+ci*wi;printf(,f(%10=%lf',e,s);printf(,nh);printf("請(qǐng)輸入要求的導(dǎo)數(shù)的值:”)；scanf(”d”,&m);printf(nr(%d)=%lfn",m,(cm+l-cm-l)/2.0);輸出結(jié)果：請(qǐng)輸入要插的值：4.563f(4.563)= 1.517932請(qǐng)輸入要求的導(dǎo)數(shù)的值：4.563f，(4.563)= 0.249350四、問題討論：在給均勻分劃的插值函數(shù)x賦值時(shí)，由于使用for循環(huán)，

32、誤將xi二i + 1寫成xi二i,導(dǎo)致運(yùn)算錯(cuò)誤。此程 序具有一定通用性，對(duì)于任意劃分的三彎矩插值法，只許改動(dòng)可。求解方程組mj時(shí)，要用到三對(duì)角方程 組的追趕法(也稱thomas算法)。變量較多，應(yīng)注意區(qū)分。求導(dǎo)時(shí)注意正負(fù)號(hào)。上機(jī)實(shí)習(xí)題四一、題目：用newton法求方程x7-28x4+14=0在(0.1,1.9)中的近似根(初始近似值取為區(qū)間端點(diǎn)，迭代6次或誤差小于0.00001 )o二、解題方法及理論依據(jù):newton迭代法是平方收斂于方程f(x)二0在區(qū)間a,b上的唯一解。，收斂速度較快，循環(huán)次數(shù)少。 方法應(yīng)用條件：i ) f (a) f (b) <0ii) f" (x)在區(qū)

33、間a, b上不變號(hào).iii) f (x)h0iv) |f (c) |/b-aw |f (c)|其中c是a, b中使min | f (a), f (b)達(dá)到的一個(gè)，則對(duì)任意時(shí)近似值x0 a, b, newton迭代過程為xk+i=* (xk)=xf(xk)/f (xk),k= 1,2,3算法：令= x7 - 28x4 +14,/(0.1) > 0,/(l .9) < 0 fx) = 7x6-1 12x3 = lxx3-16)<0 fx) = 42x5 -336x2 = 42x2(x3 -8)<0 r(1.9)./(1.9)>0故以1.9為起點(diǎn)三、1 計(jì)算程序:#in

34、elude "math, h" main ()float x, y, f, fl;scanf ("%f, &x);doy=x；f=pow(y, 7)-28*pow(y, 4)+14; fl=7*pow(y, 6)-112*pow(y, 3); x=y-f/fl;while(fabs(x-y)>=le5);printf("nthe resuit of the question is/*定義f (x)的表達(dá)式*/ /*定義f' (x)的表達(dá)式*/* newton 迭代法*/*控制誤差小于0. 00001*/%fn", x);2

35、.打印結(jié)果:請(qǐng)輸入端點(diǎn)值：1.9x=0.8454970.13.030577四、問題討論:程序較為簡(jiǎn)單。它的兒何意義為xe是f(x)在點(diǎn)xk的切線與x軸交點(diǎn)，故也稱為切線法，它是平方收斂的, 此處取xf1.9收斂性較妬 要注意判斷fz (xk)是否為零。上機(jī)實(shí)習(xí)題五一、題目:用 romberg 算法求3xx14 (5x + 7) sin x1 dx (允許誤差 £ 二0.00001)。二、解題方法及理論依據(jù)：龍貝格(romberg)方法求數(shù)值積分 rti<0)=(b-a)/2*f(a)+f(b)j tig)= (1/2) * ti<1_1)+ (b-a)/21_1*efa+

36、 (2i-l)*(b-a)/21三、1.計(jì)算程序:#include"math. h/*求 f(x)=3xxl.4(5x+7)sinx2 的值*/int a=l, b=3; double f (double x)double z;z=pow(3, x)*pow(x, 1. 4)*(5*x+7)*sin(x*x);return (z);double s (int 1)/*求 t1 中的(b-a)/2'f a+(2il)*(b-a)/21*/extern a,b; int i, m; double z=0.0;m=pow (2, 1-1);for(i=l;i<=m;i+)z+

37、=f(a+1. 0*(2*i-l)/m);z*=l. 0*(b-a)/m;return (z);main ()( extern a, b;double t2020;i nt m, n, 1=0;tl0 = (b-a)/2. 0*(f (a) +f (b);do/*龍貝格(romberg)算法*/(1+；tll=0. 5*(t1 1-11+s)；n=l-l;for (m=2; n>=0; m+, n)/*求解 ti */tm n = (pow(4, mt)*tm-l n+1-tmt n)/(pow(4, mt)t. 0);while(fabs(t0-tl-l0)>=le-5);pri

38、ntf (，znt%d 0=%f，z, 1, tl 0);2打印結(jié)果：t80=440. 536017四、問題討論：此程序較繁，計(jì)算t嚴(yán)需要復(fù)化梯形公式，還要用到richardson外推法，構(gòu)造新序列，計(jì)算新分點(diǎn)的值時(shí), 這些數(shù)值個(gè)數(shù)成倍增加。應(yīng)用給出所要求的誤差5當(dāng)口#-,+嚴(yán)|£吋控制循環(huán)。程序具有廣泛的通用性。上機(jī)實(shí)習(xí)題六一、題目：用定步長(zhǎng)四階runge-kutta法求解c dyi/dt=ldyz/dt 二 y3dy3/dt=1000-1000y2_100y3'vi (0) =0y2(0)=0宀3(0)二0h=0. 0005,打印 yi(0.025), y：(0.045), y：(0.085), yi(0. 1), (i=l,2, 3)二、解題方法及理論依據(jù)：高階方程組的runge-kutta解法c 丫辺丸汁(1/6) * (k1+2k2+2k3+k jklh*f(x“, yn)j k2=h*f(xn+h/2, yn+ki/2)k產(chǎn)h*f(xn+h/2, ymk2/2).k1二h*f(xn+h, yn+ks)適用條件：使用于那些用普通的積分方法解不了的微分方程組只要知道函數(shù)z間的關(guān)系和初值就可以不用 解出表達(dá)式而直接求解函數(shù)在要求點(diǎn)的值。三* 1計(jì)算程序:include <stdio. h>