基于高斯消元法的三对角矩阵LU分解

本文格式为Word版，下载可任意编辑——基于高斯消元法的三对角矩阵LU分解三对角与块三对角方程组课程设计

一、基于高斯消元法的三对角方程组求解

三对角矩阵是一类重要的特别矩阵，在数学计算和工程计算中有广泛应用。例如，二阶常微分方程边值问题数值求解，一维热传导方程数值求解，以及三次样条函数计算等都会涉及到三对角方程组求解。由于三对角矩阵的稀疏性质，用直接法求解三对角方程组的算法效率较高，很有实用价值。

考虑n阶三对角矩阵和n维向量

?f1???1?1??f??????122?，f??2?A=???????????????n?1n???fn?求解方程组Ax=f的高斯消元法的程序如下

functionf=triGauss(gama,alpha,bata,f)

%SolvingTriDiag(gama,alpha,bata)systemsbyGaussmethodn=length(alpha);fork=1:n-1

m=gama(k)/alpha(k);

alpha(k+1)=alpha(k+1)-m*bata(k);f(k+1)=f(k+1)-m*f(k);end

f(n)=f(n)/alpha(n);fork=n-1:-1:1

f(k)=(f(k)-bata(k)*f(k+1))/alpha(k);end

由程序知，对于n阶三对角方程组，高斯消元法只用到5n–4次乘法和除法。例1．求二阶常微分方程边值问题

??y???y?x,x?(0,1)?y(0)?0,y(1)?0?sinhx数值解，并与解析解：y(x)?x?作对比。

sinh1解：对正整数n，取h=1/(n+1)，令xj=jh，（j=0，1，2，……，n+1），yj≈y(xj)。由

yj?1?2yj?yj?1y?j??2h得差分方程

?整理，得

yj?1?2yj?yj?1h2?yj?xj

–yj-1+(2+h2)yj–yj+1=h2xj(j=1，…，n)

根据边界条件，有

y0=0，yn+1=0

形成三对角方程组

?2?h2???1??????12?h2??1???12?h2?1?x1???y1?????y?x?2???2?2??h???????????x?1??yn?1??n?1????x?2?h2???yn??n?取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，运行程序分别计算试验数据如下9192939n4.4146e-0051.1047e-0054.9106e-0062.7624e-006error数据说明，随着结点增加，数值解的误差变小。

二、迭代法求解三对角方程组

数值分析介绍的三种迭代格式用于例1中三对角方程组：1．JACOBI迭代

u(jk?1)?u(jk?1)?12(k)(k)[hx?u?ujj?1j?1]22?h1?1)(k)[h2xj?u(jk?1?uj?1]22?h2．SEIDEL迭代

3．SOR迭代

u(jk?1)?(1??)u(jk)??2?h2(k?1)(k)[hx?u?ujj?1j?1]2由于JACOBI迭代矩阵的谱半径

?(BJ)?SEIDEL迭代矩阵的谱半径

2cos?h

2?h2222)cos?h22?h?(BG?S)?(根据，的结果，取SOR迭代松驰因子为

??SOR迭代矩阵的谱半径

21?1??(BJ)2

?(BSOR)???1?1?1??(BJ)21?1??(BJ)2

试验数据如下(迭代精度要求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迭代whileerror>1.0e-8error=0;forj=2:n+1

temp=(hh*x(j)+u(j-1)+u(j+1))/rr;error=max(error,abs(temp-u(j)));v(j)=temp;end

u=[v,0];k1=k1+1;end

error1=max(abs(u-y));k1

u=u0;error=1;k2=0;%Seidel迭代whileerror>1.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);end

k2=k2+1;end

error2=max(abs(u-y));k2

u=u0;error=1;k3=0;%SOR迭代ro=2*cos(h*pi)/rr;

omiga=2/(1+sqrt(1-ro*ro));omiga1=1-omiga;whileerror>1.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);end

k3=k3+1;end

error3=max(abs(u-y));k3

[error1,error2,error3][ro,ro*ro,-omiga1]

三、块三对角方程组的迭代方法

拉普拉斯(Laplace)是法国著名数学家和天文学家。他用数学方法证明白行星的轨道大小只有周期性变化，这就是著名拉普拉斯的定理。1796年，他发表《宇宙体系论》，因研究太阳系稳定性的动力学问题被誉为法国的牛顿和天体力学之父。拉普拉斯在数学和物理学方面也有重要贡献，以他的名字命名的拉普拉斯变换和拉普拉斯方程，在科学技术的各个领域有着广泛的应用。对于二维拉普拉斯方程问题，方程的解是二元函数。对给定的边界条件，求数值解需要解块三对角方程组。

例2．正方形区域上狄利克莱边值问题

?uxx?uyy?0,0?x,y?1??u(0,y)?u(x,0)?u(x,1)?0?u(1,y)?sin?y?sh?xsin?y。sh?取正整数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)

确凿解：u(x,y)?该方程组的系数矩阵是块三对角矩阵，与三对角矩阵相比，依旧具有稀疏性。但问题的规模比较大了。对于n=4，其矩阵是16阶方程，矩阵非零元分布如下图所示

图2.块三对角矩阵非零元分布图

1．五点差分格式的Seidel迭代方法

1k?1)(k?1)(k)(k)ui(,kj?1)?[ui(?1,j?ui,j?1?ui?1,j?ui,j?1]，(i，j=0，1，…，N)

42．SOR迭代算法：

(k)ui(,kj?1)?(1??)uij??4k?1)(k?1)(k)(k)[ui(?1,j?ui,j?1?ui?1,j?ui,j?1]，(i，j=0，1，…，N)

其中，最正确松驰因子为

??2

1?sin?h6024291260.48407.3660e-005五点差分格式seidle迭代试验数据表(误差限要求：1e-008)2102202402结点数n1826062077迭代次数0.29704.1158.531CPU时间0.00236.4274e-0041.6814e-004误差超松驰迭代法，取最正确松驰因子：??2

1?sin?h五点差分格式SOR迭代试验数据表(误差限要求：1e-008)2102202402602结点数n4074139201迭代次数0.110.65604.953015.5940CPU时间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;

whileer>1e-008er=0;k=k+1;fori=2:n+1forj=2:n+1

s=(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;endendendk

times=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;

whileer>1e-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)));