Gauss列主元消去法程序设计

《Gauss列主元消去法》实验报告

实验名称：Gauss列主元消去法程序设计成绩：___________

专业班级：数学与应用数学1202班姓名：王晓阳学号：2012254010228

实验日期： 2014 年11月10日

实验报告日期： 2014年 11月10日

一．实验目的

1.学习Gauss消去法的基本思路和迭代步骤.

2.学会运用matlab编写高斯消去法和列主元消去法程序，求解线性方程组.

3.当()k

a绝对值较小时，采用高斯列主元消去法.

kk

4.培养编程与上机调试能力.

二、实验内容

用消去法解线性方程组的基本思想是用逐次消去未知数的方法把原线性方程组Ax b

=化为与其等价的三角形线性方程组，而求解三角形线性方程组可用回代的方法求解.

1.求解一般线性方程组的高斯消去法.

(1)消元过程：

设()0k kk a ≠，第i 个方程减去第k 个方程的()()/k k ik ik kk m a a =倍，(1,

,)i k n =+，得到()()11k k A x b ++=.

()()()()()()()11,,1,,k k k ij ij ik kj k k k i i ik k a a m a i j k n b b m b ++⎧=-=+⎪⎨=-⎪⎩

经过n-1次消元，可把方程组()()11A x b =化为上三角方程组()()n n A x b =.

(2)回代过程：

()()()()()1//,1,

,1n n n n nn n i i i i i ij j ii j i x b a x b a x a i n =+⎧=⎪⎛⎫⎨=-=- ⎪⎪⎝⎭⎩∑

以解如下线性方程组为例测试结果.

1212312310773264556x x x x x x x x -=⎧⎪-++=⎨⎪-+=⎩

2.列主元消去法

由高斯消去法可知，在消元过程中可能出现()0k kk a =的情况，这是消去法将无法进行，即使主元素()0k kk a ≠但很小时，用其作除数，会导致其他元素数量级的严重增长和舍入误差的扩散，最后也使得计算解不可靠.这时就需要选取主元素，假定线性方程组的系数矩阵A 是菲奇异的.

(1)消元过程：

对于1,2,,1k n =-，进行如下步骤：

1) 按列选主元，记

max pk ik k i n

a a ≤≤=

2) 交换增广阵A 的p,k 两行的元素。

A(k,j)=A(p,j) (j=k,…,n +1)

3)交换常数项b 的p,k 两行的元素。

b(k)=b(p)

4)计算消元

()()()()()()()11,,1,,k k k ij ij ik kj k k k i

i ik k a a m a i j k n b b m b ++⎧=-=+⎪⎨=-⎪⎩

(2)回代过程

()()()()()1//,1,

,1n n n n nn n i i i i i ij j ii j i x b a x b a x a i n =+⎧=⎪⎛⎫⎨=-=- ⎪⎪⎝⎭⎩∑

(3)以解如下线性方程组为例测试结果.

12120.00001122

x x x x +=⎧⎨+=⎩ 三、实验环境

MATLAB R2014a

四、实验步骤

1.高斯列主元消去法流程图：

2.程序设计：

(一)高斯消去法：

a=input('请输入系数阵：');

b=input('请输入常数项：');

. n=length(b);

A=[a,b];

x=zeros(n,1); %初始值

for k=1:n-1

for i=k+1:n %第k次消元

m(i,k)=A(i,k)/A(k,k);

for j=k+1:n

A(i,j)=A(i,j)-A(k,j)*m(i,k);

end

b(i)=b(i)-m(i,k)*b(k);

end

end

x(n)=b(n)/A(n,n); %回代

for i=n-1:-1:1;

s=0;

for j=i+1:n;

s=s+A(i,j)*x(j);

end

x(i)=(b(i)-s)/A(i,i)

end

(二)高斯列主元消去法：

a=input('请输入系数阵：');

b=input('请输入常数项：');

n=length(b);

A=[a,b];

x=zeros(n,1); %初始值

for k=1:n-1

if abs(A(k,k))<10^(-4);%判断是否选主元y=1

else

y=0;

end

if y; %选主元

for i=k+1:n;

if abs(A(i,k))>abs(A(k,k))

p=i;

else p=k;

end

end

if p~=k;

for j=k:n+1;

s=A(k,j);

A(k,j)=A(p,j);%交换系数

A(p,j)=s;

end

t=b(k);

b(k)=b(p);%交换常数项

b(p)=t;

end

end

for i=k+1:n

m(i,k)=A(i,k)/A(k,k); %第k次消元

for j=k+1:n

A(i,j)=A(i,j)-A(k,j)*m(i,k);

end

b(i)=b(i)-m(i,k)*b(k);

end

end

x(n)=b(n)/A(n,n); %回代

for i=n-1:-1:1;

s=0;

for j=i+1:n;

s=s+A(i,j)*x(j);

end

x(i)=(b(i)-s)/A(i,i)