高斯消元法和三角分解法

实验四线性方程组的直接数值解法

信息与计算科学金融崔振威201002034031

一、实验目的：

编程实现高斯消元法和三角分解法

二、实验内容：p108.15、p120.1

三、实验要求：

1、分别对所给算例利用高斯消元法和三角分解法进行数值求解

2、分析系数矩阵p108.15的算法稳定性

主程序

高斯消元法：

function [x,det,flag]=Gauss(A,b)

%求线性方程组的列主元Gauss消去法

%A为方程组的系数矩阵

%b为方程组的右端项

%x为方程组的解

%det为系数矩阵A的行列式的值

%flag为指标向量，flag为failure表示计算失败，为ok时表示计算成功[n,m]=size(A);nb=length(b);

%当方程组与列的维数不相等时停止计算，并输出出错信息

if n~=m

error('方程组与右端项列的维数不相等，错误。');

return;

end

%赋初始值，计算

flag='OK';det=1;x=zeros(n,1);

for k=1:n-1

%选主元

max1=0;

for i=k:n

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

max1=abs(A(i,k));r=i;

end

end

if max1<1e-10

flag='failure';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;det=-det;

end

%消元过程

for i=k+1:n

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

for j=k+1:n

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

end

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

end

det=det*A(k,k);

end

det=det*A(n,n);

%回代过程

if abs(A(n,n))<1e-10

flag='failure';return;

end

for k=n:-1:1

for j=k+1:n

b(k)=b(k)-A(k,j)*x(j);

end

x(k)=b(k)/A(k,k);

end

x(k)=b(k)/A(k,k);

end

三角分解法：

function [L,U,flag]=LU_Decom(A)

%求矩阵A的LU分解

%A为被分解的矩阵

%L为单位下三角阵

%U为单位上三角阵

%flag为指标向量，flag为failure表示计算失败，为ok时表示计算成功[n,m]=size(A);

%要求所分解的矩阵是方阵，否则停止计算，并输出错误信息

if n~=m

error('所分解的矩阵不是一个方阵');

return;

end

L=eye(n);U=zeros(n);flag='OK';

for k=1:n

for j=k:n

z=0;

for q=1:k-1

z=z+L(k,q)*U(q,j);

end

U(k,j)=A(k,j)-z;

end

if abs(U(k,k))

flag='failure';return;

end

for i=k+1:n

z=0;

for q=1:k-1

z=z+L(i,q)*U(q,k);

end

L(i,k)=(A(i,k)-z)/U(k,k);

end

end

P108 15

(a)

解：1、利用高斯消元，在matlab的命令窗口中输入命令，并得出结果：

>> A=[1 1/2 1/3 1/4;1/2 1/3 1/4 1/5;1/3 1/4 1/5 1/6;1/4 1/5 1/6 1/7],b=[1 0 0

0],[x,det,flag]=Gauss(hilb(4),b)

A =

1.00000000000000 0.50000000000000 0.33333333333333 0.25000000000000

0.50000000000000 0.33333333333333 0.25000000000000 0.20000000000000

0.33333333333333 0.25000000000000 0.20000000000000 0.16666666666667

0.25000000000000 0.20000000000000 0.16666666666667 0.14285714285714

b =

1 0 0 0

x =

1.0e+002 *

0.15999999999999

-1.19999999999992

2.39999999999980

-1.39999999999987