数值计算方法上机实验报告

1、数值计算方法上机实验报告数值计算方法上机实验报告实验目的：复习和巩固数值计算方法的基本数学模型，全面掌握运用计算机进行数值计算的具体过程及相关问题。利用计算机语言独立编写、调试数值计算方法程序，培养学生利用计算机和所学理论知识分析解决实际问题的能力。上机练习任务：利用计算机基本C语言编写并调试一系列数值方法计算通用程序，并能正确计算给定题目，掌握调试技能。 掌握文件使用编程技能，如文件的各类操作，数据格式设计、通用程序运行过程中文件输入输出运行方式设计等。一、 各算法的算法原理及计算机程序框图1. 列主元高斯消去法l 算法原理：高斯消去法是利用现行方程组初等变换中的一种变换，即用一个不为零的数

2、乘一个方程后加只另一个方程，使方程组变成同解的上三角方程组，然后再自下而上对上三角方程组求解。列选住院是当高斯消元到第步时，从列的以下（包括）的各元素中选出绝对值最大的，然后通过行交换将其交换到的位置上。交换系数矩阵中的两行（包括常数项），只相当于两个方程的位置交换了，因此，列选主元不影响求解的结果。l 计算机程序框图如上l 源程序：#define N 200#include stdio.h#include math.hFILE *fp1,*fp2;void LZ() int n,i,j,k=0,l; double d,t,t1; static double xN,aNN;fp1=fopen(

3、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+) aij-=aik*akj; k+; while(k=

4、0;i-) /*回代*/ t1=0; for(j=i+1;jn;j+) t1+=aij*xj; xi=ain-t1; for(i=0;in;i+) fprintf(fp2,n方程组的根为x%d=%lf,i+1,xi); fclose(fp1); fclose(fp2); main() LZ(); l 具体算例及求解结果：用列选主元法求解下列线性方程组输入3 输出结果：方程组的根为x1=6.0000001 2 -3 8 方程组的根为x2=4.000000 2 1 3 22 方程组的根为x3=2.0000003 2 1 28l 输入变量、输出变量说明：输入变量：系数矩阵元素，常向量元素输出变量：解

5、向量元素2. 杜里特尔分解法解线性方程l 算法原理：求解线性方程组时，当对进行杜里特尔分解，则等价于求解，这时可归结为利用递推计算相继求解两个三角形（系数矩阵为三角矩阵）方程组，用顺代，由求出，再利用回带，由求出。计算机程序框图：源程序：#include stdio.h#include math.hFILE *fp1,*fp2;void main()int i,j,k,N;double s,A200200,B200,x200,y200;static double L200200,U200200;fp1=fopen(a2.txt,r);fp2=fopen(b2.txt,w);fscanf(fp1

6、,%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;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矩阵为:);

7、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(fp2,%10.3f,Lij); 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

8、%d=%10.3f,i+1,xi); fclose(fp1); fclose(fp2);l 具体算例及求解结果：用杜里特尔分解法求解方程组输入数据 输出结果：U矩阵为: 2.000 3.000 4.000 0.000 -6.500 -4.000 0.000 0.000 -2.5383 L矩阵为:2 3 4 1.000 0.000 0.0003 -2 2 1.500 1.000 0.0004 2 3 2.000 0.615 1.00039 14 43 方程组解为:x1= 6.000x2= 5.000x3= 3.000l 输入变量、输出变量说明：输入变量：系数矩阵元素，常向量元素输出变量：解向量元

9、素3. 拉格朗日插值法l 算法原理：首先构造基函数，可以证明基函数满足下列条件：，对于给定个节点，次拉格朗日插值多项式由下式给出：其中由于是一个关于的次多项式，所以为关于的不高于次的代数多项式。当时，满足插值条件。l 计算机程序框图：源程序：#includestdio.h#includemath.hint n,m,i,j; float x2,x3,z1=0.0,z=0.0,t,x50,y50,c50,A50;main()FILE *fp1,*fp2; fp1=fopen(a3.txt,r); fp2=fopen(b3.txt,w); fscanf(fp1,%d,&n); for(i=0;in;

10、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; t=xi; xi=xj; xj=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+

11、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); fclose(fp2);具体算例及求解结果：对于一组数据表进行二次数值插值编程，根据下面数值表计算f(0.49)和f(0.51)0.20.40.60.8f(x)16201510输入数据： 输出结果：x=

12、 0.4900000处的函数值为:y=18.86375054 x= 0.5100000处的函数值为:y=18.36375240.2,160.4,200.6,150.8,1030.490.51l 输入变量、输出变量说明：输入变量：插值节点输出变量：插值所得到被插函数在插值点的近似值4. 曲线拟合l 算法原理：对于给定的一组数据，=1，2，寻求做次多项式使性能指标为最小。由于性能指标可以被看做关于，=0，1，的多元函数，故上述拟合多项式的构造问题可转化为多元函数的极值问题。令从而有正则方程组求解即得多项式系数。l 计算机程序框图：l 源程序：#include #include main()int

13、i,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+) Si+=pow(xj,i); for(i=0;i=m;i+) Ti=0.

14、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-) /*回代*/ for(j=i-1;j=0;j-) for

15、(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);for(i=0;iN;i+)yb+=(AiN/Aii)*pow(xb,i);fprintf(fp2,nP(%f)=%10.7f,xb,yb);fclose(f

16、p1);fclose(fp2);l 具体算例及求解结果：对于一组数据表进行二次多项式曲线拟合，根据以下数据胡二次拟合曲线求y(5)123456789101.62.83.64.95.46.87.99.210.211.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 拟合多项式为: 6 6.8 P(x)=-0.3012450+( 1.3167338)x1+(-0.0166439)x2 7 7.9 P(5.000000)= 5.866

17、3259 8 9.2 9 10.210 11.45.0 输入变量、输出变量说明：输入变量：已知数据点输出变量：拟合多项式的系数 5. 改进欧拉法l 算法原理：当取值较小时，让梯形法的迭代公式只迭代一次就结束。这样先用欧拉公式求得一个初步近似值，称之为预报值，预报值的精度不高，用它替代梯形法右端的，再直接计算得出，并称之为校正值，这时得到预报-校正公式。将预报-校正公式称为改进欧拉公式。l 计算机程序框图：l 源程序：#include stdio.h#include math.hFILE *fp1,*fp2;float func(float x,float y) float dy; dy=sqr

18、t(2*x*x+3*y*y); return (dy); /*定义函数的导*/ main()int i; float h,yp,yc,y0,x1,x2; if(fp1=fopen(a5.txt,r)=NULL) printf(cannot open this filen); 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);

19、 fprintf(fp2,节点%6.2f处的值=%10.7fn,x1+(i+1)*h,y0); fclose(fp1); fclose(fp2);l l 具体算例及求解结果：求解初值问题。取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处的值=

20、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处的值

21、=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处的值=3865

22、19.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

23、 节点 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处的值=7008656

24、96.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处的值=.0000000节点 13.00处的值=21220453376.0000000

25、节点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处的值=64250

26、2164480.0000000节点 15.20处的值=9.0000000节点 15.40处的值=1270862577664.0000000节点 15.60处的值=1787354021888.0000000节点 15.80处的值=25.0000000节点 16.00处的值=3535367438336.0000000节点 16.20处的值=4972176474112.0000000 节点 16.40处的值=6992919527424.0000000节点 16.60处的值=9834913071104.0000000 节点 16.80处的值=536.0000000 节点 17.00处的值=194533

27、54967040.0000000节点 17.20处的值=27359394660352.0000000 节点 17.40处的值=38478530215936.0000000 节点 17.60处的值=54116592123904.0000000节点 17.80处的值=76110125596672.0000000 节点 18.00处的值=1456.0000000 节点 18.20处的值=8528.0000000节点 18.40处的值=211728063266816.0000000 节点 18.60处的值=297776508305408.0000000 节点 18.80处的值=4312.0000000

28、 节点 19.00处的值=588998829408256.0000000 节点 19.20处的值=828373899149312.0000000 节点 19.40处的值=11658.0000000 节点 19.60处的值=1638514923929600.0000000 节点 19.80处的值=23808.0000000 节点 20.00处的值=32428.0000000 输入变量、输出变量说明：输入变量：处置点，区间长度，计算次数输出变量：初值问题的数值解法结果6. 四阶龙格-库塔法l 算法原理：用区间内四个不同点上的函数值的线性组合就得到四阶龙格-库塔法。四阶龙格-库塔法其中，均为待定系数

29、。类似于前面的讨论，把分别在点展开成的幂级数，代入并进行花间，然后与在点上的泰勒展开式比较，使其两式比较，使其两式右端直到的系数相等，经过复杂的数学演算可得到关于的一组特解从而得到下列常用的经典公式经典的龙格-库塔法每一步需要4次计算函数值，它具有四阶精度，即局部截断误差是。l 计算机程序框图：源程序：#include stdio.h#include math.hFILE *fp1,*fp2;float ds(float x,float y) float d; d=sqrt(x*x+y*y); return d; void main() int i; float x101,y101,h,b,k

30、1,k2,k3,k4; if(fp1=fopen(a6.txt,r)=NULL) printf(cannot open this filen); exit(0); fp2=fopen(b6.txt,w); fscanf(fp1,%f,%f,%f,%f,&h,&x0,&b,&y0); for(i=1;i=100;i+) if(x0+i*h)=b+h) xi=x0+i*h; k1=ds(xi-1,yi-1); k2=ds(xi-1+h/2,yi-1+h*k1/2); k3=ds(xi-1+h/2,yi-1+h*k2/2); k4=ds(xi,yi-1+h*k3); yi=yi-1+h*(k1+2*

31、k2+2*k3+k4)/6.0; else break;for(i=0;i=100;i+) if(x0+i*h)=b+h) xi=x0+i*h; fprintf(fp2,n%10.5f点处的值=%10.7f,xi,yi); else break; fclose(fp1); fclose(fp2);具体算例及求解结果：求解初值问题，取h=0.2,用四阶 龙格-库塔法求解下列初值问题输入数据：0.2,0,20,6输出结果： 0.00000点处的值= 6.0000000 0.20000点处的值= 7.3286018 0.40000点处的值= 8.9523811 0.60000点处的值=10.9372

32、101 0.80000点处的值=13.3631592 1.00000点处的值=16.3277931 1.20000点处的值=19.9501495 1.40000点处的值=24.3755608 1.60000点处的值=29.7815247 1.80000点处的值=36.3848495 2.00000点处的值=44.4503784 2.20000点处的值=54.3016205 2.40000点处的值=66.3337555 2.60000点处的值=81.0294876 2.80000点处的值=98.9784317 3.00000点处的值=120.9007721 3.20000点处的值=147.676

33、1780 3.40000点处的值=180.3790741 3.60000点处的值=220.3218079 3.80000点处的值=269.1072998 4.00000点处的值=328.6933289 4.20000点处的值=401.4711609 4.40000点处的值=490.3615112 4.60000点处的值=598.9317017 4.80000点处的值=731.5388794 5.00000点处的值=893.5048828 5.20000点处的值=1091.3298340 5.40000点处的值=1332.9528809 5.60000点处的值=1628.0709229 5.80

34、000点处的值=1988.5278320 6.00000点处的值=2428.7895508 6.20000点处的值=2966.5251465 6.40000点处的值=3623.3151855 6.60000点处的值=4425.5185547 6.80000点处的值=5405.3295898 7.00000点处的值=6602.0703125 7.20000点处的值=8063.7695313 7.40000点处的值=9849.0888672 7.60000点处的值=12029.6777344 7.80000点处的值=14693.0488281 8.00000点处的值=17946.0898438 8

35、.20000点处的值=21919.3554688 8.40000点处的值=26772.3007813 8.60000点处的值=32699.6894531 8.80000点处的值=39939.4023438 9.00000点处的值=48781.98828139.20000点处的值=59582.3203125 9.40000点处的值=72773.8437500 9.60000点处的值=88885.9765625 9.80000点处的值=108565.3359375 10.00000点处的值=132601.7031250 10.20000点处的值=161959.7187500 10.40000点处的

36、值=197817.5937500 10.60000点处的值=241614.4062500 10.80000点处的值=295107.8437500 11.00000点处的值=360444.7187500 11.20000点处的值=440247.1875000 11.40000点处的值=537717.9375000 11.60000点处的值=656768.6875000 11.80000点处的值=802177.2500000 12.00000点处的值=979779.3125000 12.20000点处的值=1196702.5000000 12.40000点处的值=1461652.3750000 1

37、2.60000点处的值=1785262.2500000 12.80000点处的值=2180519.2500000 13.00000点处的值=2663286.2500000 13.20000点处的值=3252937.7500000 13.40000点处的值=3973138.250000013.60000点处的值=4852791.0000000 13.80000点处的值=5927199.0000000 14.00000点处的值=7239481.0000000 14.20000点处的值=8842302.0000000 14.40000点处的值=10799988.0000000 14.60000点处的值=13191105.0000000 14.80000点处的值=16111616.0000000 15.00000点处的值=19678728.0000000 15.20000点处的值=24035598.0000000 15.40000点处的值=29357080.0000000 15.60000点处的值=35856736.