三对角与块三对角方程组课程设计.docx

三对角与块三对角方程组课程设计一、基于高斯消元法的三对角方程组求解三对角矩阵是一类重要的特殊矩阵，在数学计算和丁•稈计算中有广泛应用例如，二阶常 微分方程边值问题数值求解,一维热传导方程数值求解,以及三次样条函数计算等都会涉及到 三对角方程组求解由于三对角矩阵的稀疏性质，用肓接法求解三对角方程组的算法效率较高, 很有实用价值考虑n阶三对角矩阵和n维向量0]/1 5 角L Yn-X J求解方程纟R Ax = f的高斯消元法的程序如下function f=triGauss(gama,alpha,bata,f)%Solving TriDiag(ganici,alpha,bata)systenis by Gauss method n=le ngth(alpha);for k=l:n-lm=gama(k)/alpha(k);alpha(k+1 )=alpha(k+ l)-m*bata(k);f(k+1 )=f(k+ l)-m*f(k);endf(n )=f(n)/alpha(n);for k=n-l :-l :1f(k)=(f(k)-bata(k)*f(k+ l))/alpha(k);由程序知，对于n阶三对角方稈组，高斯消元法貝用到5n -4次乘法和除法。

例求二阶常微分方程边值问题+ y =工，x e (0,1) b(o)= o, xn = o sinh r数值解，并与解析解：作对比解：对正整数⑦取 h= \!(n + 1),令Xj=jh, ( J=o, 1, 2, ,刃+1),乃—Xx/)° 由得羞分方程整理，得根据边界条件，有 形成三对角方程组「2 +沪-1Ji・H・-12 + X-1y2兀2•■■• •• •• •■■■=h2■■■-1 2 + /-1儿一1心一1■-12 + h2儿MM取n=9，得数值解和解析解的最人误差为：4.41460-005,将数值解和解析解数据绘图如下稈序如下n=input(finput n:=f); h=l/(n+l); x=O:h:l;ux=x-sinh(x)./sinh( 1); alpha=(2+h*h)*ones(n,l); bata=-ones(n-l ,l);gama=bata; f(l :n)=h*h*x(2:n+1);uk=triGauss(gama,alpha,bata,f)i Uk=lO,uk,OJ;error=max(abs(ux-Uk)) plot(x,ux,x,Uk/ror)取不同的n,运行程序分别计算实验数据如下n9192939error4.4146 — 0051. 1047 — 0054.9106—0062.7624e-006数据表明，随看结点增加，数值解的误差变小。

二、迭代法求解三对角方程纟R数值分析介绍的三种迭代格式用于例1中三对角方稈组:1. JACOBI 迭代7 2 + h2 J J J2. SEIDEL 迭代_—ih2x2 + 03. SOR迭代=c-讪罗+ 宀厅厂+必丁 +晴］2 + /T由于JACOBI迭代矩阵的谱半径—— cos加2 + h2SEIDEL迭代矩阵的谱半径2p(BG_s) = (—^7Vf cos2 7ih2 + /7根据，的结果，取SOR迭代松驰因子为2I + - P(Bj『1-Jl-Q©)21+J1 —0(7)2SOR迭代炬阵的谱半径nJACOBISEIDELSOR迭代数误差迭代数误差迭代次数误差5880.1191 e-3470」190 e-3180」190 e-392300.4432 e-41230.4422 e-4290.4415 e-4198290」173 e-44420」137 e-4560.1106 e-43929130.5635 e-515610.4176 e-51070.2789 e-5P(Bsor) =⑴一' =实验数据如下(迭代精度要求108)数据表明，三种迭代法计算岀的数值解误弟一致但对三对角方程纟R而言，尽管迭代法程序比 较简单，但迭代法的效率不如直接法。

程序如下n=input(4nput n:=*); h=l/(n+l);x=O:h:l; y=x-sinh(x)/sinh(l); hh=h*h;rr=2+hh; uO=zeros(size(x)); u=uO;k 1 =0;error= 1.0;%-…Jacobi迭代 while error> 1.0e-8error=0;for j=2:n+lteinp=(hh*x(j)+u(j-l)+u(j+1 ))/rr; error=max(erronabs(temp-u(j))); v(j)=temp;endu=[v,0];kl=kl+l;enderror 1 =max(abs(u ・y));kl u=u0;error= 1 ;k2=0;%■…Seidel迭代 while error> 1.0e-8error=0;for j=2:n+lteiTip=u(j);u(j)=(hh*x(j)+u(j-l)+u(j+1 ))/rr; error=max(abs(teiTip-u(j)),error); endk2=k2+1;enderror2=max(abs(u-y));k2u=uO;error= 1 ；k3=0;%一一SOR迭代ro=2*cos(h*pi)/rr;omiga=2/(l +sqrt(l -ro*ro));omigal = l-omiga;while error> 1 .Oe-8error=0;for j=2:n+ltemp=u(j);u(j)=omiga 1 *temp+omiga*(hh*x(j)+u(j-1 )+u(j+1 ))/rr; error=max(abs(temp-u(j))5error);end k3=k3+1;enderror3=max(abs(u-y)); k3[error 1 ,error2,error3] [ro,ro*roromigal ]三、块三对角方程组的迭代方法拉普拉斯(Laplace)是法国著名数学家和天文学家。

他用数学方法证明了行星的轨道大小只 有周期性变化，这就是著名拉普拉斯的定理1796年，他发表《宇宙体系论》，因研究太阳系 稳定性的动力学问题被誉为法国的牛顿和天体力学之父拉普拉斯在数学和物理学方面也有重 要贡献，以他的名字命名的拉普拉斯变换和拉普拉斯方程，在科学技术的各个领域有着广泛的 应用对于二维拉普拉斯方程问题，方程的解是二元函数对给定的边界条件，求数值解需要 解块三对角方稈组例2.正方形区域上狄利克莱边值问题% + Uyy =0, 0 < X, J/ < 1< w(0,y) = w(x,0) = w(x,l) — 0"(1,尹)=sin 矽准确解：w(x, j) = ^^sin 初 osh7T取正整数⑦记治1/(卄1),对区域做网格剖分：Xj = i h (7 = 0, 1, , 〃+1 )yj =j h (j=0, 1, , 77+I)在点(七，yj)处记uy = 刃)离散拉普拉斯方稈，可得常见的的五点处差分格式⑷-1 j+ ⑷j-i - 4 Ujj + ui+]J + Uij+]= 0(j，丿=1，…，n )叫=0，%j = sin(7z/7/7) (j = l,…，乃)⑷。

0, uiyn = 0 (z = l,…，n )该方程组的系数矩阵是块三对角矩阵，与三对角矩阵相比，仍然具有稀疏性但问题的规模比 较大了对于n=4,其矩阵是16阶方程，矩阵非零元分布如下图所示510150 5 10 15图2.块三对角矩阵非零元分布图1. 五点差分格式的Seidel迭代方法p*+l）+um,j +Mi*+J » （j，_/=0, 1,…，N）2. SOR迭代算法：酱）=（1-効碣）+彳同第+邨）+凋j +喘J，（i，j = 0,其中，最佳松驰因了为2a)= 1 + sin 加说点差分为式seidle迭佯实验数爭表聲差限要求：1巳；008）“上 也 一2 1 n2 on2 I /in2结点数n21022()'402602迭代次数18260620774291CPU时间0.29704.1158.531260.4840误差0.00236.4274e-0041.6814e-0047.3660e-005超松驰迭代法，取最佳松驰因子：3 = 1 + sm 7di五点旁分?各式S?R洋代实验数甲表（•误差限要求：严-?08） I 上好 一2 1 A2 OA2 ACX1结点数n210’20~40J60J迭代次数4074139201CPU时间0.110.65604.953015.5940误差0.00236.4306e-0041.6944e-0047.6600e-005对比两种迭代方法，知SOR迭代无论从迭代次数还是从计算机运行时间比较部优于seidle迭代 法。

但两种方法的数值解误差儿乎是一样的迭代计算绘图显示如下20 40 60 80 100图3 •二维拉普拉斯方程边值问题解图形市于边界条件中，右侧边界条件非零，而其余边界条件均为零，故冇此图形图中红色区域表示对应的函数值数值大，蓝色区域部分则表明对应的函数值数值小附：MATLAB程序1 •九点羌分格式的seidel迭代法程序：n=inputfn=r);h= 1 /(n+ l)；x=O:h:l ;y=xr; u=zeros(n+2,n+2); u(:,n+2)=sin(pi*y); k=O;er=l;t=cputime;while er>le-008er=O;k=k+l;for i=2:n+lforj=2:n+ls=(u(i-l ,j)+u(i+1 ,j)+u(i,j-1 )+u(i,j+1))/4; er=max(er,abs(s-u(ij)));u(i,j)=s;endendendktimes=cputime-tpcolor(u) w=sin(pi*y)*si nh(pi*x)/si nh(pi);error=max(max(abs(w-u)))2.五点斧分格式SOR迭代程序n=inputCn=r);h=l/(n+l);x=0:h:l ;y=xf; u=zeros(n+2, n+2); u(:,n+2)=sin(pi*y); w=2/(l+sin(pi*h));wl = l-w;k=0;er=l;t=cputime;while er> 1 e-008er=0;k=k+l;for i=2:n+lfor j=2:n+ls=u(i,j); u(i,j)=w 1 *s+w*(u(i-l ,j)+u(i+1。