数值计算方法上机实验报告.docx
- 文档编号:7878049
- 上传时间:2023-01-26
- 格式:DOCX
- 页数:28
- 大小:264.84KB
数值计算方法上机实验报告.docx
《数值计算方法上机实验报告.docx》由会员分享,可在线阅读,更多相关《数值计算方法上机实验报告.docx(28页珍藏版)》请在冰豆网上搜索。
数值计算方法上机实验报告
数值计算方法上机实验报告
实验目的:
复习和巩固数值计算方法的基本数学模型,全面掌握运用计算机进行数值计算的具体过程及相关问题。
利用计算机语言独立编写、调试数值计算方法程序,培养学生利用计算机和所学理论知识分析解决实际问题的能力。
上机练习任务:
利用计算机基本C语言编写并调试一系列数值方法计算通用程序,并能正确计算给定题目,掌握调试技能。
掌握文件使用编程技能,如文件的各类操作,数据格式设计、通用程序运行过程中文件输入输出运行方式设计等。
一、各算法的算法原理及计算机程序框图
1.列主元高斯消去法
●算法原理:
高斯消去法是利用现行方程组初等变换中的一种变换,即用一个不为零的数乘一个方程后加只另一个方程,使方程组变成同解的上三角方程组,然后再自下而上对上三角方程组求解。
列选住院是当高斯消元到第
步时,从
列的
以下(包括
)的各元素中选出绝对值最大的,然后通过行交换将其交换到
的位置上。
交换系数矩阵中的两行(包括常数项),只相当于两个方程的位置交换了,因此,列选主元不影响求解的结果。
●计算机程序框图如上
●源程序:
#defineN200
#include"stdio.h"
#include"math.h"
FILE*fp1,*fp2;
voidLZ()
{
intn,i,j,k=0,l;
doubled,t,t1;
staticdoublex[N],a[N][N];
fp1=fopen("a1.txt","r");
fp2=fopen("b1.txt","w");
fscanf(fp1,"%d",&n);
for(i=0;i for(j=0;j<=n;++j) { fscanf(fp1,"%lf",&a[i][j]); } do { d=a[k][k]; l=k; i=k+1; do { if(fabs(a[i][k])>fabs(d))/*选主元*/ { d=a[i][k]; l=i; } i++; }while(i if(d==0) {printf("\n输入矩阵有误! \n");} else {/*换行*/ if(l! =k) { for(j=k;j<=n;j++) { t=a[l][j]; a[l][j]=a[k][j]; a[k][j]=t; } } } for(j=k+1;j<=n;j++)/*正消*/ a[k][j]/=a[k][k]; for(i=k+1;i for(j=k+1;j<=n;j++) a[i][j]-=a[i][k]*a[k][j]; k++; }while(k if(k! =0) { for(i=n-1;i>=0;i--)/*回代*/ { t1=0; for(j=i+1;j t1+=a[i][j]*x[j]; x[i]=a[i][n]-t1; } } for(i=0;i fprintf(fp2,"\n方程组的根为x[%d]=%lf",i+1,x[i]); fclose(fp1); fclose(fp2); } main() { LZ(); } ●具体算例及求解结果: 用列选主元法求解下列线性方程组 输入3输出结果: 方程组的根为x[1]=6.000000 12-38方程组的根为x[2]=4.000000 21322方程组的根为x[3]=2.000000 32128 ●输入变量、输出变量说明: 输入变量: 系数矩阵元素, 常向量元素 输出变量: 解向量元素 2.杜里特尔分解法解线性方程 ●算法原理: 求解线性方程组 时,当对 进行杜里特尔分解,则等价于求解 ,这时可归结为利用递推计算相继求解两个三角形(系数矩阵为三角矩阵)方程组,用顺代,由 求出 ,再利用回带,由 求出 。 计算机程序框图: 源程序: #include"stdio.h" #include"math.h" FILE*fp1,*fp2; voidmain() { inti,j,k,N; doubles,A[200][200],B[200],x[200],y[200]; staticdoubleL[200][200],U[200][200]; fp1=fopen("a2.txt","r"); fp2=fopen("b2.txt","w"); fscanf(fp1,"%d",&N); for(i=0;i {for(j=0;j fscanf(fp1,"%lf",&A[i][j]); } for(i=0;i fscanf(fp1,"%lf",&B[i]); for(i=0;i { for(j=i;j { s=0.0; for(k=0;k s+=L[i][k]*U[k][j]; U[i][j]=A[i][j]-s; } for(j=i+1;j { s=0.0; for(k=0;k s+=U[k][i]*L[j][k]; L[j][i]=(A[j][i]-s)/U[i][i]; } } for(i=0;i for(j=0;j L[i][i]=1; fprintf(fp2,"\nU矩阵为: "); for(i=0;i {fprintf(fp2,"\n"); for(j=0;j fprintf(fp2,"%10.3f",U[i][j]); } fprintf(fp2,"\nL矩阵为: "); for(i=0;i {fprintf(fp2,"\n"); for(j=0;j fprintf(fp2,"%10.3f",L[i][j]); } for(i=0;i {s=0.0; for(k=0;k s+=L[i][k]*y[k]; y[i]=B[i]-s; } for(i=N-1;i>=0;i--) {s=0.0; for(k=i+1;k s+=U[i][k]*x[k]; x[i]=(y[i]-s)/U[i][i]; } fprintf(fp2,"\n方程组解为: "); for(i=0;i fprintf(fp2,"\nx%d=%10.3f",i+1,x[i]); fclose(fp1); fclose(fp2); } ●具体算例及求解结果: 用杜里特尔分解法求解方程组 输入数据输出结果: U矩阵为: 2.0003.0004.000 0.000-6.500-4.000 0.0000.000-2.538 3L矩阵为: 2341.0000.0000.000 3-221.5001.0000.000 4232.0000.6151.000 391443方程组解为: x1=6.000 x2=5.000 x3=3.000 ●输入变量、输出变量说明: 输入变量: 系数矩阵元素, 常向量元素 输出变量: 解向量元素 3.拉格朗日插值法 ●算法原理: 首先构造基函数 ,可以证明基函数满足下列条件: , 对于给定 个节点, 次拉格朗日插值多项式由下式给出: 其中 由于 是一个关于 的 次多项式,所以 为关于 的不高于 次的代数多项式。 当 时, ,满足插值条件。 ●计算机程序框图: 源程序: #include"stdio.h" #include"math.h" intn,m,i,j; floatx2,x3,z1=0.0,z=0.0,t,x[50],y[50],c[50],A[50]; main() {FILE*fp1,*fp2; fp1=fopen("a3.txt","r"); fp2=fopen("b3.txt","w"); fscanf(fp1,"%d",&n); for(i=0;i fscanf(fp1,"%f,%f",&x[i],&y[i]); fscanf(fp1,"%d",&m); fscanf(fp1,"%f",&x2); fscanf(fp1,"%f",&x3); for(i=0;i c[i]=fabs(x[i]-x2); for(i=0;i for(j=i+1;j if(c[i]>c[j]){t=c[i];c[i]=c[j];c[j]=t; t=x[i];x[i]=x[j];x[j]=t; t=y[i];y[i]=y[j];y[j]=t;} for(i=0;i {A[i]=1.0; for(j=0;j if(i! =j)A[i]=A[i]*(x2-x[j])/(x[i]-x[j]); z=z+A[i]*y[i]; } for(i=0;i {A[i]=1.0; for(j=0;j if(i! =j)A[i]=A[i]*(x3-x[j])/(x[i]-x[j]); z1=z1+A[i]*y[i]; } fprintf(fp2,"\nx=%10.7f处的函数值为: y=%10.7f",x2,z); fprintf(fp2,"\nx=%10.7f处的函数值为: y=%10.7f",x3,z1); fclose(fp1); fclose(fp2); } 具体算例及求解结果: 对于一组数据表进行二次数值插值编程,根据下面数值表计算f(0.49) 和f(0.51) 0.2 0.4 0.6 0.8 f(x) 16 20 15 10 输入数据: 输出结果: x=0.4900000处的函数值为: y=18.8637505 4x=0.5100000处的函数值为: y=18.3637524 0.2,16 0.4,20 0.6,15 0.8,10 3 0.49 0.51 ●输入变量、输出变量说明: 输入变量: 插值节点 输出变量: 插值所得到被插函数在插值点的近似值 4.曲线拟合 ●算法原理: 对于给定的一组数据 , =1,2…, ,寻求做 次多项式 使性能指标 为最小。 由于性能指标 可以被看做关于 , =0,1,…, 的多元函数,故上述拟合多项式的构造问题可转化为多元函数的极值问题。 令 从而有正则方程组 求解即得多项式系数。 ●计算机程序框图: ●源程序: #include #include main() { inti,j,k,m,n,l,N,t,t1; doublemax,A[50][50],x[50],y[50],S[50],T[50],X[50];floatyb=0.0,xb,a1,a2,a0; FILE*fp1,*fp2; fp1=fopen("a4.txt","r"); fp2=fopen("b4.txt","w"); fscanf(fp1,"%d%d\n",&n,&m); for(i=0;i {fscanf(fp1,"%lf%lf",&x[i],&y[i]); } fscanf(fp1,"%f",&xb); for(i=0;i<=2*m;i++) {S[i]=0.0; for(j=0;j S[i]+=pow(x[j],i); } for(i=0;i<=m;i++) {T[i]=0.0; for(j=0;j T[i]+=y[j]*pow(x[j],i); } N=m+1; for(i=0;i {for(j=0;j {l=i+j; A[i][j]=S[l];} A[i][N]=T[i]; } for(i=0;i { max=fabs(A[i][i]);/*选主*/ for(j=i+1;j if(fabs(A[j][i])>max){max=fabs(A[j][i]);m=j;} if(m! =i) for(k=0;k<=N;k++) {t=A[i][k];A[i][k]=A[m][k];A[m][k]=t;} for(j=i+1;j for(k=N;k>=0;k--) A[j][k]=A[j][k]-A[i][k]*A[j][i]/A[i][i]; } for(i=N-1;i>=0;i--)/*回代*/ for(j=i-1;j>=0;j--) for(k=N;k>=0;k--) A[j][k]=A[j][k]-A[i][k]*A[j][i]/A[i][i]; fprintf(fp2,"\n解为: ");/*输出结果*/ for(i=0;i fprintf(fp2,"\na%d=%10.7lf",i,A[i][N]/A[i][i]); fprintf(fp2,"\n拟合多项式为: \n"); fprintf(fp2,"P(x)=%10.7lf",A[0][N]/A[0][0]); for(i=1;i fprintf(fp2,"+(%10.7lf)x^%d",A[i][N]/A[i][i],i); for(i=0;i yb+=(A[i][N]/A[i][i])*pow(xb,i); fprintf(fp2,"\nP(%f)=%10.7f",xb,yb); fclose(fp1); fclose(fp2); } ●具体算例及求解结果: 对于一组数据表进行二次多项式曲线拟合,根据以下数据胡二次拟合曲线 求y(5) 1 2 3 4 5 6 7 8 9 10 1.6 2.8 3.6 4.9 5.4 6.8 7.9 9.2 10.2 11.4 试用最小二乘法求二次拟合多项式 输入数据: 102输出结果: 11.6解为: 22.8a0=-0.3012450 33.6a1=1.3167338 44.9a2=-0.0166439 55.4拟合多项式为: 66.8P(x)=-0.3012450+(1.3167338)x^1+(-0.0166439)x^2 77.9P(5.000000)=5.8663259 89.2 910.2 1011.4 5.0 输入变量、输出变量说明: 输入变量: 已知数据点 输出变量: 拟合多项式的系数 5.改进欧拉法 ●算法原理: 当 取值较小时,让梯形法的迭代公式只迭代一次就结束。 这样先用欧拉公式求得一个初步近似值 ,称之为预报值,预报值的精度不高,用它替代梯形法右端的 ,再直接计算得出 ,并称之为校正值,这时得到预报-校正公式。 将预报-校正公式 称为改进欧拉公式。 ●计算机程序框图: ●源程序: #include"stdio.h" #include"math.h" FILE*fp1,*fp2; floatfunc(floatx,floaty) {floatdy; dy=sqrt(2*x*x+3*y*y); return(dy);/*定义函数的导*/ } main() {inti; floath,yp,yc,y0,x1,x2; if((fp1=fopen("a5.txt","r"))==NULL) {printf("cannotopenthisfile\n"); exit(0);} fp2=fopen("b5.txt","w"); fscanf(fp1,"%f,%f,%f,%f",&x1,&x2,&y0,&h); for(i=0;i<(x2-x1)/h;i++) {yp=y0+h*func(x1+i*h,y0); yc=y0+h*func(x1+(i+1)*h,yp); y0=0.5*(yp+yc); fprintf(fp2,"节点%6.2f处的值=%10.7f\n",x1+(i+1)*h,y0); } fclose(fp1); fclose(fp2); ●} ●具体算例及求解结果: 求解初值问题。 取h=0.2,用改进欧拉法求解下列初值问题 输入数据: 0,20,5,0.2输出结果节点0.20处的值=7.0323939 节点0.40处的值=9.8918471节点0.60处的值=13.9148121 节点0.80处的值=19.5739155节点1.00处的值=27.5336847 节点1.20处的值=38.7287178节点1.40处的值=54.4735222 节点1.60处的值=76.6169243节点1.80处的值=107.7592354 节点2.00处的值=151.5576096节点2.20处的值=213.1555939 节点2.40处的值=299.7871094节点2.60处的值=421.6261139 节点2.80处的值=592.9812927节点3.00处的值=833.9766235 节点3.20处的值=1172.9145508节点3.40处的值=1649.6000366 节点3.60处的值=2320.0151367节点3.80处的值=3262.8935547 节点4.00处的值=4588.9670410节点4.20处的值=6453.9699707 节点4.40处的值=9076.9287109节点4.60处的值=12765.8852539 节点4.80处的值=17954.0703125节点5.00处的值=25250.7871094 节点5.20处的值=35512.9648438节点5.40处的值=49945.7949219 节点5.60处的值=70244.2773438节点5.80处的值=98792.2734375 节点6.00处的值=138942.4609375节点6.20处的值=195410.0937500 节点6.40处的值=274826.7343750节点6.60处的值=386519.1406250 节点6.80处的值=543604.4218750节点7.00处的值=764530.8125000 节点7.20处的值=1075243.9062500节点7.40处的值=1512233.8750000 节点7.60处的值=2126821.1875000节点7.80处的值=2991183.0000000 节点8.00处的值=4206830.1250000节点8.20处的值=5916528.2500000 节点8.40处的值=8321065.2500000节点8.60处的值=11702831.0000000 节点8.80处的值=16458980.5000000节点9.00处的值=23148077.0000000 节点9.20处的值=32555688.0000000节点9.40处的值=45786650.0000000 节点9.60处的值=64394808.0000000节点9.80处的值=90565512.0000000 节点0.00处的值=127372260.0000000节点10.20处的值=179137632.0000000 节点10.40处的值=251940992.0000000节点10.60处的值=354332368.0000000 节点10.80处的值=498336624.0000000节点11.00处的值=700865696.0000000 节点11.20处的值=985704608.0000000节点11.40处的值=1386304896.0000000 节点11.60处的值=1949713344.0000000节点11.80处的值=2742096640.000000012.00处的值=3856512512.0000000节点12.20处的值=5423838208.0000000 节点12.40处的值=7628141056.0000000节点12.60处的值=10728294912.0000000节点12.80处的值=150********.0000000 节点13.00处的值=21220453376.0000000节点13.20处的值=29844662272.0000000节点13.40处的值=41973835776.0000000 节点13.60处的值=59032424448.0000000节点13.80处的值=83023802368.0000000节点14.00处的值=116765528064.0000000 节点14.20处的值=164220223488.0000000节点14.40处的值=230960996352.0000000节点14.60处的值=324825874432.0000000 节点14.80处的值=456838414336.0000000节点15.00处的值=642502164480.0000000节点15.20处的值=903621574656.0000000 节点15.40处的值=1270862577664.0000000节点15.60处的值=178********88.0000000节点15.80处的值=2513752817664.0000000 节点16.00处的值=3535367438336.0000000节点16.20处的值=4972176474112.0000000节点16.40处的值=6992919527424.0000000 节点16.60处的值=9834913071104.0000000节点16.80处的值=138********536.0000000节点17.00处的值=19453354967040.0000000 节点17.20处的值=27359394660352.0000000节点17.40处的值=38478530215936.0000000节点17.60处的值=54116592123904.0000000 节点17.80处的值=76110125596672.0000000节点18
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 数值 计算方法 上机 实验 报告
![提示](https://static.bdocx.com/images/bang_tan.gif)