数值分析上机实验报告七.docVIP

实验报告七 题目： 线性方程组的直接解法 摘要：在数值计算中，关于线性方程组的解法有两类：直接法和迭代法。本实验研究线性方程组的直接解法。 前言：（目的和意义） 掌握线性方程组的直接解法，如高斯消去法，LU解法，平方根法和列住元消去法等。 数学原理： 高斯消去法： 第一步先进行消元计算 第二步进行回代计算 平方根分解法： 设A为对称阵，且A的所有顺序主子式均不为零，利用A的对称性将U再分解，即 其中D为对角阵，为单位上三角阵。于是 又由分解的唯一性得，代入(1)式得到对称矩阵A的分解式 程序设计： 本实验采用Matlab的M文件编写。 对于Hilbert矩阵，用LU分解的源程序是； i=n;%%n代表阶数 h=hilb(i) [L,U]=lu(h) 平方根法: 直接在Matlab命令窗口编写 先对A矩阵作分解: L=chol(h) 检验其正确性 L*L 将L转化为下三角矩阵: L=L 利用LU分解法求解线性方程组(十阶Hilbert矩阵)，程序如下： H=[1,1/2,1/3,1/4,1/5,1/6,1/7,1/8,1/9,1/10;1/2,1/3,1/4,1/5,1/6,1/7,1/8,1/9,1/10,1/11;1/3,1/4,1/5,1/6,1/7,1/8,1/9,1/10,1/11,1/12;1/4,1/5,1/6,1/7,1/8,1/9,1/10,1/11,1/12,1/13;1/5,1/6,1/7,1/8,1/9,1/10,1/11,1/12,1/13,1/14;1/6,1/7,1/8,1/9,1/10,1/11,1/12,1/13,1/14,1/15;1/7,1/8,1/9,1/10,1/11,1/12,1/13,1/14,1/15,1/16;1/8,1/9,1/10,1/11,1/12,1/13,1/14,1/15,1/16,1/17;1/9,1/10,1/11,1/12,1/13,1/14,1/15,1/16,1/17,1/18;1/10,1/11,1/12,1/13,1/14,1/15,1/16,1/17,1/18,1/19] b=[1/2,1,1/3,1/4,0,0,1/2,0,1,1]; [L,U]=lu(H); x=U\(L\b) 直接建立求解该方程组的M文件Gauss.m如下: A=[0.4096,0.1234,0.3678,0.2943;0.2246,0.3872,0.4015,0.1129;0.3645,0.1920,0.3781,0.0643;0.1784,0.4002,0.2786,0.3927]; b=[0.4043,0.1550,0.4240,-0.2557]; [m,n]=size(A); if m~=n error; return; end if m~=size(b) error; return; end if rank(A)~=rank([A,b]) error; return; end c=n+1; A(:,c)=b; for k=1:n-1 A(k+1:n, k:c)=A(k+1:n, k:c)-(A(k+1:n,k)/ A(k,k))*A(k, k:c); end x=zeros(length(b),1); x(n)=A(n,c)/A(n,n); for k=n-1:-1:1 x(k)=(A(k,c)-A(k,k+1:n)*x(k+1:n))/A(k,k); end disp(x=); disp(x); 直接建立求解该方程的M文件Gauss_line.m,求解程序编制如下: A=[0.4096,0.1234,0.3678,0.2943;0.2246,0.3872,0.4015,0.1129;0.3645,0.1920,0.3781,0.0643;0.1784,0.4002,0.2786,0.3927]; b=[0.4043,0.1550,0.4240,-0.2557]; [m,n]=size(A); if m~=n error; return; end if m~=size(b) error; return; end if rank(A)~=rank([A,b]) error; return; end c=n+1; A(:,c)=b; for k=1:n-1 [r,m]=max(abs(A(k:n,k))); m=m+k-1; if(A(m,k)~=0) if(m~=k)