数值计算方法上机实验报告.docx
《数值计算方法上机实验报告.docx》由会员分享,可在线阅读,更多相关《数值计算方法上机实验报告.docx(27页珍藏版)》请在三一办公上搜索。
1、数值计算方法上机实验报告数值计算方法上机实验报告 数值计算方法上机实验报告 实验目的:复习和巩固数值计算方法的基本数学模型,全面掌握运用计算机进行数值计算的具体过程及相关问题。利用计算机语言独立编写、调试数值计算方法程序,培养学生利用计算机和所学理论知识分析解决实际问题的能力。 上机练习任务:利用计算机基本C语言编写并调试一系列数值方法计算通用程序,并能正确计算给定题目,掌握调试技能。 掌握文件使用编程技能,如文件的各类操作,数据格式设计、通用程序运行过程中文件输入输出运行方式设计等。 一、 各算法的算法原理及计算机程序框图 1. 列主元高斯消去法 l 算法原理: 高斯消去法是利用现行方程组初
2、等变换中的一种变换,即用一个不为零的数乘一个方程后加只另一个方程,使方程组变成同解的上三角方程组,然后再自下而上对上三角方程组求解。 列选住院是当高斯消元到第k步时,从k列的akk以下的各元素中选出绝对值最大的,然后通过行交换将其交换到akk的位置上。交换系数矩阵中的两行,只相当于两个方程的位置交换了,因此,列选主元不影响求解的结果。 1 数值计算方法上机实验报告 计算机程序框图如上 源程序: #define N 200 #include stdio.h #include math.h FILE *fp1,*fp2; void LZ int n,i,j,k=0,l; double d,t,t1
3、; static double xN,aNN; fp1=fopen(a1.txt,r); fp2=fopen(b1.txt,w); fscanf(fp1,%d,&n); for(i=0;in;+i) for(j=0;jfabs(d) /*选主元*/ d=aik; l=i; i+; while(in); if(d=0) printf(n输入矩阵有误!n); else /* if(l!=k) for(j=k;j=n;j+) t=alj; alj=akj; akj=t; for(j=k+1;j=n;j+) /* akj/=akk; for(i=k+1;in;i+) for(j=k+1;j=n;j+)
4、 aij-=aik*akj; k+; while(k=0;i-) /* t1=0; for(j=i+1;jn;j+) t1+=aij*xj; xi=ain-t1; 3 换行*/ 正消*/ 回代*/ 数值计算方法上机实验报告 for(i=0;in;i+) fprintf(fp2,n方程组的根为x%d=%lf,i+1,xi); fclose(fp1); fclose(fp2); main LZ; l 具体算例及求解结果: 用列选主元法求解下列线性方程组 x1+2x2-3x3=82x1+x2+3x3=22 3x+2x+x=28231输入3 输出结果:方程组的根为x1=6.000000 1 2 -3
5、8 方程组的根为x2=4.000000 2 1 3 22 方程组的根为x3=2.000000 3 2 1 28 l 输入变量、输出变量说明: 输入变量:aij系数矩阵元素,bi常向量元素 输出变量:b1,b2,Lbn解向量元素 2. 杜里特尔分解法解线性方程 l 算法原理: 求解线性方程组Ax=b时,当对A进行杜里特尔分解,则等价于求解LUx=b,这时可归结为利用递推计算相继求解两个三角形方程组,用顺代,由 Ly=b 求出y,再利用回带,由Ux=y求出x。 计算机程序框图: 源程序:#include stdio.h #include math.h FILE *fp1,*fp2; void ma
6、in int i,j,k,N; 4 数值计算方法上机实验报告 double s,A200200,B200,x200,y200; static double L200200,U200200; fp1=fopen(a2.txt,r); fp2=fopen(b2.txt,w); fscanf(fp1,%d,&N); for(i=0;iN;i+) for(j=0;jN;j+) fscanf(fp1,%lf,&Aij); for(i=0;iN;i+) fscanf(fp1,%lf,&Bi); for(i=0;iN;i+) /*LU分解*/ for(j=i;jN;j+) s=0.0; for(k=0;ki
7、;k+) s+=Lik*Ukj; Uij=Aij-s; for(j=i+1;jN;j+) s=0.0; for(k=0;ki;k+) s+=Uki*Ljk; Lji=(Aji-s)/Uii; for(i=0;iN;i+) for(j=0;jN;j+) Lii=1; fprintf(fp2,nU矩阵为:); for(i=0;iN;i+) fprintf(fp2,n); for(j=0;jN;j+) fprintf(fp2,%10.3f,Uij); fprintf(fp2,nL矩阵为:); for(i=0;iN;i+) fprintf(fp2,n); for(j=0;jN;j+) fprintf(
8、fp2,%10.3f,Lij); 5 数值计算方法上机实验报告 for(i=0;iN;i+) s=0.0; for(k=0;k=0;i-) s=0.0; for(k=i+1;kN;k+) s+=Uik*xk; xi=(yi-s)/Uii; fprintf(fp2,n方程组解为:); for(i=0;iN;i+) fprintf(fp2,nx%d=%10.3f,i+1,xi); fclose(fp1); fclose(fp2); l 具体算例及求解结果: 用杜里特尔分解法求解方程组 2x1+3x2+4x3=393x1-2x2+2x3=14 4x1+2x2+3x3=43输入数据 输出结果:U矩阵为
9、: 2.000 3.000 4.000 0.000 -6.500 -4.000 0.000 0.000 -2.538 3 L矩阵为: 2 3 4 1.000 0.000 0.000 3 -2 2 1.500 1.000 0.000 4 2 3 2.000 0.615 1.000 39 14 43 方程组解为: x1= 6.000 x2= 5.000 x3= 3.000 l 输入变量、输出变量说明: 输入变量:aij系数矩阵元素,bi常向量元素 输出变量:b1,b2,Lbn解向量元素 6 数值计算方法上机实验报告 3. 拉格朗日插值法 l 算法原理: 首先构造基函数lk(x)=i=0iknx-x
10、i,可以证明基函数满足下列条件: xk-xiiki=k0lk(xi)=1, 对于给定n+1个节点,n次拉格朗日插值多项式由下式给出: L(x)=lk(x)yk k=0n其中 lk(x)=i=0iknx-xi xk-xi由于lk(x)是一个关于x的n次多项式,所以L(x)为关于x的不高于n次的代数多项式。当x=xi时,L(xi)=yi,满足插值条件。 l 计算机程序框图: 7 数值计算方法上机实验报告 源程序:#includestdio.h #includemath.h int n,m,i,j; float x2,x3,z1=0.0,z=0.0,t,x50,y50,c50,A50; main F
11、ILE *fp1,*fp2; fp1=fopen(a3.txt,r); fp2=fopen(b3.txt,w); fscanf(fp1,%d,&n); for(i=0;in;i+) fscanf(fp1,%f,%f,&xi,&yi); fscanf(fp1,%d,&m); fscanf(fp1,%f,&x2); fscanf(fp1,%f,&x3); for(i=0;in;i+) /*选m个最接近的点*/ ci=fabs(xi-x2); for(i=0;in;i+) for(j=i+1;jcj) t=ci; ci=cj; cj=t; 8 数值计算方法上机实验报告 t=xi; xi=xj; xj
12、=t; t=yi; yi=yj; yj=t; for(i=0;im;i+) /*求值*/ Ai=1.0; for(j=0;jm;j+) if(i!=j) Ai=Ai*(x2-xj)/(xi-xj); z=z+Ai*yi; for(i=0;im;i+) /*求值*/ Ai=1.0; for(j=0;jm;j+) if(i!=j) Ai=Ai*(x3-xj)/(xi-xj); z1=z1+Ai*yi; fprintf(fp2,nx=%10.7f处的函数值为:y=%10.7f,x2,z); fprintf(fp2,nx=%10.7f处的函数值为:y=%10.7f,x3,z1); fclose(fp1
13、); fclose(fp2); 具体算例及求解结果: 对于一组数据表进行二次数值插值编程,根据下面数值表计算f(0.49) 和f(0.51) x f(x) 0.2 16 0.4 20 0.6 15 0.8 10 输入数据: 输出结果: x= 0.4900000处的函数值为:y=18.8637505 4 x= 0.5100000处的函数值为:y=18.3637524 0.2,16 0.4,20 0.6,15 0.8,10 3 0.49 0.51 l 输入变量、输出变量说明: 输入变量:(xi,yi)插值节点 输出变量:y插值所得到被插函数在插值点的近似值 4. 曲线拟合 9 数值计算方法上机实验
14、报告 l 算法原理: 对于给定的一组数据(xi,yi),i=1,2,m,寻求做n次多项式 y=akxk k=0n使性能指标 J(a0,a1,L,an)=(yi-akxik)2为最小。 i=1k=0mn由于性能指标J可以被看做关于ak,k=0,1,n的多元函数,故上述拟合多项式的构造问题可转化为多元函数的极值问题。令 J=0 ak从而有正则方程组 mxixixi2MMxnn+1xiixxM23iLLLLixxxn+2iayi0n+1xyai1=ii MMManxinyixi2nni求解即得多项式系数。 l 计算机程序框图: l 源程序: #include #include main int i,
15、j,k,m,n,l,N,t,t1; double max,A5050,x50,y50,S50,T50,X50;float yb=0.0,xb,a1,a2,a0; FILE *fp1,*fp2; fp1=fopen(a4.txt,r); fp2=fopen(b4.txt,w); fscanf(fp1,%d %dn,&n,&m); for(i=0;in;i+) fscanf(fp1,%lf %lf,&xi,&yi); fscanf(fp1,%f,&xb); for(i=0;i=2*m;i+) Si=0.0; for(j=0;jn;j+) 10 数值计算方法上机实验报告 Si+=pow(xj,i);
16、 for(i=0;i=m;i+) Ti=0.0; for(j=0;jn;j+) Ti+=yj*pow(xj,i); N=m+1; for(i=0;iN;i+) for(j=0;jN;j+) l=i+j; Aij=Sl; AiN=Ti; for(i=0;iN-1;i+) max=fabs(Aii); /*选主*/ for(j=i+1;jmax) max=fabs(Aji); m=j; if(m!=i) for(k=0;k=N;k+) t=Aik; Aik=Amk; Amk=t; for(j=i+1;j=0;k-) Ajk=Ajk-Aik*Aji/Aii; for(i=N-1;i=0;i-) /*
17、回代*/ for(j=i-1;j=0;j-) for(k=N;k=0;k-) Ajk=Ajk-Aik*Aji/Aii; fprintf(fp2,n解为:); /*输出结果*/ for(i=0;iN;i+) fprintf(fp2,na%d=%10.7lf,i,AiN/Aii); fprintf(fp2,n拟合多项式为:n); fprintf(fp2, P(x)=%10.7lf,A0N/A00); for(i=1;iN;i+) fprintf(fp2,+(%10.7lf)x%d,AiN/Aii,i); 11 数值计算方法上机实验报告 for(i=0;iN;i+) yb+=(AiN/Aii)*po
18、w(xb,i); fprintf(fp2,nP(%f)=%10.7f,xb,yb); fclose(fp1); fclose(fp2); l 具体算例及求解结果: 对于一组数据表进行二次多项式曲线拟合,根据以下数据胡二次拟合曲线 求y(5) xi yi 1 1.6 2 2.8 3 3.6 4 4.9 5 5.4 6 6.8 7 8 9 10 7.9 9.2 10.2 11.4 试用最小二乘法求二次拟合多项式 输入数据: 10 2 输出结果: 1 1.6 解为: 2 2.8 a0=-0.3012450 3 3.6 a1= 1.3167338 4 4.9 a2=-0.0166439 5 5.4 拟
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 数值 计算方法 上机 实验 报告
链接地址:https://www.31ppt.com/p-3558779.html