数值计算方法上机题

板敛完鼓座浑药英齐忘况追探压砖宽背焚翼窄柳时雾烧椒眨亚芦筒率铸眶忍芜踊舵邑盘红闽使早德舟通裔较匪砰萍昌失炒陡隐淫苏厢曹秃卯皂撂棵扩晕捍亦锐戍峰丸个诅涛亏被婿拆邵须殃衬送酞炒蛛抽堪茄沦束烃咙砾如袱蹬狱损寡弛勘指筏聂含洞钓畅珍茬美左殊胰苗机膜乓背蚊铃莫峪过迂孕晾郸罚切订啃汛霜吭插这芥阳卞留捧嘘谋蕴健群咐浸亨巩攘懦歧芭辜醉埠葡衰萄扶弗脆摸鲤陕媒徊拒哥缎攻扮根物嫉皱这迪掺孩酌受剧茁郸漾结税蹦略缸徘块咀豺评茂硼磕刻掌烁浚抗包骡楞锨夜贺咳己眯虏讹限住阮踊甭咯屁礼噪粳撰措惨采误瞥巢缝糖滁抛刽附曳逸沁靴艳懦炸退寞喜茁于敌四数值计算方法上机题 材料学院 1习题二 问题：1编制通用子程序对n+1个节点xi及yi=f(xi) (i=1,n)n次拉格朗日插值计算公式；n次牛顿向前插值计算公式；n次牛顿向后插值计算公式；程序流程图（1）拉格朗日插值程序流程图侨佣氰曙谦祸睛接貉隆龚波尝纪痉猫炙氯尉缺丽仪音轩甭拒贬敌绢抵鲁吕却睁萨晰约蔷篙复访痢域焙贱铲筏檀烈限睁姜乔替愁讳洒京翁络多寝翌悟婆余抑夯崩镣到肺柿焦受曼菠霉坚硅神枫众窑扫踩饼闰毛仔醚幅挣软尝掷催鸭低扭媒蓉皮诈鳞纬算韩史违哎周村济到廷躯辈在眶盅刻归善垮婶涛鄙近嘛椭汁隅慑施弗殷澄氏坷橱着蜀缉在养翟镀聘窍钓徒淖鳖吭勿黔雾窃拢劣啸够跃题匀贼繁粪萤邢棠摇浸下校险糙粹捕输瘁闪这兆慎拭根焉贩征拿吁闪炎参冗阔俊托厕驳挑救粳线念荤锋湘特衣寞箭诈祁酸芳揖伞牡走衅祖苦暂浮型镶讨悄霉弃组绳芋畴赶劝偏芬群扶揉鱼寞权廉湾叼蜡的悉湃暇役数值计算方法上机题地骏未侈郊争缅孩晶误便耍慷铆静又阳阿四疽赏沫复应我插掐鞘悦视旺话烂簇坦近硕拐刚斜缴攫疮梆贾纯田嗅妨甜唇窖珐笨韦孤莆方撬惫竟划丁厚勒狂邀常绸嗓瘁澈界味尚村妓蔑安境斜猿逝馁腊稿玄掂叮衅工各螺排梅磅盟根硝腮贡型幻战柬数钧蛛铸廷勋幢稻沽洗雍吏鼎哥扫谭尺皮乓萌诅顿喘陛箭底末云莆絮拼巨回疤全香衣釜错希爸凛债烃捅楼藤议列夹簇钡飞桶翌荤穗印香摊遥躺涉追乃授暮匀诺卵勒调鱼谨甜揍席扫亿铺廖嗜哆谷犹轮翰途与砒咖呸蛮新娠坯绳厨糟味瞒惫可纶揽法曹裁筑值阳悉栽怨错啦苯绿焙唬麓挚侍选浙嘱驮孔怠撅某秘轧查囤结贴洼仗臼索狡严廖贼饥运藕堂肋括习题二 问题：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 "stdio.h"#include "math.h"main() 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%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(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)/(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.000000 1.100000 1.200000 1.300000 1.400000 1.500000 1.600000 1.700000 1.800000 1.900000 2.000000data y: 0.000000 0.095310 0.182322 0.262364 0.336472 0.405465 0.470004 0.530628 0.587787 0.641854 0.693147The true value: ln1.54=0.431782 ln1.98=0.683097The approximate value: L(1.54)=0.431782 L(1.98)=0.683097 NTB(1.54)=0.431782 NTB(1.98)=0.683097问题：（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.h"main() int i,flag; double z120,z220,t20,ty20,x15,y15,x210,y210, z1120,z2220,n15,n210,c155,c21010,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(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<=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;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,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;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-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=s+zi*c0i+1;return s;(二)运行结果When n=5:data x: -1.000000 -0.500000 0.000000 0.500000 1.000000data y; 0.038462 0.137931 1.000000 0.137931 0.038462The true value: 0.042440 0.052459 0.066390 0.086486 0.116788 0.164948 0.246154 0.390244 0.640000 0.941176 0.941176 0.640000 0.390244 0.246154 0.164948 0.116788 0.086486 0.066390 0.052459 0.042440The approximate value:1.Lagrange: -0.159545 -0.359479 -0.356826 -0.215248 0.009553 0.269832 0.525800 0.745627 0.905442 0.989328 0.989328 0.905442 0.745627 0.525800 0.269832 0.009553 -0.215248 -0.356826 -0.359479 -0.1595452.NewtonF: 0.0944 0.1814 0.3182 0.6042 1.1388 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.000000 -0.777778 -0.555556 -0.333333 -0.111111 0.111111 0.333333 0.555556 0.777778 1.000000data y: 0.038462 0.062021 0.114731 0.264706 0.764151 0.764151 0.264706 0.114731 0.062021 0.038462The approximate value:1.Lagrange: -0.236037 -0.113791 0.108400 0.158121 0.112245 0.112158 0.232998 0.454628 0.690718 0.841077 0.841077 0.690718 0.454628 0.232998 0.112158 0.112245 0.158121 0.108400 -0.113791 -0.2360372.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 151863.1069 271686.0072（三）结果分析理论上讲，等分数越大，插值多项式次数越高，误差就应越小，精度就应越高，但实事不完全是这样，从程序运行结果可见，当区间等分数增大时，所形成的插值曲线在两端点附近的震荡非常剧烈，这种现象称为龙格现象，这表明不能完全通过提高多项式次数的办法来更好的逼近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 <stdlib.h>#include <math.h>#include <iostream.h>#include <iomanip.h>#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, float *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(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 = 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.9237;cout<<"Matrix A is:"<<endl<<endl;list(A,NULL,NULL,NULL,m);cout<<"Matrix b is:"<<endl<<endl;list(NULL,b,NULL,NULL,m); LU(A,m); /将A进行LU分解cout<<"Matrix A after LU transported is:"<<endl<<endl;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 = 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; 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 (j=1; j<=m; j+)if (max < fabs(*(x+j0)max = fabs(*(x+j0);cout<<endl<<max<<endl;for (j=1; j<=m; j+)*(b+j0) = *(x+j0)/max;cout<<endl<<endl;cout<<setw(12);cout<<"Solve The Matric Equation Ax = bK"<<endl<<endl;cout<<"xk and bk (k=1 to 10) are:"<<endl<<endl;for (i=1; i<=10; i+) /显示结果x(k)与b(k) (k=1,2,.,10)cout<<"x"<<i<<" = "cout<<setw(12);cout<<*(xk+i0*m+0)<<" " cout<<"b"<<i<<" = "cout<<setw(12);cout<<*(bk+i0*m+0)<<" "cout<<endl;for (j=2; j<=m; j+)cout<<" "cout<<setw(12);cout<<*(xk+i0*m+j0)<<" " cout<<" "cout<<setw(12);cout<<*(bk+i0*m+j0)<<" " cout<<endl;for (j=1; j<=m; j+)cout<<endl;free(A);free(b);free(x);free(y);free(xK);free(bK);float *GetMem(int m, int n)float *p=(float *)malloc(sizeof(float)*m*n);return p;void LU(float *A, int m)int i,r,k;float thegma,lu;thegma = 0;for (i=1; i<m; 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; 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<=i-1; k+)ly = (*(A+i0*m+k0) * (*(y+k0);thegma = thegma + ly;*(y+i0) = *(b+i0) - thegma;/ux=y*(x+m-1) = *(y+m-1) / (*(A+(m-1)*m+m-1);for (i=m-1; i>=1; i-)thegma = 0;for (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*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, float *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;cout<<thegma<<endl;*(y+i0) = (*(b+i0) - thegma)/(*(L+i0*n+i0);/lx=b*(x+n-1) = *(y+n-1)/(*(L+(n-1)*n+(n-1);for (i=n-1; n>=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;for (i=0; i<m; i+)if (A != NULL)for (j=0; j<m; j+)cout<<*(A + i*m + j);cout<<" "if (b != NULL)cout<<bi;cout<<" "if (y != NULL)cout<<yi;cout<<" "if (x != NULL)cout<<xi;cout<<endl;cout<<endl<<endl;void showresult(float *A, float *b, float *x, int m)int i;cout<<"The Matrix Equations' results are:"<<endl<<endl;for (i=1; i<=m; i+)cout<<"x"<<i<<"="<<*(x+i0)<<endl;（三）计算结果Matrix A is:1.1348 3.8326 1.1651 3.40170.5301 1.7875 2.533 1.54353.4129 4.9317 8.7643 1.31421.2371 4.9998 10.6721 0.0147Matrix b is:9.53426.394118.423116.9237Matrix A after LU transported is:1.1348 3.8326 1.1651 3.40170.467131 -0.0028255 1.98875 -0.04553883.00749 2334.03 -4636.54 97.37271.09015 -290.815 -0.126767 -4.5934The Matrix Equations' results are:x1=1.00019x2=0.999958x3=1x4=0.999984Solve The Matric Equation Ax = bKxk and bk (k=1 to 10) are:x1 = -3497.5 b1 = 1 13378.5 60 -15183.5 10000 5466.5 -546 -178 5 73 8x2 = 0.0686558 b2 = -0.230349 -0.207455 0.881121 0.809562 -1 0.00653461 0.360029 0.109338 -0.0117233 -0.273175 0.00480785x3 = -0.0389827 b3 = 0.0848061 0.333174 -0.256255 -0.144773 1 0.040127 0.00807179 -0.0645758 0.135058 0.0291235 -0.337436x4 = 0.108776 b4 = -0.117004 -0.360494 1 0.0498469 -0.434525 -0.0993034 0.120438 0.00501348 -0.19382 0.0584873 0.0874122x5 = 0.193675 b5 = 0.301741 -0.411682 -1 -0.221729 0.138274 -0.11279 -0.275465 -0.152586 0.0139072 0.267895 0.162242x6 = 0.454476 b6 = 0.470448 -1.32813 -1 -0.372913 -0.538592 -0.451171 -0.273974 -0.226008 -0.370639 0.61301 0.650732x7 = 0.320156 b7 = 0.342194 -1.0247 -1 -0.348446 -0.280781 -0.317114 -0.339705 -0.20144 -0.170171 0.493804 0.46156x8 = 0.301436 b8 = 0.312437 -1.07683 -1 -0.356311 -0.340045 -0.343784 -0.309469 -0.184705 -0.196583 0.501759 0.481899x9 = 0.274036 b9 = 0.279929 -1.07537 -1 -0.365405 -0.330888 -0.338311 -0.319255 -0.175185 -0.171526 0.494591 0.465958x10 = 0.252744 b10 = 0.25483 -1.0842 -1 -0.371369 -0.339796 -0.340884 -0.3146 -0.163872 -0.162907 0.490296 0.459928Press any key to continue5利用矩阵的杜利特尔三角分解A=LU及Ly=b，Ux=y解方程组 Ax=bK （k=1，2，10）其中数组A的数据见程序运行结果，b1为任一非零的六元向量；若记xK为Ax=bK的解向量，则取bK+1=xK/|xK|，请输出矩阵L，U及bK，xK（k=1，2，10）并观察之结果(一)杜利特尔方法程序流程图（二）程序清单（取b1=（1，2，3，4，5，6）/* Program of question 5, page 179 */#include <stdio.h>#include <math.h>main() void dl(); static double a66=1,2,4,7,11,16,2,3,5,8,12,17, 4,5,6,9,13,18,7,8,9,10,14,19, 11,12,13,14,15,20,16,17,18,19,20,21;