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

数值计算方法上机实验报告数值计算方法上机实验报告 数值计算方法上机实验报告 实验目的：复习和巩固数值计算方法的基本数学模型，全面掌握运用计算机进行数值计算的具体过程及相关问题。利用计算机语言独立编写、调试数值计算方法程序，培养学生利用计算机和所学理论知识分析解决实际问题的能力。 上机练习任务：利用计算机基本C语言编写并调试一系列数值方法计算通用程序，并能正确计算给定题目，掌握调试技能。 掌握文件使用编程技能，如文件的各类操作，数据格式设计、通用程序运行过程中文件输入输出运行方式设计等。 一、 各算法的算法原理及计算机程序框图 1. 列主元高斯消去法 l 算法原理： 高斯消去法是利用现行方程组初等变换中的一种变换，即用一个不为零的数乘一个方程后加只另一个方程，使方程组变成同解的上三角方程组，然后再自下而上对上三角方程组求解。 列选住院是当高斯消元到第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; static double xN,aNN; fp1=fopen("a1.txt","r"); fp2=fopen("b1.txt","w"); fscanf(fp1,"%d",&n); for(i=0;i<n;+i) for(j=0;j<=n;+j) fscanf(fp1,"%lf",&aij); 2 ll 数值计算方法上机实验报告 do d=akk; l=k; i=k+1; do if(fabs(aik)>fabs(d) /*选主元*/ d=aik; l=i; i+; while(i<n); 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;i<n;i+) for(j=k+1;j<=n;j+) aij-=aik*akj; k+; while(k<n); if(k!=0) for(i=n-1;i>=0;i-) /* t1=0; for(j=i+1;j<n;j+) t1+=aij*xj; xi=ain-t1; 3 换行*/ 正消*/ 回代*/ 数值计算方法上机实验报告 for(i=0;i<n;i+) fprintf(fp2,"n方程组的根为x%d=%lf",i+1,xi); fclose(fp1); fclose(fp2); main LZ; l 具体算例及求解结果： 用列选主元法求解下列线性方程组 ìx1+2x2-3x3=8ïí2x1+x2+3x3=22 ï3x+2x+x=2823î1输入3 输出结果：方程组的根为x1=6.000000 1 2 -3 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 main 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;i<N;i+) for(j=0;j<N;j+) fscanf(fp1,"%lf",&Aij); for(i=0;i<N;i+) fscanf(fp1,"%lf",&Bi); for(i=0;i<N;i+) /*LU分解*/ for(j=i;j<N;j+) s=0.0; for(k=0;k<i;k+) s+=Lik*Ukj; Uij=Aij-s; for(j=i+1;j<N;j+) s=0.0; for(k=0;k<i;k+) s+=Uki*Ljk; Lji=(Aji-s)/Uii; for(i=0;i<N;i+) for(j=0;j<N;j+) Lii=1; fprintf(fp2,"nU矩阵为:"); for(i=0;i<N;i+) fprintf(fp2,"n"); for(j=0;j<N;j+) fprintf(fp2,"%10.3f",Uij); fprintf(fp2,"nL矩阵为:"); for(i=0;i<N;i+) fprintf(fp2,"n"); for(j=0;j<N;j+) fprintf(fp2,"%10.3f",Lij); 5 数值计算方法上机实验报告 for(i=0;i<N;i+) s=0.0; for(k=0;k<i;k+) s+=Lik*yk; yi=Bi-s; for(i=N-1;i>=0;i-) s=0.0; for(k=i+1;k<N;k+) s+=Uik*xk; xi=(yi-s)/Uii; fprintf(fp2,"n方程组解为:"); for(i=0;i<N;i+) fprintf(fp2,"nx%d=%10.3f",i+1,xi); fclose(fp1); fclose(fp2); l 具体算例及求解结果： 用杜里特尔分解法求解方程组 ìï2x1+3x2+4x3=39í3x1-2x2+2x3=14 ïî4x1+2x2+3x3=43输入数据 输出结果：U矩阵为: 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=0i¹knx-xi，可以证明基函数满足下列条件： xk-xii¹ki=kì0lk(xi)=íî1， 对于给定n+1个节点，n次拉格朗日插值多项式由下式给出： L(x)=ålk(x)yk k=0n其中 lk(x)=Õi=0i¹knx-xi xk-xi由于lk(x)是一个关于x的n次多项式，所以L(x)为关于x的不高于n次的代数多项式。当x=xi时，L(xi)=yi，满足插值条件。 l 计算机程序框图： 7 数值计算方法上机实验报告 源程序：#include"stdio.h" #include"math.h" int 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;i<n;i+) fscanf(fp1,"%f,%f",&xi,&yi); fscanf(fp1,"%d",&m); fscanf(fp1,"%f",&x2); fscanf(fp1,"%f",&x3); for(i=0;i<n;i+) /*选m个最接近的点*/ ci=fabs(xi-x2); for(i=0;i<n;i+) for(j=i+1;j<n;j+) if(ci>cj) t=ci; ci=cj; cj=t; 8 数值计算方法上机实验报告 t=xi; xi=xj; xj=t; t=yi; yi=yj; yj=t; for(i=0;i<m;i+) /*求值*/ Ai=1.0; for(j=0;j<m;j+) if(i!=j) Ai=Ai*(x2-xj)/(xi-xj); z=z+Ai*yi; for(i=0;i<m;i+) /*求值*/ Ai=1.0; for(j=0;j<m;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) 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 数值计算方法上机实验报告 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从而有正则方程组 émxiåêêåxiåxi2êMêMêxnn+1xååiiëåxåxM23iLLLLiåxåxåxn+2iùéaùéåyiùú0úêún+1êxyaúiê1ú=êåiiú úêúêúMMMúêúêúanûêåxinyiúëåxi2núëûûni求解即得多项式系数。 l 计算机程序框图： l 源程序： #include <math.h> #include <stdio.h> main int 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;i<n;i+) fscanf(fp1,"%lf %lf",&xi,&yi); fscanf(fp1,"%f",&xb); for(i=0;i<=2*m;i+) Si=0.0; for(j=0;j<n;j+) 10 数值计算方法上机实验报告 Si+=pow(xj,i); for(i=0;i<=m;i+) Ti=0.0; for(j=0;j<n;j+) Ti+=yj*pow(xj,i); N=m+1; for(i=0;i<N;i+) for(j=0;j<N;j+) l=i+j; Aij=Sl; AiN=Ti; for(i=0;i<N-1;i+) max=fabs(Aii); /*选主*/ for(j=i+1;j<N;j+) if(fabs(Aji)>max) 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<N;j+) /*消元*/ for(k=N;k>=0;k-) Ajk=Ajk-Aik*Aji/Aii; for(i=N-1;i>=0;i-) /*回代*/ for(j=i-1;j>=0;j-) for(k=N;k>=0;k-) Ajk=Ajk-Aik*Aji/Aii; fprintf(fp2,"n解为:"); /*输出结果*/ for(i=0;i<N;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;i<N;i+) fprintf(fp2,"+(%10.7lf)x%d",AiN/Aii,i); 11 数值计算方法上机实验报告 for(i=0;i<N;i+) yb+=(AiN/Aii)*pow(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 拟合多项式为: 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，并称之为校正值，这时得到预报-校正公式。将预报-校正公式 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("cannot 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 具体算例及求解结果： 求解初值问题。取h=0.2,用改进欧拉法求解下列初值问题 ìy/=2x2+3y2üïïïïíy(0)=5ý ï0£x£20ïïïîþ输入数据：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 14 数值计算方法上机实验报告 节点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处的值=15088381952.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处的值=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.0000000 节点 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处的值=297776508305408.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,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)n2n211222ï3ïîk4=f(xn+l3h,yn+m31k1h+m32k2h+m33k3h)其中，w,l,m均为待定系数。 类似于前面的讨论，把k2,k3,k4分别在xn点展开成h的幂级数，代入yn+1并进行花间，然后与y(xn+1)在xn点上的泰勒展开式比较，使其两式比较，使其两式右端直到h4的系数相等，经过复杂的数学演算可得到关于w,l,m的一组特解 ìï1ïl1=l2=m11=m22=2ïïm21=m31=m32=0ï íl3=m33=1ï1ïw1=w4=6ïï1ïw2=w3=3î从而得到下列常用的经典公式 ìhïyn+1=yn+(k1+2k2+2k3+k4)6ïïk1=f(xn,yn)ïhïk=f(x,y+k1) í21nn+22ïïhk=f(x,y+k2)ï31nn+2ï2ïk=f(x,y+hk)n+1n3î4经典的龙格-库塔法每一步需要4次计算函数值f(x,y)，它具有四阶精度，即局部截断误差是O(h5)。 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,%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.7f",xi,yi); else break; fclose(fp1); fclose(fp2); 具体算例及求解结果： 求解初值问题，取h=0.2,用四阶 龙格-库塔法求解下列初值问题 ìy/=x2+y2üïïïïíy(0)=6ý ï0£x£20ïïïîþ输入数据：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.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.6761780 3.40000点处的值=180.3790741 3.60000点处的值=220.3218079 3.80000点处的值=269.1072998 4.00000点处的值=328.6933289 4.20000点处的值=401.4711609 4.40