数值分析报告.doc

中国矿业大学工力13-02班02130857刘志强
第PAGE\*MERGEFORMAT8页

计算方法实验报告

实验：求解线性方程组的两种方法
班级：工力13-02
姓名：刘志强
学号：02130857

实验内容
分别用列主元素法和LU分解法编程求解，并对A或b做微小改动后观察结果

1-12-10 6
10110 4
213-44 X= -2
0-11-14 5
37823 1

实验原理
列主元素法

方法说明（以4阶为例）：

第1步消元——在增广矩阵（A，b）第一列中找到绝对值最大的元素，将其所在行与第一行交换，再对（A，b）做初等行变换使原方程组转化为如下形式：

第2步消元——在增广矩阵（A，b）中的第二列中（从第二行开始）找到绝对值最大的元素，将其所在行与第二行交换，再对（A，b）做初等行变换使原方程组转化为：

第3步消元——在增广矩阵（A，b）中的第三列中（从第三行开始）找到绝对值最大的元素，将其所在行与第二行交换，再对（A，b）做初等行变换使原方程组转化为：

按x4x3x2x1的顺序回代求解出方程组的解。

LU分解法
将HYPERLINK"http://baike.baidu.com/view/3503330.htm"\t"http://baike.baidu.com/_blank"系数矩阵A转变成等价两个矩阵L和U的乘积，其中L和U分别是下三角和上HYPERLINK"http://baike.baidu.com/view/3372514.htm"\t"http://baike.baidu.com/_blank"三角矩阵。当A的所有HYPERLINK"http://baike.baidu.com/view/1643301.htm"\t"http://baike.baidu.com/_blank"顺序主子式都不为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,N
integerid_max
real*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)=A
Ab(:,N+1)=b

dok=1,N-1
elmax=dabs(Ab(k,k))
id_max=k

!找最大元素标号

doi=k+1,n
if(dabs(Ab(i,k))>elmax)then
elmax=Ab(i,k)
id_max=i
endif
enddo

!交换元素，其他不变
vtemp1=Ab(k,:)
vtemp2=Ab(id_max,:)

Ab(k,:)=vtemp2
Ab(id_max,:)=vtemp1

!消元法
doi=k+1,N

temp=Ab(i,k)/Ab(k,k)

Ab(i,:)=Ab(i,:)-temp*Ab(k,:)

enddo
enddo

Aup(:,:)=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,N
real*8A(N,N),b(N),x(N)
x(N)=b(N)/A(N,N)
!回带部分
doi=n-1,1,-1

x(i)=b(i)
doj=i+1,N
x(i)=x(i)-a(i,j)*x(j)
enddo
x(i)=x(i)/A(i,i)
enddo
endsubroutineuptri

programmain
integer,parameter::N=5 !确定N
integer::i,j
real*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,*)b
callsolve(A,b,x,N)
!输出结果到DataOut_G