数值分析5

三对角与块三对角方程组课程设计一、基于高斯消元法的三对角方程组求解三对角矩阵是一类重要的特殊矩阵，在数学计算和工程计算中有广泛应用。例如，二阶常微分方程边值问题数值求解，一维热传导方程数值求解，以及三次样条函数计算等都会涉及到三对角方程组求解。由于三对角矩阵的稀疏性质，用直接法求解三对角方程组的算法效率较高，很有实用价值。考虑n阶三对角矩阵和n维向量A=nn122111，nffff21求解方程组Ax=f的高斯消元法的程序如下functionf=triGauss(gama,alpha,bata,f)%SolvingTriDiag(gama,alpha,bata)systemsbyGaussmethodn=length(alpha);fork=1:n-1m=gama(k)/alpha(k);alpha(k+1)=alpha(k+1)-m*bata(k);f(k+1)=f(k+1)-m*f(k);endf(n)=f(n)/alpha(n);fork=n-1:-1:1f(k)=(f(k)-bata(k)*f(k+1))/alpha(k);end由程序知，对于n阶三对角方程组，高斯消元法只用到5n–4次乘法和除法。例1．求二阶常微分方程边值问题0)1(,0)0()1,0(,yyxxyy数值解，并与解析解：1sinhsinh)(xxxy作对比。解：对正整数n，取h=1/(n+1)，令xj=jh，（j=0，1，2，……，n+1），yj≈y(xj)。由2112hyyyyjjjj得差分方程jjjjjxyhyyy2112整理，得–yj-1+(2+h2)yj–yj+1=h2xj(j=1，…，n)根据边界条件，有y0=0，yn+1=0形成三对角方程组nnnnxxxxhyyyyhhhh121212122222112112112取n=9，得数值解和解析解的最大误差为：4.4146e-005，将数值解和解析解数据绘图如下图1.数值解与解析解比较程序如下n=input('inputn:=');h=1/(n+1);x=0:h:1;ux=x-sinh(x)./sinh(1);alpha=(2+h*h)*ones(n,1);bata=-ones(n-1,1);gama=bata;f(1:n)=h*h*x(2:n+1);uk=triGauss(gama,alpha,bata,f);Uk=[0,uk,0];error=max(abs(ux-Uk))plot(x,ux,x,Uk,'ro')取不同的n，运行程序分别计算实验数据如下n9192939error4.4146e-0051.1047e-0054.9106e-0062.7624e-006数据表明，随着结点增加，数值解的误差变小。二、迭代法求解三对角方程组数值分析介绍的三种迭代格式用于例1中三对角方程组：1．JACOBI迭代][21)(1)(122)1(kjkjjkjuuxhhu2．SEIDEL迭代][21)(1)1(122)1(kjkjjkjuuxhhu3．SOR迭代][2)1()(1)1(122)()1(kjkjjkjkjuuxhhuu由于JACOBI迭代矩阵的谱半径hhBJcos22)(2SEIDEL迭代矩阵的谱半径hhBSG222cos)22()(根据，的结果，取SOR迭代松驰因子为2)(112JBSOR迭代矩阵的谱半径22)(11)(111)(JJSORBBB实验数据如下(迭代精度要求10-8)nJACOBISEIDELSOR迭代数误差迭代数误差迭代次数误差5880.1191e-3470.1190e-3180.1190e-392300.4432e-41230.4422e-4290.4415e-4198290.1173e-44420.1137e-4560.1106e-43929130.5635e-515610.4176e-51070.2789e-5数据表明，三种迭代法计算出的数值解误差一致。但对三对角方程组而言，尽管迭代法程序比较简单，但迭代法的效率不如直接法。程序如下n=input('inputn:=');h=1/(n+1);x=0:h:1;y=x-sinh(x)/sinh(1);hh=h*h;rr=2+hh;u0=zeros(size(x));u=u0;k1=0;error=1.0;%----Jacobi迭代whileerror1.0e-8error=0;forj=2:n+1temp=(hh*x(j)+u(j-1)+u(j+1))/rr;error=max(error,abs(temp-u(j)));v(j)=temp;endu=[v,0];k1=k1+1;enderror1=max(abs(u-y));k1u=u0;error=1;k2=0;%----Seidel迭代whileerror1.0e-8error=0;forj=2:n+1temp=u(j);u(j)=(hh*x(j)+u(j-1)+u(j+1))/rr;error=max(abs(temp-u(j)),error);endk2=k2+1;enderror2=max(abs(u-y));k2u=u0;error=1;k3=0;%----SOR迭代ro=2*cos(h*pi)/rr;omiga=2/(1+sqrt(1-ro*ro));omiga1=1-omiga;whileerror1.0e-8error=0;forj=2:n+1temp=u(j);u(j)=omiga1*temp+omiga*(hh*x(j)+u(j-1)+u(j+1))/rr;error=max(abs(temp-u(j)),error);endk3=k3+1;enderror3=max(abs(u-y));k3[error1,error2,error3][ro,ro*ro,-omiga1]三、块三对角方程组的迭代方法拉普拉斯(Laplace)是法国著名数学家和天文学家。他用数学方法证明了行星的轨道大小只有周期性变化，这就是著名拉普拉斯的定理。1796年，他发表《宇宙体系论》，因研究太阳系稳定性的动力学问题被誉为法国的牛顿和天体力学之父。拉普拉斯在数学和物理学方面也有重要贡献，以他的名字命名的拉普拉斯变换和拉普拉斯方程，在科学技术的各个领域有着广泛的应用。对于二维拉普拉斯方程问题，方程的解是二元函数。对给定的边界条件，求数值解需要解块三对角方程组。例2．正方形区域上狄利克莱边值问题yyuxuxuyuyxuuyyxxsin),1(0)1,()0,(),0(1,0,0准确解：yshxshyxusin),(。取正整数n，记h=1/(n+1)，对区域做网格剖分：xi=ih(i=0，1，……，n+1)yj=jh(j=0，1，……，n+1)在点(xi，yj)处记uij=u(xi，yj)。离散拉普拉斯方程，可得常见的的五点处差分格式ui-1,j+ui,j-1–4uij+ui+1,j+ui,j+1=0(i，j=1，…，n)u0,j=0，un,j=sin(πj/n)(j=1，…，n)ui,0=0，ui,n=0(i=1，…，n)该方程组的系数矩阵是块三对角矩阵，与三对角矩阵相比，仍然具有稀疏性。但问题的规模比较大了。对于n=4，其矩阵是16阶方程，矩阵非零元分布如下图所示图2.块三对角矩阵非零元分布图1．五点差分格式的Seidel迭代方法][41)(1,)(,1)1(1,)1(,1)1(,kjikjikjikjikjiuuuuu，(i，j=0，1，…，N)2．SOR迭代算法：][4)1()(1,)(,1)1(1,)1(,1)()1(,kjikjikjikjikijkjiuuuuuu，(i，j=0，1，…，N)其中，最佳松驰因子为hsin12五点差分格式seidle迭代实验数据表(误差限要求：1e-008)结点数n2102202402602迭代次数18260620774291CPU时间0.29704.1158.531260.4840误差0.00236.4274e-0041.6814e-0047.3660e-005超松驰迭代法，取最佳松驰因子：hsin12五点差分格式SOR迭代实验数据表(误差限要求：1e-008)结点数n2102202402602迭代次数4074139201CPU时间0.110.65604.953015.5940误差0.00236.4306e-0041.6944e-0047.6600e-005对比两种迭代方法，知SOR迭代无论从迭代次数还是从计算机运行时间比较都优于seidle迭代法。但两种方法的数值解误差几乎是一样的。迭代计算绘图显示如下图3.二维拉普拉斯方程边值问题解图形由于边界条件中，右侧边界条件非零，而其余边界条件均为零，故有此图形。图中红色区域表示对应的函数值数值大，蓝色区域部分则表明对应的函数值数值小。附：MATLAB程序1．五点差分格式的seidel迭代法程序：n=input('n=');h=1/(n+1);x=0:h:1;y=x';u=zeros(n+2,n+2);u(:,n+2)=sin(pi*y);k=0;er=1;t=cputime;whileer1e-008er=0;k=k+1;fori=2:n+1forj=2:n+1s=(u(i-1,j)+u(i+1,j)+u(i,j-1)+u(i,j+1))/4;er=max(er,abs(s-u(i,j)));u(i,j)=s;endendendktimes=cputime-tpcolor(u)w=sin(pi*y)*sinh(pi*x)/sinh(pi);error=max(max(abs(w-u)))2．五点差分格式SOR迭代程序n=input('n=');h=1/(n+1);x=0:h:1;y=x';u=zeros(n+2,n+2);u(:,n+2)=sin(pi*y);w=2/(1+sin(pi*h));w1=1-w;k=0;er=1;t=cputime;whileer1e-008er=0;k=k+1;fori=2:n+1forj=2:n+1s=u(i,j);u(i,j)=w1*s+w*(u(i-1,j)+u(i+1,j)+u(i,j-1)+u(i,j+1))/4;er=max(er,abs(s-u(i,j)));endendendktimes=cputime-t%surf(u)%figure,pcolor(u)z=sin(pi*y)*sinh(pi*x)/sinh(pi);error=max(max(abs(z-u)))实验题1．求下列二阶常微分方程数值解xxxexuxux0),cos2(sin)()(0)(,0)0(uu该问题解析解为：u(x)=exsinx2．求正方形区域上的泊松方程数值解222yxu，0x，y1边界条件u(x，0)=0，221)1,(xxu，0x1yyusin),0(，221sin),1(yyeyu，0y1精确解：2)(21sin),(xyyeyxux