东南大学计算方法上机报告实验报告完整版.doc

实习题11. 用两种不同的顺序计算，试分析其误差的变化解：从n=1开始累加，n逐步增大，直到n=10000；从n=10000开始累加，n逐步减小，直至1。算法1的C语言程序如下：#include#includevoid main()float n=0.0;int i;for(i=1;i=1;j-)m=m+1.0/(j*j);printf(%-7f,m);printf(n);运行后结果如下：结论：4.设，已知其精确值为。1）编制按从大到小的顺序计算的程序；2）编制按从小到大的顺序计算的程序；3）按2种顺序分别计算，并指出有效位数。解：1）从大到小的C语言算法如下：#include#include#includeusing namespace std;void main()float n=0.0;int i;int N;coutPlease input NN;for(i=N;i1;i-) n=n+1.0/(i*i-1);N=N-1;printf(%-100f,n);printf(n);执行后结果为：N=2时，运行结果为：N=3时，运行结果为：N=100时，运行结果为：N=4000时，运行结果为：2）从小到大的C语言算法如下：#include#include#includeusing namespace std;void main()float n=0.0;int i;int N;coutPlease input NN;for(i=2;i=N;i+)n=n+1.0/(i*i-1);printf(%-100f,n);printf(n);执行后结果为：N=2时，运行结果为：N=3时，运行结果为：N=100时，运行结果为：N=4000时，运行结果为：结论：通过比较可知：N的值较小时两种算法的运算结果相差不大，但随着N的逐渐增大，两种算法的运行结果相差越来越大。3）用从小到大的C语言算法计算，结果如下：N=1000时，运行结果：N=10000时，运行结果：N=30000时，运行结果：用从大到小顺序的C语言算法分别计算，结果如下：N=1000时，运行结果：N=10000时，运行结果：N=30000时，运行结果：实习题21. 用牛顿法法解下列方程的根：1）；3）。解：1）给定初始值，为根的容许误差，为的容许误差，N为迭代次数的容许值。如果或迭代次数大于N，则算法失败，结束，转而执行。在或时，则输出，程序结束，否则执行。令，转向。根据大致的草图，判断出方程只有一个根，且在（-2，0），因此可以取初始值-1.000000，执行程序。C语言算法如下：#include #include#include using namespace std;#define N 100#define eps 1e-6#define eta 1e-8float Newton(float(*f)(float),float(*f1)(float),float x0)float x1,d;int k=0;dox1=x0-(*f)(x0)/(*f1)(x0);if(k+N|fabs(*f1)(x1)eps)printf(n Newton迭代发散);break;d=fabs(x1)eps&fabs(*f)(x1)eta);return x1;float f(float x)return float(x*x-exp(x);float f1(float x)return float(2*x-exp(x);void main()float x0,y0;printf(请输入迭代初值x0n);scanf(%f,&x0);printf(x(0)=%fn,x0);y0=Newton(f,f1,x0);printf(方程的根为：%fn,y0);运行结果为：由此可见，利用牛顿法求该方程的近似解为-0.703467，精确到小数点后六位，同时可以发现，迭代次数只有4次，若是用二分法进行求解，那么迭代次数将会多很多，可见，牛顿法的收敛速度还是较快的。3）2、编写一个割线法的程序，求解上述各方程。1）；解：解：1）给定初始值，为根的容许误差，为的容许误差，N为迭代次数的容许值。如果或迭代次数大于N，则算法失败，结束，转而执行。在或时，则输出，程序结束，否则执行。令，转向。根据大致的草图，判断出方程只有一个根，且在（-2，0），因此可以取初始值-1.000000，执行程序。C语言算法如下：#include #include#include using namespace std;#define N 100#define eps 1e-6#define eta 1e-8float Newton(float(*f)(float),float(*f1)(float),float x0)float x1,d;int k=0;dox1=x0-(*f)(x0)/(*f1)(x0);if(k+N|fabs(*f1)(x1)eps)printf(n Newton迭代发散);break;d=fabs(x1)eps&fabs(*f)(x1)eta);return x1;float f(float x)return float(x*x-exp(x);float f1(float x)return float(2*x-exp(x);void main()float x0,y0;printf(请输入迭代初值x0n);scanf(%f,&x0);printf(x(0)=%fn,x0);y0=Newton(f,f1,x0);printf(方程的根为：%fn,y0);运行结果为：由此可见，利用牛顿法求该方程的近似解为-0.703467，精确到小数点后六位，同时可以发现，迭代次数只有4次，若是用二分法进行求解，那么迭代次数将会多很多，可见，牛顿法的收敛速度还是较快的。实习题31、用列主元消去法解方程组：1）解：将方程用增广矩阵 A|b=(aij)n(n+1)表示,分别进行消元和回代过程、消元过程：对于选主元；如果=0，则矩阵A奇异，程序结束；否则执行；如果，则交换第k行与第行对应元素位置；消元， 对，计算，对，计算；、回代过程 若，则矩阵A奇异，程序结束；否则执行； ；对，计算；C语言程序代码为：#include#include#includeusing namespace std;void 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+1,xi);void ColPivot(float *c,int n,float x)int i,j,k,t;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!=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)/(*(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)+n)/(*(c+i*(n+1)+i);运行结果为：所以原方程组的近似解为，此解与精确解几乎一致，因此可以发现，列主元消去法可以比较精确地求得多元一次方程组的解，比一般普通方法精确度要高。2、用LU分解法解方程组Ax=b，其中 ，解：、将方程组中的A分解成，其中L为单位下三角形矩阵，U为上三角形矩阵，则方程组化解为2个方程组Ly=b，Ux=y，具体算法如下：对j=1,2,3, n计算对i=2,3, , n计算对k=2,3，n：a.对j=k,k+1, ,n计算b.对i=k+1,k+2, ,n计算，对k=2,3, ,n计算，对k=n-1,n-2, ,2,1计算、程序代码如下：#includevoid 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);DirectLU(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;i=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=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);、运行结果为：4、编写用追赶法解三对角线性方程组的程序，并解下列方程组。2），其中 ，解：应用高斯消去法求解时，消元过程可以简化，得到同解的三角方程组，进行回代消去过程为： 回代过程为 C语言程序为：#includevoid main()float x10;int i;float a1011=-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,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;void Direct(float*,int,float);Direct(a0,10,x);for(i=0;i=9;i+)printf(x%d=%fn,i,xi);void Direct(float*u,int n,float x)int i,r,k;for(r=0;r=n-1;r+)for(i=r;i=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=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);运行结果为：所以原方程组的近似解为：x1 =8.705758, x2=7.823033, x3=7.586371, x4=7.522452, x5=7.503439, x6=7.491305, x7=7.461785, x8=7.355835, x9=6.961556, x10=5.4903895、分别用雅可比迭代法与高斯-赛德尔迭代法解下列方程组：1）解：用雅可比迭代法求该方程组、算法：设方程组的系数矩阵的对角线元素，M为迭代次数所容许的最大值，为容许误差。取初值向量，令k=0.对i=1,2,n计算如果，则输出，结束；否则执行。如果，则不收敛，终止程序；否则，转。、程序代码如下：#include#include#define eps 1e-6#define max 100void Jacobi(float *a,int n,float x)int i,j,k=0;float epsilon,s;float *y=new floatn;for(i=0;in;i+)xi=0;while(1)epsilon=0;k+;for(i=0;in;i+)s=0;for(j=0;jn;j+)if(j=i)continue;s+=*(a+i*(n+1)+j)*xj;yi=(*(a+i*(n+1)+n)-s)/(*(a+i*(n+1)+i);epsilon+=fabs(yi-xi);for(i=0;in;i+)xi=yi;if(epsilon=max)printf(迭代发散);return;delete y;void main()int i;float a45=10,-1,2,0,-11,0,8,-1,3,-11,2,-1,10,0,6,-1,3,-1,11,25;float x4;Jacobi(a0,4,x);for(i=0;i4;i+)printf(x%d=%fn,i,xi);、运行结果如下：用高斯-赛德尔迭代法求该方程组：、算法：设方程组的系数矩阵的对角线元素，M为迭代次数所容许的最大值，为容许误差。取初值向量，令k=0.对i=1,2,n计算如果，则输出，结束；否则执行。如果，则不收敛，终止程序；否则，转。、程序代码如下：#include#include#define N 500void main()int i;float x4;float c45=10,-1,2,0,-11,0,8,-1,3,-11,2,-1,10,0,6,-1,3,-1,11,25;void GauseSeidel(float *,int,float);GauseSeidel(c0,4,x); for(i=0;i=3;i+)printf(x%d=%fn,i,xi);void GauseSeidel(float *a,int n,float x)int i,j,k=1;float d,dx,eps;for(i=0;i=n-1;i+)xi=0.0;while(1)eps=0;for(i=0;i=n-1;i+)d=0;for(j=0;j=n-1;j+)if(j=i)continue;d+=*(a+i*(n+1)+j)*xj;dx=(*(a+i*(n+1)+n)-d)/(*(a+i*(n+1)+i);eps+=fabs(dx-xi);xi=dx;if(epsN)printf(迭代发散n);return;k+;、运行结果如下：分析：用雅可比迭代法和高斯-赛德尔迭代法解该方程组时，当二者精度相同时，高斯-赛德尔迭代法的迭代次数比雅可比迭代法的迭代次数要少8次，大大减少了计算量。实习题4 2、按下列数据xi0.300.420.500.580.660.72yi1.044031.084621.118031.156031.198171.23223作五次插值，并求x1=0.46, x2=0.55, x3=0.60时函数的近似值。解：、输入令、对， 计算C语言程序为：#includefloat Lagrange(float x,float y,float xx,int n)int i,j;float *a,yy=0;a=new floatn;for(i=0;in;i+)ai=yi;for(j=0;j A=11.522.533.544.5 55.566.577.58;B=33.479.5122.65159.05189.15214.15238.65252.5 267.55280.5296.65301.4310.4318.15325.15;S=zeros(1,5);T=zeros(3,1); for k=1:5S(k)=sum(A.(k-1);end for k=1:3T(k)=sum(A.(k-1).*B);end C=zeros(3,3); a=zeros(3,1); for i=1:3for j=1:3C(i,j)=S(i+j-1);endend a=CT; for k=1:3disp(a(k);end、运行结果为：a0= -45.3333a1= 94.2302a2= -6.1316、拟合曲线：(2)用 指数曲线拟合将等式两边取自然对数得：lny=lna+bx记 y=lny a=lna y=a+bx得到新数据：x11.522.533.544.555.566.577.58y3.514.384.815.075.245.375.475.535.595.645.695.715.745.765.78、C语言程序代码为：#include#include#includevoid main()void ColPivot(float*,int,float);int i;float x2;float c23=15,67.5,79.29, 67.5,373.5,373.5;ColPivot(c0,2,x);for(i=0;i=1;i+)printf(x%d=%fn,i,xi);void ColPivot(float *c,int n,float x)int i,j,k,t;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!=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)/(*(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)+n)/(*(c+i*(n+1)+i);、运行结果为：即：x0=4.208903x1=0.239355、.拟合曲线为：分析：带入数据进行检验，我们可以发现抛物线与指数曲线拟合相比，抛物线更加精确，而指数曲线具有较大的误差，因为在指数曲线拟合时，将指数函数按照一次函数来进行处理的，而拟合时，拟合函数的复杂度越高越准确，所以抛物线的拟合方法更准确。实习题6用复化梯形公式和复化辛卜生公式计算积分和积分，观察n为多少时，所得近似值是六位有效数字。解：（1）用复化梯形公式求积分、源程序如下：#include#includeint n;float AutoTrap(float(*f)(float),float a,float b)int i;float x,s,h=b-a;float t1,t2=h/2*(*f)(a)+(*f)(b);n=1;dos=0;t1=t2;for(i=0;i0.5e-6);return t2;float f(float x)return sqrt(1+cos(x)*cos(x);void main()float s;double pi=3.141592653;s=AutoTrap(f,0,pi/2);printf(T(%d)=%fn,n,s);、运行结果为：、结果分析由运行结果知：T8=1.910099该运行结果表示n为8时所得近似值具有6位有效数字。（2）用复化辛卜生公式计算、源程序：#include#includefloat Simpson(float(*f)(float),float a,float b,int n)int k;float s,s1,s2=0;float h=(b-a)/n;s1=(*f)(a+h/2);for(k=1;k=n-1;k+)s1+=(*f)(a+k*h+h/2);for(k=1;k=n-1;k+)s1+=(*f)(a+k*h+h/2);s2+=(*f)(a+k*h);s=h/6*(*f)(a)+4*s1+2*s2+(*f)(b);return s;float f(float x)if (x=0) return 1; else return sin(x)/x;void main()double pi=3.141592653;int i,n=2;float s;for(i=0;i=2;i+)s=Si