计算方法实习报告

5.29.20041一．给出一个有效的算法和无效的算法计算积分y(n)=∫(x^n)/(4x+1)dx,n=0,1,2,…,10,积分限为（0，1）1．有效算法利用递推公式y(n)=-y(n-1)/4+1/(4n)，取y0=(log5)/4程序为：#includeiostream.h#includemath.hvoidmain(){doubley0,y1;y0=1/4.0*log(5.0);couty0=y0;for(intn=1;n=10;n++){y1=-1.0/4.0*y0+1.0/(4.0*n);coutyn=y1;y0=y1;if(n%3==0)coutendl;}}其结果为：y0=0.402359y1=0.14941y2=0.0876475y3=0.0614215y4=0.0471446y5=0.0382138y6=0.0321132y7=0.027686y8=0.0243285y9=0.0216957y10=0.0195761Pressanykeytocontinue2．无效算法利用递推公式y(n-1)=-4y(n)+1/n,又由广义积分中值定理可得y(n)=1/((4n+1)(4ζ+1)),ζ∈(0,1),则1/(5(n+1))y(n)1/(n+1),所以可取y(n)≈[1/(5(n+1))+1/(n+1)]/2=3/(5(n+1)程序为：#includeiostream.h#includemath.hvoidmain(){floaty9,y10;y10=3.0/55.0;couty10=y10;for(intn=9;n=0;n--){y9=-4.0*y10+1.0/(n+1);coutyn=y9;y10=y9;if(n%3==0)coutendl;}}其结果为：y10=0.0545455y9=-0.118182y8=0.583838y7=-2.21035y6=8.98427y5=-35.7704y4=143.282y3=-572.877y2=2291.84y1=-9166.86y0=36668.43．心得体会5.29.20042由有效算法与无效算法的结果可以知道：有效算法的误差的传递是逐步缩小的，而无效算法的误差的传递是逐步扩大的。因此有效算法是数值稳定的算法。二．用牛顿法求下列方程lnx+x-2=0的根。1．算法给定初值X0，ε为根的容许误差，η为|F（X）|的容许误差，N为迭代次数的容许值。(1).如果F（X0）的导数等于零或者迭代次数大于N，则算法失败，结束；否则执行2。(2).计算X1=X0-F（X0）/F’（X0）。(3).若|X1-X0|ε,或者|F（X1）|η,则输出X1，程序结束，否则执行4。(4).令X0=X1，转向1。2.程序#includeiostream.h#includemath.h#defineEPS1e-8#defineETA1e-8#defineN100floatf(floatx){returnlog(x)+x-2;}floatf1(floatx){return1+1/x;}voidmain(){floatx1,x0=1.5,d;coutx0=x0;for(inti=0;i=N;){if(f1(x0)==0)break;elsex1=x0-f(x0)/f1(x0);d=(x1-x0);if(fabs(f(x1))ETA&&fabs(d)EPS){coutx1endl;break;}elsex0=x1;i++;coutxi=x0;while(i%3==0)coutendl;}coutTherootoftheequationisx=x0endl;}5.29.200433.结果x0=1.5x1=1.55672x2=1.557151.55715Therootoftheequationisx=1.557154．心得体会通过牛顿迭代算法的编程实习，我了解到牛顿法把非线性方程线性化，加快了收敛速度。三．编写一个用牛顿前插公式计算函数值的程序，要求先输出差分表，再计算x点的函数值，并应用于下面的问题：xi2021222324yi1.301031.322221.342421.361731.38021求x=21.4时的插值多项式的值。1．算法(1).输入n=4,xi,yi(i=0,1,2,3,4).(2).计算各阶差分f00,f10,f20,f30,f40。(3).计算函数值：f00+f10*t+f20*t*(t-1)/2+f30*t*(t-1)*(t-2)/6+f40*t*(t-1)*(t-2)*(t-3)/24;2．程序#includeiostream.hfloatNewTon(floatf00,floatf10,floatf20,floatf30,floatf40,floatt){returnf00+f10*t+f20*t*(t-1)/2+f30*t*(t-1)*(t-2)/6+f40*t*(t-1)*(t-2)*(t-3)/24;}voidmain(){floatf1[4],f2[3],f3[2],f4;floatf0[5]={1.30103,1.32222,1.34242,1.36173,1.38021};floatx[5]={20,21,22,23,24};cout差分表endl;for(inti=0;i=4;i++){coutxi=x[i];}coutendl;for(i=0;i=4;i++){coutf0[i];}coutendl1阶差分;for(i=0;i=3;i++){f1[i]=f0[i+1]-f0[i];coutf1[i];5.29.20044}coutendl2阶差分;for(i=0;i=2;i++){f2[i]=f1[i+1]-f1[i];coutf2[i];}coutendl3阶差分;for(i=0;i=1;i++){f3[i]=f2[i+1]-f2[i];coutf3[i];}coutendl4阶差分;f4=f3[1]-f3[0];coutf4endl;floatt=21.4-x[0];floatN4=NewTon(f0[0],f0[1],f0[2],f0[3],f0[4],t);coutN(21.4)=N4endl;}3．结果差分表x0=20x1=21x2=22x3=23x4=241.301031.322221.342421.361731.380211阶差分0.02118990.02020.019310.01848012阶差分-0.000989914-0.000890017-0.0008299353阶差分9.98974e-0056.00815e-0054阶差分-3.98159e-005N(21.4)=3.482684．心得体会由于此题目选用的节点是等距的，所以可用等距节点的插值公式代替节点的任意分布的牛顿插值公式，这样就可以简化插值公式，同时可以避免作除法运算，而且程序算法的设计也相对比较简单。四．用改进欧拉方法（取h=0.1）求解下列微分方程初值问题：y’=(x^2+y^2),x∈[0,1.0]y(0)=01．算法解一阶常微分方程初值问题y’=f(x,y),y(x0)=y0,x∈[a,b],将区间[a,b]做n等分，取步长h=(b-a)/n.欧拉公式为：y(i+1)=y(i)+hf(x(i),y(i))梯形公式为：y(i+1)=y(i)+h[f(x(i),y(i))+f(x(i+1),y(i+1))]/2改进欧拉法，采用公式：ŷ(i+1)=y(i)+hf(x(i),y(i))y(i+1)=y(i)+h[f(x(i),y(i))+f(x(i+1),ŷ(i+1))]/22．程序#includeiostream.h#defineN105.29.20045floatf(floatx,floaty){returnx*x+y*y;}voidmain(){floatx0=0,y0=0,a=0,b=1.0;floath=(b-a)/N;coutx[0]=x0y[0]=y0endl;floatx1=x0,y1=y0;for(inti=0;i=N;){y1=y0+h*f(x0,y0);/*欧拉公式*/x1=x0+h;y1=y0+h*(f(x0,y0)+f(x1,y1))/2;/*梯形公式*/y0=y1;x0=x1;i++;coutx[i]=x1y[i]=y1endl;}}3．结果x[0]=0y[0]=0x[1]=0.1y[1]=0.0005x[2]=0.2y[2]=0.00300013x[3]=0.3y[3]=0.00950303x[4]=0.4y[4]=0.0220247x[5]=0.5y[5]=0.0426214x[6]=0.6y[6]=0.0734421x[7]=0.7y[7]=0.116817x[8]=0.8y[8]=0.175396x[9]=0.9y[9]=0.252374x[10]=1y[10]=0.35183x[11]=1.1y[11]=0.4792944．心得体会改进欧拉公式和欧拉公式相比，它们同为单步显式公式，但是前者的局部截断误差比后者的局部截断误差高一阶。欧拉公式具有一阶精度，梯形公式和改进欧拉公式具有二阶精度。五．用阿当姆斯方法（取h=0.1）求解下列微分方程初值问题（用四阶龙格-库塔公式提供出发值）：y’=0.1(x^3+y^2),x∈[0,1]y(0)=11.算法概要阿当姆斯方法是一种线性多步法，其四阶显式公式为y(i+1)=y(i)+h(55f(i)-59f(i-1)+37f(i-2)-9f(i-3))/245.29.20046四阶隐式公式为y(i+1)=y(i)+h(9f(i+1)+19f(i)-5f(i-1)+f(i-2))/24以上两个公式均具有四阶精度，将它们组成阿当姆斯预测-校正系统：ŷ(i+1)=y(i)+h(55f(i)-59f(i-1)+37f(i-2)-9f(i-3))/24f’(i+1)=f(x(i+1),ŷ(i+1))y(i+1)=y(i)+h(9f’(i+1)+19f(i)-5f(i-1)+f(i-2))/24f(i+1)=f(x(i+1),y(i+1))这是一个四步方法，计算y(i+1)时要用到y(i),y(i-1)，y(i-2),y(i-3),因此，它不是自开始的，一般借助于同阶的龙格-库塔公式为其提供出发值：y1,y2,y3.2.程序#includeiostream.h#includeconio.h#defineN10floatf(floatx,floaty){return0.1*(x*x*x+y*y*y);}voidmain(){floatyy[3],y,x,h,a,b,k1,k2,k3,k4,f1;a=0;b=1.0;y=1;x=0;yy[0]=1;h=(b-a)/N;inti;coutx[0]=xy[0]=yendl;for(i=1;i=3;i++){k1=f(x,y);k2=f(x+h/2,y+h*k1/2);k3=f(x+h/2,y+h*k2/2);k4=f(x+h,y+h*k3);y=y+h*(k1+2*k2+2*k3+k4)/6;x=x+h;yy[i]=y;coutx[i]=xy[i]=yendl;}for(i=3;i=N;i++){x=a+i*h;y=yy[3];yy[3]=yy[3]+h*(55*f(x,yy[3])-59*f(x-h,yy[2])+37*f(x-2*h,yy[1])-9*f(x-3*h,yy[0]))/24;f1