数值分析报告

计算方法实验报告
实验：求解线性方程组的两种方法班级：工力13-02
姓名：刘志强
学号：02130857
实验内容
分别用列主元素法和LU 分解法编程求解，并对A 或b 做微小改动后观察结果
1 -1
2 -1 0 6 1 0 1 1 0
4 2 1 3 -4 4
X = -2 0 -1 1 -1 4
5 3 7 8 2 3
1
实验原理
列主元素法
方法说明（以4阶为例）：
⎥⎥⎥⎥⎦
⎤⎢⎢⎢⎢⎣⎡=⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡⋅⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡n n nn n n n n b b b x x x a a a a a a a a a 21212122221
11211 第1步消元——在增广矩阵（A ，b ）第一列中找到绝对值最大的元素，将其所在行与第一行交换，再对（A ，b ）做初等行变换使原方程组转化为如下形式：
⎥⎥⎥⎥⎦
⎤⎢⎢⎢⎢⎣⎡=⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡⋅⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡*******0***0***0****4321x x x x 第2步消元——在增广矩阵（A ，b ）中的第二列中（从第二行开始）找到绝对值最大的元素，将其所在行与第二行交换，再对（A ，b ）做初等行变换使原方程组转化为：
⎥⎥⎥⎥⎦
⎤⎢⎢⎢⎢⎣⎡=⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡⋅⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡******00**00***0****4321x x x x 第3步消元——在增广矩阵（A ，b ）中的第三列中（从第三行开始）找到绝对值最大的元素，将其所在行与第二行交换，再对（A ，b ）做初等行变换使原方程组转化为：
⎥⎥⎥⎥⎦
⎤⎢⎢⎢⎢⎣⎡=⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡⋅⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡*****000**00***0****4321x x x x 按x 4 → x 3→ x 2→ x 1 的顺序回代求解出方程组的解。

LU分解法
将系数矩阵A转变成等价两个矩阵L和U的乘积，其中L和U分别是下三角和上三角矩阵。

当A的所有顺序主子式都不为0时，矩阵A可以分解为A=LU，且当L的对角元全为1时分解唯一。

其中L是下三角矩阵，U是上三角矩阵。

1.2 程序设计：
列主元素法
! A-(N,N)系数矩阵b(N)-右向量N-方程维数
subroutine solve(A,b,x,N)
implicit real*8(a-z)
integer i,k,N
integer id_max
real*8 A(N,N),b(N),x(N)
real*8 Aup(N,N),bup(N)
real*8 Ab(N,N+1)
real*8 vtemp1(N+1),vtemp2(N+1)
Ab(1:N,1:N)=A
Ab(:,N+1)=b
do k=1,N-1
elmax=dabs(Ab(k,k))
id_max=k
!找最大元素标号
if (dabs(Ab(i,k))>elmax) then elmax=Ab(i,k)
id_max=i
end if
end do
!交换元素，其他不变
vtemp1=Ab(k,:)
vtemp2=Ab(id_max,:)
Ab(k,:)=vtemp2
Ab(id_max,:)=vtemp1
!消元法
do i=k+1,N
temp=Ab(i,k)/Ab(k,k)
Ab(i,:)=Ab(i,:)-temp*Ab(k,:)
end do
end do
Aup(:,:)=Ab(1:N,1:N)
bup(:)=Ab(:,N+1)
call uptri(Aup,bup,x,n)
end subroutine solve
! 上三角方程组回带subroutine uptri(A,b,x,N) implicit real*8(a-z) integer::i,j,N
real*8 A(N,N),b(N),x(N)
x(N)=b(N)/A(N,N)
!回带部分
do i=n-1,1,-1
x(i)=b(i)
x(i)=x(i)-a(i,j)*x(j)
end do
x(i)=x(i)/A(i,i)
end do
end subroutine uptri
program main
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
call solve(A,b,x,N)
!输出结果到DataOut_Gauss
write(12,100)x
100 format(T5,'列主元消去法计算结果',/,T4,'x=',6(/F12.8)) write(*,*)'列主元消去法计算结果输出至DataOut_Gauss' end
LU分解法
! A-(N,N)系数矩阵b(N)-右向量N-方程维数
subroutine solve(A,b,x,N)
implicit real*8(a-z)
integer::N
real*8::A(N,N),b(N),x(N)
real*8::L(N,N),U(N,N)
real*8::y(N)
call Resolve(A,L,U,N)
call downtri(L,b,y,N)
call uptri(U,y,x,N)
end subroutine solve
subroutine Resolve(A,L,U,N)
!LU分解
implicit real*8(a-z)
integer::N,i,k,r
real*8::A(N,N),L(N,N),U(N,N)
U(1,:)=A(1,:)
L(:,1)=a(:,1)/U(1,1)
do k=2,N
l(k,k)=1
do j=k,n
s=0
do m=1,k-1
s=s+l(k,m)*u(m,j)
end do
u(k,j)=a(k,j)-s
end do
do i=k+1,n
s=0
do m=1,k-1
s=s+l(i,m)*u(m,k)
end do
l(i,k)=(a(i,k)-s)/u(k,k) !wrong?
end do
end do
end subroutine Resolve
subroutine uptri(A,b,x,N)
! 上三角方程组的回带方法
implicit real*8(a-z)
integer::i,j,k,N
real*8::A(N,N),b(N),x(N)
x(N)=b(N)/A(N,N)
!回带部分
do i=n-1,1,-1
x(i)=b(i)
do j=i+1,N
x(i)=x(i)-a(i,j)*x(j)
end do
x(i)=x(i)/A(i,i)
end do
end subroutine uptri
subroutine downtri(A,b,x,N)
! 下三角方程组回带
implicit real*8(a-z)
integer::i,j,N
real*8::A(N,N),b(N),x(N)
x(1)=b(1)/a(1,1)
do k=2,N
x(k)=b(k)
do i=1,k-1
x(k)=x(k)-a(k,i)*x(i)
end do
x(k)=x(k)/a(k,k)
end do
end subroutine downtri
program main
integer,parameter::N=5
real*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,*)b
call solve(A,b,x,N)
write(12,101)x
101 format(T5,'LU分解计算线性方程组计算结果',/,T4,'x=',6(/F12.8)) write(*,*)'lu分解结果输出至'
end
数据处理
在输入文件中输入下图数据
分别运行两段程序得到下图数据
对输入数据微小改动为
计算结果为下图
结果分析：
计算结果是符合方程的。

比较两种计算方法，可看出列主元素法精度更高，这是因为此方法选取了绝对值大的数做主元，避免了和比较小的数相除，从而提高了算法的数值稳定性。

Lu分解法实质上是高斯消去法，可能会出现数值不稳定的情况。

但当需要求解多个系数矩阵为A的的线性方程组时，LU分解法显得特别优越，只要对系数矩阵做一次LU分解，接三角形方程组即可。

当对系数矩阵和右端向量做微小变化时，原方程的解变化微小，此为良态方程。