数值分析实验报告.docx
《数值分析实验报告.docx》由会员分享，可在线阅读，更多相关《数值分析实验报告.docx（42页珍藏版）》请在冰豆网上搜索。
数值分析实验报告
数值分析上机实验报告
课题一线性方程组的直接算法
实验源程序
#include"stdio.h"
#include"math.h"
doubleA1[10][10]={4,2,-3,-1,2,1,0,0,0,0,
8,6,-5,-3,6,5,0,1,0,0,
4,2,-2,-1,3,2,-1,0,3,1,
0,-2,1,5,-1,3,-1,1,9,4,
-4,2,6,-1,6,7,-3,3,2,3,
8,6,-8,5,7,17,2,6,-3,5,
0,2,-1,3,-4,2,5,3,0,1,
16,10,-11,-9,17,34,2,-1,2,2,
4,6,2,-7,13,9,2,0,12,4,
0,0,-1,8,-3,-24,-8,6,3,-1,};
doubleA2[8][8]={4,2,-4,0,2,4,0,0,
2,2,-1,-2,1,3,2,0,
-4,-1,14,1,-8,-3,5,6,
0,-2,1,6,-1,-4,-3,3,
2,1,-8,-1,22,4,-10,-3,
4,3,-3,-4,4,11,1,-4,
0,2,5,-3,-10,1,14,2,
0,0,6,3,-3,-4,2,19,};
doubleA3[10][10]={4,-1,0,0,0,0,0,0,0,0,
-1,4,-1,0,0,0,0,0,0,0,
0,-1,4,-1,0,0,0,0,0,0,
0,0,-1,4,-1,0,0,0,0,0,
0,0,0,-1,4,-1,0,0,0,0,
0,0,0,0,-1,4,-1,0,0,0,
0,0,0,0,0,-1,4,-1,0,0,
0,0,0,0,0,0,-1,4,-1,0,
0,0,0,0,0,0,0,-1,4,-1,
0,0,0,0,0,0,0,0,-1,4,};
doubleB1[10]={5,12,3,2,3,46,13,38,19,-21,};
doubleB2[8]={0,-6,20,23,9,-22,-15,45,};
doubleB3[10]={7,5,-13,2,6,-12,14,-4,5,-5,};
voidlzy（doublea[10][10],doubleb[10],intn）
{
inti,j,k;
doublel,x[10],temp;
for（k=0;k{
for（j=k,i=k;j{
if（j==k）
temp=fabs（a[j][k]）;
elseif（temp{
temp=fabs（a[j][k]）;
i=j;
}
}
if（temp==0）
{
printf（"1无解!
\n"）;
return;
}
else
{
for（j=k;j{
temp=a[k][j];
a[k][j]=a[i][j];
a[i][j]=temp;
}
temp=b[k];
b[k]=b[i];
b[i]=temp;
}
for（i=k+1;i{
l=a[i][k]/a[k][k];
for（j=k;ja[i][j]=a[i][j]-l*a[k][j];
b[i]=b[i]-l*b[k];
}
}
if（a[n-1][n-1]==0）
{
printf（"2无解!
\n"）;
return;
}
x[n-1]=b[n-1]/a[n-1][n-1];
for（i=n-2;i>=0;i--）
{
temp=0;
for（j=i+1;jtemp=temp+a[i][j]*x[j];
x[i]=（b[i]-temp）/a[i][i];
}
for（i=0;i{
printf（"x%d=%lf\t",i+1,x[i]）;
printf（"\n"）;
}
}
voidpfg（doublea[8][8],doubleb[8],intn）
{
inti,k,m;
doublex[8],y[8],temp;
for（k=0;k{
temp=0;
for（m=0;mtemp=temp+pow（a[k][m],2）;
if（a[k][k]return;
a[k][k]=pow（（a[k][k]-temp）,1.0/2.0）;
for（i=k+1;i{
temp=0;
for（m=0;mtemp=temp+a[i][m]*a[k][m];
a[i][k]=（a[i][k]-temp）/a[k][k];
}
temp=0;
for（m=0;mtemp=temp+a[k][m]*y[m];
y[k]=（b[k]-temp）/a[k][k];
}
x[n-1]=y[n-1]/a[n-1][n-1];
for（k=n-2;k>=0;k--）
{
temp=0;
for（m=k+1;mtemp=temp+a[m][k]*x[m];
x[k]=（y[k]-temp）/a[k][k];
}
for（i=0;i{
printf（"x%d=%lf\t",i+1,x[i]）;
printf（"\n"）;
}
}
voidzgf（doublea[10][10],doubleb[10],intn）
{
inti;
doublea0[10],c[10],d[10],a1[10],b1[10],x[10],y[10];
for（i=0;i{
a0[i]=a[i][i];
if（ic[i]=a[i][i+1];
if（i>0）
d[i-1]=a[i][i-1];
}
a1[0]=a0[0];
for（i=0;i{
b1[i]=c[i]/a1[i];
a1[i+1]=a0[i+1]-d[i+1]*b1[i];
}
y[0]=b[0]/a1[0];
for（i=1;iy[i]=（b[i]-d[i]*y[i-1]）/a1[i];
x[n-1]=y[n-1];
for（i=n-2;i>=0;i--）
x[i]=y[i]-b1[i]*x[i+1];
for（i=0;i{
printf（"x%d=%lf\t",i+1,x[i]）;
printf（"\n"）;
}
}
voidmain（）
{
Printf（"线性方程组的解：
（Gauss列主元法）\n"）;
lzy（A1,B1,10）;
printf（"正定线性方程组的解：
（平方根法）\n"）;
pfg（A2,B2,8）;
printf（"三对角方程组的解：
（追赶法）\n"）;
zgf（A3,B3,10）;
}
实验结果：
结果分析
从运行结果可以看出，对方程组一和方程组三解得的结果和题中提供的理论值基本相符，解误差可以忽略，但方程组二采用平方根方法的到的结果与理论值相差甚远，原因是方程组的系数矩阵A的条件数Cond（A）非常大，导致系数矩阵A是一个病态矩阵，方程组是一个病态方程组。
课题三线性方程组的迭代法
实验源程序
#include
#include
usingnamespacestd;
#definemaxsize2000
//Jacobi迭代法所对应的函数
voidJacobi（doublee,intn,doubleA[][50],doubleb[]）
{
intk,j,i;
doublesum,y;
doublebestf=0;
doublex[2][10]={0,0,0,0,0,0,0,0,0,0};
for（k=0;k{bestf=0;
for（i=0;i{sum=0;
for（j=0;jsum+=A[i][j]*x[0][j];
x[1][i]=x[0][i]+（b[i]-sum）/A[i][i];
y=x[0][i]-x[1][i];
if（bestfbestf=fabs（y）;
}
for（i=0;ix[0][i]=x[1][i];
if（bestf<=e）
break;
k++;
}
for（i=0;icout<cout<cout<<"迭代了"<"<}
//GaussSeidol所对应的函数
voidGaussSeidol（doublee,intn,doubleA[][50],doubleb[]）
{
intk,j,i;
doublesum,y;
doublebestf=0;
doublex[2][10]={0,0,0,0,0,0,0,0,0,0};
for（k=0;k{bestf=0;
for（i=0;i{sum=0;
for（j=0;jsum+=A[i][j]*x[0][j];
x[1][i]=x[0][i]+（b[i]-sum）/A[i][i];
y=x[0][i]-x[1][i];
if（bestfbestf=fabs（y）;
x[0][i]=x[1][i];
}
k++;
if（bestf<=e）
break;
}
for（i=0;icout<cout<cout<<"迭代了"<"<}
//sor迭代法所对应的函数
voidsor（doublew,doublee,intn,doubleA[][50],doubleb[]）
{
intk,j,i;
doublesum,y;
doublebestf=0;
doublex[2][10]={0,0,0,0,0,0,0,0,0,0};
for（k=0;k{bestf=0;
for（i=0;i{sum=0;
for（j=0;jsum+=A[i][j]*x[0][j];
x[1][i]=x[0][i]+w*（b[i]-sum）/A[i][i];
y=x[0][i]-x[1][i];
if（bestfbestf=fabs（y）;
x[0][i]=x[1][i];
}
k++;
if（bestf<=e）
break;
}
for（i=0;icout<cout<cout<<"迭代了"<"<}
//主函数
voidmain（）
{
intn,i,j;
doubleA[50][50];
doubleb[50];
cout<<"请输入数组维数：
"<cin>>n;
cout<<"请输入数组数据：
"<for（i=0;ifor（j=0;jcin>>A[i][j];
cout<<"请输入右端项：
"<for（i=0;icin>>b[i];
doublew,e;
cout<<"请输入精度：
"<cin>>e;
cout<<"请输入sor的松弛因子：
"<cin>>w;
cout<<"Jacobi结果：
"<Jacobi（e,n,A,b）;
cout<<"Gauss-Seidol结果：
"<GaussSeidol（e,n,A,b）;
cout<<"SOR结果：
"<sor（w,e,n,A,b）;
}
运行结果：
１
２：
精度同为0.001，松弛因子为1.4的ＳＯＲ结果迭代次数比松弛因子为1.5时增加了16次。
松弛因子同为1.5，精度为0.0001的ＳＯＲ结果的迭代次数比精度为0.001时多了144
３：
结果分析
1.体会迭代法求解线性方程组，并能与消去法做以比较
消去法是一种求解的直接方法，不考虑计算过程中的舍入误差，运用此类方法经过有限次算术运算就能求出线性方程组的精确解，但由于舍入误差的存在，一般求得的也是近似解。
而迭代法与直接方法不同，即使无舍入误差，迭代法也难获得精确解，是一种逐次近似的方法，适于解大型稀疏线性方程组。
2.分别对不同精度要求，如
由迭代次数体会该迭代法的收敛快慢
由于在程序运行结果中采用的是精度e=0.001，在下面主要是计算了e=0.0001和e=0.00001并主要针对第二个方程来说明
当e=0.001时，三种不同迭代方法的迭代次数依次为1069，133，14
当e=0.0001时，三种不同迭代方法的迭代次数依次为1069，570，545
当e=0.00001时，三种不同迭代方法的迭代次数依次为1069，1006，1079
针对上述数据，可以得出随着精度要求的提高，迭代次数依次增大，收敛的越慢。
3.对方程组2，3使用SOR方法时，选取松弛因子
=0.8，0.9，1，1.1，1.2等，试看对算法收敛性的影响，并能找出你所选用的松弛因子的最佳者
由于
=0.8在程序运行结果中，已经给出了其他几组数的运行结果
对于方程组2，
=0.7时，迭代次数为18
=0.8时，迭代次数为13
=0.9时，迭代次数为14
=1时，迭代次数为133
=1.1时，迭代次数为197
=1.2时，迭代次数为228
发现w=0.8时，迭代次数最小，为13次。
对于方程组3，
=0.8时，迭代次数为9
=0.9时，迭代次数为8
=1时，迭代次数为9
=1.1时，迭代次数为9
=1.2时，迭代次数为12
发现
=0.9时，迭代次数最小，为8次。
课题六函数插值方法
5次Lagrange插值多项式
#include
#include
usingnamespacestd;
#defineN5
floatx[]={0.4,0.55,0.65,0.80,0.95,1.05};
floaty[]={0.41075,0.57815,0.69675,0.90,1.00,1.25382};
floatp（floatxx）
{
inti,k;
floatpp=0,m1,m2;
for（i=0;i<=N;i++）
{
m1=1;m2=1;
for（k=0;k<=N;k++）
if（k!
=i）{
m1*=xx-x[k];
m2*=x[i]-x[k];}
pp+=y[i]*m1/m2;
}
returnpp;
}
voidmain（）
{
floatx,y;
cout<<"请输入x,y的值：
"<cin>>x;
cin>>y;
cout<<"结果为：
"<cout<<"f（"<cout<<"f（"<}
运行结果
课题九数值积分
源程序
#include
#include
usingnamespacestd;
//问题1所对应的梯形公式
doubleT1（floata,floatb,intn）
{
doublex,sum;
doubleh;
intk=0;
h=（b-a）/n;
sum=sqrt（4-sin（a）*sin（a））;
for（k=1;k{
x=a+k*h;
sum+=2*sqrt（4-sin（x）*sin（x））;
}
sum+=sqrt（4-sin（b）*sin（b））;
sum=（h/2）*sum;
returnsum;
}
//问题1所对应的Simpson公式
doubleS1（floata,floatb,intn）
{
doublex,sum=0;
doubleh;
intk=0;
h=（b-a）/n;
sum=sqrt（4-sin（a）*sin（a））;
for（k=1;k{
x=a+（k-0.5）*h;
sum+=4*sqrt（4-sin（x）*sin（x））;
x=a+k*h;
sum+=2*sqrt（4-sin（x）*sin（x））;
}
x=a+（n-0.5）*h;
sum+=4*sqrt（4-sin（x）*sin（x））;
sum+=sqrt（4-sin（b）*sin（b））;
sum=（h/6）*sum;
returnsum;
}
//问题1所对应的Romberg公式
doubleR1（doublee,doublea,doubleb）
{
intk,j,i;
doublec,bestf=1000;
doublesum;
doublet[2][10]={0};
c=b-a;
t[0][0]=c*（1+sqrt（4-sin（b）*sin（b）））/2;
k=1;
while（bestf>e）
{
sum=0;
for（i=1;i<=pow（2,k-1）;i++）
sum+=sqrt（4-sin（a+（2*i-1）*c/pow（2,k））*sin（a+（2*i-1）*c/pow（2,k）））;
t[1][0]=0.5*t[0][0]+c*sum/pow（2,k）;
for（j=1;j<=k;j++）
t[1][j]=（pow（4,j）*t[1][j-1]-t[0][j-1]）/（pow（4,j）-1）;
bestf=t[1][j-1]-t[0][j-1];
bestf=fabs（bestf）;
for（j=0;j<=k;j++）
t[0][j]=t[1][j];
k++;
}
sum=t[1][k-2];
returnsum;
}
//问题2所对应的梯形公式
doubleT2（floata,floatb,intn）
{
doublex,sum;
doubleh;
intk=0;
h=（b-a）/n;
sum=1;
for（k=1;k{
x=a+k*h;
sum+=2*sin（x）/x;
}
sum+=sin（b）/b;
sum=（h/2）*sum;
returnsum;
}
//问题2所对应的Simpson公式
doubleS2（floata,floatb,intn）
{
doublex,sum=0;
doubleh;
intk=0;
h=（b-a）/n;
sum=1;
for（k=1;k{
x=a+（k-0.5）*h;
sum+=4*sin（x）/x;
x=a+k*h;
sum+=2*sin（x）/x;
}
x=a+（n-0.5）*h;
sum+=4*sin（x）/x;
sum+=sin（b）/b;
sum=（h/6）*sum;
returnsum;
}
//问题2所对应的Romberg公式
doubleR2（doublee,doublea,doubleb）
{
intk,j,i;
doublec,bestf=1000;
doublesum;
doublet[2][10]={0};
c=b-a;
t[0][0]=c*（1+sin（b）/（b））/2;
k=1;
while（bestf>e）
{
sum=0;
for（i=1;i<=pow（2,k-1）;i++）
sum+=sin（a+（2*i-1）*c/pow（2,k））/（a+（2*i-1）*c/pow（2,k））;
t[1][0]=0.5*t[0][0]+c*sum/pow（2,k）;
for（j=1;j<=k;j++）
t[1][j]=（pow（4,j）*t[1][j-1]-t[0][j-1]）/（pow（4,j）-1）;
bestf=t[1][j-1]-t[0][j-1];
bestf=fabs（bestf）;
for（j=0;j<=k;j++）
t[0][j]=t[1][j];
k++;
}
sum=t[1][k-2];
returnsum;
}
//问题3所对应的梯形公式
doubleT3（floata,floatb,intn）
{
doublex,sum;
doubleh;
intk=0;
h=（b-a）/n;
sum=exp（a）/（4+a*a）;
for（k=1;k{
x=a+k*h;
sum+=2*exp（x）/（4+x*x）;
}
sum+=exp（x）/（4+x*x）;
sum=（h/2）*sum;
returnsum;
}
//问题3所对应的Simpson公式
do
升级会员 