计算机数值方法实验报告

精选文库本科实验报告课程名称： 计算机数值方法 实验项目：方程求根 线性方程组的直接解法 线性方程组的迭代解法代数插值和最小二乘法拟合多项式 实验地点： 逸夫302 专业班级： 软件 学号： 学生姓名： 指导教师：田华 2013年 4 月 24 日学生姓名实验成绩实验名称 实验一 方程求根实验目的和要求（必填）熟悉使用、迭代法、牛顿法、割线法等方法对给定的方程进行根的求解。选择上述方法中的两种方法求方程：二分法f(x)=x3+4x2-10=0在1,2内的一个实根，且要求满足精度|x*-xn|0.510-5实验内容和原理（必填）函数f(x)在区间（x，y）上连续，先在区间（x，y）确定a与b，若f(a)，f(b)异号，说明在区间(a，b)内存在零点，然后求f(a+b)/2。假设F(a)0,ab， 如果f(a+b)/2=0，该点即为零点； 如果f(a+b)/20，则区间(a,(a+b)/2)内存在零点，(a+b)/2b；返回重新循环，不断接近零点。通过每次把f(x)的零点所在区间收缩一半的方法，使区间内的两个端点逐步逼近函数零点，最终求得零点近似值。主要仪器设备台式或笔记本计算机实验记录(写出实验内容中的程序代码和运行结果)(可分栏或加页)代码 1. 二分法：#include#include#includeint main() double a=1.0, b=2.0; double x,s; while(1) x=(a+b)/2; s=pow(x,3)+4*x*x-10; if (-0.000005 s & s 0.000005) break; else if(s 0) b=x; printf(%ft%fn,a,b); printf(%fn,x); printf(%fn,s); return 0;2. 割线法：#includestdio.h#includemath.hint 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.00001)break; b=c; printf(%fn,b); printf(%fn,c); 流程图；运行结果；1二分法2 割线法 实验结果和分析 两种方法均能求出方程的解，但割线法比二分法的收敛速度更快，且程序的代码更简洁。心得体会（遇到的问题和解决方法）通过实验，加深了对方程求根方法的理解，加强了实践操作能力，实现了理论和实践相结合。实验名称 实验二 线性方程组的直接求解实验目的和要求合理利用Gauss消元法、LU分解法、追赶法求解下列方程组： （n=5,10,100,）实验内容高斯消元： lik=aik/akk aij= aij- lik* akj （ k=1,2,n-1 i=k+1,k+2, ,n j=k+1,k+2, ,n+1 ）由回代过程求得原方程组的解： xn= ann+1/ ann xk=( akn+1-akj xj)/ akk追赶法： 当矩阵A为三对角矩阵，在 A 的LU 分解中, L取下三角阵, U 取单位上三角阵,这样求解方程组Ax=d 的方法称为追赶法. LU分解法：将系数矩阵A转化为A=L*U，L为单位下三角矩阵，U为普通上三角矩阵，然后通过解方程组l*y=b，u*x=y,来求解x。主要仪器设备台式或笔记本计算机实验记录(写出实验内容中的程序代码和运行结果)(可分栏或加页)高斯消元法#include#define n 3main()int i,j,k;float ann,cnn,bn,dn;for(i=0;in;i+)for(j=0;jn;j+)scanf(%f,&aij);cij=aij;scanf(%f,&bi);di=bi;for(k=0;kn;k+)bk=dk/ckk;for(i=0;in;i+)if(i=k) continue;cik=cik/ckk;for(j=k+1;jn;j+)akj=ckj/ckk;aij=cij-cik*ckj;bi=di-cik*dk;for(i=0;in;i+)di=bi;for(j=k+1;jn;j+)cij=aij;for(i=0;in;i+)printf(b%d=%fn,i,bi);LU分解法：#include #include #define L 30double aLL,bL,lLL,uLL,xL,yL;int main() int n,i,j,k,r;printf(请输入矩阵元次：n); scanf(%d,&n);printf(请输入矩阵各项：n); for(i=1;i=n;+i) for(j=1;j=n;+j) scanf(%lf,&aij); printf(请输入方程组的常数项：n); 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*xj; xi/= uii; for(i=1;i=n;+i) printf(%0.2lfn,xi); return 0;追赶法#include stdio.h#define n 5main() float an,bn,cn-1,dn,t;int i;scanf(%f%f%f,&b0,&c0,&d0);for(i=1;in-1;i+)scanf(%f%f%f%f,&ai,&bi,&ci,&di);scanf(%f%f%f,&an-1,&bn-1,&dn-1);c0=c0/b0;d0=d0/b0;for(i=1;i=0;i-) di=di-ci*di+1;for(i=0;in;i+)printf(d%d=%fn,i,di);实验结果和分析1高斯消元法2LU分解3追赶法实验分析：高斯消元法，是先消元，再回带的过程。由程序段可以发现，始终消去对角线下方的元素。从消元过程可以看出，对于n阶线性方程组，只要各步主元素不为零，经过n-1步消元，就可以得到一个等价的系数矩阵为上三角形阵的方程组，然后再利用回代过程可求得原方程组的解。LU分解法，分解矩阵为单位下三角阵L与上三角阵U的乘积，然后解方程组Ly=b，回代，解方程组Ux=y。其中的L为n阶单位下三角阵、U为上三角阵.对于追赶法，追赶法是适用于三角矩阵的线性方程组的求解的方法，并不适用于其他类型矩阵。心得体会（遇到的问题和解决方法）本次实验难度比较大，在编译时经常出现各种错误，程序代码也比较繁琐，深深感觉到自己的上机操作能力有限，应加强自己的编程能力，以后要继续努力。实验名称 实验三 线性方程组的迭代求解实验目的和要求使用雅可比迭代法或高斯-赛德尔迭代法对下列方程组进行求解。 实验内容设线性方程组 Ax=b的系数矩阵A可逆，且主对角元素a11,a22,ann均不为零，令D=diag(a11,a22,ann)并将A分解成 A=(A-D)+D从而线性方程组可写成 Dx=(D-A)x+b则有迭代公式x(k+1)=B1x(k)+f1其中，B1=I-D-1A,f1=D-1b。主要仪器设备台式或笔记本计算机实验记录(写出实验内容中的程序代码和运行结果)(可分栏或加页)雅可比迭代法#include#include#define n 3void f(float *b,float x)float yn+1=0,0,0,1;int i,j,k;dok=0;for(i=0;in+1;i+)xi=yi;for(i=0;in;i+)yi=0;for(j=0;jn+1;j+)yi+=*(b+(n+1)*i+j)*xj;for(i=0;in;i+)if(fabs(yi-xi)=0.5e-3) k+;if(k=3) break;while(1);for(i=0;in;i+)printf(y%d=%fn,i,yi);main()float bnn+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);高斯赛德尔迭代法#include iostream#include iomanipusing namespace std;int main()int i,j,k=0,m,n;double t1,t2,e1,e2=0.0; coute1;coutm;coutn;coutendl;double (*a)=new double *m;/生成二维动态数组for(i=0;i=m;i+)ai=new doublen;double (*b)=new double m;double (*x)=new double n;cout请输入系数矩阵：endl;cout-endl;for(int num1=0;num1m;num1+)for(int num2=0;num2anum1num2;coutendl;cout输入的系数矩阵为：endl;for (int num3=0;num3m;num3+)for(int num4=0;num4n;num4+)coutanum3num4 ;coutendl;cout-endl;cout请输入矩阵b：endl;cout-endl;for(int num5=0;num5bnum5;cout输入的矩阵b为：endl;for(int num6=0;num6m;num6+)coutbnum6 ; coutendl; cout-endl;for(int num7=0;num7n;num7+)xnum7=0.0000;do cout第k次迭代值：;e2=0.0;for(i=0;im;i+) double sum=0.0;for(j=0;j=0?(xi)-t1:t1-(xi);e2=(e2=t2?e2:t2);coutsetprecision(8)xi ;cout=e1&k30);cout共迭代了k次;deletea;deleteb;deletex;return 0 ;实验结果和分析1雅克比迭代2高斯赛德尔迭代 实验分析：使用这两种方法都可以求出方程的解，高斯赛德尔迭代法所需的迭代次数比雅克比迭代少，能够更早的达到精度要求。但是雅克比的时效性要比高斯赛德尔的好。心得体会（遇到的问题和解决方法）本次试验，让我对这两种方法更加理解，在编程操作上也更加熟练，今后继续努力，不断丰富自己的知识，增强操作能力。实验名称 实验四 代数插值和最小二乘法拟合实验目的和要求1使用拉格朗日插值法或牛顿插值法求解：已知f(x)在6个点的函数值如下表所示，运用插值方法，求f(0.596)的近似值。x0.400.550.650.800.901.05f(x)0.410750.578150.696750.888111.026521.25386 2给定数据点（xi ,yi），用最小二乘法拟合数据的多项式，并求平方误差。xi00.50.60.70.80.91.0yi11.751.962.192.442.713.00实验内容1设函数在区间a,b上n+1互异节点x0,x1,xn上的函数值分别为y0,y1,yn,求n次插值多项式Pn(x),满足条件Pn(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)2 建立正规方程组：（xij+k）ak=xijyi ,j=0,1,n 平方误差：I=（akxik-yi）2 对给定数据点(Xi，Yi)(i=0,1，m），在取定的函数类 中，求p(x），使误差的平方和E2最小，E2=p(Xi)-Yi2。从几何意义上讲，就是寻求与给定点 (Xi，Yi)(i=0,1，m）的距离平方和为最小的曲线y=p(x）。函数p(x）称为拟合函数或最小二乘解，求拟合函数p(x）的方法称为曲线拟合的最小二乘法。得到的两个关于a0、 a1为未知数的两个方程组，解这两个方程组得出：a0 = （Yi) / m - a1（Xi) / m a1 = mXi Yi - （Xi Yi) / mXi2 - （Xi)2 ) 即最终的拟合多项式各项系数主要仪器设备台式或笔记本计算机实验记录(写出实验内容中的程序代码和运行结果)(可分栏或加页)1#include #include #include #include void difference(float *x,float *y,int n) float *f; int k,i; f=(float *)malloc(n*sizeof(float); for(k=1;k=n;k+) f0=yk; for(i=0;ik;i+) fi+1=(fi-yi)/(xk-xi); yk=fk; return; int main() int i,n; float x20,y20,xx,yy; printf(请输入数据个数n：); scanf(%d,&n);printf(n); for(i=0;i=0;i-)yy=yy*(xx-xi)+yi; printf(n近似值为：(%f)=%fn,xx,yy); 2#include#include#define N 15double power(double &a,int n)double b=1;for(int i=0;in;i+)b*=a;return b;void Gauss();double XN,YN,sumXN,sumYN,aNN,bN,lNN,xN;int main()ofstream outdata;ifstream indata;double s;int i,j,k,n,index;coutn;coutendl;cout请输入X和Y:endl; /输入给定数据for(i=0;in;i+)coutXiXi;sumX1+=Xi;coutYiYi;sumY1+=Yi;coutendl;coutsumX1=sumX1tsumY1=sumY1endl;coutindex;coutendl;i=n;sumX0=i;for(i=2;i=2*index;i+)sumXi=0;for(j=0;jn;j+)sumXi+=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