计算数值方法实验报告--太原理工大学.doc

《计算数值方法实验报告--太原理工大学.doc》由会员分享，可在线阅读，更多相关《计算数值方法实验报告--太原理工大学.doc（25页珍藏版）》请在沃文网上搜索。

1、计算数值方法实验报告计算数值方法实验报告学院：软件学院专业：软件工程班级：软件1012班学号：2010004719姓名：乔婧峰太原理工大学学生实验报告学院名称软件学院专业班级软件1012班学号2010004719学生姓名乔婧峰实验日期20124成绩课程名称数值计算方法实验题目实验一 二分法一、课题名称方程求根：熟悉使用、迭代法、牛顿法、割线法等方法对给定的方程进行根的求解。选择上述方法中的两种方法求方程：二分法f(x)=x3+4x2-10=0在1,2内的一个实根，且要求满足精度|x*-xn|0.510-5迭代法：用迭代公式x=f(x)进行迭代计算，直到满足|x*-xn|0.510-5 为止 。

2、割线法：x=x-f(x)/g(x),其中f(x)为给定的函数，g(x)为给定函数的导数，直到满足|x*-xn|0.510-5 为止 。二、目的和意义（1）了解非线性方程求根的常见方法，如二分法、牛顿法、割线法。（2）加深对方程求根方法的认识，掌握算法。（3）会进行误差分析，并能对不同方法进行比较。三、计算公式f(x)在区间（x，y）上连续 先找到a、b属于区间（x，y），使f(a)，f(b)异号，说明在区间(a,b)内一定有零点，然后求f(a+b)/2, 现在假设f(a)0,ab 如果f(a+b)/2=0，该点就是零点， 如果f(a+b)/2a，从开始继续使用 中点函数值判断。 如果f(a+b

3、)/20，则在区间(a,(a+b)/2)内有零点，(a+b)/2a; float t, x; x=a; do x=sqrt(10-x*x*x)/4); t=a; a=x; while(fabs(a-t)0.5*1e-5); printf(x=%f,a); system(pause);割线法: #include stdafx.h#includestdio.h#includemath.h#includeiostreamusing namespace std;float main() float c,a=1.0,b=2.0; /cinab; while(1) c=b-(b*b*b+4*b*b-10)

4、*(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;六、结果讨论和分析割线法: 迭代法:分析：使用不同的方法，可以不同程度的求得方程的解，不同的方法速度不同。实验地点ZSA401指导教师李志学院名称软件学院专业班级软件1012班学号2010004719学生姓名乔婧峰实验日期2012.4成绩课程名称数值计算方法实验题目实验二 线性方程组的直接解法一、课题名称线性方程组的直接解法合理利用Gauss消元法、LU分解法、追赶法求解下列方程组： （n=5,10,100）二、目的和意义（1）了解线性方

5、程组常见的直接解法，如Guass消元法、LU分解法、追赶法。（2）加深对线性方程组求解方法的认识，掌握算法。（3）会进行误差分析，并能对不同方法进行比较。三、计算公式 高斯分解法：将原方程组化为三角形方阵的方程组：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 (k=n-1,n-2, ,2,1) LU分解法：将系数矩阵A转化为A=L*U, L为单位下三角矩阵，U为普通上三角矩阵，然后通过解方程组l

6、*y=b,u*x=y,来求解x.追赶法:用来求对角方程组；将系数矩阵A转化为A=L*U, L为普通下n-1对角矩阵，U为单位上n-1对角矩阵，然后通过解方程组l*y=b,u*x=y,来求解x.四、主要仪器设备Vc2008,hp五、结构程序设计 Gauss消元法：#include stdafx.h#includestdio.h#includeiostreamusing namespace std; float main() float a34=1,2,3,14,0,1,2,8,2,4,1,13; float x3; float sum=0; int k,i,j; for(k=0;k2;k+) f

7、or(i=k+1;i3;i+) for(j=k+1;j4;j+)aij=aij-aik/akk*akj; for(i=0;i3;i+) for(j=0;j4;j+) printf(a%d%d=%f,i,j,aij); cout=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 stdafx.h#include #include #define L 30 double a L L , b L , l L L , u L

8、L , x L , y L ; int main() int n, i, j, k, r; scanf( %d, &n ); for ( i = 1; i = n; +i ) for ( j = 1; j = n; +j ) scanf( %lf, &a i j ); for ( i = 1; i = n; +i ) scanf( %lf, &b i ); for ( i = 1; i = n; +i ) for ( j = 1; j = n; +j ) l i j =0; u i j = 0.0; for ( k = 1; k = n; +k ) for ( j = k; j = n; +j

9、 ) u k j = a k j ; for ( r = 1; r k; +r ) u k j -= l k r * u r j ; for ( i = k + 1; i = n; +i ) l i k = a i k ; for ( r = 1; r k; +r ) l i k -= l i r * u r k ; l i k /= u k k ; l k k = 1.0; for ( i = 1; i = n; +i ) y i = b i ; for ( j = 1; j 0; -i ) x i = y i ; for ( j = i + 1; j = n; +j ) x i -= u

10、i j * x j ; x i /= u i i ; for ( i = 1; i = n; +i ) printf( %0.2lfn, x i ); return 0;追赶法：#include stdafx.h#include stdio.hvoid main() 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,

11、&bi,&ci,&di); fscanf(f,%lf%lf%lf,&an,&bn,&dn); fclose(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);Zgf.txt文件中的内容是：52 1 -71 2 1 -5 1 2 1 -5 1 2 1 -5 1 2 -5六、结果讨论和分析Gauss消元法：LU分解法：追赶法：分析从消元过程可以看出，对于n阶线性方程组，只要各步主元素不为零，经过n-1步消元，就可以得到一个等价的系数矩阵

12、为上三角形阵的方程组，然后再利用回代过程可求得原方程组的解.消元过程相当于分解 A为单位下三角阵L与上三角阵U的乘积，解方程组Ly=b回代过程就是解方程组Ux=y。其中的L为n阶单位下三角阵、U为上三角阵. 在 A 的LU 分解中, L取下三角阵, U 取单位上三角阵,这样求解方程组Ax=d 的方法称为追赶法.实验地点ZSA401指导教师李志学院名称软件学院专业班级软件1012学号2010004719学生姓名乔婧峰实验日期20114成绩课程名称数值计算方法实验题目实验三 线性方程组的迭代解法一、课题名称线性方程组的迭代解法使用雅可比迭代法或高斯-赛德尔迭代法对下列方程组进行求解。 二、目的和意

13、义学习使用雅可比迭代法或高斯-赛德尔迭代法三、计算公式雅克比迭代法：设线性方程组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。四、主要仪器设备Vc2008,hp五、结构程序设计高斯-赛德尔迭代法：#include stdafx.h#include #include void main() float a33=10,-1,-2,-1,10,-2,-1,-1,5,b3=7.2,8.

14、3,4.2;float x3=0,0,0,sum1,sum2;int i,j,k,n=3;for (k=0;k10;k+) for(i=0;in;i+) sum1=0; sum2=0;for(j=0;ji;j+) sum1=sum1+aij*xj; for(j=i+1;j3;j+) sum2=sum2+aij*xj; xi=(bi-sum1-sum2)/aii; for(i=0;in;i+) printf(x%d=%f,i+1,xi); printf(n); 雅克比迭代：#include stdafx.h#include #include void main() float a33=10,-1

15、,-2,-1,10,-2,-1,-1,5,b3=7.2,8.3,4.2;float x3=0,0,0,sum1;int i,j,k,n=3;for (k=0;k10;k+) for(i=0;i3;i+) sum1=0; for(j=0;jn;j+) if(i=j) continue; sum1=sum1+aij*xj; xi=(bi-sum1)/aii; for(i=0;in;i+) printf(x%d=%f,i+1,xi);printf(n); 六、实验结果与分析：高斯-赛德尔迭代法：雅克比迭代：分析：使用高斯-赛德尔和雅克比迭代都可以求出方程组的解，但是利用高斯-赛德尔迭代法所需的迭代次

16、数比雅克比迭代少，能够更早的达到精度要求。实验地点ZSA401指导教师李志学院名称软件学院专业班级软件1012班学号2010004719学生姓名乔婧峰实验日期20114成绩课程名称数值计算方法实验题目实验四 最小二乘法拟合多项式一、课题名称（1）了解矩阵特征值与特征向量问题解法，掌握幂法。（2）加深对矩阵特征值与特征向量问题求解方法的认识，掌握算法。（3）会进行误差分析。二、目的和意义 学习使用最小二乘法拟合多项式三、计算公式幂法：由已知的非零向量x0和矩阵A的乘幂构造向量序列xn以计算矩阵A的按模最大特征值及其特征向量的方法，称为幂法。迭代公式：结果可取 四、主要仪器设备Vc2008,hp五

17、、结构程序设计 五、结果讨论和分析分析： 幂法是一种求任意矩阵A的按模最大特征值及其对应特征向量的迭代算法。该方法的最大优点是计算简单，容易在计算机上实现，对稀疏矩阵较为适合，但有时收敛速度很慢。实验地点综合楼六层606室指导教师王峥学院名称计算机科学与技术专业班级计算机学号1111111111学生姓名某某实验日期2011-6-20成绩课程名称数值计算方法实验题目实验五 代数插值一、 课题名称使用拉格朗日插值法或牛顿插值法求解：已知f(x)在6个点的函数值如下表所示，运用插值方法，求f(0.596)的近似值。x0.400.550.650.800.901.05f(x)0.410750.57815

18、0.696750.888111.026521.25386二、 目的和意义学习使用拉格朗日插值法或牛顿插值法求解三、 计算公式设函数在区间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)四、 结构程序设计#inclu

19、de#include#includetypedef struct data float x; float y;Data; /变量x和函数值y的结构Data d20; /最多二十组数据float f(int s,int t) /牛顿插值法，用以返回插商 if(t=s+1) return (dt.y-ds.y)/(dt.x-ds.x); else return (f(s+1,t)-f(s,t-1)/(dt.x-ds.x); float Newton(float x,int count) int n; while(1) coutn; if(n=count-1)/ 插值次数不得大于count1次 br

20、eak; else system(cls); float t=1.0; float y=d0.y; float yt=0.0; for(int j=1;j=n;j+) t=(x-dj-1.x)*t; yt=f(0,j)*t; y=y+yt; return y;float lagrange(float x,int count) float y=0.0; for(int k=0;kcount;k+)/这儿默认为count1次插值 float p=1.0;/初始化p for(int j=0;jcount;j+) /计算p的值 if(k=j)continue;/判定是否为同一个数 p=p*(x-dj.

21、x)/(dk.x-dj.x); y=y+p*dk.y;/求和 return y;/返回y的值void main() float x,y; int count; while(1) coutcount; if(count=20) break;/检查输入的是否合法 system(cls); /获得各组数据 for(int i=0;icount;i+) cout请输入第i+1di.x; cout请输入第i+1di.y; system(cls); coutx; while(1) int choice=3; cout请您选择使用哪种插值法计算：endl; cout (0):退出endl; cout (1)

22、:Lagrangeendl; cout (2):Newtonendl; coutchoice;/取得用户的选择项 if(choice=2) cout你选择了牛顿插值计算方法，其结果为：; y=Newton(x,count);break;/调用相应的处理函数 if(choice=1) cout你选择了拉格朗日插值计算方法，其结果为：; y=lagrange(x,count);break;/调用相应的处理函数 if(choice=0) break; system(cls); cout输入错误!endl; coutx , yendl;/输出最终结果 五、结果讨论和分析分析：拉格朗日插值的优点是插值多

23、项式特别容易建立，缺点是增加节点是原有多项式不能利用，必须重新建立，即所有基函数都要重新计算，这就造成计算量的浪费。实验地点综合楼六层606室指导教师王峥学院名称计算机科学与技术专业班级计算机学号1111111111学生姓名某某实验日期2011-6-20成绩课程名称数值计算方法实验题目实验六 最小二乘法拟合多项式一、课题名称给定数据点（xi ,yi），用最小二乘法拟合数据的多项式，并求平方误差。xi00.50.60.70.80.91.0yi11.751.962.192.442.713.00二、目的和意义1熟练运用已学计算方法求解方程组2加深对计算方法技巧，选择正确的计算方法来求解各种方程组3培

24、养使用电子计算机进行科学计算和解决问题的能力三、计算公式建立正规方程组：（xij+k）ak=xijyi ,j=0,1,n 平方误差：I=（akxik-yi）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;void main()ofstream outdata;ifstream indata;double s;int i,j,

25、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+)su

26、mYi+=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=1;i-)s=0;for(j=i+1;j=index+1;j+)s=s+aij*xj;xi=(bi-s)/aii;cout拟合系数为:; /输出拟合系数for(i=1;i=index+1;i+)coutxit;double m=0;coutendl平方误差为:;for(i=0;in;i+) double t=x1+x2*Xi-Yi;m=m+power(t,2);coutmendl;五、结果讨论和分析分析：数据拟合的具体作法是：对给定的数据（xi ,yi）（i=0,1,m），在取定的函数类中，求p(x)属于此函数类，使误差ri=p(xi)- yi (i=0,1,m)的平方和最小，即ri2=（p(xi)-yi）2=min从几何意义上讲，就是寻求与给定点（x