论文全部材料 - 图文(5)

(k)(k)?第k步：在矩阵?A,b??的第k列中选主元，如

ai(kkk)，使

ai(kkk)?maxaik(k).将

k?i?n(k)(k)i??A,b??的第k行与第k行互换，进行第k次消元.

如此经过n?1步，增广矩阵（3.2）被化成上三角形，最后由回代过程求解.容易证明，

(k)只要det(A)?　0，列主元素法就可以顺利完成，即不会出现主元素akk?0的情形.

3.1.3 高斯消元法的Matlab实现 （1）高斯消去法的Matlab程序如下：

% gauss1.m

function x = gauss1(a, b,flag)

% 用途：顺序Gauss消去法解线性方程组ax=b

% 格式：x = gauss1(a, b,flag) a为系数矩阵，b为右端列向量，flag若为0，则显 % 示中间过程，否则不显示，默认值为0，x为解向量 if nargin<3,flag=0;end n = length(b); a=[a,b]; % 消元

for k = 1 : (n-1)

a((k+1):n,(k+1):(n+1))=a((k+1):n,(k+1):(n+1)-a((k+1):n,k)/a(k,k)*a(k,(k+1):(n+1));

a((k+1):n,k)= zeros(n-k,1)

if flag==0,a, end end % 回代 x=zero(n,1); x(n)=a(n,n+1)/a(n,n); for k=n-1:-1:1

x(k,:)=(a(k,n+1)-a(k,(k+1):n)*x((k+1):n))/a(k,k); end

（2）高斯列主元消去法的Matlab程序如下：

% gauss2.m

function x = gauss2(a, b,flag)

% 用途：列主元Gauss消去法解线性方程组ax=b

% 格式：x = gauss2(a, b,flag) a为系数矩阵，b为右端列向量，flag若为0，则显

7

% 示中间过程，否则不显示，默认值为0，x为解向量 if nargin<3,flag=0;end n = length(b); a=[a,b]; for k = 1 : (n-1) % 选主元

[ap,p]=max(abs(a(k:n,k)))；p=p+k-1; if p>k

t=a(k,:)；a(k，:)=a(p，:)；a(p，:)=t； t=b(k);b(k)=b(p);b(p)=t; end % 消元

a((k+1):n,(k+1):(n+1))=a((k+1):n,(k+1):(n+1)-a((k+1):n,k)/a(k,k)*a(k,(k+1):(n+1));

a((k+1):n,k)= zeros(n-k,1)

if flag==0,a, end end % 回代 x=zeros(n,1); x(n)=a(n,n+1)/a(n,n) for k=n-1:-1:1

x(k,:)=(a(k,n+1)-a(k,(k+1):n)*x((k+1):n))/a(k,k); end

3.2 矩阵三角分解法 3.2.1 直接三角分解法

如果方程组Ax?b的系数矩阵A能分解成

Ａ＝LU,

其中, L下三角矩阵, U是上三角矩阵,这时方程组就可化为两个容易求解的三角形方程组

Ly?b,Ux?y,

先由Ly?b解出向量y，再由Ux?y解出向量x，这就是原方程组Ax?b的解.若是单位下三角阵，则称相应的分解为 Doolittle分解(也称为LU分解).

定理 3.1 矩阵A??aii?n?n有唯一的Doolittle分解的充分必要条件是A的前n?1个

8

顺序主子式Dk?0(k?1,2,???n?1).

.

设矩阵A??aii?n?n非奇异，且A的前n?1个顺序主子式都不为零.则A有Dool-ittle分解

?　　0??1　 0　　?l　　?1　　?　　　021　?, A?LU????　　?　　?　　　????l　?　　l１?n,n?1　　??n1　?由矩阵的乘法可知

A?LU,

aij?u1j,j?1,2,?,n;　　　ai1?li1u11,i?2,?,n.

当k?2,3???,n时有

akj??lktutj?ukj,j?k,k?1,?,n,t?1nk?1aik??litutu??litutk?likukk,i?k?1,k?2,?,n,t?1t?1k?1

对于Doolittle分解算法,我们知道: (1) 对k?1

uij?aij,j?1,2,???,n, ai1lii?,i?2,???,n,u11(2) 对k?2,???,n

akj??lkrurj,j?k,k?1,?,n,r?1k?1lik?aik??litutkt?1k?1

ukk,i?k?1,k?2,?,n,其中可以利用LU分解求解线性方程组的算法, 则可以先求解Ly?b，即

?y1?b1?ly?y?b?2112, ??　　????ln1y1?ln2y2???ln,n?1yn?b9

所以

?y1?b1?k?1, ??yk?bk??lktytt?1?再求解Ux?y，即

?u11x1?u12x2???u1nxn?y1?ux???ux?y?2222nn2, ?????unnxn?yn回代求解

yn?x??nunn??n. ?y?ux?jttj?t?j?1?xj?,j?n?1,n?2,?,1u?jj?3.2.2 追赶法

设n元线性方程组Ax?b的系数矩阵为非奇异矩阵，A的三对角矩阵

?a1　 c1　　0　　0　　0??d　 a　c　　?0　　　02　2　　2　　??A??0　　?　　?　　?　　?　0?,

??0　　0　?　　?　　?　cn?1???0　　0　　0　ｄ　　a?ｎｎ??这种方程组称为三对角线性方程组.

这类方程组具有许多明显的应用背景，在求微分方程数值解、三次样条函数等问题中， 都会遇到这样的线性方程组.

设A的前n?1个顺序主子式都不为零，则A有唯一的LU分解，并且A的LU分解有如下形式：

??1　0　 0　 0　 　0　　??ｕ　ｖ　　０　　０　　０　　１１??l　1 0 0 　 0??０　ｕ　　ｖ　　０　　０２２??2??A?LU??0　　?　　?　0 0??０　０　　?　　　?　　　０?,

????0　　0　　　? ? 　0０　０　　０　　　　?　ｖn?1???????０　０　　０　　０　　　ｕ?0 0 0 ln　 0???ｎ?其中

10

vi?ci,ui?ai,　　i?1,2,?,n?1;li?di/ui?1,ui?ai?livi?1,　i?2,3,?,n;三对角方程组的追赶法

(1) 向前\追\的过程

.

ui?ai,vi?ci,yi?bi,

对i?2,3,???,n计算

li?di,ui?ai?livi?1,ui?1

vi?ci,yi?bi?liyi?1,(2) 往回\赶\的过程

xn?yn, un对i?n?1,n?2,?,1，计算

xi?yi?cixi?1. ui3.2.3 直接三角分解法的Matlab实现

（1）直接三角分解法的Matlab程序如下：

% nalu.m

function [1,u]= nalu(a)

% 用途：用LU分解法解方程组ax=b

% 格式：[1,u]= nalu(a) a为可逆方阵，l为单位下三角矩阵，u为上三角矩阵

n=length(a);

u=zeros(n,n);l=eye(n,n);

u(1,:)=a(1,:)；l(2:n,1)=a(2:n,1)/u(1,1); for k=2：n

u(k,k:n)=a(k,k:n)-l(k,1:k-1)*u(1:k-1,k:n);

l(k+1:n,k)=(a(k+1:n,k)-l(k+1:n,1:k-1)*u(1:k-1,k))/u(k,k); end

（2）追赶法的Matlab程序如下：

% nalu2.m

function [ x,L,U]=thomas(a,b,c,f) n=length(b);

% 对A进行分解

11