《数值分析》上机习题

1、数值分析上机实验报告数值分析上机实验报告学 院：机械工程学院专 业：机械工程所在班级：机研1402学 号：2014021847姓 名：王浦舟指导老师：代新敏老师2015年1月5日题目1.用Newton法求方程 X7 - 28X4 + 14=0在（0.1，1.9）中的近似根（初始近似值取为区间端点，迭代6次或误差小于0.00001）。（1）算法说明：设函数在有限区间上二阶导数存在，且满足条件 1、；2、在区间上不变号；3、；4、其中是中使达到的一个。则对任意初始近似值，由Newton迭代过程 所生的迭代序列平方收敛于方程在区间上的唯一解。令 故以1.9为起点 （2）程序清单#include &q

2、uot;stdio.h"#include "math.h"main()double x1=1.9,eps=0.00001,x0; /* 输入初值x1 */do x0=x1; x1=x0-(pow(x0,7)-28*pow(x0,4)+14)/(7*pow(x0,6)-112*pow(x0,3); /* 利用算法公式迭代 */ while(fabs(x1-x0)>=eps); /* 结束条件表述 */ printf("x1=%5fn",x1);（3）运行结果题目2.已知函数值如下表：x12345f(x)00.693147181.098612

3、31.38629441.6094378x678910f(x)1.79175951.94591012.0794452.19722462.3025851f(x)f(1)=1f(10)=0.1试用三次样条插值求f(4.563)及f(4.563)的近似值。（1）算法说明a. 三次样条插值的三弯矩法 (j=1,2,N-1)原方程可转化为： LUx=d 即。记，则：=由已知得出的，再由方程：=最终的得到解方程Ax=d的算法：（2）程序清单#include <math.h>#include <stdio.h>main( )double x10=1.0,2.0,3.0,4.0,5.0,

4、6.0,7.0,8.0,9.0,10.0,y10=0.0,0.69314718,1.0986123,1.3862944,1.6094378,1.7917595,1.9459101,2.079445,2.1972246,2.3025851,xo=4.563,h9,k8,l8,d10,yl0=1.0,yl9=0.1,a9,b10=2.0,2.0,2.0,2.0,2.0,2.0,2.0,2.0,2.0,2.0,c9,p10,q10,w9,m10,f1,f2;int i;for(i=0;i<9;i+)hi=xi+1-xi;for(i=0;i<8;i+)ki=hi/(hi+hi+1);li=

5、1-ki;ai=ki;ci+1=li;a8=c0=1;d0=(6/h0)*(y1-y0)/h0-yl0);d9=(6/h8)*(yl9-(y9-y8)/h8); for(i=1;i<9;i+)di=(6/(hi-1+hi)*(yi+1-yi)/hi-(yi-yi-1)/hi-1); p0=b0;q0=d0;for(i=0;i<9;i+)wi=ai/pi;pi+1=bi+1-wi*ci;qi+1=di+1-wi*qi;m9=q9/p9; for(i=8;i>=0;i-)mi=(qi-ci*mi+1)/pi; f1=m3*pow(x4-xo,3)/(6.0*h3)+m4*pow(

6、xo-x3,3)/(6.0*h3)+(y3-m3*pow(h3,2)/6.0)*(x4-xo)/h3+(y4-m4*pow(h3,2)/6.0)*(xo-x3)/h3;f2=m4*pow(xo-x3,2)/(2.0*h3)-m3*pow(x4-xo,2)/(2.0*h3)-(y3-m3*pow(h3,2)/6.0)/h3+(y4-m4*pow(h3,2)/6.0)/h3;printf("f(4.563)=%lf f'(4.563)=%lfn",f1,f2);（3）运行结果题目3.用Romberg算法求.（1）算法说明引入记号Tm(k)称为Romberge序列,当|T

7、l(0)-Tl+1(0)|< 0.00001就停止计算: （2）程序清单#include<math.h>#include<stdio.h>main() int j=1,m,k,i; float a=1,b=3; float c,d=0,f,x; float fa=pow(3,a)*pow(a,1.4)*(5*a+7)*sin(pow(a,2);/*函数在积分两端点值*/ float fb=pow(3,b)*pow(b,1.4)*(5*b+7)*sin(pow(b,2); float t100100; /*定义二维数组存放迭代值*/ t10=(b-a)*(fb+fa

8、)/2; do for(i=1;i<=pow(2,j-1);i+) x=a+(2*i-1)*(b-a)/(pow(2,j); f=pow(3,x)*pow(x,1.4)*(5*x+7)*sin(pow(x,2); d=d+f; c=(b-a)*d/(pow(2,j-1); /*求TM+1K-1*/ t1j=(t1j-1+c)/2; for(m=1;m<=j;m+) for(k=1;k<=j-m+1;k+) tm+1k-1=(pow(4,m)*tmk-tmk-1)/(pow(4,m)-1); d=0; j+; while(abs(tm-10-tm0)>=0.00001);

9、 printf("%f",tm0);（3）运行结果题目4. 用定步长四阶Runge-Kutta求解h=0.0005，打印yi(0.025) ， yi(0.045) ， yi(0.085) ， yi(0.1) ，(i=1，2，3)（1）算法说明四阶Runge-Kutta方法：（2）程序清单#include<stdio.h>float y3=0,0,0,f34;float b43,a4=0.025,0.045,0.085,0.10;main ( )int i,m,j=0,n=1; float h=0.0005;dofor(m=0;m<3;m+)for(i=0;

10、i<=3;i+)fmi=0;f00=h;f10=h*y2;f20=h*(1000-1000*y1-100*y2);f01=h;f11=h*(y2+(1/3.0)*f00);f21=h*(1000-1000*(y1+(1/3.0)*f10)-100*(y2+(1/3.0)*f10);f02=h;f12=h*(y2+(1/3.0)*f10+f11);f22=h*(1000-1000*(y1+(1/3.0)*f20+f21)-100*(y2+(1/3.0)*f20+f21);f03=h;f13=h*(y2+f10-f11+f12);f23=h*(1000-1000*(y1+f20-f21+f2

11、2)-100*(y2+f20-f21+f22);for(i=0;i<3;i+) yi=yi+(fi0+3*fi1+3*fi2+fi3)/8;if(n=50|n=90|n=170|n=200)for(i=0;i<3;i+) bji=yi;j+;n+;while(n<=200);printf("The result is:n");for(j=0;j<4;j+) for(i=0;i<3;i+) printf("y%d(%5.4f)=%.7f,",i+1,aj,bji); printf("n");（3）运行结果题

12、目5. 用列主元消去法求解Ax=b。（1）算法说明列主元消去法：对矩阵作恰当的调整，选取绝对值尽量大的元素作为主元素。然后把矩阵化为上三角阵。再进行回代，求出方程的解。 , , , ,（2）程序清单#include <math.h>#include <stdio.h>main()double a 9=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.563

13、849,-0.784165,1.112348,3.123124,-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,-0.113584,2.189736,2.031454,4.101011,19.897918,0.431637,-3.111223,2.121314,1.784317,0.718719,1.56

14、3849,1.836742,-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.713846,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

15、.213474,4.446782,40.00001; int i,j,n,fi,p;double s,w,da,h; w=1.4; /*输入松弛因子w*/ double xi9,bi910,x9;double b9=2.1874369,33.992318,-25.173417,0.84671695,1.784317,-86.612343,1.1101230,4.719345,-5.6784392; for(i=0;i<9;i+)xii=0;for(i=0;i<9;i+)for(j=0;j<9;j+)if(i=j)biij=0;else biij=-aij/aii;for(i=

16、0;i<9;i+)bii9=bi/aii;for(n=0;n<9;n+)for(i=0;i<9;i+)s=0;for(j=0;j<9;j+)s=s+biij*xij;s=s+bii9;xii=xii*(1-w)+w*s;printf("the resule of songchi is:n");for(i=0;i<9;i+)if (i=5) printf("n"); printf("X(%d)=%-10.6f",i,xii);printf("n");printf("the re

17、sulu of liezhu is:");printf("n"); for(i=0;i<9;i+) xi=0; for(p=0;p<9;p+) fi=p; da=app; for(i=p+1;i<9;i+) /*选主元素*/ if(fabs(aip)>fabs(da) fi=i; da=aip; if(fi!=p) for(i=0;i<9;i+) /*循环行交换*/ h=api; api=afii; afii=h; h=bp;bp=bfi;bfi=h;h=app;for(j=0;j<9;j+) /*将p行除以 */apj=apj/h;bp=bp/h;for(i=0;i<9;i+) /*将p行以外各行的p列变为零*/h=aip;for(j=0;j<