中国石油大学(华东)计算方法上机实验程序

1.二分法#includeiostream.hdoublef(doublex);voidmain(){doublea,b,x,ep;a=0;b=5;ep=0.000001;while(b-a=ep){x=(a+b)/2;if(f(a)*f(x)0)a=x;elseb=x;}x=(a+b)/2;couttherootoff(x)isxendl;}doublef(doublex){doublef1;f1=x*x*x-6*x-1;returnf1;}2.牛顿迭代法#includeiostream.h#includemath.hdoublef1(doublex);doublef2(doublex);voidmain(){intflag,nmax,k;doublex,ep,x0;x=2;ep=000001;nmax=200;flag=1;k=1;while(fabs(x-x0)ep&&(knmax)){x0=x;if(fabs(f2(x0))ep){flag=0;break;}x=x0-f1(x0)/f2(x0);k=k+1;}if(k=nmax)flag=0;if(flag==0)coutthenewtonmethodisfailureendl;elsecouttherootoff(x)isxendl;}doublef1(doublex){doublef1;f1=x*x*x-3*x-1;returnf1;}doublef2(doublex){doublef1;f1=3*x*x-3;returnf1;}3.列主元消去法#includeiostream.h#includemath.hvoidmain(){staticdoublea[3][3]={{7,8,11},{5,1,-3},{1,2,3}};staticdoubleb[3]={-3,-4,1};doubleep=0.001;intip=1;intn=3;doubledmax,temp,s;intm,i,j,k;for(k=0;k=n-2;k++){dmax=fabs(a[k][k]);m=k;for(i=k+1;i=n-1;i++){if(fabs(a[i][k])dmax){dmax=fabs(a[i][k]);m=i;}}if(dmaxep){ip=-1;break;}if(m!=k){for(j=k;j=n-1;j++){temp=a[k][j];a[k][j]=a[m][j];a[m][j]=temp;}temp=b[k];b[k]=b[m];b[m]=temp;}for(i=k+1;i=n-1;i++){a[i][k]=a[i][k]/a[k][k];for(j=k+1;j=n-1;j++){a[i][j]=a[i][j]-a[i][k]*a[k][j];}b[i]=b[i]-a[i][k]*b[k];}}b[n-1]=b[n-1]/a[n-1][n-1];for(i=n-2;i=0;i--){s=0;for(j=i+1;j=n-1;j++){s=s+a[i][j]*b[j];}b[i]=(b[i]-s)/a[i][i];}if(ip==-1)coutgaussmethodisfailureendl;elsecoutthesolutionofequtionisb[0],b[1],b[2]endl;}4.LU分解法#includeiostream#includefstreamusingnamespacestd;#defineN3intmain(){floatA[N][N],U[N][N],L[N][N],X[N],Y[N],B[N];floats=0;inti,j,k;cout请输入将要分解的矩阵:endl;for(i=0;iN;i++)//输入待分解的矩阵for(j=0;jN;j++)cinA[i][j];for(i=0;iN;i++)//初始化for(j=0;jN;j++){L[i][j]=0;U[i][j]=0;if(i==j)L[i][j]=1;}for(j=0;jN;j++)//求U矩阵第一行U[0][j]=A[0][j];for(i=1;iN;i++)L[i][0]=A[i][0]/U[0][0];//求L矩阵的第一列for(k=1;kN;k++)//用k记录化解的行列数，在行和列分解都完成时再进行下一行列的分解{for(j=k;jN;j++)//本循环将完成行分解的工作{for(intu=0;u=k-1;u++){s=s+L[k][u]*U[u][j];}U[k][j]=A[k][j]-s;s=0;//每行的和不一致，一行使用完后要置0}for(i=k+1;iN;i++)//在行分解完成后进行列分解{for(intl=0;l=k-1;l++){s=s+L[i][l]*U[l][k];}L[i][k]=(A[i][k]-s)/U[k][k];s=0;//每列的和不一致，一行使用完后要置0}}coutendlLendl;for(i=0;iN;i++){//输出LUfor(j=0;jN;j++)coutL[i][j]\t;coutendl;}coutendlUendl;for(i=0;iN;i++){for(j=0;jN;j++)coutU[i][j]\t;coutendl;}cout请输入矩阵Bendl;coutendl;for(i=0;iN;i++)cinB[i];Y[0]=B[0];//求Yfor(k=1;kN;k++){for(j=0;jk;j++){s=s+L[k][j]*Y[j];}Y[k]=B[k]-s;s=0;}coutendlY:endl;//输出Yfor(i=0;iN;i++)coutY[i]\t;coutendl;X[N-1]=Y[N-1]/U[N-1][N-1];//求Xfor(k=N-2;k=0;k--){for(j=k+1;jN;j++){s=s+U[k][j]*X[j];}X[k]=(Y[k]-s)/U[k][k];s=0;}coutendlX:endl;//输出Xfor(i=0;iN;i++)coutX[i]\t;coutendl;system(pause);return0;}5.追赶法#includeiostreamusingnamespacestd;intmain(){doublea[15],b[15],c[15],d[15];doublet;inti,n;cout请输入n=;cinn;coutendl;cout请输入b[1],c[1],d[1]:endl;cinb[1]c[1]d[1];coutendl;cout输入2行到n-1行:endl;for(i=2;i=n-1;i++){cina[i]b[i]c[i]d[i];}coutendl;cout请输入a[n],b[n],d[n]:endl;cina[n]b[n]d[n];c[1]=c[1]/b[1];d[1]=d[1]/b[1];for(i=2;i=n-1;i++){t=b[i]-c[i-1]*a[i];c[i]=c[i]/t;d[i]=(d[i]-d[i-1]*a[i])/t;}d[n]=(d[n]-d[n-1]*a[n])/(b[n]-c[n-1]*a[n]);for(i=n-1;i=1;i--)d[i]=d[i]-c[i]*d[i+1];for(i=1;i=n;i++)coutendlx[i]=d[i];coutendl;return0;}6.平方根法#includestdio.h#includestdlib.h#includemath.h#defineEPS1.0e-8#defineN20doublea[N][N],b[N],x[N];intn;intzhuyuan(introw);/*选主元*/voidhangjiaohuan(introw1,introw2);/*行交换*/voidxiaoyuan(introw);/*消元*/voidhuidai();/*回代*/voidmain(){printf(请输入方程的维数n!\nn=);scanf(%d,&n);getchar();printf(\n输入%d行%d列矩阵\n,n,n);for(inti=0;in;i++){for(intj=0;jn;j++)scanf(%lf,&a[i][j]);getchar();}printf(\n输入线性方程组右端项b[%d]:,n);for(i=0;in;i++){scanf(%lf,&b[i]);}getchar();for(i=0;in-1;i++){doublerowmark=zhuyuan(i);if(rowmark==-1){printf(多解!);system(pause);return;}if(rowmark!=i){hangjiaohuan(i,rowmark);}xiaoyuan(i);}huidai();printf(\n线性方程组的解为:);for(i=0;in;i++){printf(x%d=%lf,i+1,x[i]);}printf(\n);system(pause);}intzhuyuan(introw){doubleelem=a[row][row];introwmark=row;for(inti=row+1;in;i++){if(elema[i][row]){elem=a[i][row];rowmark=i;}}if(fabs(elem)=EPS){return-1;}returnrowmark;}voidhangjiaohuan(introw1,introw2){doubletmp;tmp=b[row1];b[row1]=b[row2];b[row2]=tmp;for(intj=0;jn;j++){tmp=a[row1][j];a[row1][j]=a[row2][j];a[row2][j]=tmp;}}voidxiaoyuan(introw){for(inti=row+1;in;i++){intj=row;doubletmp=a[i][j]/a[row][row];b[i]-=tmp*b[row];for(a[i][j++]=0;jn;j++){a[i][j]-=tmp*a[row][j];}}}voidhuidai(){x[n-1]=b[n-1]/a[n-1][n-1];for(inti=n-2;i=0;i--){doublesum=0.0;for(intj=i+1;jn;j++){sum-=a[i][j]*x[j];}x[i]=(b[i]+sum)/a[i][i];}}7.G-S迭代法#includeiostream#includemath.hvoidmain(){doubleep=0.000001;intn=3,ip=1,nmax=200;doublec[3][3]={{0,0.2,0.1},{0.2,0,0.1},{0.2,0.4,0}};staticdoubled[3]={0.3,1.5,2};staticdoublex[3]={0.3,1.5,2};inti,j,k;doubleemax,s;k=0;do{emax=0;for(i=0;in;i++){s=d[i];for(j=0;jn;j++){s=s+c[i][j]*x[j];}if(fabs(s-x[i])emax)emax=fabs(s-x[i]);x[i]=s;}k=k+1;}while((em