三对角线性方程组的解法

1数值分析课程设计题目：三对角线性方程组的解法指导老师：罗婷海南大学2目录•一、实验分析‥‥‥‥‥‥‥3•二、实验内容‥‥‥‥‥‥‥4•三、实验要求‥‥‥‥‥‥‥5•四、实验过程‥‥‥‥‥‥‥6•五、结果分析‥‥‥‥‥‥‥19•六、程序代码‥‥‥‥‥‥‥21•七、课程信息‥‥‥‥‥‥‥303一、实验分析•关于线性方程组的数值解法一般有两类：直接法和迭代法。而对于系数矩阵较特殊的一类线性方程组，解法也相对有所不同。在实际问题中，例如解常微分方程边值问题，解热传导问题等，都会要求解系数矩阵为对角占优的三对角方程组。•如何求解三对角方程组就成为关键。三对角线性方程组是系数矩阵为三对角矩阵的一类方程组，其数值解法有：追赶法，Jacobi迭代法，Gauss迭代法和SOR法。4考虑线性方程组：其中A为三对角矩阵，编制程序求解该方程组二、实验内容,,nnnAxbARbR5（1）考虑不同的数值解法，对方程组进行编程，编写其通用代码；（2）对不同的程序代码用实例进行验证，取矩阵：方程初值：。（3）比较不同方法对求解三对角线性方程组的优异性。三、实验要求210001121000,012101001210000120Ab0(0,0,0,0,0,0)Tx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）。(一)追赶法fLyyUx111/bc)/(1iiiiiabcfLy111/bfy)/()(11iiiiiiiabyafyyUxnnyx1iiiixyx8•运行结果：（点击查看源程序）a=-1*ones(1,4);c=a;b=2*ones(1,5);f=[10100];[x]=ZhuiGan(a,b,c,f)‘n=5x=1.333333333333331.666666666666672.000000000000001.333333333333330.66666666666667注：其中a为对角下向量，b为对角向量，c为对角上向量，f为方程常系数。x为方程组的解。9（二）迭代法•对于，A为非奇异矩阵，除了可用追赶法计算之外，也可建立迭代法求解。迭代法的一般格式为：，其中B为迭代矩阵，对由迭代格式产生迭代序列，若，则即为方程的解。•迭代法有以下几种方法：Jacobi迭代法，Gauss迭代法和SOR法。,,nnnAxbARbR(1)()kkxBxf()kx()*limkkxx**xBxf*xAxb10Jacobi迭代法1.迭代过程：对于，设A可逆，且对于令D为A的对角元素部分，则,继而对于方程组有，所以有解，令为迭代矩阵，。则Jacobi迭代法的迭代格式为：。Axb0,1,2,,iiain()AADD()DxDAxb11()xDDAxDb11()()JBDDADLU1JfDb(1)()kkjJxBxf112.算法设计：（1）将方程组系数矩阵用表示，常数项用表示。（2）用存储初值；（3）用库函数diag,tril,triu求取对角矩阵D,下三角矩阵L,上三角矩阵U；（4）利用得到的D,L,U，求取迭代矩阵J，和f;（5）循环判断，存取循环值，用二范式迭代循环，当时，终止循环。00001.010xxnnijaA)(1)(nijaB0x1x123.运行结果：（点击查看源程序）A=[2-1000;-12-100;0-12-10;00-12-1;000-12]';B=[10100]';x0=[00000]';[x1]=Jacobi(A,B,x0)n=78x1=1.333319924553121.666646553496341.999973182439571.333313220163010.66665325788645注：其中A为方程组的系数矩阵，B为常数向量，x0为方程组的初值。迭代次数n=78，x1为程序运行后方程组的解。13Gauss迭代法1.迭代过程：Gauss迭代法的迭代过程与Jacobi迭代法的迭代过程相似，其迭代格式为其中，,。算法设计同Jacobi迭代法，只是用库函数计算时的迭代矩阵B和f不同。GkGkfxBx)()1(ULDBG1)(bLDfG1)(142.运行结果：（点击查看源程序）A=[2-1000;-12-100;0-12-10;00-12-1;000-12]';B=[10100]';x0=[00000]';n=41x1=1.333324114796941.666652838862071.999986172195401.333322962479890.66666148123994注：其中迭代次数n=41，x1为迭代结果。15SOR法1.迭代过程：•逐次超松弛迭代法，即SOR方法是利用松弛思想来解决问题的。引入松弛变量w（其中w0），得到新的迭代格式：其中且。当w=1时，SOR法就是Gauss迭代法•SOR方法是Gauss迭代法的一种修正，主要思想是：设已知及已计算是的分量•首先，用Gauss迭代法定义辅助变量：（1式），•再由与加权平均定义，即（2式），•将（1式）代入（2式）得到的SOR迭代式。wkwkfxBx)()1())1(()(1wUDwwLDBwwbfw)(kx)1(kx)1,,2,1()1(ijxkjiiijnijkjijkjijikiaxaxabx/)(111)()1()1()1(kix)(kix)1(kix)1(kix)()1()()1()1()1()()1(kikikikikikixxwxxwxwxAxb162.算法设计：在算法设计中，SOR方法的程序思想与Jacobi迭代法、Gauss迭代法的思想相近，都是基于库函数完成的，只是迭代矩阵的差异；另外，SOR方法每迭代一次主要运算量是计算一次矩阵与向量的乘法；用控制迭代终止程序。2)(2)(kkAxbr173.运行结果为：（点击查看源程序）A=[2-1000;-12-100;0-12-10;00-12-1;000-12]';B=[10100]';x0=[00000]';w=1.4w=1.40000000000000n=15x1=0.952379996155791.190475159374251.428571495546800.952381268826180.47619074994783184.算法分析：通过该迭代法程序代码，我们可以逐次检验1.0≤w≤1.9时，w取不同数值时SOR方法的迭代次数以及此时对应的解y。得到以下松弛因子与迭代次数表格：松弛因子w迭代次数n松弛因子w迭代次数n1.0411.5191.1331.6241.2261.7351.3181.8531.4151.911219五、结果分析•追赶法是求解三对角矩阵的常用方法，但从整体编程角度分析，其程序编写较迭代法复杂，但通用性较好。另外，对于追赶法和迭代法的有效数字位数来看，总体上，迭代法的有效数字位数较高。迭代法中，迭代次数由Jacobi迭代法、Gauss－Seidel迭代法到SOR迭代法逐次降低，特别是在最佳松弛因子处最低。•由于矩阵A为三个对角占优的对角矩阵，则根据P251页定理7可知，对于Jacobi和Gauss－Seidel迭代法是收敛的。另外，由于Jacobi迭代法的迭代次数n=78，而Gauss－Seidel迭代法的迭代次数n=41。所以Gauss－Seidel迭代法的收敛速度较Jacobi迭代法的收敛速度快。20（接上页）•对于SOR迭代法，由定理9可以确定：当0w2时，SOR迭代法是收敛的，其迭代次数则会因选择的松弛因子的不同，而有所不同。在以上表格中可以看出：当w=1.4时，迭代次数最小为n=15，此时的w为最佳松弛因子。同时，我们也可以验证：当w=2时，SOR迭代法不收敛。•综上，我们可以得到：当使用w=1.4的SOR迭代法时，迭代次数最小，收敛速度最快。21六、程序源代码1.追赶法程序代码：•function[x]=ZhuiGan(a,b,c,f)•%a为对角下向量；b为对角向量；c为对角上向量；f为方程常系数•formatlong•r=size(a);•m=r(2);•r=size(b);•n=r(2);•ifsize(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/b1•t=0;22（接上页）•fori=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的待定值Bi•end••%回代过程•Y(n)=(f(n)-a(n-1)*Y(n-1))/(b(n)-a(n-1)*p(n-1));•x(n)=Y(n);••fori=n-1:-1:1•x(i)=Y(i)-p(i)*x(i+1);•end•n23•追赶法运行结果：•a=-1*ones(1,4);•c=a;•b=2*ones(1,5);•f=[10100];•[x]=ZhuiGan(a,b,c,f)’•n=5•x=•1.33333333333333•1.66666666666667•2.00000000000000•1.33333333333333•0.66666666666667242.Jacobi迭代法程序代码：•function[x1]=Jacobi(A,B,x0)•%A为系数矩阵；B为常数项；x0为初值•formatlong•D=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;•end•n25•Jacobi迭代法运行结果：•A=[2-1000;-12-100;0-12-10;00-12-1;000-12]';•B=[10100]';•x0=[00000]';•[x1]=Jacobi(A,B,x0)•n=78•x1=•1.33331992455312•1.66664655349634•1.99997318243957•1.33331322016301•0.66665325788645263.Gauss迭代法程序代码：•formatlong•A=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;•end•n,x127•Gauss迭代法运行结果：•输入三对角矩阵A=[2-1000;-12-100;0-12-10;00-12-1;000-12]';•输入常数项B=[10100]';•输入方程