数值计算方法上机题(共41页).doc

《数值计算方法上机题(共41页).doc》由会员分享，可在线阅读，更多相关《数值计算方法上机题(共41页).doc（41页珍藏版）》请在淘文阁 - 分享文档赚钱的网站上搜索。

1、精选优质文档-倾情为你奉上板敛完鼓座浑药英齐忘况追探压砖宽背焚翼窄柳时雾烧椒眨亚芦筒率铸眶忍芜踊舵邑盘红闽使早德舟通裔较匪砰萍昌失炒陡隐淫苏厢曹秃卯皂撂棵扩晕捍亦锐戍峰丸个诅涛亏被婿拆邵须殃衬送酞炒蛛抽堪茄沦束烃咙砾如袱蹬狱损寡弛勘指筏聂含洞钓畅珍茬美左殊胰苗机膜乓背蚊铃莫峪过迂孕晾郸罚切订啃汛霜吭插这芥阳卞留捧嘘谋蕴健群咐浸亨巩攘懦歧芭辜醉埠葡衰萄扶弗脆摸鲤陕媒徊拒哥缎攻扮根物嫉皱这迪掺孩酌受剧茁郸漾结税蹦略缸徘块咀豺评茂硼磕刻掌烁浚抗包骡楞锨夜贺咳己眯虏讹限住阮踊甭咯屁礼噪粳撰措惨采误瞥巢缝糖滁抛刽附曳逸沁靴艳懦炸退寞喜茁于敌四数值计算方法上机题 材料学院 1习题二 问题：1编制通用子程序

2、对n+1个节点xi及yi=f(xi) (i=1,n)n次拉格朗日插值计算公式；n次牛顿向前插值计算公式；n次牛顿向后插值计算公式；程序流程图（1）拉格朗日插值程序流程图侨佣氰曙谦祸睛接貉隆龚波尝纪痉猫炙氯尉缺丽仪音轩甭拒贬敌绢抵鲁吕却睁萨晰约蔷篙复访痢域焙贱铲筏檀烈限睁姜乔替愁讳洒京翁络多寝翌悟婆余抑夯崩镣到肺柿焦受曼菠霉坚硅神枫众窑扫踩饼闰毛仔醚幅挣软尝掷催鸭低扭媒蓉皮诈鳞纬算韩史违哎周村济到廷躯辈在眶盅刻归善垮婶涛鄙近嘛椭汁隅慑施弗殷澄氏坷橱着蜀缉在养翟镀聘窍钓徒淖鳖吭勿黔雾窃拢劣啸够跃题匀贼繁粪萤邢棠摇浸下校险糙粹捕输瘁闪这兆慎拭根焉贩征拿吁闪炎参冗阔俊托厕驳挑救粳线念荤锋湘特衣寞箭诈祁

3、酸芳揖伞牡走衅祖苦暂浮型镶讨悄霉弃组绳芋畴赶劝偏芬群扶揉鱼寞权廉湾叼蜡的悉湃暇役数值计算方法上机题地骏未侈郊争缅孩晶误便耍慷铆静又阳阿四疽赏沫复应我插掐鞘悦视旺话烂簇坦近硕拐刚斜缴攫疮梆贾纯田嗅妨甜唇窖珐笨韦孤莆方撬惫竟划丁厚勒狂邀常绸嗓瘁澈界味尚村妓蔑安境斜猿逝馁腊稿玄掂叮衅工各螺排梅磅盟根硝腮贡型幻战柬数钧蛛铸廷勋幢稻沽洗雍吏鼎哥扫谭尺皮乓萌诅顿喘陛箭底末云莆絮拼巨回疤全香衣釜错希爸凛债烃捅楼藤议列夹簇钡飞桶翌荤穗印香摊遥躺涉追乃授暮匀诺卵勒调鱼谨甜揍席扫亿铺廖嗜哆谷犹轮翰途与砒咖呸蛮新娠坯绳厨糟味瞒惫可纶揽法曹裁筑值阳悉栽怨错啦苯绿焙唬麓挚侍选浙嘱驮孔怠撅某秘轧查囤结贴洼仗臼索狡严廖贼饥

4、运藕堂肋括习题二 问题：1编制通用子程序对n+1个节点xi及yi=f(xi) (i=1,n)（1） n次拉格朗日插值计算公式；（2） n次牛顿向前插值计算公式；（3） n次牛顿向后插值计算公式；（一） 程序流程图（1）拉格朗日插值程序流程图（2）牛顿向前插值程序流程图（3）牛顿向后插值程序流程图（二）源程序 见主程序清单问题：2计算（1）已知f(x)=lnx,，a，b=1，2，取h=0.1，xi=1+ih，i=0，1，10。用通用程序（1），（3）计算ln1.54及ln1.98的近似值；（一）程序清单/* program of question 2.1, page 61 */#include

5、stdio.h#include math.hmain() int i,flag=0; double z1,z2,x11,y11,t,s1,s2,z10,c1111,log(double),ntb(),L(); for(i=0;i=10;i+)xi=1+0.1*i;yi=log(xi); printf(data x:n); for(i=0;i=10;i+) flag+;printf(%11.6f,xi);if(flag%4=0)printf(n); printf(ndata y:n);flag=0; for(i=0;i=10;i+) flag+;printf(%11.6f,yi);if(flag

6、%4=0)printf(n); printf(nThe true value:n); printf( ln1.54=%f ln1.98=%fn,log(1.54),log(1.98); z1=L(x,y,10,1.54);z2=L(x,y,10,1.98);t=(1.54-x10)/0.1; s1=ntb(y,10,t,z,c);s2=ntb(y,10,t,z,c);t=(1.98-x10)/0.1; printf(The approximate value:n); printf( L(1.54)=%f L(1.98)=%fn,z1,z2); printf( NTB(1.54)=%f NTB(

7、1.98)=%fn,s1,s2); double L(double x,double y,int n,double t) int i,k; double z=0.0,s;if(n=1)z=y0; for(k=0;k=n;k+) s=1.0;for(i=0;i=n;i+) if(i!=k)s=s*(t-xi)/(xk-xi); z=z+s*yk; return z; double ntb(double y,int n,double t,double z,double c11) int i,j,sn=n;double s;z0=t; for(i=1;i=n-1;i+) zi=zi-1*(t+i)/

8、(i+1); for(i=0;i=n;i+) ci0=ysn-; for(j=1;j=n;j+) for(i=0;i=n-j;i+) cij=cij-1-ci+1j-1; s=yn; for(i=0;i=n-1;i+) s=s+zi*c0i+1; return s; (二)运行结果data x: 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 2.data y: 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.The true value: ln1.54=0. ln1.98=0.The approximate value: L(1.54)=0. L(1.98)=0.

9、 NTB(1.54)=0. NTB(1.98)=0.问题：（2）f(x)=1/(1+25x2), |x|1取等距节点n=5和n=10，用通用程序（1），（2）依次计算x=-0.95+ih（i=0，1，19，h=0.1）处f(x)的近似值，并将其结果与其真实值相比较。（一）程序清单/* program of question 2.2,page 61 */#include stdio.h#include math.hmain() int i,flag; double z120,z220,t20,ty20,x15,y15,x210,y210, z1120,z2220,n15,n210,c155,c2

10、1010,m; double L(),ntf1(),ntf2(); for(i=0;i=4;i+) x1i=-1+0.5*i;y1i=1.0/(1+25*x1i*x1i); printf(When n=5:ndata x:n); for(i=0;i=4;i+) printf(%10.6f,x1i); printf(ndata y;n); for(i=0;i=4;i+) printf(%10.6f,y1i); for(i=0;i=19;i+) ti=-0.95+i*0.1;z1i=L(x1,y1,4,ti); m=(ti-x10)/0.1;z11i=ntf1(y1,4,m,n1,c1); for

11、(i=0;i=19;i+) tyi=1.0/(1+25*ti*ti); printf(nThe true value:n);flag=0; for(i=0;i=19;i+) flag+;printf(%10.6f,tyi);if(flag%5=0) printf(n); printf(The approximate value:n);printf(1.Lagrange:n); flag=0; for(i=0;i=19;i+) flag+;printf(%10.6f,z1i);if(flag%5=0)printf(n); printf(2.NewtonF:n);flag=0; for(i=0;i

12、=19;i+) flag+;printf(%12.4f,z11i);if(flag%4=0)printf(n); for(i=0;i=9;i+) x2i=-1+(2.0/9)*i;y2i=1.0/(1+25*x2i*x2i); printf(nWhen n=10:ndata x:n); for(i=0;i=9;i+) printf(%10.6f,x2i);if(i=4)printf(n); printf(ndata y:n); for(i=0;i=9;i+) printf(%10.6f,y2i);if(i=4)printf(n); for(i=0;i=19;i+) ti=-0.95+i*0.1

13、;z2i=L(x2,y2,9,ti); m=(ti-x20)/0.1;z22i=ntf2(y2,9,m,n2,c2); printf(nThe approximate value:n);printf(1.Lagrange:n); flag=0; for(i=0;i=19;i+) flag+;printf(%10.6f,z2i);if(flag%5=0)printf(n); printf(2.NewtonF:n);flag=0; for(i=0;i=19;i+) flag+;printf(%12.4f,z22i);if(flag%4=0)printf(n); double L(double x,

14、double y,int n,double t) int i,k;double z=0.0,s;if(n=1) z=y0; for(k=0;k=n;k+)s=1.0;for(i=0;i=n;i+) if(i!=k) s=s*(t-xi)/(xk-xi);z=z+s*yk; return z; double ntf1(double y,int n,double t,double z,double c5) int i,j;double s;z0=t;for(i=1;in-1;i+) zi=zi-1*(t-i)/(i+1);for(i=0;i=n;i+) ci0=yi;for(j=1;j=n;j+)

15、for(i=0;i=n-j;i+) cij=ci+1j-1-cj-1j-1;s=y0;for(i=0;i=n-1;i+) s=s+zi*c0i+1;return s; double ntf2(double y,int n,double t,double z,double c10) int i,j;double s;z0=t;for(i=1;i=n-1;i+) zi=zi-1*(t-i)/(i+1);for(i=0;i=n;i+) ci0=yi;for(j=1;j=n;j+)for(i=0;i=n-j;i+) cij=ci+1j-1-cij-1;s=y0;for(i=0;i=n-1;i+) s=

16、s+zi*c0i+1;return s;(二)运行结果When n=5:data x: -1. -0. 0. 0. 1.data y; 0. 0. 1. 0. 0.The true value: 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.The approximate value:1.Lagrange: -0. -0. -0. -0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. -0. -0. -0. -0.2.NewtonF: 0.0944 0.1814 0.3182 0.6042 1.1

17、388 2.0216 3.3520 5.2295 7.7536 11.0236 15.1392 20.1997 26.3046 33.5535 42.0457 51.8807 63.1581 75.9772 90.4376 106.6387When n=10:data x: -1. -0. -0. -0. -0. 0. 0. 0. 0. 1.data y: 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.The approximate value:1.Lagrange: -0. -0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. -

18、0. -0.2.NewtonF: -0.2093 0.1604 0.1155 0.5232 0.8615 0.5232 0.1155 0.1604 -0.2093 3.2172 46.9321 265.7964 1021.3447 3120.8273 8160.6933 19033.4039 40641.6516 80869.2455 .1069 .0072（三）结果分析理论上讲，等分数越大，插值多项式次数越高，误差就应越小，精度就应越高，但实事不完全是这样，从程序运行结果可见，当区间等分数增大时，所形成的插值曲线在两端点附近的震荡非常剧烈，这种现象称为龙格现象，这表明不能完全通过提高多项式次

19、数的办法来更好的逼近f(x)，另外，插值多项式的次数越高，计算工作量就越大，累计误差也越大，因此，在实际计算中，应采用分段低次插值的方法。习题五 问题：1.（3）实现矩阵三角分解的杜利特尔及乔利斯方法及用此方法解Ax=b的过程。（一）流程图：输入n,A对于j=1,nu1j=a1jj=1k=j+1对于i=k,n计算lij和ukjj=i-1?打印L，Uj=i+1=（源程序见第5题主程序清单）问题：2.用列主消去法程序第1题中（1）解方程组如程序清单中所示，并比较计算结果精度（x1=x2=x3=x4=1）（一）流程图（二）源程序#include #include #include #include

20、#define i0 (i-1)#define j0 (j-1)#define r0 (r-1)#define k0 (k-1)/* GetMem获取矩阵内存空间 * 获取m乘n维的float类型内存空间 */float *GetMem(int m, int n);/* LU将矩阵A进行LU分解 * *A为数组A的地址，m为矩阵阶数 */void LU(float *A, int m);/* Solve根据LU分解后的A阵解方程Ax = b， * *A为数组A的地址，*b为数组b的地址， * *x为数组x的地址,*y为数组y的地址,m为矩阵阶数 */void Solve(float *A, f

21、loat *b, float *x, float *y, int m);/* LL将矩阵A进行LL分解 * *A为数组A的地址，m为矩阵阶数 */void LL(float *A, float *L, int n);/* SolveLL根据LL分解后的A阵解方程Ax = b， * *A为数组A的地址，*b为数组b的地址， * *x为数组x的地址,*y为数组y的地址,m为矩阵阶数 */void SolveLL(float *L, float *b, float *x, float *y, int n);/* list输出矩阵 * A为A阵地址，b为b阵地址，xK为x(k)向量，xkadd1为x(

22、k+1)向量， * * n为方阵阶数 * */void list(float *A, float *b, float *x, float *y, int m);/* showresult输出方程结果 * xK为x(k)向量，n为方阵阶数 */void showresult(float *A, float *b, float *x, int m);main()int i,j;int m = 4;int n = 4;float *A = NULL;float *b = NULL;float *y = NULL;float *x = NULL;float *xK = NULL;float *bK =

23、 NULL;A = GetMem(m,n);b = GetMem(m,1);y = GetMem(m,1);x = GetMem(m,1);A0=1.1348; A1=3.8326; A2=1.1651; A3=3.4017; /输入原始数组数据A4=0.5301; A5=1.7875; A6=2.5330; A7=1.5435;A8=3.4129; A9=4.9317; A10=8.7643; A11=1.3142;A12=1.2371; A13=4.9998; A14=10.6721; A15=0.0147;b0=9.5342; b1=6.3941; b2=18.4231; b3=16.9

24、237;coutMatrix A is:endlendl;list(A,NULL,NULL,NULL,m);coutMatrix b is:endlendl;list(NULL,b,NULL,NULL,m); LU(A,m); /将A进行LU分解coutMatrix A after LU transported is:endlendl;list(A,NULL,NULL,NULL,m); /显示A阵Solve(A,b,x,y,m); /根据LU分解后的A阵解方程Ax = bshowresult(A,b,x,m); /显示方程的解xfree(A);free(b);free(x);free(y);m

25、 = 6;n = 6;A = GetMem(m,n);b = GetMem(m,1);y = GetMem(m,1);x = GetMem(m,1);xK = GetMem(m,10);bK = GetMem(m,10);A0=1; A1=2; A2=4; A3=7; A4=11; A5=16; /输入原始数组数据A6=2; A7=3; A8=5; A9=8; A10=12; A11=17;A12=4; A13=5; A14=6; A15=9; A16=13; A17=18;A18=7; A19=8; A20=9; A21=10; A22=14; A23=19;A24=11; A25=12;

26、A26=13; A27=14; A28=15; A29=20;A30=16; A31=17; A32=18; A33=19; A34=20; A35=21;b0=1; b1=60; b2=10000; b3=-546; b4=5; b5=8;for (i=1; i=10; i+) /求x(k)与b(k) (k=1,2,.,10)for (j=1; j=m; j+)*(bk+i0*m+j0) = *(b+j0);LU(A,m);Solve(A,b,x,y,m);for (j=1; j=m; j+)*(xk+i0*m+j0) = *(x+j0);float max = 0; /|xk|inffor

27、 (j=1; j=m; j+)if (max fabs(*(x+j0)max = fabs(*(x+j0);coutendlmaxendl;for (j=1; j=m; j+)*(b+j0) = *(x+j0)/max;coutendlendl;coutsetw(12);coutSolve The Matric Equation Ax = bKendlendl;coutxk and bk (k=1 to 10) are:endlendl;for (i=1; i=10; i+) /显示结果x(k)与b(k) (k=1,2,.,10)coutxi = ;coutsetw(12);cout*(xk+

28、i0*m+0) ; coutbi = ;coutsetw(12);cout*(bk+i0*m+0) ;coutendl;for (j=2; j=m; j+)cout ;coutsetw(12);cout*(xk+i0*m+j0) ; cout ;coutsetw(12);cout*(bk+i0*m+j0) ; coutendl;for (j=1; j=m; j+)coutendl;free(A);free(b);free(x);free(y);free(xK);free(bK);float *GetMem(int m, int n)float *p=(float *)malloc(sizeof

29、(float)*m*n);return p;void LU(float *A, int m)int i,r,k;float thegma,lu;thegma = 0;for (i=1; im; i+)*(A+i*m) = *(A+i*m)/(*A);for (r=2; r=m; r+)for (i=r; i=m; i+)thegma = 0;for (k=1; k=r-1; k+)lu = (*(A+r0*m+k0) * (*(A+k0*m+i0);thegma = thegma + lu;*(A+r0*m+i0) = *(A+r0*m+i0) - thegma;for (i=r+1; i=m

30、; i+)thegma=0;for (k=1; k=r-1; k+)lu = (*(A+i0*m+k0) * (*(A+k0*m+r0);thegma = thegma + lu;*(A+i0*m+r0) = (*(A+i0*m+r0) - thegma)/(*(A+r0*m+r0);void Solve(float *A, float *b, float *x, float *y, int m)int i,k;float thegma,ly,ux;/ly=b*y = *b;for (i=2; i=m; i+)thegma = 0;for (k=1; k=1; i-)thegma = 0;fo

31、r (k=i+1; k=m; k+)ux = (*(A+i0*m+k0) * (*(x+k0);thegma = thegma + ux;*(x+i0) = (*(y+i0) - thegma)/(*(A+i0*m+i0);void LL(float *A, float *L, int n)int i,j,k;float thegma,ljk,likljk;*L=sqrt(*A);for (j=2; j=n; j+)*(L+j0*n) = *(A+j0*n)/(*L);for (j=2; j=n; j+)thegma = 0;for (k=1; k=j-1; k+)ljk = (*(L+j0*

32、n+k0) * (*(L+j0*n+k0);thegma = thegma + ljk;*(L+j0*n+j0) = sqrt(*(A+j0*n+j0) - thegma);for (i=j+1; i=n; i+)thegma = 0;for (k=1; k=j-1; k+)ljk = (*(L+i0*n+k0) * (*(L+j0*n+k0);thegma = thegma + ljk;*(L+i0*n+j0) = (*(A+i0*n+j0) - thegma)/(*(L+j0*n+j0);list(L,NULL,NULL,NULL,n);void SolveLL(float *L, flo

33、at *b, float *x, float *y, int n)int i,k;float thegma,ly,lx;/ly=b*y = *b/(*L);for (i=2; i=n; i+)thegma = 0;for (k=1; k=i-1; k+)ly = (*(L+i0*n+k0) * (*(y+k0);thegma = thegma + ly;coutthegma=1; n-)thegma = 0;for (k=i+1; k=n; k+)lx = (*(L+k0*n+i) * (*(x+k0);thegma = thegma + lx;*(x+i0) = (*(y+i0) - thegma)/(*(L+i0*n+i0);void list(float *A, float *b, float *x, float *y, int m)int i,j;f