计算方法上机实习题大作业(实验报告)

1、计算方法实验报告班级： 学号:姓名： 成绩:1舍入误差及稳定性一、实验目的(1)通过上机编程，复习巩固以前所学程序设计语言及上机操作指令；(2)通过上机计算，了解舍入误差所引起的数值不稳定性二、实验内容1000021、用两种不同的顺序计算n ,分析其误差的变化n 12、已知连分数f b0 a,利用下面的算法计算f :“ a? / b2 a3 /(. an / bn)ai 1dn bn, di bi (in 1,n 2,0) f d0写一程序，读入n,b0,6,，bn, a1,，an,计算并打印f3、给出一个有效的算法和一个无效的算法计算积分yndx0 4x 1(n0,1，,10)N 1一 1

2、3114、设SN,已知其精确值为1 3j 2 j 12 2 N N 1(1)编制按从大到小的顺序计算SN的程序(2)编制按从小到大的顺序计算SN的程序(3)按两种顺序分别计算S1000,S10000,S30000,并指出有效位数三、实验步骤、程序设计、实验结果及分析1000021、用两种不同的顺序计算n2,分析其误差的变化n 1(1)实验步骤：分别从110000和从100001两种顺序进行计算，应包含的头文件有 stdio.h和math.h(2)程序设计： a.顺序计算#include<stdio.h>#include<math.h>void main()double

3、sum=0;int n=1;while(1)sum=sum+(1/pow(n,2);if(n%1000=0)printf("sun%d=%-30f",n,sum);if(n>=10000)break;n+;printf("sum%d=%fn",n,sum);b.逆序计算#include<stdio.h>#include<math.h>void main()double sum=0;int n=10000;while(1)sum=sum+(1/pow(n,2);if(n%1000=0)printf("sum%d=%-

4、30f",n,sum);if(n<=1)break;n-；printf("sum%d=%fn",n,sum);(3)实验结果及分析：程序运行结果：a.顺序计算b.逆序计算而逆序计算误结果分析：两种不同顺序计算结果是一样的，顺序计算误差从一开始就很小，差最开始十分大，后来结果正确。2、已知连分数 fbo a，计算f ：bia2 / b2 a3 /(an/bn)(1)实验步骤：a利用 dnbn,dibi (i n 1,n 2,.,0) , f d0,计算 fdi i(2)程序设计#include<stdio.h>#include<math.h&

5、gt;void main()int i=0,n;float a1024,b1024,d1024;printf("please input n,n=");scanf("%d",&n);printf("nplease input a1 to an:n");for(i=1;i<=n;i+)printf("a%d=",i);scanf("%f",&ai);printf("nplease input b0 to bn:n");for(i=0;i<=n;i+)

6、printf("b%d=",i); scanf("%f",&bi);dn=bn；for(i=n-1;i>=0;i-)di=bi+ai+1/di+1；printf("nf=%fn",d0);(3)实验结果程序运行结果：3、给出一个有效的算法和一个无效的算法计算积分yndx0 4x 1(n 0,1,10)(1)实验步骤利用C语言编写程序，分别使用数值稳定的和数值不稳定的计算公式所建立的递推公 式进行计算。(2)程序设计#include<stdio.h>#include<math.h>main()dou

7、ble y_0=(1/4.0)*log(5),y_1;double y_2=(1.0/55.0+1.0/1，0)/2,y_3;int n=1,m=10;printf("有效算法输出结果：n");printf("y0=%-20f",y_0);while(1)y_1=1.0/(4*n)+y_0/(-4.0);printf("y%d=%-20f",n,y_1);if(n>=10) break;y_0=y_1;n+;if(n%3=0) printf("n");printf("n无效算法的输出结果：n&quo

8、t;);printf("y10=%-20f",y_2);while(1)y_3=1.0/n-4.0*y_2;printf("y%d=%-20f",m-1,y_3);if(m<=1) break;y_2=y_3;m-;if(m%2=0) printf("n");(3)实验结果及分析程序运行结果：结果分析：无效算法数值不稳定，误差造成的影响特别大N 1, , ,13114、设SN2J- ,已知其精确值为1 士，j 2 j 12 2 N N 1(1)实验步骤先编程按从大到小的顺序计算SN的程序，再编程按从小到大的顺序计算SN的程序,然

9、后按两种顺序分别计算S1000 , §0000 , S30000。(2)程序设计#include<stdio.h>main()int N;double SN30000;SN30000=(3.0/2.0-1.0/30000.0-1/30001.0)/2.0;for(N=30000;N>=2;N-)SNN-1=SNN-1.0/(N*N-1);printf(" 从 入 到 小 顺 序 计 算nSN1000=%fnSN10000=%fnSN30000=%fn",SN1000,SN10000,SN30000);SN2=(3.0/2-1.0/2.0-1/3.

10、0)/2.0;for(N=3;N<=30000;N+)SNN=SNN-1+1.0/(N*N-1);printf(" 从 小 到 大 顺 序 计 算nSN1000=%fnSN10000=%fnSN30000=%fn",SN1000,SN10000,SN30000);(3)实验结果及分析程序运行结果：结果分析：不同顺序计算所得结果是一样的。四、总结通过这次上机，学习了解了舍入误差在不同算法时对结果的影响不同，稳定的算法才能获得正确的结果。2方程求根一、实验目的(1)通过对二分法与牛顿迭代法做编程练习和上机运算，进一步体会二分法和牛顿法的不 同。(2)编写割线迭代法的程序，

11、求非线性方程的解，并与牛顿迭代法作比较。二、实验内容1、用牛顿法求下列方程的根(1) x2 ex 0(2) xex 10(3) lg x x 2 02、编写割线法程序求解第一问的方程三、实验步骤、程序设计、实验结果及分析1、牛顿法(1)实验步骤通过定义牛顿法求方程的子函数，用main函数调用子函数求根(2)程序设计#include <stdio.h>#include <math.h>typedef float (*p)(float);float ff1(float x)return x*x-exp(x);float ff2(float x)return x*exp(x)

12、-1;float ff3(float x)return log(x)+x-2;float answer(float(*p)(float)int k=2;float m=1,n=-1,x2,a,b,c;if (p=ff3)n=2;printf("x0 = %.4f, x1 = %.4f, ",m,n);while (1)if (fabs(m-n)<1e-4) break;a=p(n)*(n-m);b=p(n)-p(m);c=a/b;x2=n-c;m = n;n = x2;printf("x%d = %.4f, ",k,x2);k+;if (k%3=0

13、) printf("n");if (k%3!=0) printf("n");printf("iteration times: %d, roots: %.4fn ",k-2,n);return 0;main()printf("x*x-exp(x),n");answer(ff1);printf("x*exp(x)-1,n");answer(ff2);printf("lg(x)+x-2,n");answer(ff3);return 0;(3)实验结果及分析2、割线法(1)程序设计#

14、include<stdio.h>#include<math.h>float gexian(float,float);float f(float);main()int i,j；float x1=2.2;float x2=2,x3;scanf("%d",&i);if(i=1) printf("%f",x1);else if(i=2) printf("%f",x2);elsefor(j=3;j<=i;j+)x3=gexian(x1,x2);x1=x2;x2=x3;printf("%f"

15、;,gexian(x1,x2);float f(float x)return (x*x-exp(x);float gexian(float x1,float x2)return (x2-(f(x2)/(f(x2)-f(x1)*(x2-x1);(3)实验结果及分析X8 = 1.0000, Xl = -1.0000, x21 = -0,4621, xE3 = -fl.6757, x4 = -B.70&3, xS = -B.7034F x(61 = -fi.7g35. iterat ion t ine s ； 5 A i*oot s ； -0.7035pc*e>p<x>-l

16、 rx0 = 1 闾通皿 xl J = 1 ,B000, xL2J = -0.1135, x3 = 3.5495, x4J = -0.0807, xE51 =一目-明23. y 6 1 = 1.14E8, x7 = 0.234, 近凿=0.4545, xt91 = 0.5,?63, x13=回.E644. icCll = 0.5671, xll21 = 8,571, itevat ion t ime s ： 11. Mots ： B.S&71父)$3(2 , xtO = 1 .6068, xl = 2.0090. xC2J = 1.5906.“31-1.5S55. x4J - 1.55

17、72. xl - 1.5571, itevat ion time号2 4# i'-oots 1 1 .5571四、总结了解和学习了二分法和牛顿迭代法的思想以及程序设计的方法，比较了迭代法和牛顿 法的特点：牛顿法收敛速度较快，但对初值选取要求较高；割线法计算量少。3线性方程组数值解法一、实验目的(1)熟悉求解线性方程组的有关理论和方法；(2)会编制列主元消去法， LU分解法，雅可比及高斯-赛德尔迭代法的程序；(3)通过实际计算，进一步了解各种方法的优缺点，选择合适的数值方法。二、实验内容1、用列主元消去法解方程组2、用LU分解法解方程组三、实验步骤、程序设计、实验结果及分析1、用列主元消

18、去法解方程组(1)程序设计#include<stdio.h>#include<math.h>void ColPivot(float*,int,float);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<=n-1;j+)if(fabs(*(c+j*(n+1)+i)>(fabs(*(c+k*(n+1)+i) k=j;if(k!=i)for(j=i;j<=n;j+)(p=*(c+i*(n+1)+j);*(c+i*(

19、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<=n;t+)*(c+j*(n+1)+t)-=p*(*(c+i*(n+1)+t);for(i=n-1;i>=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);void main()(int i;float x4;float c

20、45=1,1,03421,-1,1,1,3,-1,-1,3,-3,-1,23-1,4;ColPivot(c0,4,x);for(i=0;i<=3;i+)printf("x%d=%fn",i,xi);(2)实验结果及分析(1)题(2)题2、用LU分解法解方程组(1)程序设计#include<stdio.h>void main()float x4;int i;float a45=48,-24,0,-12,4, -24,24,12,12,4, 0,6,20,2,-2, -6,6,2,16,-2；void DirectLU(float*,int,float);Di

21、rectLU(a0,4,x);for(i=0;i<=3;i+)printf("x%d=%fn",i,xi);)void DirectLU(float*u,int n,float x)(int i,r,k;for(r=0;r<=n-1;r+)(for(i=r;r<=n;i+)for(k=0;k<=r-1;k+)*(u+r*(n+1)+i)-=*(u+r*(n+1)+k)*(*(u+k*(n+1)+i);for(i=r+1;i<=n-1;i+)(for(k=0;k<=r-1;k+)*(u+i*(n+1)+r)-=*(u+i*(n+1)+k)*

22、(*(u+k*(n+1)+r);*(u+i*(n+1)+r)/=*(u+r*(n+1)+r);)for(i=n-1;i>=0;i-)(for(r=n-1;r>=i+1;r-)*(u+i*(n+1)+n)-=*(u+i*(n+1)+r)*xr;xi=*(u+i*(n+1)+n)/(*(u+i*(n+1)+i);)四、总结掌握了用列主元消去法和LU分解法求解方程组程序编写的技巧。4插值法一、实验目的(1)熟悉拉格朗日插值法多项式和牛顿插值多项式，注意其不同点；(2)掌握三次样条插值解决一些实际问题。二、实验内容1、按所给数据做二次插值，并求给定点的函数值2、按所给数据做五次插值，并求给

23、定点的函数值3、牛顿前插公式计算函数值三、实验步骤、程序设计、实验结果及分析1、二次插值(1)程序设计#include<stdio.h>float Lagrange(float x,float y,float xx,int n)/n 为(n+1)次插值;int i,j；float *a,yy=0;a=new floatn; for(i=0;i<=n-1;i+) ai=yi；for(j=0;j<=n-1;j+)if(j!=i)ai*=(xx-xj)/(xi-xj);yy+=ai; delete a; return yy; void main() float x5=-3.0

24、,-1.0,1.0,2.0,3.0; float y5=1.0,1.5,2.0,2.0,1.0; float xx1=-2,xx2=0,xx3=2.75,yy1,yy2,yy3;yy1=Lagrange(x,y,xx1,3);yy2=Lagrange(x,y,xx2,3);yy3=Lagrange(x,y,xx3,3);printf("x1=%-20f,y1=%fn",xx1,yy1);printf("x2=%-20f,y2=%fn",xx2,yy2);printf("x3=%-20f,y3=%fn",xx3,yy3);(2)实验结果

25、2、五次插值(1)程序设计#include<stdio.h>float Lagrange(float x,float y,float xx,int n)/n 为(n+1)次插值;int i,j;float *a,yy=0;a=new floatn; for(i=0;i<=n-1;i+)ai=yi；for(j=0;j<=n-1;j+)if(j!=i)ai*=(xx-xj)/(xi-xj);yy+=ai;delete a;return yy;void main()float x6=0.30,0.42Q50Q58Q66Q72;float y6=1.04403,1.08462,

26、1.11803,1.15603,1.19817,1.23223;float xx1=0.46,xx2=0.55,xx3=0.60,yy1,yy2,yy3;yy1=Lagrange(x,y,xx1,6);yy2=Lagrange(x,y,xx2,6);yy3=Lagrange(x,y,xx3,6);printf("x1=%-20f,y1=%fn",xx1,yy1);printf("x2=%-20f,y2=%fn",xx2,yy2);printf("x3=%-20f,y3=%fn",xx3,yy3);(2)实验结果3、牛顿前插公式计算函数

27、值(1)程序设计#include<stdio.h>#define N 3void Differencefloat y,float f44,int n) int k,i;f00=y0;fi0=yi;f20=y2;f30=y3；for(k=1;k<=n;k+)for(i=0;i<=(N-k);i+)fik=fi+1k-1-fik-1;return;void main()int i,k=1;float a,b=1,m=21.4,t=1.4,f44=0;float x5=20,21,22,23,24;float y5=1.30103,1.32222,1.34242,1.3617

28、3,1.38021;Difference(y,f,N);a=f00;for(i=1;i<=N;i+)k=k*i;b=b*(t-i+1);a=a+b*f0i/k;printf("x(k)n");for (i=0;i<=4;i+)printf( "%-20f",xi);printf("ny(k)n");for (i=0;i<=4;i+)printf("%-20f",yi);for(k=1;k<=3;k+)printf("nF(%d)n ",k);for(i=0;i<=(

29、3-k);i+)printf("%-20f",fik);printf ("n");printf("f(%f)=%-20f",m,a);printf ("n");(2)实验结果K(k>W.眄000021.0000022.独0鲍鲍23,0萌萌01,WAWA“k，L.3010301,3222261.3424201.3617301.386210Ml S.Q21170目.雕麒胭F<2) -0.000?©-0,000890F(3> O.MLEn f<21.400000>=i.350413

30、 请按任意键继续- - - H.01931B四、总结学习了插值法，学会了利用插值法编程求多项式的解，可以求解很多问题，让求解多 项式解变得非常简单。5曲线拟合一、实验目的(1)了解最小二乘法的基本原理，通过计算机解决实际问题；(2) 了解超定方程组的最小二乘解法。二、实验内容1、分别用抛物线y a bx cx2和指数曲线yaebx拟合所给数据，并比较这两个拟合函数的优劣。22、按所给实验数据，用形如 y a bx的抛物线进行最小二乘拟合。三、程序设计、结果分析 2 bx1、分别用抛物线y a bx cx和指数曲线y ae拟合所给数据a.抛物线(1)程序设计：#include<stdio.

31、h>#include<math.h>void main()int i;float a3;float x15=1,1.5,2,2.5,335,4,4.5,5,5.5,6,6.5,775,8;floaty15=33.4,79.50,122.65,159.05,189.15,214.15,238.65,252.50,267.55,280.50,296.65,301.40,310. 40,318.15,325.15;void Approx(float口,float,int,int,float);Approx(x,y,15,2,a);for(i=0;i<=2;i+)printf(

32、"a%d=%fn",i,ai);void Approx(float x口,float y口,int m,int n,float a口)int i,j,t;float *c=new float(n+1)*(n+2);float power(int,float);void ColPivot(float *,int,float);for(i=0;i<=n;i+)for(j=0;j<=n;j+)*(c+i*(n+2)+j)=0;for(t=0;t<=m-1;t+)*(c+i*(n+2)+j)+=power(i+j,xt);*(c+i*(n+2)+n+1)=0;fo

33、r(j=0;j<=m-1;j+)*(c+i*(n+2)+n+1)+=yj*power(i,xj);ColPivot(c,n+1,a);delete c;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<=n-1;j+)if(fabs(*(c+j*(n+1)+i)>(fabs(*(c+k*(n+1)+i)k=j;if(k!=i)for(j=i;j<=n;j+)p=*(c+i*(n+1)+j);*(c+i*(n+1)+j)=*(c+k*

34、(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<=n;t+)*(c+j*(n+1)+t)-=p*(*(c+i*(n+1)+t);for(i=n-1;i>=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);float power(int i,float v)float a=1;while(i-)a*=v;r

35、eturn a;(2)实验结果aE01=-45.333195 all 3=94.330103 aL2=-fc.13157请按任意键继续. . .2、最小二乘拟合(1)程序设计#include<stdio.h>#include<math.h>void main()int i,n;float a2;float x15=1,1.5,2,2.5,335,4,4.5,5,5.5,6,6.5,775,8,z15;floaty15=33.4,79.50,122.65,159.05,189.15,214.15,238.65,252.50,267.55,280.50,296.65,301

36、.40,310.40,318.15,325.15;for(n=0;n<=14;n+)增加了数组 z;zn=log(yn/xn);void Approx(float,float,int,int,float);Approx(x,z,15,1,a);变成一次拟合；/for(i=0;i<=1;i+)/printf("a%d=%fn",i,ai);printf("a=exp(a0)=%fn",exp(a0);printf("b=-a1=%fn",-a1);void Approx(float x,float y,int m,int n

37、,float a)int i,j,t;float *c=new float(n+1)*(n+2);float power(int,float);void ColPivot(float *,int,float);for(i=0;i<=n;i+)for(j=0;j<=n;j+)*(c+i*(n+2)+j)=0;for(t=0;t<=m-1;t+)*(c+i*(n+2)+j)+=power(i+j,xt);*(c+i*(n+2)+n+1)=0;for(j=0;j<=m-1;j+)*(c+i*(n+2)+n+1)+=yj*power(i,xj);ColPivot(c,n+1,a

38、);delete c;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<=n-1;j+)if(fabs(*(c+j*(n+1)+i)>(fabs(*(c+k*(n+1)+i) k=j;if(k!=i)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)

39、/(*(c+i*(n+1)+i);for(t=i;t<=n;t+)*(c+j*(n+1)+t)-=p*(*(c+i*(n+1)+t);)for(i=n-1;i>=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);)float power(int i,float v)(float a=1;while(i-)a*=v;return a;)(2)实验结果四、总结通过曲线拟合，最小二乘法的基本原理的学习， 我学会了利用计算机解决现实实际问 题中

40、的曲线拟合。6数值积分一、实验目的（1）通过实际计算体会各种方法的精确度;(2)会编写用龙贝格算法求定积分的程序。二、实验内容编写复化柯特斯求积分公式，并计算例题1例题2,观察n为多少时有6位有效数字。三、程序设计、实验结果及分析(1)程序设计#include<stdio.h>#include<math.h>float Cotes(float(*f)(float),float a,float b,int n)int k;float c,c1=0,c2,c3,c4;float h=(b-a)/n;c2=(*f)(a+h/4);c3=(*f)(a+h/2);c4=(*f)(

41、a+3*h/4);for(k=1;k<=n-1;k+)c1+=(*f)(a+k*h);c2+=(*f)(a+k*h+h/4);c3+=(*f)(a+k*h+h/2);c4+=(*f)(a+k*h+3*h/4);c=h/90*(7*(*f)(a)+(*f)(b)+14*c1+32*c2+12*c3+32*c4);return c;float f(float x)return 1/sqrt(1+x*x*x);void main()int i,n=4;float c;for(i=0;i<=4;i+)c=Cotes(f,0,1,n);printf("C(%d)=%fn"

42、,n,c);n*=2;#include<stdio.h>#include<math.h>float Cotes(float(*f)(float),float a,float b,int n)int k;float c,c1=0,c2,c3,c4;float h=(b-a)/n;c2=(*f)(a+h/4);c3=(*f)(a+h/2);c4=(*f)(a+3*h/4);for(k=1;k<=n-1;k+)c1+=(*f)(a+k*h);c2+=(*f)(a+k*h+h/4);c3+=(*f)(a+k*h+h/2);c4+=(*f)(a+k*h+3*h/4);c=h

43、/90*(7*(*f)(a)+(*f)(b)+14*c1+32*c2+12*c3+32*c4); return c;float f(float x)/ return 1/sqrt(1+x*x*x);if (x=0)return 1;else return sin(x)/x;void main()int i,n=4;float c;for(i=0;i<=4;i+)/c=Cotes(f,0,1,n);c=Cotes(f,0,5,n);printf("C(%d)=%fn",n,c);n*=2;(2)实验结果及分析四、总结学习了复化辛卜生公式， 自适应梯形公式，龙贝格算法，运

44、用求解定积分并控制精度 的方法。7常微分方程数值解法一、实验目的(1)熟悉求解常微分方程初值问题的有关方法和理论，主要是改进欧拉公式，四阶龙格库塔法和阿当姆斯方法；(2)编制上述方法计算机程序，包括求解微分方程组的计算程序；(3)针对实习题编制程序，并上机计算其所需要的结果；(4)体会各种解法的功能，优缺点及适用场合，会选取适当的求解方法。二、实验内容1、分别用改进欧拉法与四阶龙格-库塔公式(取h 0.1)求解下列微分方程初值问题2、用四阶龙格-库塔公式(取h 0.1)解下列微分方程组初值问题三、实验步骤、程序设计、实验结果及分析1、分别用改进欧拉法与四阶龙格-库塔公式(取h 0.1)求解下列

45、微分方程初值问题a.改进欧拉法(1)程序设计#include<stdio.h>void ModEuler(float(*f)(float,float),float x0,float y0,float xn,int n)int i;float yp,yc,x=x0,y=y0,h=(xn-x0)/n;printf("x0=%fty0=%fn",x,y);for(i=1;i<=n;i+)yp=y+h*(*f)(x,y);x=x0+i*h;yc=y+h*(*f)(x,yp);y=(yp+yc)/2;printf("x%d=%fty%d=%fn"

46、,i,x,i,y);float f(float x,float y)/return x*x+y*y;/题(1)/return 1/(1+y*y);/题(2)return y-2*x/y;/题(3)void main()float xn=1.0,x0=0,y0=1;ModEuler(f,x0,y0,xn,10);(2)实验结果及分析(1)(3)b.四阶龙格-库塔公式(1)程序设计#include<stdio.h>void Runge_Kutta(float(*f)(float x,float y),float afloat b,float y0,int N) float x=a,y=y0,K1,K2K3K4;float h=(b-a)/N;int i;printf("x0=%fty0=%fn",x,y);for(i=1;i<=N;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