数值分析 课程设计.docx

1、数值分析 课程设计数值分析 课程设计求解线性方程组作者姓名：孙怡丰学号：200930980221指导教师 ：张昕 副教授学院名称理 学 院专 业 名 称统计学提交日期2011年6月 10 日一、 问题的提出分别用SOR方法和高斯消元的LU分解算法（lii=1, i=1,n）求解给定的线性方程组AX=B, 以感受迭代法和直接法的不同特点。二、 实验内容1. 自定义函数 SOR(A, B, w, MAXN, TOL)，以实现SOR方法求解线性方程组AX=B，其中A系数矩阵；B常数列向量；w松弛因子；MAXN迭代的最大次数TOL达到的精度上限返回值有以下四种可能：a) -2：SOR方法不收敛；（不收

2、敛的依据为的某个分量值超出区间-108, 108。）b) -1：矩阵有一列全为0；c) 0：算法经过MAXN次迭代还未收敛；d) k：SOR方法经k次迭代收敛，求得方程组的解向量X记录下来.2. 自定义函数Direct(A, B)，以实现高斯LU分解的方法求解线性方程组AX=B，其中A系数矩阵；B常数列向量；返回值有两种可能：a) “LU decompsition failed.”：分解过程中U的对角线元素至少一个为0；b) X：分解过程中3. 分别使用步骤1中定义的函数SOR(A, B, w, MAXN, TOL)和步骤2中定义的函数Direct(A, B)进行测试，记录返回值及X值（算法收

3、敛或有效的情形, 保留4位小数）:(1) 测试1：MAXN =1000，TOL =10-9，w分别取1， 1.05， 1.1， 1.2， 1.3， 1.6， 1.95；(2) 测试2：MAXN =1000，TOL =10-9，w=1；(3) 测试3：MAXN =1000，TOL =10-9，w=1.2；(4) 测试4：MAXN =1000，TOL =10-9，w=1， 1.1， 1.3， 1.8；(5) 测试5：: n阶Hilbert矩阵定义为取n=3, MAXN =1000，TOL =10-9，w=1, 1.3, 1.6, 1.9;测试6：A为4阶Hilbert矩阵， MAXN =10000

4、，TOL =10-6，w=1, 1.3, 1.6, 1.8, 1.9.三、 实验结果及分析SOR方法测试1：W=1W=1.05W=1.1W=1.2W=1.3W=1.6W=1.95测试2：测试3：测试4：w=1w=1.1w=1.3w=1.8测试5：w=1w=1.3w=1.6w=1.9测试6：w=1w=1.3w=1.6w=1.8w=1.9高斯LU方法测试1：测试2：测试3：测试4：测试5：测试6：四、 关于本设计的体会虽然在多数实验中，两种方法的结果是一样的，但是在通过松弛度w的调整中，SOR可能会得出结果，但是在精确度上一般会差点。高斯在数字的精确度上好点但是，其使用起来并不能全部算出结果，适用

5、的范围相对较小。所以，两种方法各有优点，在做题目时应结合两种方法，对结果比较后才可以做出正确的判断。五、 参考文献数值分析（第三版） 北京理工大学出版社六、 附录原程序：SOR:#include stdlib.h#include stdio.h#include conio.h#include string.h#include math.h#define N 100/初始的一个大矩阵/float cre_sch(int n,float *w,float aNN,float bN,float lwNN,float fwN) int i,j,k; float tmp1NN,tmp2NN,revNN,

6、tmp; for(i=0;in;i+) for(j=0;jn;j+) if(j=i) tmp1ij=aij; tmp2ij=(1-*w)*aij; revij=1; else if(ji) tmp1ij=*w*aij; tmp2ij=0; revij=0; else tmp1ij=0; tmp2ij=-*w*aij; revij=0; for(j=0;jn;j+) for(i=0;in;i+) for(k=0;k=j;k+) if(ij) continue; if(i=j) revik*=1/tmp1jj; continue; revik+=revjk*(-tmp1ij); for(i=0;i

7、n;i+) for(j=0;jn;j+) tmp=0.0; for(k=0;kn;k+) tmp+=revik*tmp2kj; lwij=tmp; for(i=0;in;i+) tmp=0.0; for(k=0;kn;k+) tmp+=*w*revik*bk; fwi=tmp; float Table(int n,float aNN,float bN,float *w) int i,j,k; float lwNN,fwN; printf(Please input the matrix A by row!n); for(i=0;in;i+) printf(Row %d:,i+1); k=0; f

8、or(j=0;jn;j+) scanf(%f,&aij); if(aij=0) k=k+1; if(k=n) printf(the matrix A has a line of 0); exit(0); printf(Please input the vector b:); for(i=0;in;i+) scanf(%f,&bi); printf(Input w:); scanf(%f,w); cre_sch(n,w,a,b,lw,fw); printf(nThe matrix A and vector b:n); for(i=0;in;i+) for(j=0;jn;j+) printf(%1

9、0.5f,aij); printf(%10.5f,bi); printf(n); printf(nThe SOR iterative scheme(matrix Lw & vector fw):n); for(i=0;in;i+) for(j=0;jn;j+) printf(%10.5f,lwij); printf(%10.5f,fwi); printf(n); float init_vec(int n,float xN) int i; printf(nPlease input the initial iteration vector x:); for(i=0;in;i+) scanf(%f,

10、&xi); printf(nThe initial iteration vector x:n); for(i=0;in;i+) printf(%10.5f,xi); printf(n);float sor(int n,float aNN,float bN,float xN,float w) int i,j,k; float tmp1,tmp2,x2N; for(k=0;k+) for(i=0;in;i+) x2i=xi; for(i=0;in;i+) tmp1=0.0; tmp2=0.0; for(j=0;ji;j+) tmp1+=aij*xj; for(j=i;jn;j+) tmp2+=ai

11、j*x2j; xi=x2i+w*(bi-tmp1-tmp2)/aii; for(i=0,j=0;in;i+) if(fabs(x2i-xi)0.000001) j+; /精度要求,可改变/ if(j=n) printf(nThis SOR iterative scheme is convergent!n); printf(Number of iterations: %d,k+1); break; if(k=10000)/最大迭代次数，可改变/ printf(This SOR iterative scheme may be not convergent!); break; printf(nThe

12、 results:n); for(i=0;in;i+) printf(%12.7f,xi); int main() int n; float xN,aNN,bN,w; printf(Input n:); scanf(%d,&n); Table(n,a,b,&w); init_vec(n,x); sor(n,a,b,x,w); 高斯LU：/* LU分解法 */#define AL 100 /局阵最大行#define ROW 100 /局阵最大列#include#include#include#includevoid main() int n;/阶 double AALROW;/系数矩阵 doub

13、le bAL;/右端项 int i,j; char flag=n;/选择标志 /. 文件读入 . /*coutflag; if(flag=y|flag=Y) coutfile; coutendl; ifstream fin(file); cout输入阶: n; coutnendl; cout输入系数矩阵A: endl; for(i=0;in;i+) for(j=0;jAij; coutAij ; coutendl; cout输入右端项b: endl; for(i=0;ibi; coutbi ; coutendl; coutendl; */ /. 手工输入 ./ else cout输入阶: n;

14、 cout输入系数矩阵A: endl; for(i=0;in;i+) for(j=0;jAij; cout输入右端项b: endl; for(i=0;ibi; coutendl;/ /.计算L矩阵和U矩阵(都放到A矩阵里). int k;/ int m; double sum_1=0; for(k=0;kn-1;k+)/循环/L矩阵算法 for(i=k+1;in;i+) sum_1=0; for(int m=0;mk;m+) sum_1+=Aim*Amk; Aik=(Aik-sum_1)/Akk; /U矩阵算法 for(i=k+1;in;i+) sum_1=0; for(int m=0;mk+

15、1;m+) sum_1+=Ak+1m*Ami; Ak+1i-=sum_1; /.判断矩阵对角线上是否为0. int wr; for(int num=0;numn;num+) if (Anumnum=0) wr=1+wr; if(wr=n)coutnLU decompsition failed. ;exit (0);elsecoutn分解过程中对角线上没有0; /. LUX=b . /第一步:UX=X_*,算LX_*=b,得到Y(即bi) for(i=0;i0;j-) sum_1+=Aij-1*bj-1; bi=bi-sum_1; /第二步:UY=b,得出答案bi for(i=n-1;i=0;i-) sum_1=0; for(j=i+1;jn;j+) sum_1=sum_1+Aij*bj; bi=(bi-sum_1)/Aii; coutn答案输出: ; for(i=0;in;i+) coutbi ; coutendl;/L矩阵和U矩阵都仍旧存放在A矩阵里/* for(i=0;in;i+) for(j=0;jn;j+) coutAij ; coutendl;*/ 七、 教师评价