计算方法实验课程.docx
- 文档编号:3285592
- 上传时间:2022-11-21
- 格式:DOCX
- 页数:33
- 大小:224.18KB
计算方法实验课程.docx
《计算方法实验课程.docx》由会员分享,可在线阅读,更多相关《计算方法实验课程.docx(33页珍藏版)》请在冰豆网上搜索。
计算方法实验课程
《计算方法》实验课程指导书
徐中宇
2012年3月
目录
第一章插值法1
§1拉格朗日插值1
§2牛顿插值5
第二章线性方程组的解法8
§1高斯消去法8
§2 列主元消去法12
§3 线性方程组的迭代解法16
第三章方程求根21
§1二分法21
§2牛顿迭代法24
§3埃特金(Aitken)迭代法26
第一章插值法
目的与要求:
1.掌握不同的输入、输出语句,注意节约内存方法。
2.熟悉拉格朗日插值和牛顿插值公式,并体会它们不同的特点。
§1拉格朗日插值
1.方法概要:
拉格朗日n次插值多项式
,其中
可用双重循环来实现;内循环为连乘;外循环为连加。
用数组表下标变量,给出
。
2.程序流程图
图1-1拉格朗日插值法程序流程图
3.程序及例
例1.已知函数表:
xi0.561600.562800.564010.56521
yi0.827410.826590.825770.82495
用三次拉格朗日插值多项式求x=0.5635时的函数值。
程序为
#include"stdafx.h"
#include
constintn=3;
constdoublexx=0.5635;
voidmain()
{
doublex[n+1]={0.56160,0.56280,0.56401,0.56521};
doubley[n+1]={0.82741,0.82659,0.82577,0.82495};
doubleL=0,Li=1;
for(inti=0;i<=n;i++)
{
Li=1;
for(intj=0;j<=n;j++)
{
if(j!
=i)
{
Li=Li*(xx-x[j])/(x[i]-x[j]);
}
}
L+=Li*y[i];
}
cout<<"Theresultis:
"< } 例2.已知函数表,用拉格朗日插值多项式求0.5,0.7,0.85三点处的函数值。 程序的结构: 在程序的说明部分将已知数据分别存入一维数组x和y。 在主函数main中提示人工键入插值点的数目m,继而在m次循环计算中分别提示键入插值点,计算插值,并输出结果。 程序中的主要常量,变量及函数: n为插值基点的最大下标。 x[n+1],y[n+1]分别存放插值基点及其函数值。 M为插值点个数,X存放插值点的值。 程序清单: #include"stdafx.h" #include #definen4/*插值基点的最大下标*/ voidmain() { doublex[n+1]={0.4,0.55,0.8,0.9,1};/*插值基点的值*/ doubley[n+1]={0.41075,0.57815,0.88811,1.02652,1.17520}; intm,k,i,j; floatX,L,P; printf("\nPleaseenterm="); scanf("%d",&m);/*键入插值点的个数*/ for(k=1;k<=m;k++) { printf("\nPleaseenterX%d=",k); scanf("%f",&X); P=0.0; for(i=0;i<=n;i++) { L=1.0; for(j=0;j<=n;j++) { if(j! =i) { L=L*(X-x[j])/(x[i]-x[j]); } } P=P+y[i]*L; } printf("P(%f)=%f\n",X,P); } } 计算结果: (当提示“Pleaseenterm=”时,键入插值点个数3,再分别根据提示X1=或X2=,X3=依次键入插值点0.5,0.7,0.85),输出的结果分别为 P(0.500000)=0.521090P(0.700000)=0.758589P(0.850000)=0.956119 §2牛顿插值 1.方法概要: 牛顿n次插值多项式 Pn(x)=Pn-1(x)+an(x-x0)(x-x1)…(x-xn-1) 在计算机上计算均差表时,可用一维数祖x存放基点值,用二维数组F存放和计算由个阶均差构成的下三角矩阵。 计算公式如下: Fi0=f(xi),i=0,1,…n Fij=(Fi,j-1-Fi-1,j-1)/(xi-xi-j),j=1,2,…ni=j,j+1,…n 2.程序及例 例3.已知函数表同例2,试根据上面公式的算法计算均差表存入矩阵F,并用牛顿均差插值多项式求解上题。 注意到 ,那么 计算 的流程为: (1) 。 (2)对 ,计算 。 程序的结构: 在程序开始时将初始数据分别存入一维数组x,y。 主函数main首先调用自定义函数junca计算均差表,再提示人工键入插值点个数m,继而在m次循环中分别提示键入插值点X,并计算相应的插值,输出值结果。 程序中的主要常量,变量及函数: n,x[],y[],m,X均同上例。 Junca(doublex[],doubley[],doubleF[][n+1])是自定义函数,计算均差表并打印输出,结果存于F中,P用于计算牛顿均差插值多项式的值。 #include"stdafx.h" #include constintn=4; /*函数junca实现均差表的计算,将均差表存入数组F,并在屏幕上打印输出*/ voidjunca(floatx[],floaty[],floatF[][n+1]) { inti,j; for(i=0;i<=n;i++) { F[i][0]=y[i]; } for(j=1;j<=n;j++) { for(i=j;i<=n;i++) { F[i][j]=(F[i][j-1]-F[i-1][j-1])/(x[i]-x[i-j]); } } printf("\n%12s%12s\n","Xi","F(Xi)"); for(j=1;j<=38;j++) printf("--"); printf("\n"); for(i=0;i<=n;i++) { printf("%12f",x[i]); for(j=0;j<=i;j++) printf("%12f",F[i][j]); printf("\n"); } for(j=1;j<=38;j++) printf("--"); printf("\n"); } voidmain() { floatx[n+1]={0.4,0.55,0.8,0.9,1};/*插值基点的值*/ floaty[n+1]={0.41075,0.57815,0.88811,1.02652,1.17520}; inti,m,k; floatp,t,F[n+1][n+1]; floatX; junca(x,y,F);/*调用函数junca计算均差表*/ printf("\nPleaseenterm="); scanf("%d",&m);/*键入插值点的个数*/ for(i=1;i<=m;i++) { scanf("%f",&X);/*键入值点*/ p=F[n][n]; for(k=n-1;k>=0;k--) p=p*(X-x[k])+F[k][k]; /*计算牛顿均差值多项式*/ printf("P(%f)=%f\n",X,p); } } 计算结果: Xi F(Xi) 1阶 2阶 3阶 4阶 0.400000 0.410750 0.550000 0.578150 1.116000 0.800000 0.888110 1.239840 0.309600 0.900000 1.026520 1.384100 0.412171 0.205143 1.000000 1.175200 1.486800 0.513500 0.225175 0.033386 (当提示Pleaseenterm=时,键入插值点个数3,再分别根据提示X1=或X2,X3=依次键入插值点0.5,0.7,0.85)输出的结果为 P[0.500000]=0.521090,P[0.700000]=0.758589,P[0.850000]=0.956119 练习题 1.已知函数表: xi0.40.550.650.80.9 yi0.410750.578150.696750.888111.02652 分别用四次拉格朗日插值和牛顿插值,求x=0.596时的函数值。 第二章线性方程组的解法 本章的目的与要求: 1.熟悉求解线性方程组的有关理论与方法。 2.会编制列主元消去法、LU分解法与两种迭代法的计算程序,上机计算求得结果, 3.通过上机实习对计算数据的分析,进一步了解各种方法的功能,优缺点与适用范围,体会如何针对不同问题适当选择方法。 §1高斯消去法 1.方法概要 解线性方程组 (1) 为统一计算公式与简化计算程序,把方程组的增广矩阵记作A (1)=(A,b) (2) 解题分为消元与回代两个过程: 消元过程; 第1步,对第1列作消元,将 消为零。 第2步,对第2列作消元,将 消为零。 设 ≠0(k=1,2,3,4….,n-1)。 一般地,对于第k步消元,采用如下公式: (3) (4) (j=k,k+1,…,n+1) 直到作完n-1步消元,相应的方程组化为与原来方程组等价的方程组: (6) 回代过程: 回代过程是解方程组 ,采用如下公式 (7) 2.程序流程图(见图2—1) 流程图说明 (1)首先输入方程组的阶数n及方程组的增广矩阵A。 A为一个n×(n+1)的二维数组: 所有运算的中间结果与最后结果全部放A中相应的位置。 (2)本程序为了教学的需要,先屏幕输出原始增广矩阵,然后逐次输出每次消元后的矩阵以及每次所用的消元因子,最后输出解向量。 我们可以从中看出消元过程的具体情况。 在实用中无此必要,可以省略。 (3)这是顺序取主元的消去法。 为了检查 是否为零,设置了一个检查口,一旦遇到 ,便使其输出一个记号,以示顺序取主元的方法失败. 3.程序及例 例1.求解方程组 2.5x1+2.3x2-5.1x3=3.7 5.3x1+9.6x2+1.5x3=3.8 8.lx1+1.7x2-4.3x3=5.5 (答: x=(0.406705,0.237818,-0.419325)T) #include"stdafx.h" #include #include #definen3/*方程组的阶数*/ doubleaa[n][n+1]={{2.5,2.3,-5.1,3.7},{5.3,9.6,1.5,3.8},{8.1,1.7,-4.3,5.5}}; /*增广矩阵的初始数据*/ intgauss(doublea[][n+2],doublex[]) { inti,j,k; floatc; for(k=1;k<=n-1;k++)/*消元过程*/ { for(i=k+1;i<=n;i++) { if(fabs(a[k][k])<1e-12) { printf("\ndet=0.fail! \n"); return(0); } c=a[i][k]/a[k][k]; for(j=k+1;j<=n+1;j++) { a[i][j]=a[i][j]-c*a[k][j]; } } } for(k=n;k>=1;k--)/*回代过程*/ { x[k]=a[k][n+1]; for(j=k+1;j<=n;j++) { x[k]=x[k]-a[k][j]*x[j]; } x[k]=x[k]/a[k][k]; } return (1); } voidmain() { inti,j,det; doublea[n+1][n+2],x[n+1]; for(i=1;i<=n;i++) for(j=1;j<=n+1;j++)/*用a[1][1]~a[n][n+1]存放增广矩阵*/ a[i][j]=aa[i-1][j-1]; det=gauss(a,x);/*调用函数gauss求解方程组,并获取返回标志值*/ if(det! =0) { for(i=1;i<=n;i++) printf("\nx[%d]=%f",i,x[i]); printf("\n--------------------------\n"); } } 是 图2-1Gauss消去法程序流程图 §2 列主元消去法 1.方法概要 作第K步消元时,从A(K)的子矩阵(见§1中式(3))的第一列中选取 |,以 作为主元素,第 行称为主行。 将主行与第k行交换,然后再作这一步的消元。 2.计算框图 图2-2列主元消去法框图 3.程序及例 例2.用列主元消去法求解 程序的结构: 主函数main调用函数gaussl进行列主元高斯消去法计算,根据返回的函数值为1或0而输出计算结果或失败信息。 程序中的主要常量、变量及函数: n为常量,指定方程组的阶数。 数组aa[n][n+1]存放方程组增广矩阵的原始数据并保留不变。 数组a[n+1][n+2]、x[n+1]分别用于计算增广和方程组的解。 intgauss1(a,x)是自定义函数,用列主元高斯消去法求解n阶线性方程,当系数矩阵奇异使求解失败时,输出信息“det=0.fail! ”并返回函数值0;否则求解成功,返回函数值1。 Intdet在主函数里接受gaussl的返回值,当det不为0时主函数打印出方程组的解。 #include"stdafx.h" #include #include #definen3/*方程组的阶数*/ staticdoubleaa[n][n+1]={{1,2,-1,3},{1,-1,5,0},{4,1,-2,2}};/*增广矩阵的初始数据*/ intgauss1(doublea[][n+2],doublex[]) { inti,j,k,r; doublec; for(k=1;k<=n-1;k++)/*消元过程*/ { r=k; for(i=k;i<=n;i++)/*选列主元*/ { if(fabs(a[i][k])>fabs(a[r][k])) { r=i; } if(fabs(a[r][k])<1e-12) { printf("\ndet=0.fail! \n"); return(0); } if(r! =k) for(j=k;j<=n+1;j++)/*交换k、r两行*/ { c=a[k][j]; a[k][j]=a[r][j]; a[r][j]=c; } } for(i=k+1;i<=n;i++)/*进行消元计算*/ { c=a[i][k]/a[k][k]; for(j=k+1;j<=n+1;j++) { a[i][j]=a[i][j]-c*a[k][j]; } } } if(fabs(a[n][n])<1e-12) { printf("|ndet=0.fail! \n"); return(0); } for(k=n;k>=1;k--)/*回代过程*/ { x[k]=a[k][n+1]; for(j=k+1;j<=n;j++) { x[k]=x[k]-a[k][j]*x[j]; } x[k]=x[k]/a[k][k]; } return (1); } voidmain() { inti,j,det; doublea[n+1][n+2],x[n+1]; for(i=1;i<=n;i++) for(j=1;j<=n+1;j++)/*用a[1][1]~a[n][n+1]存放增广矩阵*/ a[i][j]=aa[i-1][j-1]; det=gauss1(a,x);/*调用函数gauss1求解方程组,并获取返回标志值*/ if(det! =0) { for(i=1;i<=n;i++) { printf("\nx[%d]=%f",i,x[i]); } printf("\n--------------------------\n"); } } 计算结果: x1=0.250000x2=1.500000x3=0.250000 §3 线性方程组的迭代解法 1.雅可比迭代法 例3.根据雅可比迭代法编制程序求解方程组 。 其中 的计算流程为: (1) ; (2)对 做: ① ②若 则 。 程序中的主要常量,变量: n,eps是常量,表示方程组的阶数及指定的精度。 aa[n][n],bb[n]分别存放方程组的系数矩阵及常数项。 intNO是最大迭代次数,当迭代达到NO次以上时程序输出失败信息,并结束。 x[n+1],y[n+1]分别用于存放相继两次的迭代向量。 用变量max计算。 程序清单: //Definestheentrypointfortheconsoleapplication. // #include"stdafx.h" #include #include #definen3/*方程组的阶数*/ #defineeps0.5e-4/*给定精度要求*/ staticdoubleaa[n][n]={{10,-1,-2},{-1,10,-2},{-1,-1,5}}; staticdoublebb[n]={7.2,8.3,4.2}; voidmain() { intk=0,NO,i,j; doubled,s,max; doublea[n+1][n+1],b[n+1],x[n+1],y[n+1]; for(i=1;i<=n;i++) { for(j=1;j<=n;j++) a[i][j]=aa[i-1][j-1]; b[i]=bb[i-1]; } printf("\nPleaseenterNo: "); scanf("%d",&NO);/*键入最大迭代次数NO*/ for(i=1;i<=n;i++) y[i]=0; do{ k++; for(i=1;i<=n;i++) x[i]=y[i]; if(k>=NO) { printf("\nFail! "); break; }/*当k>=NO时输出失败信息,结束迭代*/ max=0.0; for(i=1;i<=n;i++) { s=0.0; for(j=1;j<=n;j++) if(j! =i) { s=s+a[i][j]*x[j]; } y[i]=(b[i]-s)/a[i][i]; d=fabs(y[i]-x[i]); if(max } }while(max>=eps); if(max { printf("k=%d\n",k); for(i=1;i<=n;i++) printf("x[%d]=%f\n",i,x[i]); } } 运行结果: 产生输入提示PleaseenterNO: 后,键入适当大的正整数,可行到结果数据为 k=11 x[1]=1.099979 x[2]=1.199979 x[3]=1.299975 2.Gauss-Seidel迭代法 例4.根据高斯-赛德尔迭代法编制程序求解上题的方程组。 程序中的主要常量及变量: n,eps,aa[][],bb[],a[][],b[],x[],NO均同上例。 每次迭代计算时,用变量s临时存放x[i]的旧值。 用d计算x[i]与其旧值之差最大绝对值。 程序清单: //Definestheentrypointfortheconsoleapplication. #include"stdafx.h" #include #include #definen3 #defineeps0.5e-4 staticdoubleaa[n][n]={{10,-1,-2},{-1,10,-2},{-1,-1,5}}; staticdoublebb[n]={7.2,8.3,4.2}; voidmain() { intk=0,NO,i,j; doubled,sum,s; doublea[n+1][n+1],b[n+1],x[n+1]; for(i=1;i<=n;i++) { for(j=1;j<=n;j++) a[i][j]=aa[i-1][j-1]; b[i]=bb[i-1]; } printf("\nPleaseenterNO: "); scanf("%d",&NO); for(i=1;i<=n;i++) { x[i]=0; } do{ k++; d=0.0; if(k>=NO) { printf("\nFail! "); break; } for(i=1;i<=n;i++) { sum=0.0; s=x[i];/*s临时存放x[i]的旧值*/ for(j=1;j<=n;j++) if(j! =i) { sum=sum+a[i][j]*x[j]; } x[i]=(b[i]-sum)/a[i][i];/*计算x[i]的新值*
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 计算方法 实验 课程