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

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

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

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

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

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

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

%Solving TriDiag(gama,alpha,bata)systems by Gauss method n=length(alpha); for k=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); for k=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 = h2 xj ( 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('input n:='); 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，运行程序分别计算实验数据如下 9 19 29 39 n 4.4146e-005 1.1047e-005 4.9106e-006 2.7624e-006 error 数据表明，随着结点增加，数值解的误差变小。

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

数值分析介绍的三种迭代格式用于例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?h 22?h?(BG?S)?(根据，的结果，取SOR迭代松驰因子为

??SOR迭代矩阵的谱半径

21?1??(BJ)2

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

实验数据如下(迭代精度要求10-8) n JACOBI SEIDEL SOR 迭代数 误差 迭代数 误差 迭代次数 误差 5 88 0.1191 e-3 47 0.1190 e-3 18 0.1190 e-3 9 230 0.4432 e-4 123 0.4422 e-4 29 0.4415 e-4 19 829 0.1173 e-4 442 0.1137 e-4 56 0.1106 e-4 39 2913 0.5635 e-5 1561 0.4176 e-5 107 0.2789 e-5 数据表明，三种迭代法计算出的数值解误差一致。但对三对角方程组而言，尽管迭代法程序比较简单，但迭代法的效率不如直接法。 程序如下

n=input('input n:='); 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迭代 while error>1.0e-8 error=0; for j=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迭代 while error>1.0e-8 error=0; for j=2:n+1 temp=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; while error>1.0e-8 error=0; for j=2:n+1 temp=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 = i h ( i = 0，1，……，n+1 ) yj = j h ( j = 0，1，……，n+1)

在点(xi，yj ) 处记uij = u(xi，yj)。离散拉普拉斯方程，可得常见的的五点处差分格式

ui-1, j + ui, j-1 – 4 uij + 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?h602 4291 260.4840 7.3660e-005 五点差分格式seidle迭代实验数据表(误差限要求：1e-008) 2102 202 402 结点数 n 182 606 2077 迭代次数 0.2970 4.11 58.531 CPU时间 0.0023 6.4274e-004 1.6814e-004 误差 超松驰迭代法，取最佳松驰因子：??2

1?sin?h五点差分格式SOR迭代实验数据表(误差限要求：1e-008) 2102 202 402 602 结点数 n 40 74 139 201 迭代次数 0.11 0.6560 4.9530 15.5940 CPU时间 0.0023 6.4306e-004 1.6944e-004 7.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;

while er>1e-008 er=0;k=k+1; for i=2:n+1 for j=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; end end end k

times=cputime-t pcolor(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;

while er>1e-008 er=0;k=k+1; for i=2:n+1 for j=2:n+1 s=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))); end end end k

times=cputime-t %surf(u)

%figure,pcolor(u)

z=sin(pi*y)*sinh(pi*x)/sinh(pi);

error=max(max(abs(z-u)))

实验题

1．求下列二阶常微分方程数值解

?u??(x)?u(x)?ex(sinx?2cosx),u(0)?0,u(?)?0 x

该问题解析解为：u(x)= e sin x

2．求正方形区域上的泊松方程数值解

0?x??

?2u?x2?y2，0< x，y < 1

边界条件

u(x，0)=0，u(x,1)?12x，0< x < 1 212y，0< y < 1 2u(0,y)?sin?y，u(1,y)?e?sin?y?精确解：u(x,y)?e

?xsin?y?1(xy)2 2

本文来源：https://www.bwwdw.com/article/e93r.html