贵大数值分析上机实习报告

1、数值分析上机实习报告学院：土建学院 专业:学号: 姓名:上机实习题一、题目：已知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 (并打印b)。2. 用超松弛法求解bx=b (取松弛因子3二1.4,対°)二0,迭代9次)。3. 用列主元素消去法求解bx=bo二、解题方法的理论依据:1、用householder变换的理论依据(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 的转置*yr/a®(6 ) qr=yr-kr*ur(7 ) ar=a(r-l)-(qr*ur 的转置+ur*qr 的转置)r=l, 2,，n22、用超松弛法求解其基本思想：在高斯方法己求出x<m),的基础上，组合新的序列，从而加快收敛速度。其算式：bii -1 x z -1 + bii -xi +

6、 bii +1 x / +1 = bi却“匹心勿一1 一処也x0l +如 biibiibii<xi二 w文i + x0ix0i = xi其中3是超松弛因子，当3>1时，可以加快收敛速度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);严以下是消去法*/ 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.打印结果：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方法的基础上，已经求出x的第m解，第m-1基础上，经过重新组合 得新的序列，而在此新序列收敛速度加快。上机实习题二用对分法求b（即a）的小于23

17、且最接近23的特征值的近似值（误差小于0.1）在用反幕法求b的该特征值的更精 确近似值及相应的特征向量.二. 程序算法1. 对分法:对于对称的三对角方阵，其特征值必落在卜ii b ii , ii b ii冲，其中ii b ii是其任意范数.通过同号数的计 算,可知在0, ii b ii冲特征根的个数，在不断对分，就可得特征根所在的区间，如在ik,uk中有一根入k,且lk,uk很小, 就可取lk+uk/2作为ak的近似值.2. 反幕法:基于乘幕法屮求不太容易,因此不直接用a-1xk=xk+1,而使用axk+1=xk.用求解线性方程组求出 xk+1,再标准化，就可得到绝对值最小的特征值与相应的特征

18、向量,这就是反代法（也称反幕法）.三. 程序清单1. 对分法：#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;/*对分法所求得的特征值*/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);四. 运行结果1. 对分法求强特征值为：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对分法简单可靠，数值稳定性较高，对于求少量几个特征值特别适宜.但收敛速度较慢.2.反幕法是结合对分法使用的，所以近似值较恰当，其收敛速度是惊人的只耍迭代两次就可得到较 满意的结果但在运用中需把一般矩阵化为hessebberg阵才可计算.上机实习题三一、题目：己知函数值如下表: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)的近似值。二、解题方法的理论依据：任意划分的三弯矩插值法以及方程组解法中的三对角阵追赶算法。应用三次样条插值法能够对函数产生很好的逼近效果。而追赶算法又具有计算量少、方法简单、算法稳定的特 点。方法应用条件：适用于求复杂函数在给定区间内某一点的函数值，给出函数f(x)在区间a,b中的n个插值点，并 且给出函数在区间端点处的值。三* 1计算程序:#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(”请输入要插的值:”);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("请输入要求的导数的值:”)；scanf(”d”,&m);printf(nr(%d)=%lfn",m,(cm+l-cm-l)/2.0);输出结果：请输入要插的值：4.563f(4.563)= 1.517932请输入要求的导数的值：4.563f，(4.563)= 0.249350四、问题讨论：在给均匀分划的插值函数x赋值时，由于使用for循环，

32、误将xi二i + 1写成xi二i,导致运算错误。此程 序具有一定通用性，对于任意划分的三弯矩插值法，只许改动可。求解方程组mj时，要用到三对角方程 组的追赶法(也称thomas算法)。变量较多，应注意区分。求导时注意正负号。上机实习题四一、题目：用newton法求方程x7-28x4+14=0在(0.1,1.9)中的近似根(初始近似值取为区间端点，迭代6次或误差小于0.00001 )o二、解题方法及理论依据:newton迭代法是平方收敛于方程f(x)二0在区间a,b上的唯一解。，收敛速度较快，循环次数少。 方法应用条件：i ) f (a) f (b) <0ii) f" (x)在区

33、间a, b上不变号.iii) f (x)h0iv) |f (c) |/b-aw |f (c)|其中c是a, b中使min | f (a), f (b)达到的一个，则对任意时近似值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为起点三、1 计算程序:#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)的表达式*/ /*定义f' (x)的表达式*/* newton 迭代法*/*控制误差小于0. 00001*/%fn", x);2

35、.打印结果:请输入端点值：1.9x=0.8454970.13.030577四、问题讨论:程序较为简单。它的儿何意义为xe是f(x)在点xk的切线与x轴交点，故也称为切线法，它是平方收敛的, 此处取xf1.9收敛性较妬 要注意判断fz (xk)是否为零。上机实习题五一、题目:用 romberg 算法求3xx14 (5x + 7) sin x1 dx (允许误差 £ 二0.00001)。二、解题方法及理论依据：龙贝格(romberg)方法求数值积分 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.计算程序:#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打印结果：t80=440. 536017四、问题讨论：此程序较繁，计算t严需要复化梯形公式，还要用到richardson外推法，构造新序列，计算新分点的值时, 这些数值个数成倍增加。应用给出所要求的误差5当口#-,+严|£吋控制循环。程序具有广泛的通用性。上机实习题六一、题目：用定步长四阶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)二、解题方法及理论依据：高阶方程组的runge-kutta解法c 丫辺丸汁(1/6) * (k1+2k2+2k3+k jklh*f(x“, yn)j k2=h*f(xn+h/2, yn+ki/2)k产h*f(xn+h/2, ymk2/2).k1二h*f(xn+h, yn+ks)适用条件：使用于那些用普通的积分方法解不了的微分方程组只要知道函数z间的关系和初值就可以不用 解出表达式而直接求解函数在要求点的值。三* 1计算程序:include <stdio. h>