matlab求解代数方程组解析

第三讲 Matlab 求解代数方程组

理论介绍：直接法+迭代法，简单介绍相关知识和应用条件及注意事项 软件求解：各种求解程序

讨论如下表示含有n 个未知数、由n 个方程构成的线性方程组：

1111221121122222

1122n n n n n n nn n n

a x a x a x

b a x a x a x b a x a x a x b +++=⎧⎪+++=⎪⎨⎪

⎪+++=⎩ （1）

一、直接法 1.高斯消元法：

高斯消元法的基本原理： 在（1）中设110,a ≠将第一行乘以1

11

,k a a -

加到第(2,3,,),k k n = 得： (1)(1)(1)(1)11112211(2)(1)(2)

22112

(2)(2)(2)22n n n n n nn n n a x a x a x b a x a x b a x a x b ⎧+++=⎪++=⎪⎨⎪

⎪++=⎩

（2）

其中(1)

(1)1111,.k k a

a b b ==再设(2)22

0,a ≠将（2）式的第二行乘以(2)2

(2)22

,(3,,)k a k n a -= 加到第k 行，如此进行下去最终得到：

(1)(1)(1)(1)11112211(2)(1)(2)

22112(1)(1)(1)

1,111,1()()

n n n n n n n n n n n n n n n n nn n n a x a x a x b a x a x b a x a x b a x b --------⎧+++=⎪++=⎪

⎪

⎨⎪+=⎪

⎪=⎩

（3） 从（3）式最后一个方程解出n x ，代入它上面的一个方程解出1n x -，并如此进行

下去，即可依次将121,,,,n n x x x x - 全部解出，这样在()

0(1,2,,)k kk a k n ≠= 的假设

下，由上而下的消元由下而上的回代，构成了方程组的高斯消元法. 高斯消元法的矩阵表示：

若记11(),(,,),(,,)T T ij n n n n A a x x x b b b ⨯=== ，则（1）式可表为.Ax b =于是高斯

消元法的过程可用矩阵表示为：

121121.n n M M M Ax M M M b --=

其中：

(1)21

(1)

11

1(1)

1

(1)11

111n a a M a a ⎛⎫ ⎪ ⎪- ⎪=

⎪ ⎪ ⎪ ⎪- ⎪⎝⎭ (2)32

(2)222(2)

2(2)

221

111n a a M a a ⎛⎫

⎪

⎪ ⎪

-

⎪

=

⎪ ⎪ ⎪

⎪- ⎪⎝⎭

高斯消元法的Matlab 程序： %顺序gauss 消去法,gauss 函数 function[A,u]=gauss(a,n) for k=1:n-1

%消去过程 for i=k+1:n for j=k+1:n+1

%如果a(k,k)=0,则不能削去 if abs(a(k,k))>1e-6 %计算第k 步的增广矩阵 a(i,j)=a(i,j)-a(i,k)/a(k,k)*a(k,j); else

%a(k,k)=0,顺序gauss 消去失败 disp (‘顺序gauss 消去失败‘)； pause; exit; end end end end

%回代过程 x(n)=a(n,n+1)/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)=(a(i,n+1)-s)/a(i,i); end

%返回gauss 消去后的增广矩阵 A=triu(a); %返回方程组的解 u=x ;

练习和分析与思考： 用高斯消元法解方程组：

12345124512345124512452471523814

476192536x x x x x x x x x x x x x x x x x x x x x x ++++=⎧⎪+++=⎪⎪

++++=⎨⎪+++=⎪

+++=⎪⎩

2.列主元素消元法

在高斯消元法中进行到第k 步时，不论()

k ik a 是否为0，都按列选择()

||(,,)k ik a i k n = 中最大的一个，称为列主元，将列主元所在行与第k 行交换再

按高斯消元法进行下去称为列主元素消元法。

列主元素消元法的matlab 程序 %列主元guass 消去函数 function[A,u]=gauss(a,n) %消去过程 for k=1:n-1 %选主元

c=0;

for q=k:n

if abs(a(q,k))>c

c=a(q,k);

l=q;

end

end

%如果主元为0，则矩阵A不可逆

if abs(c)<1e-10

disp(‘error’);

pause;

exit

end

%如果l不等于k，则交换第l行和第k行if l~=k

for q=k:n+1

temp=a(k,q);

a(k,q)=a(l,q);

a(l,q)=temp;

end

end

%计算第k步的元素值

for i=k+1:n

for j=k+1:n

a(i,j)=a(i,j)-a(i,k)/a(k,k)*a(k,j);

end

end

end

%回代过程

x(n)=a(n,n+1)/a(n,n);