计算数值实验报告(太原理工大学)

《计算数值实验报告(太原理工大学)》由会员分享，可在线阅读，更多相关《计算数值实验报告(太原理工大学)（29页珍藏版）》请在装配图网上搜索。

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

2、迭代法：用迭代公式x=f(x)进行迭代计算，直到满足|x*-xn|0.510-5 为止 。二分法：设f（x）在a,b上连续，且f（a1）*f（x1）0,记（a2，b2）=（x1，b1）带入计算式进行计算 直到 |x*-xn|0.510-5 为止 。二、目的和意义（1）了解非线性方程求根的常见方法，如二分法、迭代法、牛顿法、割线法。（2）加深对方程求根方法的认识，掌握算法。会进行误差分析，并能对不同方法进行比较。三、计算公式（1）迭代法 1）.首先对给定的计算公式进行变形使其能够迭代或者找出相应迭代速度较快的式子。 2）.带入求好的式子到循环中去比如：（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)/20,则在区间（(a+b)/2，b)内有零点，反之在（a，(a+b)/2）内有零点 带入1）中继续。四、主要仪器设备Vc+ 9.0 C-free CodeBlocks五、结构程序设计迭代法: #include #include main() int i; double xn15,y,x1,x2,m ; printf(请输入x1,x2的值：n ); scanf(%lf%lf,&x1,&x2); printf

4、(请输入精度要求：n ); scanf(%lf,&m); printf( n xnn);i=0; do xn0=(x1+x2)/2 ; xni+1= sqrt(10/(4+xni); /迭代printf(%5d %5lfn,i,xni); y= fabs(xni+1-xni) ; i+; if(ym)break;while(1); 二分法：#include #include main() int m,n,o,p; double a,b,l; printf(请输入x3, x2, x的系数和常数p：n);scanf(%d%d%d%d,&m,&n,&o,&p); /1 4 0 -10 printf(

5、请输入x1,x2:n);scanf(%lf%lf,&a,&b); /1 2 printf(请输入精度要求:n);scanf(%lf,&l); /0.5x105 printf( n an bn xn f(xn)n); double x,fx; int i=1; do x=(b+a)/2; fx=m*x*x*x+n*x*x+o*x+p; printf(%5d %5f %5f %5f %5fn,i,a,b,x,fx); i+; if(fx=0) break; if(fx0) b=x; else if(fx0) a=x; if(b-a)l) break; /进行计算并返值 while(1);六、结果讨

6、论和分析 二分法： 迭代法： 分析讨论：使用不同的方法，可以不同程度的求得方程的解，不同的方法速度不同，求得的结果也稍有区别，当然和要求精度也有关系。刚开始的时候用数组对二分法进行求解，发现 循环到第二次 就无法实现值的传递，于是换了另外一种方法代替了数组。流程图：实验地点综合楼五层506室指导教师王峥太原理工大学学生实验报告学院名称计算机科学与技术专业班级计科 1002学号2010001414学生姓名 xxx实验日期2012.6.5成绩课程名称计算数值方法实验题目实验二 线性方程组的直接解法一、课题名称线性方程组的直接解法合理利用Gauss消元法、LU分解法、追赶法求解下列方程组： （n=5

7、,10,100）二、目的和意义（1）了解线性方程组常见的直接解法，如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为单位下三

8、角矩阵，U为普通上三角矩阵，然后通过解方程组l*y=b,u*x=y,来求解x.追赶法:用来求对角方程组；将系数矩阵A转化为A=L*U, L为普通下n-1对角矩阵，U为单位上n-1对角矩阵，然后通过解方程组l*y=b,u*x=y,来求解x.四、主要仪器设备Vc+ 9.0 C-free CodeBlocks五、结构程序设计 Gauss消元法： #include #include using namespace std;int main()int n,i,j,k;double a100100,b100,o;cout输入未知数个数:n;cout输入数列:endl; for (i=1;i=n;i+)fo

9、r (j=1;jaij;for (i=1;i=n;i+)for (j=i+1;j1e-7)o=aii/aji;for (k=i;k0;i-)bi=ain+1/aii;for (j=i-1;j0;j-)ajn+1=ajn+1-bi*aji;cout解得:endl;for (i=1;i=n;i+) coutbiendl;/system(pause);return 0;列主元素消元法： #include #include #define N 20 using namespace std; void load(); float aNN; int m; int main() int i,j; int c

10、,k,n,p,r; float xN,lNN,s,d; coutm; coutendl; cout请按顺序输入增广矩阵a:endl; load(); for(i=0;im;i+) for(j=i;jfabs(aii)?j:i; /*找列最大元素*/ for(n=0;nm+1;n+) s=ain; ain=acn; acn=s; /*将列最大数防在对角线上*/ for(p=0;pm+1;p+) coutaipt; coutendl; for(k=i+1;km;k+) lki=aki/aii; for(r=i;r=0;i-) d=0; for(j=i+1;jm;j+) d=d+aij*xj; xi

11、=(aim-d)/aii; /*求解*/ cout该方程组的解为:endl; for(i=0;im;i+) coutxi=xit; return 0; void load() int i,j; for(i=0;im;i+) for(j=0;jaij; LU分解法：#include void solve(float l100,float u100,float b,float x,int n) int i,j; float t,s1,s2; float y100; for(i=1;i=n;i+) /* 第一次回代过程开始 */ s1=0; for(j=1;j=1;i-) /* 第二次回代过程开始*

12、/ s2=0; for(j=n;ji;j-) t=-uij; s2=s2+t*xj; xi=(yi+s2)/uii; void main() float a100100,l100100,u100100,x100,b100; int i,j,n,r,k; float s1,s2; for(i=1;i=99;i+)/*将所有的数组置零，同时将L矩阵的对角值设为1*/ for(j=1;j=99;j+) lij=0,uij=0; if(j=i) lij=1; printf (输入方程组的个数 n:n);/*输入方程组的个数*/ scanf(%d,&n); printf (读取原矩阵 A（x的系数）:n

13、);/*读取原矩阵A*/ for(i=1;i=n;i+) for(j=1;j=n;j+) scanf(%f,&aij); printf (读取列矩阵 B（y的值）:n);/*读取列矩阵B*/ for(i=1;i=n;i+) scanf(%f,&bi); for(r=1;r=n;r+)/*求解矩阵L和U*/ for(i=r;i=n;i+) s1=0; for(k=1;k=r-1;k+) s1=s1+lrk*uki; uri=ari-s1; for(i=r+1;i=n;i+) s2=0; for(k=1;k=r-1;k+) s2=s2+lik*ukr; lir=(air-s2)/urr; prin

14、tf(输出矩阵 L:n);/输出矩阵Lfor(i=1;i=n;i+) for(j=1;j=n;j+) printf(%7.3f ,lij); printf(n); printf(输出矩阵 U:n);/输出矩阵U for(i=1;i=n;i+) for(j=1;j=n;j+) printf(%7.3f ,uij); printf(n); solve(l,u,b,x,n); printf(解为:n); for(i=1;i=n;i+) printf(x%d=%fn,i,xi); 追赶法： #include #define N 3 main() double A33,b3; printf(请按顺序输入

15、x的系数:n); int a,c; for(a=0;a3;a+) for(c=0;c3;c+) scanf(%lf,&Aac); printf(请按顺序输入y的值:n); int k; for(k=0;k3;k+)scanf(%lf,&bk); int i; A01=A01/A00; for(i=1;i2;i+) Aii+1=Aii+1/(Aii-Aii-1*Ai-1i); for(i=1;i3;i+) Aii=Aii-Aii-1*Ai-1i; b0=b0/A00; for(i=1;i=0;i-) bi=bi-Aii+1*bi+1; for(i=0;i3;i+) printf(x%d=%.6l

16、fn,i,di); 六、结果讨论和分析Gauss消元法： 列主元素消元法： LU分解法： 追赶法： 分析讨论从消元过程可以看出，对于n阶线性方程组，只要各步主元素不为零，经过n-1步消元，就可以得到一个等价的系数矩阵为上三角形阵的方程组，然后再利用回代过程可求得原方程组的解. 由于列主元素法相似且优于完全主元素法 所以省略了后者。消元过程相当于分解 A为单位下三角阵L与上三角阵U的乘积，解方程组Ly=b回代过程就是解方程组Ux=y。其中的L为n阶单位下三角阵、U为上三角阵. 在 A 的LU 分解中, L取下三角阵, U 取单位上三角阵,这样求解方程组Ax=d 的方法称为追赶法。另外是追赶法和其

17、他方法求同一方程结果不一样，我多次修改源程序，也不知道原因。再就是追赶法有很大的局限性 还待改良。流程图：实验地点综合楼五层506室指导教师王峥太原理工大学学生实验报告学院名称计算机科学与技术专业班级计科 1002学号2010001414学生姓名 xxx实验日期2012.6.5成绩课程名称计算数值方法实验题目实验三 线性方程组的迭代解法一、课题名称线性方程组的迭代解法使用雅可比迭代法或高斯-赛德尔迭代法对下列方程组进行求解。 二、目的和意义学习使用雅可比迭代法或高斯-赛德尔迭代法三、计算公式雅克比迭代法：设线性方程组Ax=b的系数矩阵A可逆且主对角元素a11,a22,ann均不为零，令D=di

18、ag(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。四、主要仪器设备 Vc+ 9.0 C-free CodeBlocks五、结构程序设计雅克比迭代法：#include #include main() int i; double x120 ,x220,x320; double x10, x20, x30; printf(请输入x1,x2,x3的初值：n);scanf(%lf%lf%lf,&x10,&x20, &x30); printf( n x1n x2n x3

19、n n); for(i=0;i18;i+) x10=x10; x20=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); 六、实验结果与分析：雅克比迭代法：分析讨论： 其实，这两个迭代法是之前迭代法的升级，多了几个迭代式子而已，而且两者相差不大比较简单，所以选择了雅克比迭代法进行求解，但是没有与另一种方法 高斯赛德尔迭代法进行实质性的比较。流程图：实验地点综合楼五层506室

20、指导教师王峥太原理工大学学生实验报告学院名称计算机科学与技术专业班级计科 1002学号2010001414学生姓名 xxx实验日期2012.6.10成绩课程名称计算数值方法实验题目实验四 矩阵特征值与特征向量问题一、课题名称使用幂法求A模为最大的特征值及其相应的特征向量。二、目的和意义（1）了解矩阵特征值与特征向量问题解法，掌握幂法。（2）加深对矩阵特征值与特征向量问题求解方法的认识，掌握算法。三、计算公式幂法：由已知的非零向量x0和矩阵A的乘幂构造向量序列xn以计算矩阵A的按模最大特征值及其特征向量的方法，称为幂法。迭代公式：结果可取 四、主要仪器设备Vc+ 9.0 C-free CodeB

21、locks五、结构程序设计源代码： #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;ifabs(Max)Max=xi; return Max; void PowerMethod(float *A) float UN,VN,r1,r2,temp; int i,j,k=0; while(kKM) k+; for(i=0;iN;i+) temp=0; for(j=0;jN;j+)temp+=*(A+i*N+j)*Uj;

22、 Vi=temp; for(i=0;iN;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=%fn,r2); for(i=0;iN;i+)printf(y%d=%fn,i+1,Ui);void main() float ANN=2,-1,0,-1,2,-1,0,-1,2 ;float UN; /A的值 U0=1; U1=1; U2=1;/x0的值 PowerMethod(A0);六、结果讨论和分析分析讨论由于该程序将A矩

23、阵和x0的值编写在程序中，所以要想修改成其他矩阵时比较麻烦，所以也有一定的局限性。 幂法是一种求任意矩阵A的按模最大特征值及其对应特征向量的迭代算法。该方法的最大优点是计算简单，容易在计算机上实现，对稀疏矩阵较为适合，但有时收敛速度很慢。流程图：实验地点综合楼五层506室指导教师王峥太原理工大学学生实验报告学院名称计算机科学与技术专业班级计科 1002学号2010001414学生姓名 xxx实验日期2012.6.10成绩课程名称计算数值方法实验题目 实验五 代数插值一、课题名称（使用拉格朗日插值法或牛顿插值法求解：已知f(x)在6个点的函数值如下表所示，运用插值方法，求f(0.596)的近似值

24、。x0.400.550.650.800.901.05f(x)0.410750.578150.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

25、,1,n再由插值多项式的唯一性，得Pn(x)Ln(x)四、主要仪器设备Vc+ 9.0 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; /coutn; /for(int i=0;i=n;i+) / /*cout请输入xi的值：xi; cout请输入yi的值：yi; Mi0=xi; Mi1=yi; */用二维保存所有数据 M00=0.40;M01=0.41075;M10=0

26、.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( 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+1endl; coutxx; for(int i=0;in;i+) X*=xx-Mi0; N+=Mi+1i+2*X; P=M01+N; cout其函数值：y=Pendl; /

27、 六、结果讨论和分析分析讨论拉格朗日插值的优点是插值多项式特别容易建立，缺点是增加节点是原有多项式不能利用，必须重新建立，即所有基函数都要重新计算，这就造成计算量的浪费。所以该程序 选择了牛顿法，由于在输入数据的时候比较麻烦 所以将数据存在了程序里面，只能在程序里修改数据。 另外由于精度问题计算结果有一定的误差。流程图：实验地点综合楼五层506室指导教师王峥太原理工大学学生实验报告学院名称计算机科学与技术专业班级计科 1002学号2010001414学生姓名 xxx实验日期2012.6.10成绩课程名称计算数值方法实验题目 实验六 最小二乘法拟合多项式一、课题名称给定数据点（xi ,yi），用

28、最小二乘法拟合数据的多项式，并求平方误差。xi00.50.60.70.80.91.0yi11.751.962.192.442.713.00二、目的和意义1熟练运用已学计算方法求解方程组2加深对计算方法技巧，选择正确的计算方法来求解各种方程组3培养使用电子计算机进行科学计算和解决问题的能力三、计算公式建立正规方程组：（xij+k）ak=xijyi ,j=0,1,n 平方误差：I=（akxik-yi）2四、主要仪器设备Vc+ 9.0 C-free CodeBlocks五、结构程序设计源代码： #include#include#define N 15double power(double &a,in

29、t 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,k,n,index;/coutn;n=7; coutendl;/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

30、=2.71;X6=1.0;Y6=3.00; /绑定数据 /可以解绑由下for循环 输入任何数据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+)s

31、umYi+=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;

32、i+)coutxit;coutendl;cout拟合曲线方程为:endl; couty(x)=x1;for(int i=1;iindex+1;i+)cout + xi;for(int j=0;ji;j+)cout*x;coutendl;double m=0;coutendl平方误差为:;for(i=0;in;i+) double t=x1+x2*Xi-Yi;m=m+power(t,2);coutmendlendl;cout特别说明：为了方便endl;cout 数据被存在程序中无需endl; cout 输入修改在代码里进行！endl;六、结果讨论和分析分析讨论 数据拟合的具体作法是：对给定的数据（xi ,yi）（i=0,1,m），在取定的函数类中，求p(x)属于此函数类，使误差ri=p(xi)- yi (i=0,1,m)的平方和最小，即ri2=（p(xi)-yi）2=min从几何意义上讲，就是寻求与给定点（xi ,yi）（i=0,1,m）的距离平方和为最小的曲线y=p(x)。说明：上述程序中 数据已被事先存入，可在源代码中任意修改想要的数据。流程图：实验地点综合楼五层506室指导教师王峥