列主元高斯消去法和列主元三角分解法解线性方程.docx

列主元高斯消去法和列主元三角分解法解线性方程.docx

《列主元高斯消去法和列主元三角分解法解线性方程.docx》由会员分享，可在线阅读，更多相关《列主元高斯消去法和列主元三角分解法解线性方程.docx（15页珍藏版）》请在冰豆网上搜索。

列主元高斯消去法和列主元三角分解法解线性方程

计算方法实验报告1

【课题名称】

用列主元高斯消去法和列主元三角分解法解线性方程

【目的和意义】

高斯消去法是一个古老的求解线性方程组的方法，但由它改进得到的选主元的高斯消去法则是目前计算机上常用的解低阶稠密矩阵方程组的有效方法。

用高斯消去法解线性方程组的基本思想时用矩阵行的初等变换将系数矩阵A约化为具有简单形式的矩阵（上三角矩阵、单位矩阵等），而三角形方程组则可以直接回带求解

用高斯消去法解线性方程组

（其中A∈Rn×n）的计算量为：

乘除法运算步骤为

，加减运算步骤为

。

相比之下，传统的克莱姆法则则较为繁琐，如求解20阶线性方程组，克莱姆法则大约要

次乘法，而用高斯消去法只需要3060次乘除法。

在高斯消去法运算的过程中，如果出现abs（A（i,i））等于零或过小的情况，则会导致矩阵元素数量级严重增长和舍入误差的扩散，使得最后的计算结果不可靠，所以目前计算机上常用的解低阶稠密矩阵方程的快速有效的方法时列主元高斯消去法，从而使计算结果更加精确。

2、列主元三角分解法

高斯消去法的消去过程，实质上是将A分解为两个三角矩阵的乘积A=LU，并求解Ly=b的过程。

回带过程就是求解上三角方程组Ux=y。

所以在实际的运算中，矩阵L和U可以直接计算出，而不需要任何中间步骤，从而在计算过程中将高斯消去法的步骤进行了进一步的简略，大大提高了运算速度，这就是三角分解法

采用选主元的方式与列主元高斯消去法一样，也是为了避免除数过小，从而保证了计算的精确度

【计算公式】

1、列主元高斯消去法

设有线性方程组Ax=b，其中设A为非奇异矩阵。

方程组的增广矩阵为

第1步（k=1）：

首先在A的第一列中选取绝对值最大的元素

，作为第一步的主元素:

然后交换（A，b）的第1行与第l行元素，再进行消元计算。

设列主元素消去法已经完成第1步到第k-1步的按列选主元，交换两行，消元计算得到与原方程组等价的方程组A（k）x=b（k）

第k步计算如下：

对于k=1，2，…，n-1

（1）按列选主元：

即确定t使

（2）如果t≠k，则交换[A，b]第t行与第k行元素。

（3）消元计算

消元乘数mik满足:

（4）回代求解

2、

列主元三角分解法

对方程组的增广矩阵经过k-1步分解后，可变成如下形式：

第k步分解，为了避免用绝对值很小的数

作除数，引进量

，于是有

=

。

如果，则将矩阵的第t行与第k行元素互换，将（i，j）位置的新元素仍记为

或

，然后再做第k步分解，这时

【列主元高斯消去法程序流程图】

【列主元高斯消去法Matlab主程序】

functionx=gauss1（A,b,c）%列主元法高斯消去法解线性方程Ax=b

if（length（A）~=length（b））%判断输入的方程组是否有误

disp（'输入方程有误！

'）

return;

end

disp（'原方程为AX=b：

'）%显示方程组

A

b

disp（'------------------------'）

n=length（A）;

fork=1:

n-1%找列主元

[p,q]=max（abs（A（k:

n,k）））;%找出第k列中的最大值，其下标为[p,q]

q=q+k-1;%q在A（k:

n,k）中的行号转换为在A中的行号

ifabs（p）

disp（'列元素太小，det（A）≈0'）;

break;

elseifq>k

temp1=A（k,:

）;%列主元所在行不是当前行，将当前行与列主

A（k,:

）=A（q,:

）;元所在行交换（包括b）

A（q,:

）=temp1;

temp2=b（k,:

）;

b（k,:

）=b（q,:

）;

b（q,:

）=temp2;

end

%消元

fori=k+1:

n

m（i,k）=A（i,k）/A（k,k）;%A（k,k）将A（i,k）消为0所乘系数

A（i,k:

n）=A（i,k:

n）-m（i,k）*A（k,k:

n）;%第i行消元处理

b（i）=b（i）-m（i,k）*b（k）;%b消元处理

end

end

disp（'消元后所得到的上三角阵是'）

A%显示消元后的系数矩阵

b（n）=b（n）/A（n,n）;%回代求解

fori=n-1:

-1:

1

b（i）=（b（i）-sum（A（i,i+1:

n）*b（i+1:

n）））/A（i,i）;

end

clearx;

disp（'AX=b的解x是'）

x=b;

【调用函数解题】

【列主元三角分解法程序流程图】

【列主元三角分解法Matlab主程序】

①自己编的程序：

functionx=PLU（A,b,eps）%定义函数列主元三角分解法函数

if（length（A）~=length（b））%判断输入的方程组是否有误

disp（'输入方程有误！

'）

return;

end

disp（'原方程为AX=b：

'）%显示方程组

A

b

disp（'------------------------'）

n=length（A）;

A=[Ab];%将A与b合并，得到增广矩阵

forr=1:

n

ifr==1

fori=1:

n

[cd]=max（abs（A（:

1）））;%选取最大列向量，并做行交换

ifc<=eps%最大值小于e，主元太小，程序结束

break;

else

end

d=d+1-1;

p=A（1,:

）;

A（1,:

）=A（d,:

）;

A（d,:

）=p;

A（1,i）=A（1,i）;

end

A（1,2:

n）=A（1,2:

n）;

A（2:

n,1）=A（2:

n,1）/A（1,1）;%求u（1,i）

else

u（r,r）=A（r,r）-A（r,1:

r-1）*A（1:

r-1,r）;%按照方程求取u（r,i）

ifabs（u（r,r））<=eps%如果u（r,r）小于e，则交换行

p=A（r,:

）;

A（r,:

）=A（r+1,:

）;

A（r+1,:

）=p;

else

end

fori=r:

n

A（r,i）=A（r,i）-A（r,1:

r-1）*A（1:

r-1,i）;%根据公式求解，并把结果存在矩阵A中

end

fori=r+1:

n

A（i,r）=（A（i,r）-A（i,1:

r-1）*A（1:

r-1,r））/A（r,r）;%根据公式求解，并把结果存在矩阵A中

end

end

end

y

（1）=A（1,n+1）;

fori=2:

n

h=0;

fork=1:

i-1

h=h+A（i,k）*y（k）;

end

y（i）=A（i,n+1）-h;%根据公式求解y（i）

end

x（n）=y（n）/A（n,n）;

fori=n-1:

-1:

1

h=0;

fork=i+1:

n

h=h+A（i,k）*x（k）;

end

x（i）=（y（i）-h）/A（i,i）;%根据公式求解x（i）

end

A

disp（'AX=b的解x是'）

x=x';%输出方程的解

②可直接得到P,L,U并解出方程解的的程序（查阅资料得子函数PLU1，其作用是将矩阵A分解成L乘以U的形式。

PLU2为调用PLU1解题的程序，是自己编的）

（Ⅰ）.

function[l,u,p]=PLU1（A）%定义子函数，其功能为列主元三角分解系数矩阵A

[m,n]=size（A）;%判断系数矩阵是否为方阵

ifm~=n

error（'矩阵不是方阵'）

return

end

ifdet（A）==0%判断系数矩阵能否被三角分解

error（'矩阵不能被三角分解'）

end

u=A;p=eye（m）;l=eye（m）;%将系数矩阵三角分解，分别求出P,L,U

fori=1:

m

forj=i:

m

t（j）=u（j,i）;

fork=1:

i-1

t（j）=t（j）-u（j,k）*u（k,i）;

end

end

a=i;b=abs（t（i））;

forj=i+1:

m

ifb

b=abs（t（j））;

a=j;

end

end

ifa~=i

forj=1:

m

c=u（i,j）;

u（i,j）=u（a,j）;

u（a,j）=c;

end

forj=1:

m

c=p（i,j）;

p（i,j）=p（a,j）;

p（a,j）=c;

end

c=t（a）;

t（a）=t（i）;

t（i）=c;

end

u（i,i）=t（i）;

forj=i+1:

m

u（j,i）=t（j）/t（i）;

end

forj=i+1:

m

fork=1:

i-1

u（i,j）=u（i,j）-u（i,k）*u（k,j）;

end

end

end

l=tril（u,-1）+eye（m）;

u=triu（u,0）

（Ⅱ）.

functionx=PLU2（A,b）%定义列主元三角分解法的函数

[l,u,p]=PLU1（A）%调用PLU分解系数矩阵A

m=length（A）;%由于A左乘p，故b也要左乘p

v=b;

forq=1:

m

b（q）=sum（p（q,1:

m）*v（1:

m,1））;

end

b

（1）=b

（1）%求解方程Ly=b

fori=2:

1:

m

b（i）=（b（i）-sum（l（i,1:

i-1）*b（1:

i-1）））;

end

b（m）=b（m）/u（m,m）;%求解方程Ux=y

fori=m-1:

-1:

1

b（i）=（b（i）-sum（u（i,i+1:

m）*b（i+1:

m）））/u（i,i）;

end

clearx;

disp（'AX=b的解x是'）

x=b;

【调用函数解题】

①

②

【编程疑难】

这是第一次用matlab编程，对matlab的语句还不是非常熟悉，因此在编程过程中，出现了许多错误提示。

并且此次编程的两种方法对矩阵的运算也比较复杂。

问题主要集中在循环控制中，循环次数多了一次或者缺少了一次，导致数据错误，一些基本的编程语句在语法上也会由于生疏而产生许多问题，但是语句的错误由于系统会提示，比较容易进行修改，数据计算过程中的一些逻辑错误，比如循环变量的控制，这些系统不会提示错误，需要我们细心去发现错误，不断修正，调试。