数值实验报告.docx

1、数值实验报告本科实验报告课程名称： 计算机数值方法 实验项目： 方程求根、线性方程组的直接求解、线性方程组的迭代求解、最小二乘法拟合和代数插值 实验地点： 专业班级： 学号： 学生姓名： 指导教师： 2014年 5月 日学生姓名实验成绩实验名称 实验一 方程求根实验目的和要求（必填）（1）了解非线性方程求根的常见方法，如二分法、牛顿法、割线法。（2）加深对方程求根方法的认识，掌握算法。（3）会进行误差分析，并能对不同方法进行比较实验内容和原理（必填）方程求根：熟悉使用二分法、牛顿法、割线法等方法对给定的方程进行根的求解。选择上述方法中的两种方法求方程：f(x)=x3+4x2-10=0在1,2内

2、的一个实根，且要求满足精度|x*-xn|0.510-5二分法原理：f(x)在区间（x，y）上连续先找到a、b属于区间（x，y），使f(a)，f(b)异号，说明在区间(a,b)内一定有零点，然后求f(a+b)/2, 现在假设f(a)0,ab1 如果f(a+b)/2=0，该点就是零点，如果f(a+b)/2a，从开始继续使用割线法原理：主要仪器设备惠普 ProBook 6470b实验记录(写出实验内容中的程序代码和运行结果)(可分栏或加页)二分法：#include iostreamusing namespace std;#includestdio.h#include math.hint i;doub

3、le a20,b20;double x20;double A,y;int main() a0=1.0,b0=2.0;for(i=0;i11;i+)printf( ai bi xiy n); xi=(bi+ai)/2; y=xi*xi*xi+4*xi*xi-10; A=ai*ai*ai+4*ai*ai-10; if(y*A)0) ai+1=ai,bi+1=xi; else ai+1=xi,bi+1=bi; printf(%f, %f, %f, %fn,ai,bi,xi,y); 割线法#includestdio.h#includemath.h#includeiostreamusing namesp

4、ace std;float main() float c,a=1.0,b=2.0; while(1) c=b-(b*b*b+4*b*b-10)*(b-a)/(b*b*b+4*b*b-(a*a*a+4*a*a); if(fabs(b-c)0.5*0.000001) break; b=c; coutc;实验结果和分析心得体会（遇到的问题和解决方法）使用不同的方法，可以不同程度的求得方程的解，不同的方法速度不同，求得的结果也稍有区别，当然和要求精度也有关系。刚开始的时候用数组对二分法进行求解，发现循环到第二次就无法实现值的传递，于是换了另外一种方法代替了数组。实验名称 实验二 线性方程组的直接求解实

5、验目的和要求（1）了解线性方程组常见的直接解法，如Guass消元法、LU分解法、追赶法。（2）加深对线性方程组求解方法的认识，掌握算法。（3）会进行误差分析，并能对不同方法进行比较。实验内容合理选择利用Gauss消元法、LU分解法、追赶法求解下列方程高斯分解法：将原方程组化为三角形方阵的方程组：lik=aik/akkaij=aij-lik*akjk=1,2,n-1i=k+1,k+2,nj=k+1,k+2,n+1由回代过程求得原方程组的解：xn=ann+1/annxk=(akn+1-akjxj)/akk(k=n-1,n-2,2,1)LU分解法：将系数矩阵A转化为A=L*U,L为单位下三角矩阵，U

6、为普通上三角矩阵，然后通过解方程组l*y=b,u*x=y,来求解x.追赶法:用来求对角方程组：将系数矩阵A转化为A=L*U,L为普通下n-1对角矩阵，U为单位上n-1对角矩阵，然后通过解方程组l*y=b,u*x=y,来求解x.主要仪器设备惠普 ProBook 6470b实验记录(写出实验内容中的程序代码和运行结果)(可分栏或加页)高斯消元法：#includestdio.h#includeiostreamusingnamespacestd;floatmain()floata34=1,2,3,14,0,1,2,8,2,4,1,13;floatx3;floatsum=0;intk,i,j;for(k

7、=0;k2;k+) for(i=k+1;i3;i+)for(j=k+1;j4;j+)aij=aij-aik/akk*akj;for(i=0;i3;i+)for(j=0;j=0;k-)sum=0;for(j=k+1;j3;j+) sum+=akj*xj;xk=(ak3-sum)/akk;for(i=0;i3;i+)printf(x%d=%f,i+1,xi);LU分解法:#include#include#defineL30doubleaLL,bL,lLL,uLL,xL,yL;intmain()intn,i,j,k,r;scanf(%d,&n); for(i=1;i=n;+i)for(j=1;j=n

8、;+j)scanf(%lf,&aij);for(i=1;i=n;+i)scanf(%lf,&bi);for(i=1;i=n;+i)for(j=1;j=n;+j)lij=0;uij=0.0;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*

9、xj;xi/=uii;for(i=1;i=n;+i)printf(%0.2lfn,xi);return 0追赶法：#include stdio.hmain() FILE *f; double a15,b15,c15,d15; double t; int i,n; /*/ f=fopen(zgf.txt,r); fscanf(f,%d,&n); fscanf(f,%lf%lf%lf,&b1,&c1,&d1); for(i=2;i=n-1;i+) fscanf(f,%lf%lf%lf%lf,&ai,&bi,&ci,&di); fscanf(f,%lf%lf%lf,&an,&bn,&dn); fcl

10、ose(f); /*/ c1=c1/b1; d1=d1/b1; for(i=2;i=1;i-) di=di-ci*di+1; printf(n*n); for(i=1;i=n;i+) printf(d%2d=%lfn,i,di);实验结果和分析高斯消元法结果：LU分解法结果：追赶法结果：心得体会（遇到的问题和解决方法）从消元过程可以看出，对于n阶线性方程组，只要各步主元素不为零，经过n-1步消元，就可以得到一个等价的系数矩阵为上三角形阵的方程组，然后再利用回代过程可求得原方程组的解.由于列主元素法相似且优于完全主元素法所以省略了后者。消元过程相当于分解A为单位下三角阵L与上三角阵U的乘积，解方

11、程组Ly=b回代过程就是解方程组Ux=y。其中的L为n阶单位下三角阵、U为上三角阵.在A的LU分解中,L取下三角阵,U取单位上三角阵,这样求解方程组Ax=d的方法称为追赶法。另外是追赶法和其他方法求同一方程结果不一样，我多次修改源程序，也不知道原因。再就是追赶法有很大的局限性还待改良实验名称 实验三 线性方程组的迭代求解实验目的和要求(1)了解线性方程组的迭代求解方法，雅可比迭代法或高斯-赛德尔迭代法等典型方法。(2)帮助学生全面消化已学的相关课程内容，深刻理解计算数值方法课程的内涵，培养使用电子计算机进行科学计算和解决问题的能力。(3)进行基本技能训练和巩固。使学生得到选择算法、编写程序、分

12、析数值结果、写数值试验报告、课堂讨论等环节的综合训练实验内容使用雅可比迭代法或高斯-赛德尔迭代法对下列方程组进行求解。 实验原理：把矩阵A分解矩阵N和P，其中N为非奇异矩阵，于是，Nx=Px+b即x=Bx+f据此，可以建立迭代公式x(k+1)=Bx（k）+1;若序列X（k）收敛,limx(k)=x*显然有x*=Bx*+f;即,极限x*便是所求方程组的解。主要仪器设备惠普 ProBook 6470b实验记录(写出实验内容中的程序代码和运行结果)(可分栏或加页)#include#includemain()inti;doublex120,x220,x320;doublex10,x20,x30;pri

13、ntf(请输入x1,x2,x3的初值：n);scanf(%lf%lf%lf,&x10,&x20,&x30);printf(nx1nx2nx3nn);for(i=0;i18;i+) x10=x10;x2 0=x20;x30=x30;x1i+1=0.1*x2i+0.2*x3i+0.72;x2i+1=0.1*x1i+0.2*x3i+0.83; x3i+1=0.2*x1i+0.2*x2i+0.84; printf(%5d%5lf%5lf%5lfn,i,x1i,x2i,x3i); 实验结果和分析心得体会（遇到的问题和解决方法）使用高斯-赛德尔和雅克比迭代都可以求出方程组的解，但是利用高斯-赛德尔迭代法所

14、需的迭代次数比雅克比迭代少，能够更早的达到精度要求。产生误差的原因，舍入误差及数的存储空间是有限的，因此出现误差。雅可比迭代法显然其速度不如高斯-赛德尔迭代法好。这是由其算法的设计本身来决定的。同时高斯-赛德尔算法由于都涉及了两次循环且同时未引入新的变量，因此相比而言高斯赛德尔迭代法较优。实验名称 实验四 代数插值和最小二乘法拟合实验目的和要求（1）了解矩阵特征值与特征向量问题解法，掌握幂法。（2）加深对矩阵特征值与特征向量问题求解方法的认识，掌握算法。（3）熟练运用已学计算方法求解方程组（4）加深对计算方法技巧，选择正确的计算方法来求解各种方程组（5）培养使用电子计算机进行科学计算和解决问题

15、的能力实验内容代数插值：用拉格朗日插值法或牛顿插值法求解：已知f(x)在6个点的函数值如下表所示，运用插值方法，求f(0.596)的近似值。x0.400.550.650.800.901.05f(x)0.410750.578150.696750.888111.026521.25386最小二乘法拟合多项式：给定数据点（xi ,yi），用最小二乘法拟合数据的多项式，并求平方误差。xi00.50.60.70.80.91.0yi11.751.962.192.442.713.00原理：设函数在区间a,b上n+1互异节点x0,x1,xn上的函数值分别为y0,y1,yn,求n次插值多项式Pn(x),满足条件P

16、n(xj)=yj,j=0,1,n令Ln(x)=y0l0(x)+y1l1(x)+ynln(x)=yili(x)其中l0(x),l1(x),ln(x)为以x0,x1,xn为节点的n次插值基函数，则Ln(x)是一次数不超过n的多项式，且满足Ln(xj)=yj,L=0,1,n再由插值多项式的唯一性，得Pn(x)Ln(x)最小二乘法原理：主要仪器设备惠普 ProBook 6470b实验记录(写出实验内容中的程序代码和运行结果)(可分栏或加页)#include#includemain() intI;charL;doubleM100100;doublex100,y100;doubleX=1,xx=0,w=1

17、,N=0,P,R=1;intn=5;/coutn;/for(inti=0;i=n;i+)/*cout请输入xi的值：xi;cout请输入yi的值：yi;Mi0=xi;Mi1=yi;*/ M00=0.40;M01=0.41075;M10=0.55;M11=0.57815;M20=0.65;M21=0.69675;M30=0.80;M31=0.88811;M40=0.90;M41=1.02652;M50=1.05;M51=1.25386;for(intj=2;j=n+1;j+) for(inti=1;i=n;i+)Mij=(Mij-1-Mi-1j-1)/(Mi0-Mi-j+10);for(inti

18、=1;i=n;i+)cout其i阶均差为：Mii+1endl;coutxx;for(inti=0;in;i+)X*=xx-Mi0;N+=Mi+1i+2*X;P=M01+N;cout其函数值：y=Pendl;/#include#include#defineN15doublepower(double&a,intn)doubleb=1;for(inti=0;in;i+)b*=a;returnb;voidGauss();doubleXN,YN,sumXN,sumYN,aNN,bN,lNN,xN;voidmain()doubles;inti,j,k,n,index;/coutn;n=7;coutendl

19、;/cout请输入X和Y:endl;/输入给定数据X0=0.0;Y0=1.00;X1=0.5;Y1=1.75; X2=0.6;Y2=1.96;X3=0.7;Y3=2.19;X4=0.8;Y4=2.44;X5=0.9;Y5=2.71;X6=1.0;Y6=3.00;/绑定数据/可以解绑由下for循环输入任何数据for(i=0;in;i+) sumX1+=Xi;sumY1+=YicoutsumX1=sumX1tsumY1=sumY1endl;coutindex;coutendl;i=n;sumX0=i;for(i=2;i=2*index;i+)sumXi=0;for(j=0;jn;j+)sumXi+

20、=power(Xj,i);coutsumXi=sumXiendl;for(i=2;i=index+1;i+)sumYi=0;for(j=0;jn;j+)sumYi+=power(Xj,i-1)*Yj;coutsumYi=sumYiendl;for(i=1;i=index+1;i+)/建立正规方程组for(j=1;j=index+1;j+)aij=sumXi+j-2;bi=sumYi;k=1;/用高斯消元法解方程组 dofor(j=k+1;j=index+1;j+)ljk=ajk/akk;for(i=k+1;i=index+1;i+)for(j=k+1;j=index+1;j+)aij=aij-lik*akj; bi=bi-lik*bk;if(k=index+1)break;k+;while(1);xindex+1=bindex+1/aindex+1ind