2904203015.李亚韦.实验三.doc

三角线性方程组的数值计算 2904203015 李亚韦一、 实验描述 求解三角线性方程组，可以通过已知函数从只有一个未知函数的式子开始，依次求解出每个未知量即采用消元法和LU分解法求解；也可以通过迭代近似的方式求解出，而迭代的方法大致有雅克比迭代、高斯-赛德尔迭代等方法。二、 实验内容 假设d1=d2=dn=12，b1=b2=bn=5;c1=c2=.=cn=a1=a2=an=-2;的情况下分别用高斯消元法，LU分解法，雅可比迭代法和高斯-赛德尔迭代法求解A矩阵为：B矩阵为：B=5;5;5;.;5;5;5;高斯消元法:根据题意，构造上述的A矩阵和B矩阵，然后通过高斯消元的过程求解出最后的变量值。流程图：开始i=1,j=1消去i行j列的元i9得消去i行后的矩阵j9输出消元后的矩阵结束YYNN流程图 1实验程序：n=input(input a n=: );A=zeros(n,n);A(1,1)=12;A(1,2)=-2;A(n,n-1)=-2;A(n,n)=12;for i=2:n-1 A(i,i)=12; A(i,i-1)=-2; A(i,i+1)=-2;end %构造A矩阵B=ones(n,1)*5; %构造B矩阵R=A B; %构造增广矩阵Rfor i=1:n-1 if R(i,i)=0 break end for j=i+1:n m=R(j,i)/R(i,i); R(j,i:n+1)=R(j,i:n+1)-m*R(i,i:n+1); endend %构造三角矩阵x=zeros(n,1);x(n)=B(n)/A(n,n);for k=n-1:-1:1 x(k)=(B(k)-A(k,k+1:n)*x(k+1:n)/A(k,k);endx(1:n) %求解X向量用LU分解求解:LU分解为非直接求解，设A为nn矩阵，由一般线性方程组如果存在置换矩阵P，使得可得以下关系式我们有替换矩阵对方程组求得Y然后通过会带求得X。通过以上分析，先通过上三角求解法求得L，U矩阵，然后求得Y ,最后即可得到X。实验流程图：开始I=1,j=1消去i行，j列的元i9得到消去i行的矩阵j9得到L,U矩阵得到Y向量输出X向量结束YYNN流程图 2实验程序：n=input(input a n=: );A=zeros(n,n);A(1,1)=12;A(1,2)=-2;A(n,n-1)=-2;A(n,n)=12;for i=2:n-1 A(i,i)=12; A(i,i-1)=-2; A(i,i+1)=-2;end %构造A矩阵B=ones(n,1)*5; %构造B矩阵for i=1:n-1 if A(i,i)=0 break end for j=i+1:n m=A(j,i)/A(i,i); A(j,i)=m; A(j,i+1:n)=A(j,i+1:n)-m*A(i,i+1:n); endend x=zeros(n,1);y=zeros(n,1);y(1)=B(n);for k=2:n y(k)=B(n-k+1)-A(k,1:k-1)*y(1:k-1); %求y向量的解endx(n)=y(n)/A(n,n);for k=n-1:-1:1 x(k)=(y(k)-A(k,k+1:n)*x(k+1:n)/A(k,k);endx(1:n) %求x向量的解雅可比迭代法：向构造A矩阵，然后根据雅可比迭代过程：求解出向量x的值。实验流程图：开始n,i=1,j=1得矩阵Ain,jn输出结束YN流程图 3实验程序：n=input(input a n=: );A=zeros(n,n);A(1,1)=12;A(1,2)=-2;A(n,n-1)=-2;A(n,n)=12;for i=2:n-1 A(i,i)=12; A(i,i-1)=-2; A(i,i+1)=-2;end %构造A矩阵B=ones(n,1)*5; %构造B矩阵x=zeros(n,1); for i=1:n for j=1:n x(j)=(B(j)-A(j,1:j-1,j+1:n)*x(1:j-1,j+1:n)/A(j,j); endendx %求解向量x 高斯赛德尔迭代：先求得A矩阵，然后通过高速-赛德尔迭代式：其中为矩阵A,B中的对应值，然后即可得到x向量的值。实验流程图：开始n,i=1,j=1得矩阵Ain,jnYK=1输入k输出x(1)K=n输出x(n)输出其余x结束流程图 4实验程序：n=input(input a n=: );A=zeros(n,n);A(1,1)=12;A(1,2)=-2;A(n,n-1)=-2;A(n,n)=12;for i=2:n-1 A(i,i)=12; A(i,i-1)=-2; A(i,i+1)=-2;end %构造A矩阵B=ones(n,1)*5; %构造B矩阵x=zeros(n,1); for j=1:n if j=1 x(1)=(B(1)-A(1,2:n)*x(2:n)/A(1,1); elseif j=n x(n)=(B(n)-A(n,1:n-1)*(x(1:n-1)/A(n,n); else x(j)=(B(j)-A(j,1:j-1)*x(1:j-1)-A(j,j+1:n)*x(j+1:n)/A(j,j); end end x %求解X向量当输入对应的n即可得到对应个数的x。三、 实验结果与分析当n为3时，所得x的结果如下：表 1X1X2X30.51470.58820.5147当n为5时，所得x的结果如下：表 2X1X2X3X4X50.51770.60610.61870.60610.5177当n为10时，所得x的结果如下：表 3X1X2X3X4X50.51780.60660.62180.62450.6249X6X7X8X9X100.62490.62450.62180.60660.5178当n的值为20时，所得结果如下：X1X2X3X4X50.51780.60660.62180.62450.6249X6X7X8X9X100.62500.62500.62500.62500.6250X11X12X13X14X150.62500.62500.62500.62500.6250X16X17X18X19X200.62490.62450.62180.60660.5178 根据结