Gauss列主元消去法

贵州师范大学数学与计算机科学学院学生实验报告
课程名称： 数值分析 班级： 实验日期： 年 月 日 学 号： 姓名： 指导教师： 实验成绩：
一、实验名称
实验五:线性方程组的数值解法
二、实验目的及要求
1. 让学生掌握用列主元gauss 消去法、超松弛迭代法求解线性方程组.
2. 培养Matlab 编程与上机调试能力.
三、实验环境
每人一台计算机，要求安装Windows XP 操作系统，Microsoft office2003、MATLAB6.5(或7.0).
四、实验内容
1. 编制逐次超松弛迭代(SOR 迭代)函数(子程序),并用于求解方程组
⎪⎪⎩⎪⎪⎨⎧=-++=+-+=++-=+++-1
41
414144321432143214321x x x x x x x x x x x x x x x x
取初始向量T x )1,1,1,1()0(=,迭代控制条件为 5)1()(102
1||||--⨯≤
-k k x x 请绘制出迭代次数与松弛因子关系的函数曲线,给出最佳松弛因子.SOR 迭代的收敛速度是否一定比Gauss-Seidel 迭代快?
2. 编制列主元 Gauss 消去法函数(子程序),并用于解 ⎪⎩⎪⎨⎧=++-=-+-=+-615318153312321
321321x x x x x x x x x
要求输出方程组的解和消元后的增广矩阵. 注：题2必须写实验报告
五、算法描述及实验步骤
Gauss 消去法：
功能 解方程组b Ax = .
输入 n ,n n ij a A ⨯=)(,T n b b b b ),,,(21 =.
输出 方程组的解T n x x x x ),,,(21 =或失败信息.
步1 对1,,2,1-=n k 执行步2→步4 .
步2 调选列主元模块 .
步3 若0=kk a ，则=x “消去法失败”，结束 .
步4 对n k k i ,,2,1 ++=执行步5→步6 .
步5 对n k k j ,,2,1 ++=执行ij kj kk ik ij a a a a a +⨯-⇐/ .
步6 i k kk ik i b b a a b +⨯-⇐/ .
步7 nn n n a b x /⇐ .
步8 对1,,2,1 --=n n i 执行ii n i j j ij i i a x a b x /)(1∑+=-
⇐ .
步9 输出T n x x x x ),,,(21 = .
选列主元模块：
功能 选列主元 .
输入 n k k i b n k k j i a i ij ,,1,,;,,1,,, +=+= .
输出 n k k i b n k k j i a i ij ,,1,,;,,1,,, +=+= .
步1 kk a m ⇐；k l ⇐ .
步2 对n k k i ,,2,1 ++=执行若m a ik >则ik a m ⇐；i l ⇐ .
步3 若k l ≠，则交换kj a 和lj a ，n k k j ,,1, +=；交换k b 和l b .
步4 返回主模块 .
六、调试过程及实验结果
>> A=[12,-3,3;-18,3,-1;1,1,1];
>> b=[15;-15;6];
>> x=Gauss1(A,b)
Ab =
-18.0000 3.0000 -1.0000 -15.0000
0 1.1667 0.9444 5.1667
0 0 3.1429 9.4286 index = 1
x = 1.0000 2.0000 3.0000
七、总结
由于数)1(-k kk
a 在Gauss 消去法中有着突出的作用，第k 步消元时，要用)1(-k kk a 作除数，如果)1(-k kk a =0消元会失败，即使主元)1(-k kk a ≠0，但很小时，舍入误差也会使计算结果面目全非，避免这种缺陷的基本方法就是选主元。

通过选主元，就可避免绝对值小的数作除数，从而避免舍入误差的恶性增长，使得Gauss 列主元消去法是解中小规模的线性方程组和某些大型稀疏线性方程组的有效方法。

八、附录（源程序清单）
function [x,index]=Gauss1(A,b)
[n,m]=size(A);x=zeros(n,1); index=1
for k=1:n-1
a_max=0;
for i=k:n
if abs(A(i,k))>a_max
a_max=abs(A(i,k));r=i;
end
end
if a_max<1e-10
index=0;return ;
end
if r>k
for j=k:n
z=A(k,j);A(k,j)=A(r,j);A(r,j)=z;
end
z=b(k);b(k)=b(r);b(r)=z;
end
for i=k+1:n
m=A(i,k)/A(k,k);
for j=k:n
A(i,j)=A(i,j)-m*A(k,j);
end
b(i)=b(i)-m*b(k);
end
end
if abs(A(n,n))==0
index=0;return ;
end
Ab=[A,b]
x(n)=b(n)/A(n,n);
for i=n-1:-1:1
for j=i+1:n
b(i)=b(i)-A(i,j)*x(j);
end
x(i)=b(i)/A(i,i);
end。