数值分析实验报告.docx
- 文档编号:9312463
- 上传时间:2023-02-04
- 格式:DOCX
- 页数:42
- 大小:767.95KB
数值分析实验报告.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;j a[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;j temp=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;m temp=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;m temp=temp+a[i][m]*a[k][m]; a[i][k]=(a[i][k]-temp)/a[k][k]; } temp=0; for(m=0;m temp=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;m temp=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(i c[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;i y[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;j sum+=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(bestf bestf=fabs(y); } for(i=0;i x[0][i]=x[1][i]; if(bestf<=e) break; k++; } for(i=0;i cout< 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;j sum+=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(bestf bestf=fabs(y); x[0][i]=x[1][i]; } k++; if(bestf<=e) break; } for(i=0;i cout< 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;j sum+=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(bestf bestf=fabs(y); x[0][i]=x[1][i]; } k++; if(bestf<=e) break; } for(i=0;i cout< cout< cout<<"迭代了"< "< } //主函数 voidmain() { intn,i,j; doubleA[50][50]; doubleb[50]; cout<<"请输入数组维数: "< cin>>n; cout<<"请输入数组数据: "< for(i=0;i for(j=0;j cin>>A[i][j]; cout<<"请输入右端项: "< for(i=0;i cin>>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); } 运行结果: 1 2: 精度同为0.001,松弛因子为1.4的SOR结果迭代次数比松弛因子为1.5时增加了16次。 松弛因子同为1.5,精度为0.0001的SOR结果的迭代次数比精度为0.001时多了144 3: 结果分析 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
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 数值 分析 实验 报告