高等教育东北大学数值分析实验报告.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
intk=1,k1=1,k2=1,k3=1;
floatx1,x2,x3;
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）;
}
main（）
{
floatx,x0;
printf（"输入x0="）;
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\n",k1）;
x0=x;
x2=two（x0）;
printf（"第二个公式迭代结果:
\n"）;
while（fabs（x0-x2）>1e-5）
{
printf（"x2=%6.5f\n",x2）;
x0=x2;
x2=two（x0）;
k2++;
}
printf（"k2=%i\n",k2）;
x0=x;
x3=three（x0）;
printf（"第三个公式迭代结果:
\n"）;
while（fabs（x0-x3）>1e-5）
{
printf（"x3=%6.5f\n",x3）;
x0=x3;
x3=three（x0）;
k3++;
}
printf（"x3=%6.5f\n",x3）;
printf（"k3=%i\n",k3）;
scanf（"%"）;
}
实验结果：
四、程序运行结果讨论和分析：
对于第一种迭代格式，收敛区间[-8.2-0.4]，在该收敛区间内迭代收敛于-1.53209，只能求得方程的一个根；
对于第二种迭代格式，收敛区间[-1.51.8]，在该收敛区间内迭代收敛于-0.34730，同样只能求得方程的一个根；
对于第三种迭代格式，收敛区间[-0.3+∞），在该收敛区间内迭代收敛于1.87937，只能求得方程的一个根；
由以上结果很容易发现，初值的选取对迭代敛散性有很大影响。
以第一种迭代格式为例，当初值大于等于-0.3时，迭代格式发散；当初值小于等于-8.3时，迭代格式也发散；只有初值在-0.3和-8.3之间时，迭代格式才收敛于—1.53209。
其他迭代格式也有这样的性质，即收敛于某个数值区间，超出这个区间迭代格式就是发散的，这就是所谓迭代格式的收敛性。
对于不同迭代格式在不同区间具有不同的敛散性的原因，我认为可以从一下两方面理解：
1、迭代法是一种逐次逼近法，其基本思想是将隐式方程归结为一组显式的计算公式，就是说，迭代过程实质上是个逐步显式化的过程。
2、我们可以用几何图像来更好地理解迭代过程。
由图可知，在某些区间选取的初始值随着迭代次数的增加会越来越逼近精确值，即收敛于精确值，而在另外一些区间选取的初始值随着迭代次数的增加却离精确值越来越远，即不会收敛于一个确定值。
课题二线性方程组的直接算法
一、问题提出
给出下列几个不同类型的线性方程组，请用适当算法计算其解。
1、设线性方程组
=
x
=（1,-1,0,1,2,0,3,1,-1,2）
2、设对称正定阵系数阵线方程组
=
x
=（1,-1,0,2,1,-1,0,2）
3、三对角形线性方程组
=
x
=（2,1,-3,0,1,-2,3,0,1,-1）
二、要求
1、对上述三个方程组分别利用Gauss顺序消去法与Gauss列主元消去法；平方根法与改进平方根法；追赶法求解（选择其一）；
2、应用结构程序设计编出通用程序；
3、比较计算结果，分析数值解误差的原因；
4、尽可能利用相应模块输出系数矩阵的三角分解式。
三、目的和意义
1、通过该课题的实验，体会模块化结构程序设计方法的优点；
2、运用所学的计算方法，解决各类线性方程组的直接算法；
3、提高分析和解决问题的能力，做到学以致用；
4、通过三对角形线性方程组的解法，体会稀疏线性方程组解法的特点。
程序代码：
1.Gsuss列主元消去法
#include
#include
#include
usingnamespacestd;
intmain（）
{
intn,i,j,k,m;
cout<<"输入维数:
";
cin>>n;
float**A=newdouble*[（n+1）];
for（i=1;i<=n;i++）
A[i]=newdouble[n+1];
float*b=newdouble[n+1];
float*x=newdouble[n+1];
floatl;
floattemp1,temp2,temp3;
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{
temp1=abs（A[k][k]）;
m=k;
for（i=k;i<=n;i++）//找最大值的列主元
{
if（temp1}
if（temp1==0）cout<<"nouniquesolution!
"<if（m!
=k）//换行
{
for（j=1;j<=n;j++）
{
temp2=A[k][j];
A[k][j]=A[m][j];
A[m][j]=temp2;
}
temp3=b[k];
b[k]=b[m];
b[m]=temp3;
}
for（i=k+1;i<=n;i++）//消元
{
l=A[i][k]/A[k][k];
for（j=k+1;j<=n;j++）
{
A[i][j]=A[i][j]-l*A[k][j];
}
b[i]=b[i]-l*b[k];
}
}
if（A[n][n]==0）
{
cout<<"nouniquesolution!
"<}
x[n]=b[n]/A[n][n];//回代求解
for（i=n-1;i>=1;i--）
{
floatsum=0;
for（j=i+1;j<=10;j++）
sum=sum+A[i][j]*x[j];
x[i]=（b[i]-sum）/A[i][i];
}
cout<<"输出结果向量x[]:
"<for（i=1;i<=10;i++）cout<system（"pause"）;
return0;
}
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<system（"pause"）;
return0;
}
3.追赶法
#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;
}
四、程序运行结果分析
在方法的选择上存在一定的误差，可以选一些更准确的方法求解；程序中对变量的类型设定，若设成double型，结果可以更精确；计算机在做运算时，会根据需要对中间结果进行舍入，这也会对最终结果有影响；
课题三线性方程组的迭代法
一、问题提出
1、设线性方程组
=
x
=（1,-1,0,1,2,0,3,1,-1,2）
J迭代算法：
程序设计流程图：
源程序代码：
#include
#include
#include
voidmain（）
{
floata[50][51],x1[50],x2[50],temp=0,fnum=0;
inti,j,m,n,e,bk=0;
printf（"使用Jacobi迭代法求解方程组:
\n"）;
printf（"输入方程组的元：
\nn="）;
scanf（"%d",&n）;
for（i=1;ix1[i]=0;
printf（"输入方程组的系数矩阵:
\n"）;
for（i=1;i{
j=1;
while（j{
scanf（"%f",&a[i][j]）;
j++;
}
}
printf（"输入方程组的常数项:
\n"）;
for（i=1;i{
scanf（"%f",&a[i][n+1]）;
}
printf（"\n"）;
printf（"请输入迭代次数:
\n"）;
scanf（"%d",&m）;
printf（"请输入迭代精度:
\n"）;
scanf（"%d",&e）;
while（m!
=0）
{
for（i=1;i{
for（j=1;j{
if（j!
=i）
temp=a[i][j]*x1[j]+temp;
}
x2[i]=（a[i][n+1]-temp）/a[i][i];
temp=0;
}
for（i=1;i{
fnum=float（fabs（x1[i]-x2[i]））;
if（fnum>temp）temp=fnum;
}
if（temp<=pow（10,-4））bk=1;
for（i=1;ix1[i]=x2[i];
m--;
}
printf（"原方程组的解为:
\n"）;
for（i=1;i{
if（（x1[i]-x2[i]）<=e||（x2[i]-x1[i]）<=e）
{
printf（"x%d=%7.4f",i,x1[i]）;
}
}
}
运行结果：
2、设对称正定阵系数阵线方程组
=
x
=（1,-1,0,2,1,-1,0,2）
GS迭代算法：
#include
#include
#include
constintm=11;
voidmain（）
{
intchoice=1;
while（choice==1）
{
doublea[m][m],b[m],e,x[m],y[m],w,se,max;
intn,i,j,N,k;
cout<<"Gauss-Seidol迭代法"<cout<<"请输入方程的个数：
";
cin>>n;
for（i=1;i<=n;i++）
{
cout<<"请输入第"<";
for（j=1;j<=n;j++）
cin>>a[i][j];
}
cout<<"请输入各个方程等号右边的常数项:
\n";
for（i=1;i<=n;i++）
{
cin>>b[i];
}
cout<<"请输入最大迭代次数：
";
cin>>N;
cout<<"请输入最大偏差：
";
cin>>e;
for（i=1;i<=n;i++）
{
x[i]=0;
y[i]=x[i];
}
k=0;
while（k!
=N）
{
k++;
for（i=1;i<=n;i++）
{
w=0;
for（j=1;j<=n;j++）
{
if（j!
=i）
w=w+a[i][j]*y[j];
}
y[i]=（b[i]-w）/double（a[i][i]）;
}
max=fabs（x[1]-y[1]）;
for（i=1;i<=n;i++）
{
se=fabs（x[i]-y[i]）;
if（se>max）
max=se;
}
if（max{
cout<for（i=1;i<=n;i++）
cout<<"x"<break;
}
for（i=1;i<=n;i++）
{
x[i]=y[i];
}
}
if（k==N）
cout<<"迭代失败！
！
"<choice=0;
}
}
3、三对角形线性方程组
=
x
=（2,1,-3,0,1,-2,3,0,1,-1）
SOR方法：
#include
#include
#include
/**********定义全局变量**********/
float**a;/*存放A矩阵*/
float*b;/*存放b矩阵*/
float*x;/*存放x矩阵*/
floatp;/*精确度*/
floatw;/*松弛因子*/
intn;/*未知数个数*/
intc;/*最大迭代次数*/
intk=1;/*实际迭代次数*/
/**********SOR迭代法**********/
voidSOR（floatxk[]）
{
inti,j;
floatt=0.0;
floattt=0.0;
float*xl;
xl=（float*）malloc（sizeof（float）*（n+1））;
for（i=1;i{
t=0.0;
tt=0.0;
for（j=1;jt=t+a[i][j]*xl[j];
for（j=i;jtt=tt+a[i][j]*xk[j];
xl[i]=xk[i]+w*（b[i]-t-tt）/a[i][i];
}
t=0.0;
for（i=1;i{
tt=fabs（xl[i]-xk[i]）;
tt=tt*tt;
t+=tt;
}
t=sqrt（t）;
for（i=1;ixk[i]=xl[i];
if（k+1>c）
{
if（t<=p）
printf（"\nReachthegivenprecision!
\n"）;
else
printf（"\noverthemaximalcount!
\n"）;
printf（"\nCountnumberis%d\n",k）;
}
else
if（t>p）
{
k++;
SOR（xk）;
}
else
{
printf（"\nReachthegivenprecision!
\n"）;
printf（"\nCountnumberis%d\n",k）;
}
}
/**********程序*****开始**********/
voidmain（）
{
inti,j;
printf（"SOR方法\n"）;
printf（"请输入方程个数:
\n"）;
scanf（"%d",&n）;
a=（float**）malloc（sizeof（float）*（n+1））;
for（i=0;ia[i]=（float*）malloc（sizeof（float）*（n+1））;
printf（"请输入三对角矩阵:
\n"）;
for（i=1;ifor（j=1;jscanf（"%f",&a[i][j]）;
for（i=1;ifor（j=1;jb=（float*）malloc（sizeof（float）*（n+1））;
printf（"请输入等号右边的值:
\n"）;
for（i=1;iscanf（"%f",&b[i]）;
x=（float*）malloc（sizeof（float）*（n+1））;
printf（"请输入初始的x:
"）;
for（i=1;iscanf（"%f",&x[i]）;
printf（"请输入精确度:
"）;
scanf（"%f",&p）;
printf（"请输入迭代次数:
"）;
sca
升级会员 