数值分析实验报告

1、 数值分析实验报告计1310 陈万全学号：132852实验一一实验内容：用两种不同的顺序计算，分析其误差的变化。二实验目的：1. 通过上机编程，复习巩固以前所学程序设计语言；2. 通过上机计算，了解舍入误差所引起的数值不稳定性。3. 通过上机计算，了解运算次序对计算结果的影响，从而尽量避免大数吃小数的现象。三算法描述 取float类型数据从大到和从小到大的累加。然后输入结果四源程序：/* 实验一 */#include using namespace std;float shunxi()float sum=0;for(int i=1;i0;i-)sum=sum+1.0/float(i*i);re

2、turn sum;int main()cout.precision(10);cout第一个顺序方法的值：shunxi()endl;cout第二个倒序方法的值：daoxu()endl;return 0;五、程序结果六实验体会 这个程序还是比较简单的，仅仅需要一个循环就可以完成，只要细心就不会有错。从这个实验里我也看到了计算机毕竟不是万能的，由于机器字长的问题对于一些特殊的计算会出错即大数吃小数这样的计算，这样会让人们在日常的生产生活中产生一些意想不到的问题。不过只要稍微对计算方法做一下改动就可以改正这种错误改变算术的顺序。至于从哪个方向计算更好，那就要具体问题具体分析了.对于程序上的解决了，co

3、ut输出的小数保留的位数。实验二一实验内容：1、拉格朗日插值按下列数据x-3.0-1.01.02.03.0y1.01.52.02.01.0作二次插值，并求x=-2，x=0，x=2.75时的函数近似值2牛顿插值按下列数据x0.300.420.500.580.660.72y1.044031.084621.118031.156031.198171.23223作五次插值，并求x=0.46，x=0.55，x=0.60时的函数近似值.二实验目的： 1. 掌握牛顿插值法。2. 掌握拉格朗日插值法。三算法描述开始时判断输入的点是否在区间内，不在则退出在则继续根据值选择节点 带入公式计算四 源程序：朗格朗日二次

4、插值#include using namespace std;double inser(double x)double sum1,sum2,sum3,sum;double xi5=-3.0,-1.0,1.0,2.0,3.0;double yi5=1.0,1.5,2.0,2.0,1.0;int i;int a,b,c;for (i=0;ixii&xxii+1|xxi4)a=i;b=(i+1)%5;c=(i+2)%5;break;sum1=yia*(x-xib)*(x-xic)/(xia-xib)*(xia-xic);sum2=yib*(x-xia)*(x-xic)/(xib-xia)*(xib-

5、xic);sum3=yic*(x-xib)*(x-xia)/(xic-xib)*(xic-xia);sum=sum1+sum2+sum3;return sum;int main()cout朗格朗日二次插值-2.0的值：inser(-2.0)endl;cout朗格朗日二次插值0.0的值：inser(0.0)endl;cout朗格朗日二次插值2.75的值：inser(2.75)endl;return 0;牛顿插值#include using namespace std;#define N 6double fiN;int chazhi();double niudun(double x)double

6、xiN=0.3,0.42,0.5,0.58,0.66,0.72;double yiN=1.04403,1.08462,1.11803,1.5603,1.19817,1.23223;double sum=0,sumxN;chazhi();sumx0=1;for (int j=1;jN;j+)sumxj=(x-xij-1)*sumxj-1;for (int i=0;i0;j-)for (i=0;ij;i+)yii=double(yii+1-yii)/double(xiflag+i-xii);:fiflag=yi0;flag+;return 0;int main()cout牛顿插值插入0.46，可得

7、值：niudun(0.46)endl; cout牛顿插值插入0.55，可得值：niudun(0.55)endl;cout牛顿插值插入0.60，可得值：niudun(0.60)endl;return 0;五、程序结果拉格朗日插值牛顿插值五 实验体会这个实验的目的主要是掌握如何在计算机上运用拉格朗日插值法和牛顿插值法来解决问题。对于拉格朗日插值，已经给出了五个节点的函数值，可以在一定的范围内使用拉格朗日插值法，所以对于-5到1上的X我利用前三个节点，而对于1到5上X我利用后三个节点，这样算出来值的误差会比较小。在牛顿插值上，题中已经给出了六个节点的函数值，利用这六个节点可以算出五阶差商，然后很容易

8、就能够得出结果。其中四阶差商为负数，但是差别不大，五阶差商又变回了正数，所以可以用到五阶差商，这样结果也会更准确。实验三一实验内容分别用复化梯形公式和复化辛卜生公式计算f(x)=sin(x)/x的积分，并与准确值比较判断精度。二实验目的1掌握复化梯形公式。2掌握复化辛卜生公式。三算法描述1.复化梯形公式代数形式求解定积分，为每个区间内的梯形公式的累加，所以除端点外的点都累加两次2复化新普生公式与梯形公式形式上差别不大，但是由于使用的是三个点一个区间所以只可以求解Sn4所以流程类似。四源程序/*复化梯形公式*#include using namespace std;#include #defin

9、e A 0#define B 1#define N 8double Yi(double xk)double s;if(xk=0)return 0.0;elses=double(sin(xk)/double(xk);return s;double fx() double sumx=0;double ki;for (int k=0;kB)cout输入的A小于B，错误endl;return 1;double sumx=fx(),sumy,sum;sumy=sumx-Yi(A)+Yi(B);sum=(double(B-A)/double(N)*(sumx+sumy)*0.5;coutsinx/x在A到

10、B区间的N段分值为:sumendl;return 0; / */*#include using namespace std;#include #define A 0#define B 1#define N 100double Yi(double xk)double s;if(xk=0)return 0.0;elses=double(sin(xk)/double(xk);return s;double fx() double sumx=0;double ki;for (int k=0;kN;k+)ki=(double(B-A)/double(N)*k+A+(double(B-A)/double(

11、2*N);sumx+=Yi(ki);return sumx;double fy() double sumy=0;double ki;for (int k=0;kB)cout输入的A小于B，错误endl;return 1; double sumx=fx(),sumy=fy(),sum; double h=abs(A-B); sum=(Yi(A)+4*sumx+2*sumy+Yi(B)*h/double(6*N); coutsinx/x在A到B区间的N段分值为:sumendl; return 0;*/复化辛卜生公式六实验体会要求函数sin(x)/x的积分，程序是通过两种不同的复化求积，通过计算机计

12、算，检验，更加理解了复化梯形公式和复化辛卜生公式。如果要想求其它区间的积分可以直接输入。复化梯形和复化辛卜生公式都是将区间分解成为若干等份，然后依次求积，编程应注意算法的严密性。实验四一实验内容 用改进欧拉方法解初值问题y=x+y; y(0)=1。0x1,取步长h=0.1计算，并与准确值 y=-x-1-2ex相比较。（p141习题第二题2）二实验目的1 了解掌握欧拉方法。三算法描述按照欧拉格式进行带入迭代计算流程图如下四源程序/用改进欧拉方法解初值问题y=x+y; y(0)=1。0x1,取步长h=0.1计算，/并与准确值 y=-x-1-2ex相比较。（p141习题第二题2）#include #

13、include #include using namespace std;#define A 0#define B 1double h;double implicit (int i,double yi)double yp;yp=yi+h*(double (i*h+A)+yi);return yp;double expli (int i,double yi,double yp)double yc;yc=yi+h*(double(i+1)*h+A)+yp);return yc;int main ()int n;cout请输入步长h。如：0.1h;n=int(double(B-A)/double(h

14、);double *y=new doublen+1,*yl=new doublen+1;cout请输入y(0)为：y0;yl0=y0;double ypp,ycc;int i;for (i=0;in;i+)ypp=implicit(i,yli);ycc=expli(i,yli,ypp);yli+1=(ypp+ycc)/double(2); /计算改进欧拉值for (i=0;i=n;i+)yi=-(i*h+A)+1-2*exp(i*h+A); /计算真实值coutyln-1 ylnendl;cout改进欧拉公式 准确值 误差endl;for (i=0;i=n;i+)cout.precision(

15、7);coutyli yi yli-yiendl;delete y;delete yl;return 0;六 实验体会课本上给出了明确的过程了，所以我就按照那个写了这个程序。我认为实验目的主要还是通过编程更好地理解改进欧拉的方法，课本上给出了它的两种形式，第一种是用于做题的，而第二个就是这个程序上的。熟悉了这两种方法就能够很好地掌握改进欧拉方法了。实验五一 实验内容分别用下列方法求f(x)=x3-3x-1=0在x0=2附近的根。根的准确值为x*=1.87938524，要求准确到四位有效数字，并对比各种算法的计算量。二分法；（2）简单迭代法；（3）牛顿迭代法二 实验目的1 了解二分法2 熟悉简单

16、迭代法3 掌握牛顿迭代法三 算法描述 这三个都是迭代算法，以上一部的值带入下一步求出近似解。四 源程序#include #include using namespace std;#define A 1#define B 3double ErFenFa()double b=B,a=A,ab=(b+a)/double(2);double yb=b*b*b-double(3)*b-double(1), ya=a*a*a-double(3)*a-double(1), yab;if(yb*ya=0)cout输入A B 区间有错0.0005)ab=double(b+a)/double(2);yab=ab*

17、ab*ab-double(3)*ab-double(1);if(yab0)b=ab;if (yab0)a=ab;if (yab=0)break;cout.precision(7);cout由二分法计算可得x=ab 其误差为:ab-1.879385240.0005|N20)x=xk;xk=pow(double(3)*x+double(1),1.0/3);cout.precision(7);cout由简单迭代法计算可得x=xk 其误差为:xk-1.879385240.0005|N20)x=xk;xk=x-(x*x*x-double(3)*x-double(1)/(double(3)*x*x-dou

18、ble(3);cout.precision(7);cout由牛顿法计算可得x=xk 其误差为:xk-1.87938524endl;return xk;int main()double x1,x2,x3;x1=ErFenFa();x2=DieDaiFa();x3=NewDunfa();return 0;五 实验结果六实验体会二分法实现比较简单，很容易写出算法来，简单迭代法按公式进行迭代，但好几次没弄清楚变量，由实验结果也可看出牛顿法的迭代次数少、精确度高。实验六一 实验内容：分别用高斯列主元消去法和直接三角分解法（LU分解）求方程组的解系数矩阵：10 7 8 7 常向量：10 7 5 6 5 8

19、 8 6 10 9 6 7 5 9 10 7精确解为：（-60，102，-27，16）二 实验目的：1.了解高斯列主元消去法的算法步骤。2.掌握直接三角分解法。三 算法描述列主元消去每次选择一个最大列元素然后交换到第一行进行消元计算然后队进行消元后的结果回带得到方程的解。四 源程序：/高斯法#define t(i,j) (i)-1)*n+(j)-1#include #include int main(void)int n,i,j,k; /n方程个数 i，j，k控制循环变量double *matrix; /方程系数矩阵double *c; /方程常向量double *sol; /方程解向量dou

20、ble t; /临时变量/输入nprintf(请输入方程个数:n);scanf(%d,&n);/申请内存matrix = (double*)malloc(sizeof(double) * n * n);c = (double*)malloc(sizeof(double) * n);sol = (double*)malloc(sizeof(double) * n);/输入数据for(i = 1; i = n; i+)printf(请输入第%d行系数:n,i);for(j = 1; j =n; j+)scanf(%lf,&matrixt(i,j);printf(请输入常向量:n,i);for(i

21、= 1; i = n; i+)scanf(%lf,&ct(1,i);/求梯形矩阵for(k = 1; k = n; k+)/akk项等于1t = matrixt(k,k);for(j = k; j = n; j+)matrixt(k,j) = matrixt(k,j) / t;ct(1,k) = ct(1,k) / t;/aik项等于0for(i = 1 + k; i = n; i+)t = -matrixt(i,k);for(j = k; j =1 ; i-) solt(1,i) = ct(1,i);for(j = i + 1; j =n; j+)solt(1,i) = solt(1,i) - matrixt(i,j)*solt(1,j);/打印梯形矩阵printf(梯形矩阵为n);for(i = 1; i = n; i+)for(j = 1; j =n; j+)printf(%5g ,matrixt(i,j);printf(%10gn,ct(1,i);/打印方程的解printf(nn线性方程组解为n);for(i = 1; i = n; i+)printf(%5g,solt(1,i);printf(