三对角线性方程组的解法ppt课件.ppt
1,数值分析课程设计 题目：三对角线性方程组的解法 指导老师： 罗 婷,海 南 大 学,2,目 录,一、实验分析 3二、实验内容 4三、实验要求 5四、实验过程 6五、结果分析 19六、程序代码 21七、课程信息 30,3,一、实验分析,关于线性方程组的数值解法一般有两类：直接法和迭代法。而对于系数矩阵较特殊的一类线性方程组，解法也相对有所不同。在实际问题中，例如解常微分方程边值问题，解热传导问题等，都会要求解系数矩阵为对角占优的三对角方程组。如何求解三对角方程组就成为关键。三对角线性方程组是系数矩阵为三对角矩阵的一类方程组，其数值解法有：追赶法，Jacobi迭代法，Gauss迭代法和SOR法。,4,考虑线性方程组：其中A为三对角矩阵，编制程序求解该方程组,二、实验内容,5,（1）考虑不同的数值解法，对方程组进行编程，编写其通用代码；（2）对不同的程序代码用实例进行验证，取矩阵： 方程初值 ： 。 （3）比较不同方法对求解三对角线性方程组的优异性。,三、实验要求,6,四、实验过程,解法,追赶法,Jacobi迭代法,Gauss迭代法,SOR迭代法,7,编程思想： 矩阵的三对角线性方程组是由矩阵的直接三角分解法来推到的来的，求解等价于解两个三角形方程组： （1） ，求y；（2） ，求x 从而得到解三对角线方程组的追赶法公式。 1. 计算的递推公式： （i=2,3,n-1）； 2. 解 （i=2,3,n）； 3. 解 （i=n-1，n-2，2,1）。,(一)追赶法,8,运行结果：（点击查看源程序） a=-1*ones(1,4); c=a; b=2*ones(1,5); f=1 0 1 0 0; x=ZhuiGan(a,b,c,f) n =5 x = 1.33333333333333 1.66666666666667 2.00000000000000 1.33333333333333 0.66666666666667注：其中a为对角下向量，b为对角向量，c为对角上向量，f为方程常 系数。x为方程组的解。,9,（二）迭代法,对于 ，A为非奇异矩阵，除了可用追赶法计算之外，也可建立迭代法求解。迭代法的一般格式为： ，其中B为迭代矩阵，对由迭代格式产生迭代序列 ，若 ，则 即 为 方程 的解。迭代法有以下几种方法：Jacobi迭代法，Gauss迭代法和SOR法。,10,Jacobi迭代法,1.迭代过程： 对于 ，设A可逆，且对于 令D为A的对角元素部分，则 ,继而对于方程组有 ，所以有解 ，令 为迭代矩阵， 。则Jacobi迭代法的迭代格式为： 。,11,2.算法设计：（1）将方程组系数矩阵用 表示，常数项用 表示。（2）用 存储初值；（3）用库函数diag,tril,triu求取对角矩阵D,下三角矩阵L,上三角矩阵U；（4）利用得到的D,L,U，求取迭代矩阵J，和f; （5）循环判断， 存取循环值，用二范式迭代循环，当 时，终止循环。,12,3.运行结果：（点击查看源程序） A=2 -1 0 0 0 ;-1 2 -1 0 0 ;0 -1 2 -1 0 ;0 0 -1 2 -1 ;0 0 0 -1 2 ; B=1 0 1 0 0; x0=0 0 0 0 0; x1=Jacobi(A,B,x0) n =78 x1 = 1.33331992455312 1.66664655349634 1.99997318243957 1.33331322016301 0.66665325788645注：其中A为方程组的系数矩阵，B为常数向量，x0为方程组的初值。 迭代次数n=78，x1为程序运行后方程组的解。,13,Gauss迭代法,1.迭代过程： Gauss迭代法的迭代过程与Jacobi迭代法的迭代过程相似，其迭代格式为 其中， , 。 算法设计同Jacobi迭代法，只是用库函数计算时的迭代矩阵B和f不同。,14,2.运行结果：（点击查看源程序） A=2 -1 0 0 0 ;-1 2 -1 0 0 ;0 -1 2 -1 0 ;0 0 -1 2 -1 ;0 0 0 -1 2 ; B=1 0 1 0 0; x0=0 0 0 0 0; n =41 x1 = 1.33332411479694 1.66665283886207 1.99998617219540 1.33332296247989 0.66666148123994注：其中迭代次数n=41，x1为迭代结果。,15,SOR法,1.迭代过程：逐次超松弛迭代法，即SOR方法是利用松弛思想来解决问题的。引入松弛变量w（其中w0），得到新的迭代格式： 其中 且 。当w=1时，SOR法就是Gauss迭代法SOR方法是Gauss迭代法的一种修正，主要思想是：设已知 及已计算是 的分量 首先，用Gauss迭代法定义辅助变量 ： （1式），再由 与 加权平均定义 ，即 （2式），将（1式）代入（2式）得到 的SOR迭代式。,16,2.算法设计： 在算法设计中，SOR方法的程序思想与Jacobi迭代法、Gauss迭代法的思想相近，都是基于库函数完成的，只是迭代矩阵的差异；另外，SOR方法每迭代一次主要运算量是计算一次矩阵与向量的乘法； 用 控制迭代终止程序。,17,3.运行结果为：（点击查看源程序） A=2 -1 0 0 0 ;-1 2 -1 0 0 ;0 -1 2 -1 0 ;0 0 -1 2 -1 ;0 0 0 -1 2 ; B=1 0 1 0 0; x0=0 0 0 0 0; w=1.4 w =1.40000000000000 n =15 x1 = 0.95237999615579 1.19047515937425 1.42857149554680 0.95238126882618 0.47619074994783,18,4.算法分析： 通过该迭代法程序代码，我们可以逐次检验1.0w1.9时，w取不同数值时SOR方法的迭代次数以及此时对应的解y。得到以下松弛因子与迭代次数表格：,19,五、结果分析,追赶法是求解三对角矩阵的常用方法，但从整体编程角度分析，其程序编写较迭代法复杂，但通用性较好。另外，对于追赶法和迭代法的有效数字位数来看，总体上，迭代法的有效数字位数较高。迭代法中，迭代次数由Jacobi迭代法、Gauss Seidel迭代法到SOR迭代法逐次降低，特别是在最佳松弛因子处最低。由于矩阵A为三个对角占优的对角矩阵，则根据P251页定理7可知，对于Jacobi和Gauss Seidel迭代法是收敛的。另外，由于Jacobi迭代法的迭代次数n=78，而Gauss Seidel迭代法的迭代次数n=41。所以Gauss Seidel迭代法的收敛速度较Jacobi迭代法的收敛速度快。,20,（接上页）对于SOR迭代法，由定理9可以确定：当0=2时，SOR迭代法不收敛。综上，我们可以得到：当使用w=1.4的SOR迭代法时，迭代次数最小，收敛速度最快。,21,六、程序源代码,1.追赶法程序代码：function x=ZhuiGan(a,b,c,f) %a为对角下向量；b为对角向量；c为对角上向量；f为方程常系数format longr=size(a);m=r(2);r=size(b);n=r(2);if size(a)=size(b)|m=n-1|size(b)=size(f) error(变量不匹配，检查变量匹配情况);end p=ones(1,m); Y=ones(1,n); x=Y; p(1)=a(1)/b(1); % c1/b1 Y(1)=f(1)/b(1); % f1/b1t=0;,22,（接上页）for i=2:m t=b(i)-a(i-1)*p(i-1); %p(i) 求出下三角矩阵U 的待定值 Ri p(i)=c(i)/t; Y(i)=(f(i)-a(i-1)*Y(i-1)/t; % Y(i)求出上三角矩阵L的待定值 Biend % 回代过程Y(n)=(f(n)-a(n-1)*Y(n-1)/(b(n)-a(n-1)*p(n-1);x(n)=Y(n); for i=n-1:-1:1 x(i)=Y(i)-p(i)*x(i+1); end n,23,追赶法运行结果：a=-1*ones(1,4);c=a; b=2*ones(1,5);f=1 0 1 0 0;x=ZhuiGan(a,b,c,f)n =5x = 1.33333333333333 1.66666666666667 2.00000000000000 1.33333333333333 0.66666666666667,24,2.Jacobi迭代法程序代码：function x1=Jacobi(A,B,x0)%A为系数矩阵；B为常数项；x0为初值format longD=diag(diag(A); %对角线矩阵 L=-tril(A,-1); %下三角矩阵U=-triu(A,1); %上三角矩阵J=(inv(D)*(L+U); %迭代矩阵f=(inv(D)*B;x1=J*x0+f;n=1;while(norm(x1-x0)=0.00001) x0=x1;x1=J*x0+f;n=n+1;endn,25,Jacobi迭代法运行结果：A=2 -1 0 0 0 ;-1 2 -1 0 0 ;0 -1 2 -1 0 ;0 0 -1 2 -1 ;0 0 0 -1 2 ;B=1 0 1 0 0;x0=0 0 0 0 0;x1=Jacobi(A,B,x0)n =78x1 = 1.33331992455312 1.66664655349634 1.99997318243957 1.33331322016301 0.66665325788645,26,3.Gauss迭代法程序代码：format longA=input(输入三对角矩阵A=);B=input(输入常数项B=);x0=input(输入方程初值x0=); D=diag(diag(A); L=-tril(A,-1); U=-triu(A,1); G=(inv(D-L)*U; f=(inv(D-L)*B;x1=G*x0+f;n=1;while(norm(x1-x0)=0.00001) x0=x1;x1=G*x0+f;n=n+1;endn ,x1,27,Gauss迭代法运行结果：输入三对角矩阵A=2 -1 0 0 0 ;-1 2 -1 0 0 ;0 -1 2 -1 0 ;0 0 -1 2 -1 ;0 0 0 -1 2 ;输入常数项B=1 0 1 0 0;输入方程初值x0=0 0 0 0 0;n =41x1 = 1.33332411479694 1.66665283886207 1.99998617219540 1.33332296247989 0.66666148123994,28,4.SOR迭代法程序代码：format longA=input(输入三对角矩阵A=);B=input(输入常数项B=);x0=input(输入方程初值x0=); D=diag(diag(A);L=-tril(A,-1);U=-triu(A,1);w=input(输入松弛因子w=)S=(inv(D-w*L)*(1-w)*D+w*U);f=(inv(D-w*L)*B;x1=S*x0+f;n=1;while(norm(x1-x0)=0.00001) x0=x1;x1=S*x0+f;n=n+1;endn ,x1,29,SOR迭代法运行结果：输入三对角矩阵A=2 -1 0 0 0 ;-1 2 -1 0 0 ;0 -1 2 -1 0 ;0 0 -1 2 -1 ;0 0 0 -1 2 ;输入常数项B=1 0 1 0 0;输入方程初值x0=0 0 0 0 0;输入松弛因子w=1.4w =1.40000000000000n =15x1 = 0.95237999615579 1.19047515937425 1.42857149554680 0.95238126882618 0.47619074994783,30,七、课程信息,题 目：三对角线性方程组的解法 组 别： 第二组 组 员： （见下页） 年 级： 2007 级 学 院：信息科学技术学院 系 别： 数学系 专 业： 信息与计算科学专业 指导教师： 罗 婷 完成日期： 2010年 6 月 4 日,31,32,
