《数值计算方法上机实验报告.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 拟
19、合多项式为: 6 6.8 P(x)=-0.3012450+( 1.3167338)x1+(-0.0166439)x2 7 7.9 P(5.000000)= 5.8663259 8 9.2 9 10.2 10 11.4 5.0 输入变量、输出变量说明： 输入变量：(xi,yi)已知数据点 输出变量：ai拟合多项式的系数 5. 改进欧拉法 l 算法原理： 当h取值较小时，让梯形法的迭代公式只迭代一次就结束。这样先用欧拉公式求得一个初步近似值yn+1(0)，称之为预报值，预报值的精度不高，用它替代梯形法右端的yn+1，再直接计算得出yn+1，并称之为校正值，这时得到预报-校正公式。将预报-校正公式
20、12 数值计算方法上机实验报告 yn+1(0)=yn+hf(xn,yn) h(0)f(xn,yn)+f(xn+1,yn+1)yn+1=yn+2称为改进欧拉公式。 l 计算机程序框图： l 源程序： #include stdio.h #include math.h FILE *fp1,*fp2; float func(float x,float y) float dy; dy=sqrt(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(cann
21、ot open this filen); exit(0); 13 定义函数的导*/ 数值计算方法上机实验报告 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.7fn,x1+(i+1)*h,y0); fclose(fp1); fclose(fp2); l l 具体算例及求解结果： 求解初
22、值问题。取h=0.2,用改进欧拉法求解下列初值问题 y/=2x2+3y2y(0)=5 0x20输入数据：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处的值
23、=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 节点
24、 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.8125
25、000 节点 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.0000
26、000 节点 9.40处的值=45786650.0000000 节点 9.60处的值=64394808.0000000 节点 9.80处的值=90565512.0000000 节点0.00处的值=127372260.0000000 节点 10.20处的值=179137632.0000000 14 数值计算方法上机实验报告 节点10.40处的值=251940992.0000000 节点10.60处的值=354332368.0000000 节点 10.80处的值=498336624.0000000 节点11.00处的值=700865696.0000000 节点 11.20处的值=985704608
27、.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处的值=15088381952.0000000 节点 13.00处的值=21220453376.0000000节点13.20处的值=298446
28、62272.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.000000
29、0节点 15.20处的值=903621574656.0000000 节点 15.40处的值=1270862577664.0000000节点 15.60处的值=1787354021888.0000000节点 15.80处的值=2513752817664.0000000 节点 16.00处的值=3535367438336.0000000节点 16.20处的值=4972176474112.0000000 节点 16.40处的值=6992919527424.0000000 节点 16.60处的值=9834913071104.0000000 节点 16.80处的值=13831921729536.0000
30、000 节点 17.00处的值=19453354967040.0000000 节点 17.20处的值=27359394660352.0000000 节点 17.40处的值=38478530215936.0000000 节点 17.60处的值=54116592123904.0000000 节点 17.80处的值=76110125596672.0000000 节点 18.00处的值=107042052243456.0000000 节点 18.20处的值=150545029398528.0000000 节点 18.40处的值=211728063266816.0000000 节点 18.60处的值=2
31、97776508305408.0000000 节点 18.80处的值=418795919245312.0000000 节点 19.00处的值=588998829408256.0000000 节点 19.20处的值=828373899149312.0000000 节点 19.40处的值=1165033501360128.0000000 节点 19.60处的值=1638514923929600.0000000 节点 19.80处的值=2304424236023808.0000000 节点 20.00处的值=3240965678563328.0000000 输入变量、输出变量说明： 输入变量：(x0
32、,y0)处置点，h区间长度，N计算次数 输出变量：(x1,y1)初值问题的数值解法结果 6. 四阶龙格-库塔法 l 算法原理： 用区间xk,xk+1内四个不同点上的函数值的线性组合就得到四阶龙格-库塔法。 15 数值计算方法上机实验报告 四阶龙格-库塔法 yn+1=yn+h(w1k1+w2k2+w3k3+w4k4)k=f(x,y)1nn k2=f(xn+l1h,yn+m11k1h)k=f(x+lh,y+mkh+mkh)n2n2112223k4=f(xn+l3h,yn+m31k1h+m32k2h+m33k3h)其中，w,l,m均为待定系数。 类似于前面的讨论，把k2,k3,k4分别在xn点展开成
33、h的幂级数，代入yn+1并进行花间，然后与y(xn+1)在xn点上的泰勒展开式比较，使其两式比较，使其两式右端直到h4的系数相等，经过复杂的数学演算可得到关于w,l,m的一组特解 1l1=l2=m11=m22=2m21=m31=m32=0 l3=m33=11w1=w4=61w2=w3=3从而得到下列常用的经典公式 hyn+1=yn+(k1+2k2+2k3+k4)6k1=f(xn,yn)hk=f(x,y+k1) 21nn+22hk=f(x,y+k2)31nn+22k=f(x,y+hk)n+1n34经典的龙格-库塔法每一步需要4次计算函数值f(x,y)，它具有四阶精度，即局部截断误差是O(h5)。
34、 l 计算机程序框图： 16 数值计算方法上机实验报告 源程序：#include stdio.h #include math.h FILE *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,k1,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
35、,%f,%f,&h,&x0,&b,&y0); 17 数值计算方法上机实验报告 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*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.7
36、f,xi,yi); else break; fclose(fp1); fclose(fp2); 具体算例及求解结果： 求解初值问题，取h=0.2,用四阶 龙格-库塔法求解下列初值问题 y/=x2+y2y(0)=6 0x20输入数据：0.2,0,20,6 输出结果： 0.00000点处的值= 6.0000000 0.20000点处的值= 7.3286018 0.40000点处的值= 8.9523811 0.60000点处的值=10.9372101 0.80000点处的值=13.3631592 1.00000点处的值=16.3277931 1.20000点处的值=19.9501495 1.4000
37、0点处的值=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.6761780 3.40000点处的值=180.3790741 3.60000点处的值=220.3218079 3.80000点处的值=269.1072998 4.00000点处的值=328.6933289 4.20000点处的值=401.4711609 4.40
