三角线性方程组的数值计算.docx

三角线性方程组的数值计算1. 实验描述 有回代的高斯消元法：如果A是NN非奇异矩阵，则存在线性方程组UX=Y与线性方程组AX=B等价，这里U是上三角矩阵，并且ukk0。当构造出U和Y后，可以用回代法求解UX=Y，并且得到方程组的解X。三角分解法：设线性方程组AX=B的系数矩阵A存在三角分解（如果非奇异矩阵A可以表示为下三角矩阵L和上三角矩阵U的乘积：A=LU 则A存在一个三角分解），则线性方程组可以表示为 LUX=B 而方程组的解可以通过定义Y=UX并求解下面的俩个方程组得到。首先对方程组LY=B求解Y，然后对方程组UX=Y求解X。 雅可比迭代：设矩阵A具有严格对角优势，则AX=B有唯一解X=P利用迭代式（xj（k+1)=bj-aj1xk-ajj-1xj-1k-ajj+1xjj+1k-ajNxN(k)ajj (其中j=1，2，N)）可产生一个向量序列Pk，而且对于任意初始向量P0，向量序列都将收敛到P0。 高斯-赛德尔迭代：xj（k+1)=bj-aj1xk+1-ajj-1xj-1k+1-ajj+1xjj+1k-ajNxN(k)ajj2. 实验内容设有如下三角线性方程组，而且系数矩阵具有严格对角优势：d1x1+c1x1 =b1a1+d2x2+c2x2 =b2 a2x2+d3x3+c3x4 =b3 . . . . . . . . . . . . aN-2xN-2+dN-1xN-1+cN-1xN=bN-1 aN-1xN-1+dNxN=bN 假设d1=d2=dn=12，b1=b2=bn=5;c1=c2=.=cn=a1=a2=an=-2; 1.用高斯消元法计算；2.采用LU分解来计算；3.采用雅可比迭代法计算；4.采用高斯-赛德尔迭代法计算。3. 实验结果及分析， 高斯消元法：分析：先做增广矩阵AugAB 通过高斯消元法：m=Aug(k,k-1)/Aug(k-1,k-1); Aug(k,k:N+1)=Aug(k,k:N+1)-m*Aug(k-1,k:N+1);从而得到上三角矩阵 向上利用X(k)=(Aug(k,N+1)-Aug(k,k+1:N)*X(k+1:N)/Aug(k,k)回代得到所求答案：input N:=6X = 0.5178 0.6065 0.6213 0.6213 0.6065 0.5178， 三角分解法：分析：利用指令L,U,P=lu(A)得出PA的分解矩阵L，Uinput N:=6L = 1.0000 0 0 0 0 0 -0.1667 1.0000 0 0 0 0 0 -0.1714 1.0000 0 0 0 0 0 -0.1716 1.0000 0 0 0 0 0 -0.1716 1.0000 0 0 0 0 0 -0.1716 1.0000U = 12.0000 -2.0000 0 0 0 0 0 11.6667 -2.0000 0 0 0 0 0 11.6571 -2.0000 0 0 0 0 0 11.6569 -2.0000 0 0 0 0 0 11.6569 -2.0000 0 0 0 0 0 11.6569 再构造N矩阵Y，令LYPB，再做增广矩阵AugLPB 利用向下回代公式Y(k)=(Aug(k,N+1)-Aug(k,k-1:-1:1)*Y(k-1:-1:1)/Aug(k,k)求出Y的值：Y = 5.0000 5.8333 6.0000 6.0294 6.0345 6.0354再令UXY，构造增广矩阵AugUY 利用向上回代公式X(k)=(Aug(k,N+1)-Aug(k,k+1:N)*X(k+1:N)/Aug(k,k)从而可以得到答案：X = 0.5178 0.6065 0.6213 0.6213 0.6065 0.5178 0.5178，雅可比迭代法：分析：利用雅可比迭代式： X(j)=(B(j)-A(j,1:j-1,j+1:N)*P(1:j-1,j+1:N)/A(j,j)迭代次数，精度为。X矩阵的初始值为P，P为N的零矩阵从而可以求的答案：input N:=6X = 0.5177 0.6065 0.6213 0.6213 0.6065 0.5177，高斯德尔赛迭代法：分析：考虑雅克比迭代法收敛速度有限 可以运用新的迭代方式增快收敛速度X(j)=(B(j)-A(j,1:j-1)*X(1:j-1)-A(j,j+1:N)*P(j+1:N)/A(j,j);令X的初始矩阵为P，P为N的零矩阵迭代次数，精度为。从而求的答案为：input N:=6X = 0.5177 0.6065 0.6213 0.6213 0.6065 0.51784. 结论四种方法所求答案基本相等，在N（未知数个数）很小时，程序运行速度基本相等。体现不出优缺点。在N很大时法3、法4优势明显，尤其是法4。代码：1，高斯消元法：N=input(input N:=); %输入N的值A=zeros(N,N); %定义A为NN的矩阵a=-2;c=-2;d=12;b=5;A(1,1)=d;A(1,2)=c;A(N,N-1)=a;A(N,N)=d; for k=2:N-1 A(k,k-1)=a; A(k,k)=d; A(k,k+1)=c; End %构造系数矩阵AB=zeros(N,1); %定义B为N的矩阵for k=1:N B(k,1)=b; %构造矩阵Bend Aug=A B; %Aug为A，B的增广矩阵 for k=2:N m=Aug(k,k-1)/Aug(k-1,k-1); Aug(k,k:N+1)=Aug(k,k:N+1)-m*Aug(k-1,k:N+1); End %向上消元得到Aug为一个下三角矩阵 X=zeros(N,1); %定义X为N的矩阵 X(N)=Aug(N,N+1)/Aug(N,N); %求出X第一行的值 for k=N-1:-1:1 X(k)=(Aug(k,N+1)-Aug(k,k+1:N)*X(k+1:N)/Aug(k,k); End X %向上回代的到X其余行的值2，三角分解法：N=input(input N:=); %输入N的值A=zeros(N,N); %定义A为NN的矩阵a=-2;c=-2;d=12;b=5;A(1,1)=d;A(1,2)=c;A(N,N-1)=a;A(N,N)=d;for k=2:N-1 A(k,k-1)=a; A(k,k)=d; A(k,k+1)=c;End %构造系数矩阵AB=zeros(N,1); %定义B为N的矩阵for k=1:N B(k,1)=b;end %构造矩阵BL,U,P=lu(A) %得出P-1*A的分解矩阵L，UAug=L P*B; % Aug为L，P*B的增广矩阵Y=zeros(N,1); %定义Y为N的矩阵 Y(1)=Aug(1,N+1)/Aug(1,1); %求出Y的第一行的值for k=2:N Y(k)=(Aug(k,N+1)-Aug(k,k-1:-1:1)*Y(k-1:-1:1)/Aug(k,k);End%向下回代求出Y其余行的值 Y Aug=U Y; % Aug为U，Y的增广矩阵 X=zeros(N,1); %定义X为N的矩阵 X(N)=Aug(N,N+1)/Aug(N,N); %求出X的第一行的值 for k=N-1:-1:1 X(k)=(Aug(k,N+1)-Aug(k,k+1:N)*X(k+1:N)/Aug(k,k); end%向上回代求出X其余行的值 X3，雅可比迭代法：clearN=input(input N:=);%输入N的值A=zeros(N,N); %定义A为NN的矩阵a=-2;c=-2;d=12;b=5;A(1,1)=d;A(1,2)=c;A(N,N-1)=a;A(N,N)=d;for k=2:N-1 A(k,k-1)=a; A(k,k)=d; A(k,k+1)=c;end %构造系数矩阵AB=zeros(N,1); %定义B为N的矩阵for k=1:N B(k,1)=b;end %构造矩阵B P=zeros(N,1);%定义P为N的矩阵为X的初始值N=length(B);for k=1:30%进行迭代次数 for j=1:N X(j)=(B(j)-A(j,1:j-1,j+1:N)*P(1:j-1,j+1:N)/A(j,j); End%利用雅可比迭代式求出行矩阵X的值 err=abs(norm(X-P);%得出相对误差 relerr=err/(norm(X)+eps); %得出绝对误差 P=X; if(err10-4|(relerr10-4) %当误差小与精度时停止迭代 break endendX=X 4，高斯-赛德尔迭代法：clearN=input(input N:=); %输入N的值A=zeros(N,N); %定义A为NN的矩阵a=-2;c=-2;d=12;b=5;A(1,1)=d;A(1,2)=c;A(N,N-1)=a;A(N,N)=d;for k=2:N-1 A(k,k-1)=a; A(k,k)=d; A(k,k+1)=c;end%构造系数矩阵AB=zeros(N,1); %定义B为N的矩阵for k=1:N B(k,1)=b;end %构造矩阵BP=zeros(N,1);%定义P为N的矩阵为X的初始值for k=1:30 for j=1:N if j=1 X(1)=(B(1)-A(1,2:N)*P(2:N)/A(1,1); %求出X（）的值 elseif j=N X(N)=(B(N)-A(N,1:N-1)*(X(1:N-1)/A(N,N); %求出X（N）的值 el