计算方法实习报告

一．给出一种有效的算法和无效的算法计算积分y(n)=∫(x^n)/(4x+1)dx,n=0,1,2,…,10,积分限为（0，1）有效算法运用递推公式y(n)=-y(n-1)/4+1/(4n)，取y0=(log5)/4程序为：#include<iostream.h>#include<math.h>voidmain(){doubley0,y1; y0=1/4.0*log(5.0); cout<<"y0="<<y0<<""; for(intn=1;n<=10;n++) { y1=-1.0/4.0*y0+1.0/(4.0*n); cout<<"y"<<n<<"="<<y1<<"";y0=y1; if(n%3==0)cout<<endl; }}其成果为：y0=0.402359y1=0.14941y2=0.0876475y3=0.0614215y4=0.0471446y5=0.0382138y6=0.0321132y7=0.027686y8=0.0243285y9=0.0216957y10=0.0195761Pressanykeytocontinue无效算法运用递推公式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)程序为：#include<iostream.h>#include<math.h>voidmain(){floaty9,y10; y10=3.0/55.0; cout<<"y10="<<y10<<""; for(intn=9;n>=0;n--) { y9=-4.0*y10+1.0/(n+1); cout<<"y"<<n<<"="<<y9<<"";y10=y9; if(n%3==0)cout<<endl; }}其成果为：y10=0.0545455y9=-0.118182y8=0.583838y7=-2.21035y6=8.98427y5=-35.7704y4=143.282y3=-572.877y2=2291.84y1=-9166.86y0=36668.4心得体会由有效算法与无效算法的成果可以懂得：有效算法的误差的传递是逐渐缩小的，而无效算法的误差的传递是逐渐扩大的。因此有效算法是数值稳定的算法。二．用牛顿法求下列方程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.程序#include<iostream.h>#include<math.h>#defineEPS1e-8#defineETA1e-8#defineN100floatf(floatx){ returnlog(x)+x-2;}floatf1(floatx){ return1+1/x;}voidmain(){ floatx1,x0=1.5,d; cout<<"x0="<<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) { cout<<x1<<endl; break; } else x0=x1; i++; cout<<"x"<<i<<"="<<x0<<""; while(i%3==0)cout<<endl; } cout<<"Therootoftheequationisx="<<x0<<endl;}3.成果x0=1.5x1=1.55672x2=1.557151.55715Therootoftheequationisx=1.557154．心得体会通过牛顿迭代算法的编程实习，我理解到牛顿法把非线性方程线性化，加紧了收敛速度。编写一种用牛顿前插公式计算函数值的程序，规定先输出差分表，再计算x点的函数值，并应用于下面的问题：xi2021222324yi1.301031.322221.342421.361731.38021求x=21.4时的插值多项式的值。算法(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;程序#include<iostream.h>floatNewTon(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++) { cout<<"x"<<i<<"="<<x[i]; } cout<<endl<<"";for(i=0;i<=4;i++) { cout<<""<<f0[i]; } cout<<endl<<"1阶差分"<<""; for(i=0;i<=3;i++) { f1[i]=f0[i+1]-f0[i]; cout<<f1[i]<<""; } cout<<endl<<"2阶差分"<<"";for(i=0;i<=2;i++) { f2[i]=f1[i+1]-f1[i]; cout<<f2[i]<<""; } cout<<endl<<"3阶差分"<<"";for(i=0;i<=1;i++) { f3[i]=f2[i+1]-f2[i]; cout<<f3[i]<<""; } cout<<endl<<"4阶差分"<<""; f4=f3[1]-f3[0]; cout<<f4<<endl; floatt=21.4-x[0]; floatN4=NewTon(f0[0],f0[1],f0[2],f0[3],f0[4],t); cout<<"N(21.4)="<<N4<<endl;}成果差分表x0=20x1=21x2=22x3=23x4=241.301031.322221.342421.361731.380211阶差分0.02118990.02020.019310.01848012阶差分-0.-0.-0.3阶差分9.98974e-0056.00815e-0054阶差分-3.98159e-005N(21.4)=3.48268心得体会由于此题目选用的节点是等距的，因此可用等距节点的插值公式替代节点的任意分布的牛顿插值公式，这样就可以简化插值公式，同步可以防止作除法运算，并且程序算法的设计也相对比较简朴。用改善欧拉措施（取h=0.1）求解下列微分方程初值问题：y’=(x^2+y^2),x∈[0,1.0]y(0)=0算法解一阶常微分方程初值问题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))]/2程序#include<iostream.h>#defineN10floatf(floatx,floaty){ returnx*x+y*y;}voidmain(){ floatx0=0,y0=0,a=0,b=1.0; floath=(b-a)/N; cout<<"x[0]="<<x0<<""<<"y[0]="<<y0<<endl; 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++; cout<<"x["<<i<<"]="<<x1<<""<<"y["<<i<<"]="<<y1<<endl; }}成果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.479294心得体会改善欧拉公式和欧拉公式相比，它们同为单步显式公式，不过前者的局部截断误差比后者的局部截断误差高一阶。欧拉公式具有一阶精度，梯形公式和改善欧拉公式具有二阶精度。五．用阿当姆斯措施（取h=0.1）求解下列微分方程初值问题（用四阶龙格-库塔公式提供出发值）：y’=0.1(x^3+y^2),x∈[0,1]y(0)=1算法概要阿当姆斯措施是一种线性多步法，其四阶显式公式为y(i+1)=y(i)+h(55f(i)-59f(i-1)+37f(i-2)-9f(i-3))/24四阶隐式公式为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.程序#include<iostream.h>#include<conio.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; cout<<"x[0]="<<x<<""<<"y[0]="<<y<<endl; 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; cout<<"x["<<i<<"]="<<x<<""<<"y["<<i<<"]="<<y<<endl; } 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=f(x+h,yy[3]); yy[3]=y+h*(9*f1+19*f(x,yy[3])-5*f(x-h,yy[2])+f(x-2*h,yy[1]))/24; yy[0]=yy[1]