资源描述:
《数值分析课程设计--松弛迭代法中松弛因子》由会员上传分享,免费在线阅读,更多相关内容在学术论文-天天文库。
1、安徽建筑大学数值分析设计报告书题目松弛迭代法中松弛因子院系数理系专业信息与计算科学班级信息②班学号12207210220姓名穆海山时间2013-12-10~2013-12-23指导教师刘华勇题目:选用Jacobi迭代法、Gauss-Seidel迭代法和超松弛迭代法求解下面的方程组(考虑等于150)=考虑初值的变化和松弛因子的变化收敛效果的影响;对上述方程组还可以采用哪些方法求解?选择其中一些方法编程上机求解上述方程组,说明最适合的是什么方法;将计算结果进行比较分析,谈谈你对这些方法的看法。一、摘要本课程设计用matlab就线性方程组数值方法,Jacobi迭代法,Gauss-Seidel迭代
2、法,超松弛法对所设计的问题进行求解,并编写程序在Matlab中实现,在文章中对各种迭代法进行了收敛性分析。接着用几种不同方法对线性方程组进行求解及结果分析,最后对此次课程设计进行了总结。关键词:线性方程组,迭代,Matlab,结果分析二、设计目的用熟悉的计算机语言编程上机求解线性方程组。三、理论基础对方程组做等价变换如:令,则则,我们可以构造序列若同时:所以,序列收敛,与初值的选取无关则转化为矩阵形式(1)令或者故迭代过程(1)化为(2)等价线性方程组为称(2)式为解线性方程组(1)的Jacobi迭代法(J法)迭代矩阵考虑迭代式(2)即将上式改为(3)(4)超松弛迭代记则可以看作在前一步上
3、加一个修正量。若在修正量前乘以一个因子有对Gauss-Seidel迭代格式迭代矩阵四、程序代码及运算结果1.利用Jacobi迭代法求解:编制名为majacobifun.m的文件,内容如下:function[x,n]=jacobifun(A,b,x0,eps)%jacobifun为编写在雅可比迭代函数%A为线性方程组的系数矩阵%b为线性方程组的常数向量%x0为迭代初始向量%eps为解的精度%x为线性方程组的解%n为求出所需精度的解实际的迭代步数D=diag(diag(A));B=D(D-A);f=Db;x=B*x0+f;n=0;ifmax(abs(eig(B)))>=1disp('War
4、ning:不收敛!');return;endwhilenorm(x-x0)>=epsx0=x;x=B*x0+f;n=n+1;end2.利用Gauss-Seidel迭代法求解:编制名为G_Seidelifun.m的文件,内容如下:nction[x,n]=G_Seidelifun(A,b,x0,eps)%G_Seidelifun为编写在高斯-赛德尔迭代函数%A为线性方程组的系数矩阵%b为线性方程组的常数向量%x0为迭代初始向量%eps为解的精度%x为线性方程组的解%n为求出所需精度的解实际的迭代步数D=diag(diag(A));L=-tril(A,-1);U=-triu(A,1);G=(D-
5、L)U;f=(D-L)b;x=G*x0+f;n=0;ifmax(abs(eig(G)))>=1disp('Warning:不收敛!');return;endwhilenorm(x-x0)>=epsx0=x;x=G*x0+f;n=n+1;end3.利用SOR迭代法求解:编制名为SORfun.m的文件,内容如下:function[x,n]=SORfun(A,b,x0,w,eps)%SORfun为编写在松弛迭代函数%A为线性方程组的系数矩阵%b为线性方程组的常数向量%x0为迭代初始向量%w为松弛因子%eps为解的精度%x为线性方程组的解%n为求出所需精度的解实际的迭代步数if(w<=0
6、
7、w
8、>=2)errorreturn;endD=diag(diag(A));L=-tril(A,-1);U=-triu(A,1);S=(D-w*L)((1-w)*D+w*U);f=w*((D-w*L)b);x=S*x0+f;n=0;ifmax(abs(eig(S)))>=1disp('Warning:不收敛!');return;endwhilenorm(x-x0)>=epsx0=x;x=S*x0+f;n=n+1;end在Matlab命令窗口输入:>>clearw=1.5;eps=1e-6;n=150;A=diag(8*ones(1,n-1),-1)+diag(6*ones(1,n))+dia
9、g(ones(1,n-1),1);b(1)=10.5;b(n)=21;fori=2:n-1b(i)=22.5;endb=b';x0=zeros(n,1);[x1,n]=jacobifun(A,b,x0,eps);[x2,n]=G_Seidelifun(A,b,x0,eps);[x3,n]=SORfun(A,b,x0,w,eps);plot(x1);holdonplot(x2,'r');holdonplot(x3,'g*'