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

三对角与块三对角方程组课程设计 一、基于高斯消元法的三对角方程组求解 三对角矩阵是一类重要的特殊矩阵，在数学计算和工程计算中有广泛应用。例如，二阶常微分方程边值问题数值求解，一维热传导方程数值求解，以及三次样条函数计算等都会涉及到三对角方程组求解。由于三对角矩阵的稀疏性质，用直接法求解三对角方程组的算法效率较高，很有实用价值。 考虑n阶三对角矩阵和n维向量 A =， 求解方程组 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．求二阶常微分方程边值问题 数值解，并与解析解：作对比。 解：对正整数n，取h= 1/(n + 1)，令xj= jh，（ j = 0，1，2，……，n+1），yj ≈ y( xj)。由 得差分方程 整理，得 – yj-1 + (2 + h2 ) yj – yj+1 = h2 xj ( j = 1，…，n) 根据边界条件，有 y0 = 0，yn+1 = 0 形成三对角方程组 取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，运行程序分别计算实验数据如下 n 9 19 29 39 error 4.4146e-005 1.1047e-005 4.9106e-006 2.7624e-006 数据表明，随着结点增加，数值解的误差变小。 二、迭代法求解三对角方程组 数值分析介绍的三种迭代格式用于例1中三对角方程组： JACOBI迭代 SEIDEL迭代 SOR迭代 由于JACOBI迭代矩阵的谱半径 SEIDEL迭代矩阵的谱半径 根据，的结果，取SOR迭代松驰因子为 SOR迭代矩阵的谱半径 实验数据如下(迭代精度要求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 error1.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 error1.0e-8 error=0; for j=2:n+1 temp=u(