数值分析报告

中国矿业大学工力13-02班02130857刘志强第1页共11页计算方法实验报告实验：求解线性方程组的两种方法班级：工力13-02姓名：刘志强学号：02130857中国矿业大学工力13-02班02130857刘志强第2页共11页实验内容分别用列主元素法和LU分解法编程求解，并对A或b做微小改动后观察结果1-12-106101104213-44X=-20-11-145378231实验原理列主元素法方法说明（以4阶为例）：nnnnnnnnbbbxxxaaaaaaaaa2121212222111211第1步消元——在增广矩阵（A，b）第一列中找到绝对值最大的元素，将其所在行与第一行交换，再对（A，b）做初等行变换使原方程组转化为如下形式：*******0***0***0****4321xxxx第2步消元——在增广矩阵（A，b）中的第二列中（从第二行开始）找到绝对值最大的元素，将其所在行与第二行交换，再对（A，b）做初等行变换使原方程组转化为：******00**00***0****4321xxxx第3步消元——在增广矩阵（A，b）中的第三列中（从第三行开始）找到绝对值最大的元素，将其所在行与第二行交换，再对（A，b）做初等行变换使原方程组转化为：*****000**00***0****4321xxxx按x4x3x2x1的顺序回代求解出方程组的解。中国矿业大学工力13-02班02130857刘志强第3页共11页LU分解法将系数矩阵A转变成等价两个矩阵L和U的乘积，其中L和U分别是下三角和上三角矩阵。当A的所有顺序主子式都不为0时，矩阵A可以分解为A=LU，且当L的对角元全为1时分解唯一。其中L是下三角矩阵，U是上三角矩阵。1.2程序设计：列主元素法!A-(N,N)系数矩阵b(N)-右向量N-方程维数subroutinesolve(A,b,x,N)implicitreal*8(a-z)integeri,k,Nintegerid_maxreal*8A(N,N),b(N),x(N)real*8Aup(N,N),bup(N)real*8Ab(N,N+1)real*8vtemp1(N+1),vtemp2(N+1)Ab(1:N,1:N)=AAb(:,N+1)=bdok=1,N-1elmax=dabs(Ab(k,k))id_max=k!找最大元素标号中国矿业大学工力13-02班02130857刘志强第4页共11页doi=k+1,nif(dabs(Ab(i,k))elmax)thenelmax=Ab(i,k)id_max=iendifenddo!交换元素，其他不变vtemp1=Ab(k,:)vtemp2=Ab(id_max,:)Ab(k,:)=vtemp2Ab(id_max,:)=vtemp1!消元法doi=k+1,Ntemp=Ab(i,k)/Ab(k,k)Ab(i,:)=Ab(i,:)-temp*Ab(k,:)enddoenddoAup(:,:)=Ab(1:N,1:N)bup(:)=Ab(:,N+1)calluptri(Aup,bup,x,n)endsubroutinesolve!上三角方程组回带subroutineuptri(A,b,x,N)implicitreal*8(a-z)integer::i,j,Nreal*8A(N,N),b(N),x(N)x(N)=b(N)/A(N,N)!回带部分doi=n-1,1,-1x(i)=b(i)中国矿业大学工力13-02班02130857刘志强第5页共11页doj=i+1,Nx(i)=x(i)-a(i,j)*x(j)enddox(i)=x(i)/A(i,i)enddoendsubroutineuptriprogrammaininteger,parameter::N=5!确定Ninteger::i,jreal*8::A(N,N),b(N),x(N)open(unit=11,file='DataIn.txt')open(unit=12,file='DataOut_Gauss.txt')!从DataIn读取数据read(11,*)((A(i,j),j=1,N),i=1,N)read(11,*)bcallsolve(A,b,x,N)!输出结果到DataOut_Gausswrite(12,100)x100format(T5,'列主元消去法计算结果',/,T4,'x=',6(/F12.8))write(*,*)'列主元消去法计算结果输出至DataOut_Gauss'endLU分解法!A-(N,N)系数矩阵b(N)-右向量N-方程维数subroutinesolve(A,b,x,N)implicitreal*8(a-z)integer::Nreal*8::A(N,N),b(N),x(N)real*8::L(N,N),U(N,N)real*8::y(N)callResolve(A,L,U,N)calldowntri(L,b,y,N)calluptri(U,y,x,N)endsubroutinesolvesubroutineResolve(A,L,U,N)!LU分解implicitreal*8(a-z)integer::N,i,k,r中国矿业大学工力13-02班02130857刘志强第6页共11页real*8::A(N,N),L(N,N),U(N,N)U(1,:)=A(1,:)L(:,1)=a(:,1)/U(1,1)dok=2,Nl(k,k)=1doj=k,ns=0dom=1,k-1s=s+l(k,m)*u(m,j)enddou(k,j)=a(k,j)-senddodoi=k+1,ns=0dom=1,k-1s=s+l(i,m)*u(m,k)enddol(i,k)=(a(i,k)-s)/u(k,k)!wrong?enddoenddoendsubroutineResolvesubroutineuptri(A,b,x,N)!上三角方程组的回带方法implicitreal*8(a-z)integer::i,j,k,Nreal*8::A(N,N),b(N),x(N)x(N)=b(N)/A(N,N)!回带部分doi=n-1,1,-1x(i)=b(i)doj=i+1,Nx(i)=x(i)-a(i,j)*x(j)enddox(i)=x(i)/A(i,i)enddoendsubroutineuptri中国矿业大学工力13-02班02130857刘志强第7页共11页subroutinedowntri(A,b,x,N)!下三角方程组回带implicitreal*8(a-z)integer::i,j,Nreal*8::A(N,N),b(N),x(N)x(1)=b(1)/a(1,1)dok=2,Nx(k)=b(k)doi=1,k-1x(k)=x(k)-a(k,i)*x(i)enddox(k)=x(k)/a(k,k)enddoendsubroutinedowntriprogrammaininteger,parameter::N=5real*8::A(n,n),L(N,N),U(N,N)real*8::b(N),x(N)open(unit=11,file='DataIn.txt')open(unit=12,file='DataOut_Lu.txt')read(11,*)((A(i,j),j=1,N),i=1,N)read(11,*)bcallsolve(A,b,x,N)write(12,101)x101format(T5,'LU分解计算线性方程组计算结果',/,T4,'x=',6(/F12.8))write(*,*)'lu分解结果输出至'end数据处理在输入文件中输入下图数据中国矿业大学工力13-02班02130857刘志强第8页共11页分别运行两段程序得到下图数据对输入数据微小改动为计算结果为下图中国矿业大学工力13-02班02130857刘志强第9页共11页结果分析：计算结果是符合方程的。比较两种计算方法，可看出列主元素法精度更高，这是因为此方法选取了绝对值大的数做主元，避免了和比较小的数相除，从而提高了算法的数值稳定性。Lu分解法实质上是高斯消去法，可能会出现数值不稳定的情况。但当需要求解多个系数矩阵为A的的线性方程组时，LU分解法显得特别优越，只要对系数矩阵做一次LU分解，接三角形方程组即可。当对系数矩阵和右端向量做微小变化时，原方程的解变化微小，此为良态方程。