东北大学数值分析实验报告.docx
《东北大学数值分析实验报告.docx》由会员分享，可在线阅读，更多相关《东北大学数值分析实验报告.docx（49页珍藏版）》请在冰点文库上搜索。
东北大学数值分析实验报告
实验一迭代格式的比较
一、问题提出
设方程
f（x）=x
-3x–1=0有三个实根x
=1.8793,x
=-0.34727,x
=-1.53209现采用下面三种不同计算格式，求f（x）=0的根x
或x
1、x=
2、x=
3、x=
二、要求
1、编制一个程序进行运算，最后打印出每种迭代格式的敛散情况；
2、用事后误差估计
来控制迭代次数，并且打印出迭代的次数；
3、初始值的选取对迭代收敛有何影响；
4、分析迭代收敛和发散的原因。
三、目的和意义
1、通过实验进一步了解方程求根的算法；
2、认识选择计算格式的重要性；
3、掌握迭代算法和精度控制；
4、明确迭代收敛性与初值选取的关系。
四、程序设计流程图
五、源程序代码
#include
#include
floatone（floatx0）
{
x1=（3*x0+1）/（x0*x0）；
return（x1）；
}
floattwo（floatx0）
{
x2=（pow（x0,3）-1）/3；
return（x2）；
}
floatthree（floatx0）
{
x3=pow（3*x0+1,0.33333）；
return（x3）；
}
voidmain（）
{
floatx,x0;
floatx1,x2,x3；
intk1,k2,k3;
printf（pleaseinputx0="）；
scanf（"%f"，&x）；
x0=x;
x1=one（x0）;
printf（"第一个公式迭代结果:
\n"）;
while（fabs（x0-x1）>1e-5）
{
printf（"x1=%6.5f\n",x1）;
x0=x1;
x1=one（x0）;
k1++;
}
printf（"x1=%6.5f\n",x1）;
printf（"k1=%i",k1）;
x0=x;
x2=two（x0）;
printf（"第二个公式迭代结果为:
\n"）;
while（fabs（x0-x2）>1e-5）
{
printf（"x2=%6.5f\n",x2）;
x0=x;
x1=two（x0）;
k2++;
}
printf（"x2=%6.5f\n",x2）;
printf（"k2=%i",k2）;
x0=x;
x3=three（x0）;
printf（"第三个公式迭代结果为:
\n"）;
while（fabs（x0-x3）>1e-5）
{
printf（"x3=%6.5f\n",x3）;
x0=x;
x1=three（x0）;
k3++;
}
printf（"x3=%6.5f\n",x3）;
printf（"k3=%i",k3）;
}
运行结果如下：
结果分析:
初值对迭代结果影响解析：
实验二线性方程组的直接算法
一、问题提出
给出下列几个不同类型的线性方程组，请用适当算法计算其解。
1、设线性方程组
=
x
=（1,-1,0,1,2,0,3,1,-1,2）
Gsuss列主元消去法
#include
#include
#include
#defineMAX100
typedefstruct{
introw,col;
floatMAT[MAX][MAX];
floatSolution[MAX];
}Matrix;
voidGauss（Matrix*M）;
voidMBack（Matrix*M）;
voidMSave（Matrix*M）;
voidMInput（Matrix*M）;
voidMOutput（Matrix*M）;
voidSolution（Matrix*M）;
voidMSort（Matrix*M,intn）;
voidmain（）
{
printf（"列主元方法如下：
\n"）;
MatrixMat;
MInput（&Mat）;
MSave（&Mat）;
Gauss（&Mat）;
MSave（&Mat）;
if（Mat.row==Mat.col-1）{
MBack（&Mat）;
Solution（&Mat）;
}
printf（"Pressanykeytohalt..."）;
getch（）;
}
voidMInput（Matrix*M）
{
inti,j;
printf（"输入行数:
"）;scanf（"%d",&M->row）;
printf（"输入列数:
"）;scanf（"%d",&M->col）;
for（i=0;irow;i++）{
printf（"第%d行:
",i+1）;
for（j=0;jcol;j++）{
scanf（"%f",&M->MAT[i][j]）;
}
}
for（i=0;irow;i++）
M->Solution[i]=0;
}
voidMOutput（Matrix*M）
{
inti,j;
printf（"MATRIX:
\n"）;
for（i=0;irow;i++）{
for（j=0;jcol;j++）
printf（"%10.3f",M->MAT[i][j]）;
printf（"\n"）;
}
printf（"---END----\n"）;
}
voidGauss（Matrix*M）
{
inti,j,k;
floattemp;
for（i=0;irow-1;i++）{
MSort（M,i）;MOutput（M）;
for（j=i+1;jrow;j++）{
temp=M->MAT[j][i];
for（k=0;kcol;k++）
if（temp!
=0）{
M->MAT[j][k]/=temp;
M->MAT[j][k]*=M->MAT[i][i];
M->MAT[j][k]-=M->MAT[i][k];
}
}
}
MOutput（M）;
}
voidMSort（Matrix*M,intn）
{
inti,j,k;
floattemp[MAX];
for（i=n;irow-1;i++）{
for（j=n;jrow-i-1;j++）{
if（fabs（M->MAT[j][n]）MAT[j+1][n]））{
for（k=0;kcol;k++）{
temp[k]=M->MAT[j+1][k];
M->MAT[j+1][k]=M->MAT[j][k];
M->MAT[j][k]=temp[k];
}
}
}
}
}
voidMBack（Matrix*M）
{
inti,j;
floatsum;
M->Solution[M->row-1]=M->MAT[M->row-1][M->col-1]
M->MAT[M->row-1][M->row-1];
for（i=M->row-2;i>=0;i--）{
sum=M->MAT[i][M->col-1];
for（j=i+1;jrow;j++）
sum-=M->MAT[i][j]*M->Solution[j];
M->Solution[i]=sum/M->MAT[i][i];
}
}
voidSolution（Matrix*M）
{
inti;
printf（"Solution:
\n"）;
for（i=0;irow;i++）
printf（"X[%d]=%f\n",i+1,M->Solution[i]）;
printf（"\n---END---\n"）;
}
voidMSave（Matrix*M）
{
inti,j;
FILE*eryar;
eryar=fopen（"Matrix.txt","a"）;
fprintf（eryar,"--------BEGIN--------\n"）;
for（i=0;irow;i++）{
for（j=0;jcol;j++）
fprintf（eryar,"%10.3f",M->MAT[i][j]）;
fprintf（eryar,"\n"）;
}
fclose（eryar）;
}
2、设对称正定阵系数阵线方程组
=
x
=（1,-1,0,2,1,-1,0,2）
平方根法:
源程序:
#include
#include
#include
usingnamespacestd;
intmain（）
{
intn,i,j,k,m;
cout<<"输入维数:
";
cin>>n;
double**A=newdouble*[（n+1）];
for（i=1;i<=n;i++）
A[i]=newdouble[n+1];
double*b=newdouble[n+1];
double*x=newdouble[n+1];
double*y=newdouble[n+1];
cout<<"输入系数对称正定矩阵A[][]:
"<for（i=1;i<=n;i++）
for（j=1;j<=n;j++）
cin>>A[i][j];
cout<<"输入向量b[]:
";
for（i=1;i<=n;i++）
cin>>b[i];
cout<for（k=1;k<=n;k++）
{
doublesum=0;
for（m=1;m<=k-1;m++）
{
sum=sum+pow（A[k][m],2.0）;
}
sum=A[k][k]-sum;
A[k][k]=sqrt（sum）;
for（i=k+1;i<=n;i++）
{
doubletemp1=0;
for（m=1;m<=k-1;m++）
{
temp1=temp1+A[i][m]*A[k][m];
}
temp1=A[i][k]-temp1;
A[i][k]=temp1/A[k][k];
}
doubletemp2=0;
for（m=1;m<=k-1;m++）
{
temp2=temp2+A[k][m]*y[m];
}
y[k]=（b[k]-temp2）/A[k][k];
}
x[8]=y[8]/A[8][8];
for（k=n-1;k>=1;k--）
{
doubletemp3=0;
for（m=k+1;m<=n;m++）
{
temp3=temp3+A[m][k]*x[m];
}
x[k]=（y[k]-temp3）/A[k][k];
}
cout<<"输出结果向量x[]:
"<for（i=1;i<=n;i++）cout<return0;
}
3、三对角形线性方程组
=
x
=（2,1,-3,0,1,-2,3,0,1,-1）
追赶法:
源程序
#include
#include
#include
usingnamespacestd;
intmain（）
{
intn,i;
cout<<"输入系数矩阵的维数:
";
cin>>n;
double*a=newdouble[n+1];
double*c=newdouble[n+1];
double*d=newdouble[n+1];
double*b=newdouble[n+1];
double*x=newdouble[n+1];
double*y=newdouble[n+1];
cout<<"输入系数矩阵A[]数据：
"<for（i=1;i<=n;i++）cin>>a[i];
for（i=1;i<=n;i++）cin>>c[i];
for（i=1;i<=n;i++）cin>>d[i];
cout<<"输入b[]:
"<for（i=1;i<=n;i++）cin>>b[i];
for（i=1;i<=n-1;i++）
{
c[i]=c[i]/a[i];
a[i+1]=a[i+1]-d[i+1]*c[i];
}
cout<<"输出解向量a[]:
"<for（i=1;i<=n;i++）cout<cout<<"输出解向量c[]:
"<for（i=1;i<=n;i++）cout<y[1]=b[1]/a[1];
for（i=2;i<=n;i++）
{
y[i]=（b[i]-d[i]*y[i-1]）/a[i];
}
cout<<"输出解向量y[]:
"<for（i=1;i<=n;i++）cout<x[n]=y[n];
for（i=n-1;i>=1;i--）
{
x[i]=y[i]-c[i]*x[i+1];
}
cout<<"输出解向量x[]:
"<for（i=1;i<=n;i++）cout<system（"pause"）;
return0;
}
运行结果
课题三线性方程组的迭代法
一、问题提出
对课题二所列目的和意义的线性方程组，试分别选用Jacobi迭代法，Gauss-Seidol迭代法和SOR方法计算其解。
二、要求
1、体会迭代法求解线性方程组，并能与消去法做以比较；
2、分别对不同精度要求，如
由迭代次数体会该迭代法的收敛快慢；
3、对方程组2，3使用SOR方法时，选取松弛因子
=0.8，0.9，1，1.1，1.2等，试看对算法收敛性的影响，并能找出你所选用的松弛因子的最佳者；
4、给出各种算法的设计程序和计算结果。
三、目的和意义
1、通过上机计算体会迭代法求解线性方程组的特点，并能和消去法比较；
2、运用所学的迭代法算法，解决各类线性方程组，编出算法程序；
3、体会上机计算时，终止步骤
<
或k>（予给的迭代次数），对迭代法敛散性的意义；
4、体会初始解x
，松弛因子的选取，对计算结果的影响。
1．Jacobi迭代法
源程序:
#include
#include
#include
usingnamespacestd;
intmain（）
{
intn,N;
inti,j,k;
doublee;
cout<<"输入维数n:
";
cin>>n;
cout<<"输入最大迭代次数N:
";
cin>>N;
cout<<"输入精度e:
";
cin>>e;
double**A=newdouble*[（n+1）];
for（i=1;i<=n;i++）
A[i]=newdouble[n+1];
double*b=newdouble[n+1];
double*x0=newdouble[n+1];
double*x=newdouble[n+1];
cout<<"输入系数矩阵A[][]:
"<for（i=1;i<=n;i++）
{
for（j=1;j<=n;j++）
{
cin>>A[i][j];
}
}
cout<<"输入右端项b[]:
"<for（i=1;i<=n;i++）cin>>b[i];
cout<<"初始化解向量x0[]:
"<for（i=1;i<=n;i++）cin>>x0[i];
for（k=1;k<=N;k++）
{
for（i=1;i<=n;i++）
{
doublesum=0;
for（j=1;j<=n;j++）
{
sum=sum+A[i][j]*x0[j];
}
x[i]=x0[i]+（b[i]-sum）/A[i][i];
}
doublemax=0;
for（i=1;i<=n;i++）
{
doublem;
m=fabs（x[i]-x0[i]）;
if（m>max）max=m;
}
if（max<=e）
{
cout<<"输出迭代次k:
"<cout<<"迭代后的解x[]:
"<for（i=1;i<=n;i++）cout<system（"pause"）;
return0;
}
else
{
for（i=1;i<=n;i++）
{
x0[i]=x[i];
}
}
}
cout<<"已达到最大迭代次数:
"<cout<<"迭代后的解x[]:
"<for（i=1;i<=n;i++）cout<system（"pause"）;
return0;
}
2．Gauss-Seidol迭代法
源程序
#include
#include
#include
usingnamespacestd;
intmain（）
{
intn,N;
inti,j,k;
doublee;
cout<<"输入维数n:
";
cin>>n;
cout<<"输入最大迭代次数N:
";
cin>>N;
cout<<"输入精度e:
";
cin>>e;
double**A=newdouble*[（n+1）];
for（i=1;i<=n;i++）
A[i]=newdouble[n+1];
double*b=newdouble[n+1];
double*x0=newdouble[n+1];
double*x=newdouble[n+1];
cout<<"输入系数矩阵A[][]:
"<for（i=1;i<=n;i++）
{
for（j=1;j<=n;j++）
{
cin>>A[i][j];
}
}
cout<<"输入右端项b[]:
"<for（i=1;i<=n;i++）cin>>b[i];
cout<<"初始化解向量x0[]:
"<for（i=1;i<=n;i++）cin>>x0[i];
for（k=1;k<=N;k++）
{
for（i=1;i<=n;i++）
{
doublesum=0;
for（j=1;j<=n;j++）
{
if（jelse
sum=sum+A[i][j]*x0[j];
}
x[i]=x0[i]+（b[i]-sum）/A[i][i];
}
doublemax=0;
for（i=1;i<=n;i++）
{
doublem;
m=fabs（x[i]-x0[i]）;
if（m>max）max=m;
}
if（max<=e）
{
cout<<"输出迭代次k:
"<cout<<"迭代后的解x[]:
"<for（i=1;i<=n;i++）cout<system（"pause"）;
return0;
}
else
{
for（i=1;i<=n;i++）
{
x0[i]=x[i];
}
}
}
cout<<"已达到最大迭代次数:
"<cout<<"迭代后的解x[]:
"<for（i=1;i<=n;i++）cout<system（"pause"）;
return0;
}
3．SOR方法
源程序
#include
#include
#include
#include
usingnamespacestd;
intmain（）
{
intn;
cout<<"输入维数:
";
cin>>n;
vector>va;
vectorvd;
cout<<"输入系数矩阵:
"<for（inti=0;ifor（intj=0;j<=n;++j）{
doubledtemp;
cin>>dtemp;
vd.push_back（dtemp）;
}
va.push_back（vd）;
vd.clear（）;
}
cout<<"输入初始解向量:
";
vectorvtemp;
for（inti=0;idoubledtemp;
cin>>dtemp;
vd.push_back（dtemp）;
}
cout<<"输入精度:
";
doubled;
cin>>d;
cout<<"输入最大迭代次数:
";
intmcount;
cin>>mcount;
cout<<"输入松弛因子:
";
doublexd;
cin>>xd;
vector>:
:
iteratoria,iaend;
ia=va.begin（）;
boolflag=true;
intcount=0;
while（flag）{
doubledmax=0.0;
vtemp.assign（vd.begin（）,vd.end（））;
for（inti=0;idoubledtemp=0.0;
for（intj=0;j
升级会员 