计算数值实验报告

1、本科实验报告课程名称：计算数值方法实验地点：综合楼五层506室专业班级：计科 1002 学号：J4学生姓名： xxx指导教师：王峥2012 年 6 月 20太原理工大学学生实验报告学院名称计算机科学与技术专业班级计科1002学号14学生姓名xxx实验日期成绩课程名称计算数值方法实验题目实验一方程求根、课题名称方程求根：熟悉使用、迭代法、牛顿法、割线法等方法对给定的方程进行根的求解。选择上述方法中的两种方法求方程：二分法f(x)=x 3+4x2-10=0在1,2内的一个实根，且要求满足精度|x*-xn|< X 10-5迭代法：用迭代公式 x=f(x)进行迭代计算，直到满足|x*-xn|&l

2、t;Xl05为止。二分法：设 f (x)在a,b上连续，且 f (al) *f (x1) <0,记(a2, b2) = (x1, bl) 带入计算式进行计算直到|x *-xn|< X10-5为止。二、目的和意义(1) 了解非线性方程求根的常见方法，如二分法、迭代法、牛顿法、割线法。(2)加深对方程求根方法的认识，掌握算法。会进行误差分析，并能对不同方法进行 比较。三、计算公式(1)迭代法1 ).首先对给定的计算公式进行变形使其能够迭代或者找出相应迭代速度较快的 式子。2 ).带入求好的式子到循环中去比如：xk 1(xk)(k 0,1, 2,L)(2)二分法：f(x)在区间(x, y

3、)上连续1) .先找到a、b属于区间(x, y),使f(a) , f(b)异号，说明在区间(a,b)内一 定有零点，然后求f(a+b)/2,2) .如果f(a+b)/2=0 ，该点就是零点，如果f(a+b)/2<0, 则在区间(a+b)/2 , b)内有零点,反之在(a, (a+b)/2 )内 有零点带入1)中继续。四、主要仪器设备Vc+ C-free CodeBlocks五、结构程序设计迭代法：#include<>#include<>main()int i;double xn15,y,x1,x2,m ;printf("请输入 x1,x2 的值：n&qu

4、ot;);scanf("%lf%lf",&x1,&x2);printf("请输入精度要求：n");scanf("%lf",&m);printf(" n xn'n");i=0;doxn0=(x1+x2)/2 ;xni+1= sqrt(10/(4+xni);0.3 10 1559.1431Xi59.174215Xi123Xi1401285.2916.13012x246.7887210x2x2211311.2952x314836X34x31211x421261120X41x1721x251

5、21 Xn 1512 Xn5追赶法:用来求对角方程组；将系数矩阵A转化为A=L*U, L为普通下n-1对角矩阵，U为单位上n-1对角矩阵，然后通过程组l*y=b,u*x=y,来求解x.、主要仪器设备c+ C-free CodeBlocks、结构程序设计Gauss消元法：#include <iostream>#include <cmath>using namespace std;int main()int n,i,j,k;double a100100,b100,o;cout<<"输入未知数个数:"<<endl;cin>&g

6、t;n;cout<<"输入数列:"<<endl;for (i=1;i<=n;i+)for (j=1;j<=n+1;j+)cin>>a皿;for (i=1;i<=n;i+)for (j=i+1;j<=n;j+)if (fabs(aji)>1e-7)o=aii/a皿i;for (k=i;k<=n+1;k+) ajk=ajk*o-aik;for (i=n;i>0;i-)bi=a皿n+1/aii;for (j=i-1;j>0;j-)a皿n+1=ajn+1-bi*a皿i;cout<<&qu

7、ot;解得:"<<endl;for (i=1;i<=n;i+) cout<<bi<<endl;lfn",i,di);、结果讨论和分析消元法:列主元素消元法:解法:追赶法:讨论从消元过程可以看出，对于n阶线性方程组，只要各步主元素不为零，经过 n-1步消元，就可以得到一个等价的系数 三角形阵的方程组，然后再利用回代过程可求得原方程组的解.由于列主元素法相似且优于完全主元素法所以省即消元过程相当于分解A为单位下三角阵L与上三角阵U的乘积，解方程组Ly=b回代过程就是解方程组Ux=y。其中白 单位下三角阵、U为上三角阵.在A的LU分解中,

8、L取下三角阵,U取单位上三角阵,这样求解方程组Ax=d的方左 法。另外是追赶法和其他方法求同一方程结果不一样，我多次修改源程序，也不知道原因。再就是追赶法有很大的局 改良。图：实验地点综合楼五层506室指导教师王峥太原理工大学学生实验报告学院名称计算机科学与技术专业班级计科1002学号14学生姓名xxx实验日期成绩课程名称计算数值方法实验题目实验三线性方程组的迭代解法、课题名称线性方程组的迭代解法使用雅可比迭代法或高斯-赛德尔迭代法对下列方程组进行求解。10x1 x2 2x37.2xi 10X2 2x3 8.3x1 x2 5x3 4.2、目的和意义学习使用雅可比迭代法或高斯-赛德尔迭代法雅克比

9、迭代法：设线性方程组Ax=b的系数矩阵A可逆且主对角元素aii,a22,ann均不为零，令D=diag（ai1,a 22,，a nn）并将A分解成A=(A-D)+D从而线性方程组可写成DX=(D-A) x+b则有迭代公式x(k+1)=Bx(k)+fi其中，Bi=I-D1A,f 尸D1b。四、主要仪器设备Vc+ C-free CodeBlocks五、结构程序设计雅克比迭代法：#include<>#include<>main()int i;double x120 ,x220,x320;double x10, x20, x30;printf("请输入 x1,x2,x

10、3 的初值：n");scanf("%lf%lf%lf",&x10,&x20, &x30);printf(" n x1nx2nx3n n");for(i=0;i<18;i+)x10=x10;x20=x20;x30=x30;x1i+1=*x2i+*x3i+;x2i+1=*x1i+*x3i+;x3i+1=*x1i+*x2i+;printf("%5d %5lf %5lf %5lfn",i,x1i,x2i,x3i); 六、实验结果与分析：雅克比迭代法:分析讨论：其实，这两个迭代法是之前迭代法的升级，多了几

11、个迭代式子而已，而且两者相差不大比较简单，所以选择了雅克比迭代法进行求解，但是没有与另一种方法 高斯赛德尔迭代法进行实质性的比较。流程图：实验地点综合楼五层506室指导教师王峥太原理工大学学生实验报告学院名称计算机科学与技术专业班级计科1002学号14学生姓名xxx实验日期成绩课程名称计算数值方法实验题目实验四 矩阵特征值与特征向量问题一课题名称使用幕法求A模为最大的特征值及其相应的特征向量。210A 1 21012二、目的和意义(1) 了解矩阵特征值与特征向量问题解法，掌握幕法。(2)加深对矩阵特征值与特征向量问题求解方法的认识，掌握算法。三、计算公式幕法：由已知的非零向量x0和矩阵A的乘幕

12、构造向量序列xn以计算矩阵A的按模最大特征值及其特征向量的方法，称为幕法。yk m kxk迭代公式：Ax k 1max( yk), k 1,2， ykmk结果可取1 mk1 yk或 1 xk四、主要仪器设备Vc+ C-free CodeBlocks五、结构程序设计源代码：#include<>#include<>#define N 3#define eps 1e-6#define KM 30float MaxValue(float x,int n)float Max=x0;int i;for (i=1;i<n;i+)if(fabs(xi)>fabs(Max)M

13、ax=xi;return Max;void PowerMethod(float *A)float UN,VNj1,r2,temp;int i,j,k=0;while(k<KM)k+;for(i=0;i<N;i+)temp=0;for(j=0;j<N;j+)temp+=*(A+i*N+j)*Uj;Vi=temp;for(i=0;i<N;i+)Ui=Vi/MaxValue(V,N);if(k=1)r1=MaxValue(V,N);else r2=MaxValue(V,N);if(fabs(r2-r1)<eps)break;r1=r2;printf("r=%f

14、n",r2);for(i=0;i<N;i+)printf("y%d=%fn",i+1,Ui);void main()float ANN尸2,-1,0,-12-1,0,-1,2 ;float UN; /A的值U0=1;U1=1;U2=1;/x0 的值PowerMethod(A0);六、结果讨论和分析分析讨论由于该程序将A矩阵和x0的值编写在程序中，所以要想修改成其他矩阵时 比较麻烦，所以也有f的局限性。幕法是一种求任意矩阵A的按模最大特征值及其对应特征向量的迭代算法。该方法的最大优点是计算简单，容易在计 算机上实现，对稀疏矩阵较为适合，但有时收敛速度很慢。流程

15、图：实验地点综合楼五层506室指导教师王峥太原理工大学学生实验报告学院名称计算机科学与技术专业班级计科1002学号14学生姓名xxx实验日期成绩课程名称计算数值方法实验题目实验五代数插值一课题名称(使用拉格朗日插值法或牛顿插值法求解：已知 f(x)在6个点的函数值如下表所示，运用 插值方法，求f的近似值。Xf(x)、目的和意义学习使用拉格朗日插值法或牛顿插值法求解三、计算公式设函数在区间a,b上n+1互异节点Xo,X 1,X n上的函数值分别为yo,y 1,y n,求n次插值多项式Pn(x),满足条件Pn(Xj)=yj, j=0,1,n令Ln(x)=y ol o(x)+y il i(x)+ +

16、ynl n(x)= /y/i(x)其中l o(x),l l(x),，l n(x)为以X0,X 1，,X n为节点的n次插值基函 数，则Ln(x)是一次数不超过n的多项式，且满足Ln(x j)=yj, L=0,1,，n再由插值多项式的唯一性，得Pn(x)三 Ln(x)四、主要仪器设备Vc+ C-free CodeBlocks五、结构程序设计#include <>#include <>main()int I;char L;double M100100;double x100,y100;double X=1,xx=0,w=1,N=0,P,R=1;int n=5;/cout&l

17、t;<”请输入所求均差阶数："/求所有阶差cin>>n;/for(int i=0;i<=n;i+)/*cout<<"请输入 x"<<i<<"的值："<<endl;cin>>xi;cout<<"请输入 y"<<i<<"的值："<<endl;cin>>yi;Mi0=xi;Mi1=yi; */用二维保存所有数据M00=;M01=;M10=;M11=;M20=;M21=;

18、M30=;M31=;M40=;M41=;M50=;M51=;for( int j=2;j<=n+1;j+)for(int i=1;i<=n;i+)Mij=(Mij-1-Mi-1j-1)/(Mi0-Mi-j+10);) for(int i=1;i<=n;i+)cout<<"其"<<i<<”阶均差为："<<Mii+1卜<endl;cout<<"请输入x的值：x="cin>>xx;for(int i=0;i<n;i+)X*=xx-Mi0;N+=Mi+

19、1i+2*X;P=M1+N;cout<<"其函数值：y="<<P<<endl;/ )六、结果讨论和分析分析讨论拉格朗日插值的优点是插值多项式特别容易建立，缺点是增加节点是原有多项式不能利用，必须重新建立，即所有基函数都要重新计算，这就造成计算量的浪费。所以该程序 选择了牛顿法，由于在输入数据的时候比较麻烦所以将数据存在了程序里面，只能在程序里修改数据。另外由于精度问题计算结果有一定的误差。流程图：实验地点综合楼五层506室指导教师王峥太原理工大学学生实验报告学院名称计算机科学与技术学生姓名课程名称xxx计算数值方法专业班级计科1002学号

20、14实验日期成绩实验题目实验六最小一乘法拟合多项式、课题名称给定数据点（Xi , yi）,用最小二乘法拟合数据的多项式，并求平方误差。Xi0yi1、目的和意义1 .熟练运用已学计算方法求解方程组2 .加深对计算方法技巧，选择正确的计算方法来求解各种方程组3 .培养使用电子计算机进行科学计算和解决问题的能力 三、计算公式建立正规方程组：.(汇 Xij+k) ak=Exijyi , j=o,i,n平方误差：I =.(. akXi -y i)四、主要仪器设备Vc+ C-free CodeBlocks五、结构程序设计源代码：#include<>#include<>#define

21、 N 15double power(double &a,int n)double b=1;for(int i=0;i<n;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,k,n,index;/cout<<"请输入已知点的个数 n="/cin>>n;n=7;cout<<endl;/cout<<"

22、请输入X和Y:"<<endl; /输入给定数据X0=;Y0=;X1=;Y1=;X2=;Y2=;X3=;Y3=;X4=;Y4=;X5=;Y5=;X6=;Y6=; /绑定数据可以解绑由下for循环输入任何数据for(i=0;i<n;i+)/cout<<"X"<<i<<"="/cin>>Xi;sumX1+=Xi;/cout<<"Y"<<i<<"="/cin>>Yi;sumY1+=Yi;/cout&l

23、t;<endl;cout<<"sumX1="<<sumX1<<"t"<<"sumY1="<<sumY1<<endl;cout<<"请输入拟合次数 index="cin>>index;cout<<endl;i=n;sumX0=i;for(i=2;i<=2*index;i+)(sumXi=0;for(j=0;j<n;j+)sumXi+=power(Xj,i);cout<<"

24、sumX"<<i<<"="<<sumXi<<endl;for(i=2;i<=index+1;i+)(sumYi=0;for(j=0;j<n;j+)sumYi+=power(Xj,i-1)*Yj;cout<<"sumY"<<i<<"="<<sumYi<<endl;for(i=1;i<=index+1;i+)/建立正规方程组for(j=1;j<=index+1;j+)aij=sumXi+j-2;bi

25、=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+1index+1;for(i=index;i>=1;i-)s=0;for(j=i+1;j<=index+1;j+)s=s+aij*xj;xi=(bi-s)/aii;cout<&l