列主元高斯消去 LU分解 迭代法

数学与信息科学学院
实验报告
课程名称：《数值计算方法》
实验名称：数值积分
实验类型：验证性■综合性□设计性□
实验室名称：数学实验室
班级学号： 090721
学生姓名：
任课教师（教师签名）：
成绩：
实验日期：2012-4-27
一、实验目的及题目
熟悉线性方程组求解原理，运用列主消元高斯消去法，LU 分解法及Jacobi 迭代法与Gauss-Seidel 迭代法丢出线性方程的解 实验题目：
1、用列主元消去法解方程组：
⎪⎪⎩⎪⎪⎨⎧=-++--=+--=+-+=++4
323331243432143214
321421x x x x x x x x x x x x x x x
2、用LU 分解法界方程组b Ax =，其中
⎪⎪⎪⎪
⎪⎭
⎫ ⎝⎛--=⎪⎪⎪⎪⎪⎭⎫ ⎝⎛----=2244,1626622060121224121202448b A
3、分别用jacobi 迭代法和gauss-seidel 迭代法求解方程组：
⎪⎪⎩⎪⎪⎨⎧=+-+-=+--=+--=+-25113610211381121043213214
32321x x x x x x x x x x x x x
二、实验原理、程序框图、程序代码等
实验原理：
1、列主元高斯消去原理：
每次消去中，选取每列的绝对值最大的元素作为消去对象并作为主元素。

然后换行使之变到主元位子上，在进行消元计算。

设)()(k k b x A =，确定第k 列主元所在位置k i ，在交换k i 行和k 行后，在进行消元。

根据矩阵理论，交换k i 和k 两个方程的位置，列主元素的消去过程相当于对交换后的新矩阵进行消元，即
)1(,+=k k i k k A A I L k
同时，右端向量)(k b 变化为 )1(,+=k k i k k b b I L k
2、LU 分解原理：
直接三角分解：先将A 分解为上三角矩阵U 和下三角矩阵L （A=LU ），则原为题等价于求解两个方程组：y Ux b Ly ==, 3、迭代法
（1）jacobi 迭代法（又称简单迭代法） 考虑n 阶线性代数方程组
⎪⎪⎩⎪⎪⎨⎧=+++=+++=+++n
n nn n n n n n n b x a x a x a b x a x a x a b x a x a x a (22112)
22221211
1212111
其矩阵形式为
b Ax =
设该方程组的系数矩阵A 非奇异且),.....,2,1(0n i a ii =≠，可将A 分解为：
U L D A --=
其中],.....,,[2211nn a a a diag D =；
-=L ⎥⎥⎥⎥⎥⎥⎦
⎤
⎢⎢⎢⎢⎢⎢⎣⎡000
1
1
1
3131
21
n n n a a a a a a
，U=－⎥⎥
⎥
⎥⎥⎥⎦⎤⎢⎢
⎢⎢⎢⎢⎣⎡-00
00,1223
11312
n n n n a a a a a a
然后化为如下等价形式
b D x U L D x 11)(--++=
简记为 J J f x B x += 选取初始向量[
]T
n
x x x x )0()
0(2
)0(1)0(, =,通过上式可得到线性代数方程组的迭代格式
)2,1,0(,)()1( =+=+k f x B x J k J k
其分量形式为 )2,1,0)(,2,1(),(11
1
)
1( ==-=∑≠=+k n i x a b a x n
j i k j ij i ii k i 即
⎪⎪⎪
⎪⎩
⎪⎪⎪
⎪⎨⎧----=----=----=--+++)(1)(1)(1)
(11,)(22)(11)1()
(2)(323)(121222
)1(2)(1)(313)(212111
)1(1
k n n n k n k n n nn k n k n
n k k k k n n k k k x a x a x a b a x x a x a x a b a x x a x a x a b a x
以上计算过程称迭代法，矩阵J B 陈为jacobi 迭代法的迭代矩阵。

（2）Gauss-Seidel 迭代法 将jacobi 迭代格式改为如下形式 )2,1,0)(,2,1(),(1111
)()()
1( ==--=∑∑-=+=+k n i x a x a b a x i i n
i j k j ij k j ij i ii k i
其矩阵形式为
b L D Ux L D x k k 1)(1)1()()(--+-+-=
于是得Gauss-Seidel 迭代公式为
)2,1,0(,)()( =+=k f x B x G k G k 称G B 为Gauss-Seidel 迭代法的迭代矩阵。

程序代码： 第一题程序代码： function GEpiv(A,b) [m,n]=size(A); nb=n+1;Ab=[A b]; for i=1:m-1
[pivot,p]=max(abs(Ab(i:n,j))); ip=p+i-1; if ip~=i
Ab([i ip],:)=Ab([ip i],:); end
pivot=Ab(i,i);
for k=i+1:m
Ab(k,i:nb)=Ab(k,i:nb)-(Ab(k,i)/pivot)*Ab(i,i:nb); end
end
x=zeros(n,1);
x(n)=Ab(n,nb)/Ab(n,n);
for i=n-1:-1:1
x(i)(Ab(i,nb)-Ab(i,i+1:n)*x(i+1:n,1))/Ab(i,i);
end
for k=1:n
fprintf('x[%d]=%f\n',x(k));
end
运行结果为
输入
A=[1 1 0 3;2 1 -1 1;3 -1 -1 3;-1 2 3 -1];
b=[4;1;-3;4];
得结果为
x =
-1.3333
2.3333
-0.3333
1.0000
第二题代码：
function x=lufj(A,b) %lu分解求值函数
[l u p]=lu(A);
y=lx(l,b);
x=us(u,y);
function y=lx(A,b) %求系数矩阵为下三角方程函数[n,m]=size(A);
y(1)=b(1)/A(1,1);
for i=2:n
s=0;
for k=1:i-1
s=s+A(i,k)*y(k);
end
y(i)=(b(i)-s)/A(i,i);
end
end
function x=us(A,b) %求系数矩阵为上三角方程函数[n,m]=size(A);
x(n)=b(n)/A(n,n);
for i=0:n-2
s=0;
for k=0:i
s=s+A(n-i-1,n-k)*x(n-k);
end
x(n-i-1)=(b(n-i-1)-s)/A(n-i-1,n-i-1);
end
end
输入矩阵
A=[48 -24 0 12;-24 4 12 12;0 6 20 2;-6 6 2 16]; b=[4;4;-2;-2];
结果为
0.5212 1.0055 -0.3757 -0.2597
第三题代码：
（1） Jacobi迭代：
function x = agui_jacobi(a,b)
n=length(b);
N=100;
e=1e-4;
x0=zeros(n,1);
x=x0;
x0=x+2*e;
k=0;
d=diag(diag(a));
l=-triu(a,-1);
u=-triu(a,1);
while norm(x0-x,inf)>e&k<N
k=k+1;
x0=x;
x=inv(d)*(1+u)*x+inv(d)*b;
k
disp(x')
end
if k==N
end
（2）Gauss-seidel迭代
function x=gau(A,b,x)
[d u l]=fenjie(A);
B=inv(d-l)*u;
for i=1:9
y=x;
x=B*y+f;
end
function [d u l]=fenjie(A)
[n m]=size(A);
for i=1:n
for j=1:n
d(i,j)=0;
end
end
for i=1:n
d(i,i)=A(i,i);
end
for i=1:n
for j=1:i
l(i,j)=-A(i,j);
u(j,i)=-A(j,i);
end
for k=i:n
l(i,k)=0;
u(k,i)=0;
end
end
（1）
输入矩阵：a=[10 -1 2 0;0 8 -1 3;2 -1 10 0; -1 3 -1 11]; >> b=[-11;-11;6;25];
>> x = agui_jacobi(a,b)
x=
-1.4674 -2.3587 0.6576 2.8424
（2）A=[10 -1 2 0;0 8 -1 3;2 -1 10 0;-1 3 -1 11];
b=[-11;-11;6;25];
x=[-1.3;-2.1;0.55;2.66];
gau(A,b,x)
结果为
x=
-1.4674 -2.3587 0.6576 2.8424
四、实验中存在的问题及解决方案
命令需多次修改，如Gauss迭代中，m文件中因end的不匹配使得程序无法运行处结果，则需细心调试，去除多余的end，才能使程序正常运行。

五、心得体会
在线性代数中，从理论上解决了线性方程组的求解问题，但常用的高斯消元法有时会导致结果产生严重的偏差，因此需要研究数值解法。

Jacobi迭代与Gauss-Seidel有一个共同特点：新的近似解是已知近似解的线性函数，并且新的近似解只与已知近似解有关。

编写程序是一个繁杂的过程，需要不断调试，并且善于发现程序中的细小问题，才能够正常运行处程序。