东南大学计算方法与实习实验报告

1、计算方法与实习实验报告学院:学号:姓名： 完成日期:实习题一4、设，已知其精确值为12(32-1N-1N+1)。1）编制按从大到小的顺序计算Sn的程序；2）编制按从小到大的顺序计算Sn的程序；3）按两种顺序分别计算S1000，S10000，S30000，并指出有效位数。l 实验代码C语言程序如下：#include #include using namespace std;int main()float Sn=0;int N;cinN;for(float j=2;j=N;j+)Sn=1/(j*j-1)+Sn;cout从小到大计算的结果为Sn=2;j-)Sn=1/(j*j-1)+Sn;cout从大

2、到小计算的结果为Snendl;return 0;l 运行窗口实习题二1、用牛顿法求下列方程的根：1）l 实验代码C语言程序代码如下：#include #include #define N 100#define eps 1e-6#define eta 1e-8using namespace std;float Newton(float f(float),float fl(float),float x0)float x1,d;int k=0;dox1=x0-f(x0)/fl(x0);if(k+N|fabs(fl(x1)eps)cout发散endl;break;d=fabs(x1)1?x1-x0:(

3、x1-x0)/x1;x0=x1;coutx=x0eps&fabs(f(x1)eta);return x1;float f(float x)return x+log10(x)-2;float fl(float x)return 1+1/x;void main()float x0,y0;cinx0;y0=Newton(f,fl,x0);cout方程的根为y0endl;l 运行窗口实习题三1、用列主元消去法解方程组：1） l 实验代码C语言程序代码如下:#include#includeusing namespace std;void ColPivot(float *c,int n,float x)i

4、nt i,j,t,k;float p;for(i=0;i=n-2;i+)k=i;for(j=i+1;j(fabs(*(c+k*(n+1)+i) k=j;if(k!=j)for(j=i;j=n;j+)p=*(c+i*(n+1)+j);*(c+i*(n+1)+j)=*(c+k*(n+1)+j);*(c+k*(n+1)+j)=p;for(j=i+1;j=n-1;j+)p=(*(c+j*(n+1)+i)/(*(c+i*(n+1)+i);for(t=i;t=0;i-)for(j=n-1;j=i+1;j-)(*(c+i*(n+1)+n)-=xj*(*(c+i*(n+1)+j);xi=*(c+i*(n+1)

5、+n)/(*(c+i*(n+1)+i);int main()void ColPivot(float*,int,float);int i;float x4;float c45=1,1,0,3,4, 2,1,-1,1,1, 3,-1,-1,3,-3, -1,2,3,-1,4;ColPivot(c0,4,x);for(i=0;i=3;i+)printf(x%d=%fn,i,xi);return 0;l 运行窗口4、编写用追赶法解三对角线性方程组的程序，并解下列方程组：2），其中A10x10= -4 1 1 -4 1 1 -4 1 . . . . . . 1 -4 1 1 -4b= -27 -15 -

6、15l 实验代码C语言程序如下：#include#includeusing namespace std;void ColPivot(float *c,int n,float x)int i,j,t,k;float p;for(i=0;i=n-2;i+)k=i;for(j=i+1;j(fabs(*(c+k*(n+1)+i) k=j;if(k!=j)for(j=i;j=n;j+)p=*(c+i*(n+1)+j);*(c+i*(n+1)+j)=*(c+k*(n+1)+j);*(c+k*(n+1)+j)=p;for(j=i+1;j=n-1;j+)p=(*(c+j*(n+1)+i)/(*(c+i*(n+

7、1)+i);for(t=i;t=0;i-)for(j=n-1;j=i+1;j-)(*(c+i*(n+1)+n)-=xj*(*(c+i*(n+1)+j);xi=*(c+i*(n+1)+n)/(*(c+i*(n+1)+i);int main()void ColPivot(float*,int,float);int i;float x10;float c1011=-4,1,0,0,0,0,0,0,0,0,-27, 1,-4,1,0,0,0,0,0,0,0,-15, 0,1,-4,1,0,0,0,0,0,0,-15, 0,0,1,-4,1,0,0,0,0,0,-15, 0,0,0,1,-4,1,0,0

8、,0,0,-15, 0,0,0,0,1,-4,1,0,0,0,-15, 0,0,0,0,0,1,-4,1,0,0,-15, 0,0,0,0,0,0,1,-4,1,0,-15, 0,0,0,0,0,0,0,1,-4,1,-15, 0,0,0,0,0,0,0,0,1,-4,-15;ColPivot(c0,10,x);for(i=0;i=9;i+)printf(x%d=%fn,i,xi);return 0;l 运行窗口实习题四2、按下列数据Xi0.300.420.500.580.660.72Yi1.044031.084621.118031.156031.198171.23223作5次插值，并求X1=

9、0.46，X2=0.55，X3=0.60时的函数近似值。l 实验代码C语言程序如下：#include#include using namespace std;#define N 5void Difference(float x,float y,int n)float *f=new floatn+1;int k,i;for(k=1;k=n;k+)f0=yk;for(i=0;i=0;i-)a=a*(varx-xi)+yi;for(i=N-1;i=0;i-)b=b*(vary-xi)+yi;for(i=N-1;i=0;i-)c=c*(varz-xi)+yi;printf(Nn(%f)=%fn,var

10、x,a);printf(Nn(%f)=%fn,vary,b);printf(Nn(%f)=%fn,varz,c);return 0;l 运行窗口实习题六1、 用复化梯形公式和复化辛卜生公式计算积分I1（f）=dx。观察n为多少时，所得近似值具有6位有效数字。l 实验代码C语言程序如下：#include#include#include using namespace std;int main() double f1=0,f2=0,p=0,q=0,e=0.,a,b,h,x1,x2,x3; int i,j=1,n1=1,n2=1; cout计算结果如下：endl; cout以下为用复化梯形公式的计算

11、结果endl; coutsetw(15)f1setw(15)前后项差值endl; for(i=0;i+) a=0; b=3.14159/2; h=(b-a)/n1; for(j=1;j=n1;j+) x1=a+(j-1)*h; x2=a+j*h; f1+=h/2*(sqrt(1+cos(x1)*cos(x1)+sqrt(1+cos(x2)*cos(x2); if(fabs(f1-p)0.)break; coutsetw(15)f1setw(15)fabs(f1-p)endl; if(fabs(f1-p)0.)break; p=f1; n1+; f1=0; coutsetprecision(7)

12、 用 复 化 梯 形 公 式 计 算 得 f1=f1, 迭 代 次 数 为 ： n=n1endlendl; cout以下为用复化辛卜生公式的计算结果endl; coutsetw(15)f1setw(15)前后项差值endl; for(i=0;i+) a=0; b=3.14159/2; h=(b-a)/n2; for(j=1;j=n2;j+) x1=a+(j-1)*h; x2=a+j*h; x3=a+(j-0.5)*h; f2+=h/6*(sqrt(1+cos(x1)*cos(x1)+4*sqrt(1+cos(x3)*cos(x3)+sqrt(1+cos(x2)*cos(x2); if(fabs

13、(f2-q)0.)break; coutsetw(15)f2setw(15)fabs(f2-q)endl; if(fabs(f2-q)0.)break; q=f2; n2+; f2=0; coutsetprecision(7) 用 复 化 辛 卜 生 公 式 计 算 得 f1=f2, 迭 代 次 数 为 ： n=n2endl; return 0; l 运行窗口实习题七2、用四阶龙格-库塔公式（取h=0.1）解下列微分方程组初值问题：（1） l 实验代码C语言程序如下：#includevoid Runge_Kutta(float(*f1)(float x,float y1,float y2),f

14、loat(*f2)(float x,float y1,float y2),float a,float b,float y0,int N)float x=a,K1,K2,K3,K4,L1,L2,L3,L4,y1=y0,y2=y0;float h=(b-a)/N;int i;printf(x0=%fty10=%fty20=%fn,x,y1,y2);for(i=1;i=N;i+)K1=(*f1)(x,y1,y2);L1=(*f2)(x,y1,y2);K2=(*f1)(x+h/2,y1+h*K1/2,y2+h*L1/2);L2=(*f2)(x+h/2,y1+h*K1/2,y2+h*L1/2);K3=(*f1)(x+h/2,y1+h*K2/2,y2+h*L2/2);L3=(*f2)(x+h/2,y1+h*K2/2,y2+h*L2/2);K4=(*f1)(x+h,y1+h*K3,y2+h*L3);L4=(*f2)(x+h,y1+h*K3,y2+h*L3);y1=y1+h*(K1+2*K2+2*K3+K4)/6;y2=y2+h*(L1+2*L2+2*L3+L4)/6;x=a+i*h;printf(x%d=%fty1%d=%fty2%d=%fn,i,x,i,y1,i,y2);float f1(float x,float y1,float y2)return 1