计算机数值方法实验报告

1、TAIWAN UNIVERSITV OF TECHNOLOGY本科实验报告课程名称:计算机数值方法实验项目:方程求根线性方程组的直接解法线性方程组的迭代解法代数插值和最小二乘法拟合多项式实验地点:逸夫302专业班级:软件学号:学生姓名:指导教师:田华2013 年 4 月 24 日精选文库学生姓名实验成绩实验名称实验一方程求根28实验目的和要求(必填)熟悉使用、迭代法、牛顿法、割线法等方法对给定的方程进行根的求解。选择上述方法中的两种方法求方程：二分法f(x)=x3+4x2-10=0在1,2内的一个实根，且要求满足精度|x*-xn|v0.5 W5实验内容和原理(必填)函数f(x)在区间(x，y)

2、上连续，先在区间(x，y)确定a与b，若f(a) ，f(b) 异号，说明在区间(a，b)内存在零点，然后求假设 F(a)0,ab 如果 f(a+b)/2=0 如果 f(a+b)/20f(a+b)/2，该点即为零点； ，则区间(a+b)/2，则区间(a,(a+b)/2)返回重新循环，不断接近零点。通过每次把,b)内存在零点，(a+b)/2 a;内存在零点，(a+b)/2 b;f(x)的零点所在区间收缩一半的方法，使区间内的两个端点逐步逼近函数零点，最终求得零点近似值。主要仪器设备台式或笔记本计算机实验记录(写出实验内容中的程序代码和运行结果)(可分栏或加页) 代码1.二分法:#in clude#

3、i nclude#in cludeint mai n()double a=1.0, b=2.0;double x,s;while(1) x=(a+b)/2;s=po w(x,3)+4*x*x-10;if (-0.000005 s & s 0.000005)break;else if(s 0)b=x;prin tf(%ft%fn,a,b);prin tf(%fn,x);prin tf(%fn,s);return 0;2.割线法：#i ncludestdio.h#i ncludemath.hint mai n()float c,a=1.0,b=2.0;while(1) c=b-(b*b*b+4*b

4、*b-10)*(b-a)/(b*b*b+4*b*b-(a*a*a+4*a*a); if(fabs(b-c)0.5*0.00001)break;b=c;prin tf(%fn,b);prin tf(%fn,c);流程图;运行结果；1二分法F：V1.2500001.2500001.3125001.3437501.3593751.5000001.3750001.3750001.3750001.3750001.3593751.3632811.363281 1.3642581.3647461.3649901.3651121.3651731.3652041.3652191.3652271.3652271.

5、3652291.3652301.3652300,0000011.367188 1.3671881.3652341.3652341.3652341.3652341.3652341.3652341.3652341.3652341.3652341.3652311.3652311.3652311 3872561.360731 1366160 1 36503S 1 365270 1.3652ZZ 1.355232 1 3652302割线法实验结果和分析两种方法均能求出方程的解，但割线法比二分法的收敛速度更快，且程序的代码更简洁。心得体会（遇到的问题和解决方法）通过实验，加深了对方程求根方法的理解，加强了

6、实践操作能力，实现了理论和实 践相结合。实验名称实验二线性方程组的直接求解实验目的和要求合理利用Gauss消元法、LU分解法、追赶法求解下列方程组:123X114012X28241X3130.3 10I 1559.1431X159.175.2916.13012X246.7811.2952X311211X424215X1287210X274836X371261120X4321X17121X25(n=5,10,100,121Xn 1512Xn5)实验内容高斯消元：l ik =aik /a kkaij = aij - l ik *akj (k=1,2,n-1i=k+1,k+2,n j=k+1,k+2

7、,，n+1 )由回代过程求得原方程组的解:Xn = a nn+l/ a nn X k = ( a kn+1-刀 akj X j )/ a kk追赶法:当矩阵A为三对角矩阵，在A的LU分解中，L取下三角阵，U取单位上三角阵，这样求解方程组Ax=d的方法称为追赶法.LU分解法:将系数矩阵A转化为A=L*U，L为单位下三角矩阵，U为普通上三角矩阵, 然后通过解方程组l*y=b，u*x=y,来求解主要仪器设备台式或笔记本计算机实验记录（写出实验内容中的程序代码和运行结果）（可分栏或加页）高斯消元法#in clude #defi ne n 3 mai n() float ann,cnn,bn,dn;fo

8、r(i=0;i n ;i+) for(j=0;j n;j+) sca nf(%f, &aij);cij=aij; sca nf(%f,&bi);di=bi; for(k=0;k n;k+) bk=dk/ckk;for(i=0;i n ;i+) if(i=k) con ti nue;cik=cik/ckk;for(j=k+1;j vn ;j+) akj=ckj/ckk;aij=cij-cik*ckj; bi=di-cik*dk; for(i=0;i vn ;i+) di=bi;for(j=k+1;j vn ;j+) cij=aij; for(i=0;i vn ;i+) prin tf(b%d=%

9、fn,i,bi);LU分解法:#in clude vstdio.h#in clude vmath.h#defi ne L 30double aLL,bL,lLL,uLL,xL,yL; int mai n()intn ,i,j,k,r;printf(”请输入矩阵元次：n); sca nf(%d,&n);printf(”请输入矩阵各项：n”)； 血0=1;)=n;+i)for(j=1;j=n;+j)scan f(%lf,&aij);printf（”请输入方程组的常数项：n）;for(i=1;i=n;+i)scan f(%lf, &bi);for(i=1;i=n;+i)for(j=1;j=n;+j)

10、lij=0;uij=O.O;for(k=1;k=n;+k)for(j=k;j=n;+j)ukj=akj; for(r=1;rk;+r)ukj-=lkr*urj; for(i=k+1;i=n ;+i)lik=aik; for(r=1;rk;+r)lik-=lir*urk;lik/= ukk; lkk=1.0;for(i=1;i=n;+i)yi = bi; for(j=1;j0;-i)xi = yi;for(j=i+1;j=n ;+j) xi-=uij*xj; xi/= uii;for(i=1;i=n;+i)prin tf(%0.2lfn,xi); return 0;追赶法#in elude st

11、dio.h #defi ne n 5main ()float an ,b n,c n-1,d n,t;int i;scan f(%f%f%f,&b0, &c0, &d0);for(i=1;i n-1;i+)scan f(%f%f%f%f,&ai, &bi,&ci, &di);scan f(%f%f%f,&an-1,&b n-1，& d n-1); c0=c0/b0;d0=d0/b0;for(i=1;i=0;i-) di=di-ci*di+1;for(i=0;i n;i+)prin tf(d%d=%fn,i,di);实验结果和分析1高斯消元法2LU分解3追赶法实验分析:高斯消元法，是先消元，再回

12、带的过程。由程序段可以发现，始终消去 对角线下方的元素。从消元过程可以看出，对于n阶线性方程组，只要各 步主元素不为零，经过n-1步消元，就可以得到一个等价的系数矩阵为上三角形阵的方程组，然后再利用回代过程可求得原方程组的解。LU分解法，分解矩阵为单位下三角阵L与上三角阵U的乘积，然后解方程组Ly=b，回代，解方程组Ux=y。其中的L为n阶单位下三角阵、U 为上三角阵.对于追赶法，追赶法是适用于三角矩阵的线性方程组的求解的方法，并 不适用于其他类型矩阵。心得体会（遇到的问题和解决方法）本次实验难度比较大，在编译时经常出现各种错误，程序代码也比较繁琐，深深感 觉到自己的上机操作能力有限，应加强自

13、己的编程能力，以后要继续努力。实验名称实验三线性方程组的迭代求解实验目的和要求使用雅可比迭代法或高斯-赛德尔迭代法对下列方程组进行求解。10期一 込一2也 =7.2* -尤1 +10乞一2壬-E_3X -工一5主二4一2匚*实验内容设线性方程组Ax=b的系数矩阵A可逆，且主对角元素a11,a22,ann均不为零，令D=diag( a11,a 22, * ,a nn)并将A分解成A=(A-D)+D从而线性方程组可写成Dx=(D-A)x+b则有迭代公式x(k+1)=B1x(k)+f1其中，B1=I-D-1A,f1=D-1b。主要仪器设备台式或笔记本计算机 实验记录（写出实验内容中的程序代码和运行结

14、果）（可分栏或加页）雅可比迭代法#in clude #in clude #defi ne n 3 void f(float *b,float x)float yn+1=0,0,0,1;dok=0;for(i=0;i vn+1;i+) xi=yi； for(i=0;i vn ;i+) yi=0;for(j=0;j vn+1;j+) yi+=*(b+( n+1)*i+j)*xj; for(i=0;i vn ;i+) if(fabs(yi-xi)v=0.5e-3) k+;if(k=3) break;while(1);for(i=0;i vn ;i+) prin tf(y%d=%fn,i,yi); m

15、ai n() float bn n+1=0,0.1,0.2,0.72,0.1,0,0.2,0.83,0.2,0.2,0,0.84;float xn+1=0,0,0,1;f(b0,x);高斯一赛德尔迭代法#in elude iostream#in elude ioma nipusing n ames pace std;int mai n()int i,j,k=O,m ,n;double t1,t2,e1,e2=0.0;coutvv请输入精度e：;cin e1;coutvv请输入系数矩阵行数: cinm;coutvv请输入系数矩阵列数：;cinn;coutvve ndl;double (*a)=n

16、ew double *m;/ 生成二维动态数组for(i=0;iv=m;i+)ai二new double n;double (*b)二new double m;double (*x)=new double n;cout请输入系数矩阵：endl;coute ndl;for(i nt nu m1二 0;nu m1 m;nu m1+)for(i nt nu m2二 0;num2 a nu m1 nu m2;coute ndl;cout输入的系数矩阵为：endl;for (int nu m3二 0;nu m3 m;nu m3+)for(i nt num4=0;num4n;num 4+)couta nu

17、 m3 num4coute ndl;coute ndl;cout请输入矩阵 b: endl;coute ndl;for(i nt nu m5二 O;nu m5bnu m5;cout输入的矩阵b为：endl;for(i nt nu m6二 O;nu m6 m;nu m6+)coutb nu m6Hcoute ndl;coute ndl;for(i nt nu m7二 0;num7n;nu m7+)x num7=0.0000;docout第k次迭代值：;e2=0.0;for(i=0;im;i+)double sum=0.0;for(j=0;j=0?(xi)-t1:t1-(xi);He2=(e2=t

18、2?e2:t2);coutset precisi on( 8)xicout=e1 &k30);coutvv共迭代了 vvkvv次;deletea;deleteb;deletex;return 0 ;实验结果和分析1雅克比迭代2高斯一一赛德尔迭代实验分析：使用这两种方法都可以求出方程的解，高斯赛德尔迭代法所需的迭代次数比雅克比迭代少，能够更早的达到精度要求。但是雅克比的时效性要比高斯赛德尔的好。心得体会（遇到的问题和解决方法）本次试验，让我对这两种方法更加理解，在编程操作上也更加熟练，今后继续努力, 不断丰富自己的知识，增强操作能力。实验名称实验四代数插值和最小二乘法拟合实验目的和要求1使用拉格

19、朗日插值法或牛顿插值法求解：已知f(x)在6个点的函数值如下表所示，运用插值方法，求f(0.596)的近似值。X0.400.550.650.800.901.05f(x)0.410750.578150.696750.888111.026521.253862给定数据点(Xi ,yi),用最小二乘法拟合数据的多项式，并求平方误差。Xi00.50.60.70.80.91.0yi11.751.962.192.442.713.00实验内容1设函数在区间a,b上n+1互异节点Xo, Xi,Xn上的函数值分别为yo, yi,，yn,求n次插值多项式Pn(x),满足条件Pn(Xj)=yj, j=0,1,，nLn

20、(x)=y ol o(x)+y il i(x)+ +ynln(x)= E yili(x)其中 l o(x),l1(X),l n(x)为以 X0,X 1,值基函数，,x n为节点的n次插则Ln(x)是一次数不超过n的多项式，且满足Ln(Xj)=yj, L=0,1,，n再由插值多项式的唯一性，得R(x)三 Ln(x)2 建立正规方程组:j+kjE(E Xi ) ak=EXi yi , j=0,1,n平方误差:Ik2I =E(E akXi -y i)对给定数据点(Xi，Yi)(i=0,1，m)，在取定的函数类 中，求p(x) ，使误差的平方和E2最小，E2=E p(Xi)-Yi2。从几何意义上讲，就

21、是寻求与给定点 (Xi , Yi)(i=0,1 ,，m)的距离平方 和为最小的曲线y=p(x)。函数p(x)称为拟合函数或最小二乘解，求拟合函数p(x) 的方法称为曲线拟合的最小二乘法。得到的两个关于aO、a1为未知数的两个方程组，解这两个方程组得出:aO =(刀 Yi) / m - a1 (刀 Xi) / ma1 = m 刀 Xi Yi(刀 Xi 刀 Yi) / m 刀Xi2(刀 Xi)2 )即最终的拟合多项式各项系数主要仪器设备台式或笔记本计算机实验记录(写出实验内容中的程序代码和运行结果)(可分栏或加页)1#in elude #in elude #in elude #in elude v

22、oid differe nce(float *x,float *y,i nt n)float *f;int k,i;f=(float *)malloe( n*sizeof(float);for(k=1;k=n ;k+)fO=yk；for(i=O;ik;i+) fi+1=(fi-yi)/(xk-xi);yk=fk;return;int mai n()int i,n;float x20,y20,xx,yy;printf(请输入数据个数n: ”);scan f(%d,&n);prin tf(n);for(i=0;i=0;i-)yy=yy*(xx-xi)+yi;prin tf(n近似值为：(%f)=%

23、fn,xx,yy);2#in clude#in clude#defi ne N 15double po wer(double &a,i nt n)double b=1;for(i nt i=0;i n;i+)b*=a;return b;void Gauss();double XN,YN,sumXN,sumYN,aNN,bN,lNN,xN;int mai n()ofstream outdata;ifstream in data;double s;int i,j,k, n,in dex;coutn;coute ndl;cout请输入 X 和 Y:endl; for(i=0;i n;i+)/输入给定

24、数据coutXiXi; sumX1+=Xi;coutYiYi; sumY1+=Yi; coute ndl; coutsumX1=sumX1tsumY1=sumY1e ndl; coutin dex;coute ndl;i=n;sumX0=i;for(i=2;i=2* in dex;i+)sumXi=0;for(j=0;j n;j+)sumXi+=p ower(Xj,i);coutsumXi=sumXie ndl;for(i=2;i=in dex+1;i+)sumYi=0;for(j=0;j n;j+)sumYi+=po wer(Xj,i-1)*Yj;coutsumYi=sumYie ndl;for(i=1;i=in dex+1;i+)/建立正规方程组for(j=1;j=in dex+1;j+) aij=sumXi+j-2; bi=sumYi;k=1;/用高斯消元法解方程组dofor(j=k+1;j=in dex+1;j+) ljk=ajk/akk;for(i=k+1;i=in dex+1;i+) for(j=k+1;j=1;i-)S=0；for(j=i+1；j=in dex+1;j+) s=s+aij*xj;xi=(bi-s)/