当前位置:首页 >> 工学 >>

数值分析第四次实习报告


学生学号

实验课成绩

武汉理工大学
学 生 实 验 报 告 书

实验课程名称: 实验课程名称: 数 值 分 析(

第四次实习)

开 课 名 称:计算机科学与技术学院 指导老师姓名: 指导老师姓名: 熊盛武 学 生 姓 名: 学生专业班级: 学生专业班级:软件工程 0803 班

2009 —— 2010 学年

第 一 学期

一、计算实习的内容: 计算实习的内容:
实习 1
给定矩阵 A 与向量 b

0? ?n ? ? ? n ?1 n ? ? M O M? ? ? 3 0? ? 2 ? 1 2 n? ? ? ? ? ? ? A=

?1 ? ? ? ? 0? ? 0? ? ? ?M ? ? 0? ? ? ? ? b= ? 0 ?

(1)求 A 的三角分解(不必输出) ; (2)利用 A 的三角分解解下列方程组: ① Ax=b ②A2x=b ③A3x=b 对第③题分析一下,如果先求 M=A2,再解 Mx=b,有何缺点?

实习 2

追赶法的优点

分别对 n=5,100,300 解下列方程组 Ax=b

1 ?2 ? ? ? 2 1 ?1 ? ? O O O ? ? ? 1 2 1? ? ? ? 1 2? A= ?

?? ? ?? ? ? b= ? ?

7? ? 5? 0? ? 5?

再用现有的三角分解法程序解此方程组,比较一下两者的效率,哪种方法机时少,可求解的 矩阵阶数高?

实习 3

Jacobi 迭代法与 Gauss-Seidel 迭代法的收敛性与收敛速度

研究用 Jacobi 迭代法与 Gauss-Seide 迭代法解下列方程组 Ax=b 的收敛性, 通过上机计算, 验证分析是否正确,并观察右端项对迭代收敛是否有影响,比较两法的收敛速度。

? 6 2 ?1? ?? 3 ? ? ? ? 1 4 ? 2? ? 2 ?? 3 1 4 ? ? ? b1= ? 4 (1)A= ? ? 1 0.8 0.8 ? ? ? ? 0.8 1 0.8 ? ? 0.8 0.8 1 ? ? (2) A= ?
? 1 3? ? ? (3) A= ? ?7 1 ?

? ? ? ? ?

? 100 ? ? ? ? ? 200 ? ? 345 ? ? b2= ? ? 5 ? ? ? ? 0 ? ? ?10 ? ? b2= ?

?3? ? ? ? 2? ?1 ? b1= ? ?

? 4? ? ? b= ? ?6 ?

实习 4

松弛因子对超松弛迭代法收敛速度的影响

编写一个超松弛迭代法解方程组 Ax=b 的计算机程序,其中

1 ? ?4 ? ? ? ?4 1 ?1 ? ? ? O O O ? ? 1 ?4 1? ? ? ? 1 ?4 ? A= ?
x(k ) ? x ( k ?1)
?6

? ?3 ? ? ? ? ?2 ? ? M ? ? ? ? ?2 ? ? ? b= ? ?3 ?

分别对不同的阶数(例如 n=10,n=100)w=1.1,1.2,……,1.9 进行迭代,记录近似解 <e= 10 时所用的迭代次数 k,观察松弛因子对收敛速度的影响。

实习 5

用欧拉公式求解

二、计算实习: 计算实习:
(一)三角分解法程序源代码
#include<iostream> using namespace std; int main() { const int MAX_N=20; static double a[MAX_N][MAX_N],b[MAX_N],x[MAX_N],y[MAX_N]; static double u[MAX_N][MAX_N],l[MAX_N][MAX_N]; int i,j,k,n,r; cout<<"请输入矩阵阶数:"; cin>>n; if(n>MAX_N) { cout<<"The input n is larger than MAX_N,please redefine the MAX_N"<<endl; return 1; } if(n<=0) { cout<<"please input n between 1 and"<<MAX_N; return 1; } //输入 a[i][j],b[i] cout<<"now input the matrix a[i][j],i,j=0..."<<n-1<<" :"; for(i=0;i<n;i++) {

for(j=0;j<n;j++) cin>>a[i][j]; } cout<<"now input the matrix b[i],i=0..."<<n-1<<" :"; for(i=0;i<n;i++) cin>>b[i]; //分解过程 for(i=0;i<n;i++) l[i][i]=1;//l 矩阵对角元素为 1 for(k=0;k<n;k++) { for(j=0;j<n-1;j++)//计算 u 矩阵 { u[k][j]=a[k][j]; for(r=0;r<=k-1;r++) u[k][j]-=l[k][r]*u[r][j]; } for(i=k+1;k<n;i++)//计算 l 矩阵 { l[i][k]=a[i][k]; for(r=0;r<=k-1;r++) l[i][k]-=l[i][r]*u[r][k]; l[i][k]/=u[k][k]; } } //ly=b for(i=0;i<n;i++) { y[i]=b[i]; for(j=0;j<=i-1;j++) y[i]-=l[i][r]*y[j]; } //ux=y for(i=n-1;i>=0;i--) { x[i]=y[i]; for(j=i+1;j<n;j++) x[i]-=u[i][j]*x[j]; x[i]/=u[i][i]; } cout<<"方程组的解为:";//输出 x for(i=0;i<n;i++) cout<<x[i]<<" "; return 0; }

(二)追赶法程序源代码
#include<iostream> using namespace std; int main() { const int MAX_N=300; int a[MAX_N],b[MAX_N],c[MAX_N],d[MAX_N]; float l[MAX_N],r[MAX_N],y[MAX_N],x[MAX_N]; int i,n; cout<<"请输入矩阵阶数:"; cin>>n; if(n>MAX_N) { cout<<"The input n is larger than MAX_N,please redefine the MAX_N"<<endl; return 1; } if(n<=0) { cout<<"please input n between 1 and"<<MAX_N; return 1; } //输入 a[n],b[n],c[n],d[n] cout<<"now input the matrix a[n],n=2..."<<n<<" :"; for(i=2;i<=n;i++) cin>>a[i]; cout<<"now input the matrix b[n],n=1..."<<n<<" :"; for(i=1;i<=n;i++) cin>>b[i]; cout<<"now input the matrix c[n],n=1..."<<n-1<<" :"; for(i=1;i<n;i++) cin>>c[i]; cout<<"now input the matrix d[n],n=1..."<<n<<" :"; for(i=1;i<=n;i++) cin>>d[i]; //追的过程 r[0]=0; y[0]=0; a[1]=0; c[n]=0; for(i=1;i<=n;i++) { l[i]=b[i]-a[i]*r[i-1]; r[i]=c[i]/l[i];

y[i]=(d[i]-a[i]*y[i-1])/l[i]; } //赶的过程 cout<<endl; x[n]=y[n]; for(i=n-1;i>=1;i--) { x[i]=y[i]-r[i]*x[i+1];\ } cout<<"方程组的解为:"; for(i=1;i<=n;i++) { cout<<x[i]<<" "; } cout<<endl; return 0; } 计算得解为:-2 -3 -3 -3 -2

(三)雅可比迭代法和高斯赛德尔迭代法程序源代码
#include <iostream> #include <cmath> using namespace std; int a,b,m; double *x0; void Jacobi(double **c,double *d,int n,double eps); void Gauss(double **c,double *d,int n,double eps); void main() { int n; double **A,*B; double e; cout<<"请选择求方程组的迭代方法!雅可比选 0,高斯-赛德尔选 1!"<<endl; cin>>n; cout<<"输入方程组的变量的个数以及方程的个数!"<<endl; cin>>a>>b; A=new double*[b]; for(int i=0;i<b;i++) { A[i]=new double[a]; } B=new double[b]; x0=new double[a]; cout<<"输入每个方程组的变量的系数以及方程右端的值!"<<endl;

for(int k=0;k<b;k++) { for(int j=0;j<a;j++) { cin>>A[k][j]; } cin>>B[k]; } cout<<"输入方程组迭代的次数及所要求的精度!"<<endl; cin>>m>>e; cout<<"输入方程组迭代的初值!"<<endl; for(int j=0;j<a;j++) { cin>>x0[j]; } switch (n) { case 0:Jacobi(A,B,m,e); break; case 1:Gauss(A,B,m,e); break; default:cout<<"你没有选择求解方程组的一种方法!!"<<endl; break; } } void Jacobi(double **c,double *d,int n,double eps) { int k,i; double *y = new double[a],*x=new double[a],s,temp=0.0; k=1; while(1) { temp = 0.0; for(i=0;i<a;i++) { s=0.0; for(int j=0;j<a;j++) { if(j!=i) { s+=c[i][j]*x0[j]; } }

s=(d[i]-s)/c[i][i]; y[i]=s; if(fabs(x0[i]-s)>temp) { temp=fabs(x0[i]-s); } } if(temp<eps) { cout<<"迭代成功!迭代结果为:"<<endl; for(i=0;i<a;i++) { cout<<"y["<<i<<"] ="<<y[i]<<endl; } break; } if(k==m) { cout<<"迭代失败!!"<<endl; break; } k+=1; for(i=0;i<a;i++) { x0[i]=y[i]; } } } void Gauss(double **c,double *d,int n,double eps) { int k,i; double *y=new double[a],*x=new double[a],s,temp=0.0; for(i=0;i<a;i++) { x[i]=x0[i]; y[i]=x[i]; } k=1; while(1) { temp=0.0; for(i=0;i<a;i++) { s=0.0;

for(int j=0;j<a;j++) { if(j!=i) { s+=c[i][j]*y[j]; } } s=(d[i]-s)/c[i][i]; y[i]=s; if(fabs(x[i]-s)>temp) { temp=fabs(x[i]-s); } } if(temp<eps) { cout<<"迭代成功!迭代结果为:"<<endl; for(i=0;i<a;i++) { cout<<"y["<<i<<"] ="<<y[i]<<endl; } break; } if(k==m) { cout<<"迭代失败!!"<<endl; break; } k+=1; for(i=0;i<a;i++) { x[i]=y[i]; } } }

用雅可比迭代法,程序运行结果如下: 用雅可比迭代法,程序运行结果如下:

用高斯赛德尔迭代法程序运行结果如下: 用高斯赛德尔迭代法程序运行结果如下

通过运行结果可以看出高斯赛德尔迭代法比雅可比迭代法所求结果精度高,收敛性好



(四)超松弛迭代法程序源代码
#include<iostream> using namespace std; #include<math.h> #define MAX_N 20 #define MAXREPT 100 #define epsilon 0.00001 int main() { int n; int i,j,k; double err,w;

static double a[MAX_N][MAX_N],b[MAX_N][MAX_N],c[MAX_N],g[MAX_N]; static double x[MAX_N],nx[MAX_N]; cout<<"input n value (dim of Ax=c):";//输入方程的维数 cin>>n; if(n>MAX_N) { cout<<"The input n is larger than MAX_N,please redefine the MAX_N"<<endl; return 1; } if(n<=0) { cout<<"please input n between 1 and"<<MAX_N; return 1; } //输入a[i][j],c[i] cout<<"now input the matrix a[i][j],i,j=0..."<<n-1<<" :"; for(i=0;i<n;i++) { for(j=0;j<n;j++) cin>>a[i][j]; } cout<<"now input the matrix c[i],i=0..."<<n-1<<" :"; for(i=0;i<n;i++) cin>>c[i]; cout<<"now input the w value:"); cin>>w; if(w<=1||w>=2) { cout<<"w must between 1 and 2."<<endl; return 1; } for(i=0;i<n;i++)//形成x_{k+1}=bx_{k}+g迭代矩阵b for(j=0;j<n;j++) { b[i][j]=-a[i][j]/a[i][i]; g[i]=c[i]/a[i][i];//为了简化程序,假设a[i][i]!=0 //否则要附加对a[i][i]的处理 } for(i=0;j<MAXREPT;i++) { for(j=0;j<n;j++) nx[j]=g[j]; for(j=0;j<n;j++) {

for(k=0;k<n;k++) { if(j==k)continue; nx[j]+=b[j][k]*nx[k];//迭代 for(k=j+1;k<n;k++) nx[j]+=b[j][k]*x[k]; nx[j]=(1-w)*x[j]+w*nx[j]; } } err=0; for(j=0;j<n;j++) if(err<fabs(nx[j]-x[j]))err=fabs(nx[j]-x[j]); for(j=0;j<n;j++) x[j]=nx[j]; if(err<epsilon) { cout<<"the solve are:"; for(i=0;i<n;i++) cout<<x[i]<<" "; cout<<endl; return 0; } } cout<<"after "<<MAXREPT<<"repeat,no result...";//输出 return 1; } ω =1.2 时 解为:0232179 0.0662066 0.494194 ω =1.4 时 解为:023218 0.0662066 0.494194 ω =1.6 时 解为:023218 0.0662074 0.494193 由以上 3 组结果可知, ω 越小,精度越大

(五)欧拉公式程序源代码
#include<iostream.h> #include<math.h> #define f(x,y) (x+y) int main() { int m; int i; double a,b,y0; double xn,yn,xn1,yn1,yn1b;

double h; cout<<"inpiut the begin and end of x:"; cin>>a; cin>>b; cout<<"input the y value:"; cin>>y0; cout<<"input m value:"; cin>>m; if(m<=0) { cout<<"please input a number lager than 1."; return 1; } h=(b-a)/m; xn=a;yn=y0; for(i=1;i<=m;i++) { xn1=xn+h; yn1b=yn+h*f(xn,yn); yn1=yn+h/2*(f(xn,yn)+f(xn1,yn1b)); cout<<"x"<<i<<" "<<xn1<<" "<<"y"<<i<<" "<<yn1; xn=xn1;yn=yn1; } return 0; }

调试结果: 调试结果:

x1=0.2 x2=0.4 x3=0.6 x4=0.8 x5=1

y1=1.24 y2=1.5768 y3=2.031 y4=2.63067 y5=3.40542

第四次实验小结及体会
1、任何非奇异矩阵都可分解为三角形式。 2、如果需要重复地求解系数矩阵相同,而右端常数项不同的线性方程组时,三角分解法是很 有效的 3、当系数矩阵是对角占优矩阵时,用追赶法求线性代数方程组是非常简单的。 4、高斯迭代法公式比较复杂,但在一般情况下比雅可比迭代法收敛快。 5、高斯迭代法是异步迭代法每次迭代时都会用到刚才已经迭代出来的数值,因此迭代对收敛 情况有所改善。 6、迭代法是否收敛,与迭代矩阵密切相关,在用迭代法进行迭代计算时应先判断系数矩阵的 谱半径的绝对值与 1 的大小关系,从而可以决定迭代法的收敛性。 7、超松弛迭代法的迭代矩阵与松弛因子 w 有关。 8、欧拉公式精度低但其数值解法给我们很好的启示。


赞助商链接
相关文章:
数值分析第二次实习报告
数值分析第次实习报告_理学_高等教育_教育专区。实验报告 一、实习目的:(1)通过编程计算实践,体会和理解 Lagrange 插值公式、Newton 插值公式、 分段插值公式和 ...
数值分析实习报告
数值分析实习报告 - 数值分析实习报告 数值分析实习报告 题 目: MATLAB 算法实习 学生姓名: 指导教师: 学学院: 号: 刘圣军 2016 年 12 月 数值分析实...
数值分析实习报告
数值分析实习报告 - 数值分析实习报告 数值分析实习报告 题 目: MATLAB 算法实习 学生姓名: 指导教师: 学学院: 号: 刘圣军 2016 年 12 月 数值分析实...
数值分析实习报告
暂无评价|0人阅读|0次下载|举报文档数值分析实习报告_其它_高等教育_教育专区。《数值分析》编程实习报告 姓学 名号 Agunk·D·Ragon 院(所) 年专级业 XD 二...
数值分析实习报告
暂无评价|0人阅读|0次下载|举报文档数值分析实习报告_工学_高等教育_教育专区。...8 . 实验三:经典的四阶Runge-Kutta法′ 3 用经典的四阶Runge-Kutta法解初值...
西南交大数值分析上机实习报告
西南交大数值分析上机实习报告_理学_高等教育_教育专区。研究生课程上机实习作业 ...计算过程便于改变步长, 缺点是计算量较大,每前进一步需要计算四次函数值 f 。...
数值分析上机实习报告
[j]); } } } 第 17 页 数值分析上机实习报告 2012 年 12 月 (3)结果截屏 (4) MATLAB 上机程序 unction [x]=mgauss2(A,b,flag) if nargin<3,flag...
2014级硕士研究生数值分析上机实习报告(答案)
(A,x,b,1e-10); } -9- 哈尔滨工业大学(威海)实验报告纸 2014 级硕士研究生数值分析上机实习 (第四次) 姓名: 学号: 学院: 2 实习题目: 分别用复化...
数值分析上机实习报告
4. 附录 4.1 问题一的 matlab 程序代码 clc;clear; syms x; a=input('按...('次多项式的方差为:');disp(S(i)); end end 数值分析上机实习报告要求 ...
数值分析上机实习报告
数值分析上机实习报告 姓名: 学号: 专业: 联系电话: 目 录 序言......[2,4]之间的任意 3 个数\n"); for(i=1;i<4;i++) { printf("第%d 个数:...
更多相关标签: