数值分析课程设计数值分析课程设计 题目：题目：三对角线性方程组的解法三对角线性方程组的解法 指导老师：指导老师： 罗罗 婷婷 海海 南南 大大 学学.目目 录录一、实验分析 3二、实验内容 4三、实验要求 5四、实验过程 6五、结果分析 19六、程序代码 21七、课程信息 30.一、实验分析关于线性方程组的数值解法一般有两类：直接法和迭代法。而对于系数矩阵较特殊的一类线性方程组，解法也相对有所不同。在实际问题中，例如解常微分方程边值问题，解热传导问题等，都会要求解系数矩阵为对角占优的三对角方程组。如何求解三对角方程组就成为关键。三对角线性方程组是系数矩阵为三对角矩阵的一类方程组，其数值解法有：追赶法，Jacobi迭代法，Gauss迭代法和SOR法。.考虑线性方程组：其中A为三对角矩阵，编制程序求解该方程组二、实验内容.（1）考虑不同的数值解法，对方程组进行编程，编写其通用代码；）考虑不同的数值解法，对方程组进行编程，编写其通用代码；（2）对不同的程序代码用实例进行验证，取矩阵：）对不同的程序代码用实例进行验证，取矩阵： 方程初值方程初值 ： 。 （3）比较不同方法对求解三对角线性方程组的优异性。）比较不同方法对求解三对角线性方程组的优异性。三、实验要求.四、实验过程解法解法追赶法追赶法Jacobi迭代法迭代法 Gauss迭代法迭代法 SOR迭代法迭代法 .编程思想编程思想： 矩阵的三对角线性方程组是由矩阵的直接三角分解法来推到的来矩阵的三对角线性方程组是由矩阵的直接三角分解法来推到的来的，求解等价于解两个三角形方程组：的，求解等价于解两个三角形方程组： （1） ，求，求y；（；（2） ，求，求x 从而得到解三对角线方程组的追赶法公式。从而得到解三对角线方程组的追赶法公式。 1. 计算的递推公式：计算的递推公式： （i=2,3,n-1）；）； 2. 解解 （i=2,3,n）；）； 3. 解解 （i=n-1，n-2，2,1）。）。(一)追赶法.运行结果：（运行结果：（点击查看源程序点击查看源程序） 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为方程组的解。.（二）迭代法迭代法 对于 ，A为非奇异矩阵，除了可用追赶法计算之外，也可建立迭代法求解。迭代法的一般格式为： ，其中B为迭代矩阵，对由迭代格式产生迭代序列 ，若 ，则 即 为 方程 的解。迭代法有以下几种方法：Jacobi迭代法，Gauss迭代法和SOR法。.Jacobi迭代法迭代法 1.迭代过程：迭代过程： 对于 ，设A可逆，且对于 令D为A的对角元素部分，则 ,继而对于方程组有 ，所以有解 ，令 为迭代矩阵， 。则Jacobi迭代法的迭代格式为： 。.2.算法设计：算法设计：（1）将方程组系数矩阵用）将方程组系数矩阵用 表示，常数项表示，常数项用用 表示。表示。（2）用）用 存储初值；存储初值；（3）用库函数）用库函数diag,tril,triu求取对角矩阵求取对角矩阵D,下三角下三角矩阵矩阵L,上三角矩阵上三角矩阵U；（4）利用得到的）利用得到的D,L,U，求取迭代矩阵，求取迭代矩阵J，和，和f; （5）循环判断，）循环判断， 存取循环值，用二范式迭代循存取循环值，用二范式迭代循环，当环，当 时，终止循环。时，终止循环。.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为程序运行后方程组的解。.Gauss迭代法迭代法 1.迭代过程：迭代过程： Gauss迭代法的迭代过程与Jacobi迭代法的迭代过程相似，其迭代格式为 其中， , 。 算法设计同Jacobi迭代法，只是用库函数计算时的迭代矩阵B和f不同。.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为迭代结果。.SOR法法 1.迭代过程：迭代过程：逐次超松弛迭代法，即逐次超松弛迭代法，即SOR方法是利用松弛思想来解决问题的。方法是利用松弛思想来解决问题的。引入松弛变量引入松弛变量w（其中（其中w0），得到新的迭代格式：），得到新的迭代格式： 其中其中 且且 。当。当w=1时，时，SOR法就是法就是Gauss迭代法迭代法SOR方法是方法是Gauss迭代法的一种修正，主要思想是：设已知迭代法的一种修正，主要思想是：设已知 及已计算是及已计算是 的分量的分量 首先，用首先，用Gauss迭代法定义辅助变量迭代法定义辅助变量 ： （1式），式），再由再由 与与 加权平均定义加权平均定义 ，即，即 （2式），式），将（将（1式）代入（式）代入（2式）得到式）得到 的的SOR迭代式。迭代式。 .2.算法设计：算法设计： 在算法设计中，SOR方法的程序思想与Jacobi迭代法、Gauss迭代法的思想相近，都是基于库函数完成的，只是迭代矩阵的差异；另外，SOR方法每迭代一次主要运算量是计算一次矩阵与向量的乘法； 用 控制迭代终止程序。.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.4.算法分析：算法分析： 通过该迭代法程序代码，我们可以逐次检验1.0w1.9时，w取不同数值时SOR方法的迭代次数以及此时对应的解y。得到以下松弛因子与迭代次数表格： 松弛因子w 迭代次数n松弛因子w 迭代次数n1.0411.5191.1331.6241.2261.7351.3181.8531.4151.9112.五、结果分析追赶法是求解三对角矩阵的常用方法，但从整体编程角度分析，其程序编写较迭代法复杂，但通用性较好。另外，对于追赶法和迭代法的有效数字位数来看，总体上，迭代法的有效数字位数较高。迭代法中，迭代次数由Jacobi迭代法、Gauss Seidel迭代法到SOR迭代法逐次降低，特别是在最佳松弛因子处最低。由于矩阵A为三个对角占优的对角矩阵，则根据P251页定理7可知，对于Jacobi和Gauss Seidel迭代法是收敛的。另外，由于Jacobi迭代法的迭代次数n=78，而Gauss Seidel迭代法的迭代次数n=41。所以Gauss Seidel迭代法的收敛速度较Jacobi迭代法的收敛速度快。.（接上页）对于SOR迭代法，由定理9可以确定：当0w=2时，SOR迭代法不收敛。综上，我们可以得到：当使用w=1.4的SOR迭代法时，迭代次数最小，收敛速度最快。.六、程序源代码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;.（接上页）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.追赶法运行结果：追赶法运行结果：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.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.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.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 .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.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 .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.七、课程信息题 目：三对角线性方程组的解法 组 别： 第二组 组 员： （见下页） 年 级： 2007 级 学 院：信息科学技术学院 系 别： 数学系 专 业： 信息与计算科学专业 指导教师： 罗 婷 完成日期： 2010年 6 月 4 日 .学号姓名学号姓名20070744012黄顺勇20070744018李婉20070744013黄永刚20070744020李玉娇20070744014李冰姬20070744021林志德20070744015李河20070744023吕怀申20070744016李亮亮20070744024潘发东20070744017李少华.