解线性方程组的直接解法

解线性方程组的直接解法

一、实验目的及要求

关于线性方程组的数值解法一般分为两大类：直接法与迭代法。直接法是在没有舍入误差的情况下，通过有限步运算来求方程组解的方法。通过本次试验的学习，应该掌握各种直接法，如：高斯列主元消去法，LU分解法和平方根法等算法的基本思想和原理，了解它们各自的优缺点及适用范围。

二、相关理论知识

求解线性方程组的直接方法有以下几种：

1、利用左除运算符直接求解

线性方程组为b

x\

=即可。

A

Ax=，则输入b

2、列主元的高斯消元法

程序流程图：

输入系数矩阵A，向量b，输出线性方程组的解x。

根据矩阵的秩判断是否有解，若无解停止；否则，顺序进行；

对于1

p

:1-

=n

选择第p列中最大元，并且交换行；

消元计算；

回代求解。（此部分可以参看课本第150页相关算法）

3、利用矩阵的分解求解线性方程组

（1）LU分解

调用matlab中的函数lu即可，调用格式如下：

[L,U]=lu（A）

注意：L往往不是一个下三角，但是可以经过行的变换化为单位下三角。

（2）平方根法

调用matlab 中的函数chol 即可，调用格式如下：

R=chol （A ）

输出的是一个上三角矩阵R ，使得R R A T =。

三、研究、解答以下问题

问题1、先将矩阵A 进行楚列斯基分解，然后解方程组b Ax =(即利用平方根法求解线性方程组，直接调用函数)：

⎪⎪⎪⎪⎪⎭⎫ ⎝⎛--------=19631699723723312312A ，⎪⎪⎪⎪⎪⎭

⎫ ⎝⎛-=71636b 解答：

程序：

A=[12 -3 2 1;-3 23 -7 -3;2 -7 99 -6;1 -3 -6 19];

R=chol(A)

b=[6 3 -16 7]';

y=inv(R')*b %y=R'\b

x=inv(R)*y %x=R\y

结果：

R =3.4641 -0.8660 0.5774 0.2887

0 4.7170 -1.3780 -0.5830

0 0 9.8371 -0.7085

0 0 0 4.2514

y =1.7321

0.9540

-1.5945

1.3940

x =0.5463

0.2023

-0.1385

0.3279

问题 2、先将矩阵A 进行LU 分解，然后解方程组b Ax =（直接调用函数）：

⎪⎪⎪⎪⎪⎪⎪⎪⎭⎫ ⎝⎛----=8162517623158765211331056897031354376231A ，⎪⎪⎪⎪⎪⎪⎪⎪⎭

⎫ ⎝⎛-=715513252b

解答：

程序：

A=[1/3 -2 76 3/4 5;3 1/sqrt(3) 0 -7 89;56 0 -1 3 13;21 65 -7 8 15;23 76 51 62 81];

b=[2/sqrt(5);-2;3;51;5/sqrt(71)];

[L,U]=lu(A)

y=inv(L)*b

x=inv(U)*y

结果：

L = 0.0060 -0.0263 1.0000 0 0

0.0536 0.0076 -0.0044 0.1747 1.0000

1.0000 0 0 0 0

0.3750 0.8553 -0.6540 1.0000 0

0.4107 1.0000 0 0 0

U =56.0000 0 -1.0000 3.0000 13.0000

0 76.0000 51.4107 60.7679 75.6607

0 0 77.3589 2.3313 6.9137

0 0 0 -43.5728 -50.0631

0 0 0 0 96.5050

y =3.0000

-0.6388

0.8598

50.9836

-11.0590

x =0.1367

0.9004

0.0526

-1.0384

-0.1146

问题3、利用列主元的高斯消去法，求解下列方程组：

⎪⎪⎩⎪⎪⎨⎧=+--=--+=-+-=+-+0

10020

1010051

1.030520

001.0204321432143214321x x x x x x x x x x x x x x x x

解答：

程序：

function [RA,RB,n,X]=liezhu(A,b)

B=[A b];n=length(b);RA=rank(A);

RB=rank(B);zhica=RB-RA;

if zhica>0

disp('Çë×¢Ò⣺RA~=RB£¬ËùÒÔ´Ë·½³Ì×éÎ޽⡣')

return

end

if RA==RB

if RA==n

disp('Çë×¢Ò⣺ÒòΪRA=RB=n,ËùÒÔ´Ë·½³Ì×éÓÐΨһ½â¡£')

X=zeros(n,1);C=zeros(1,n+1);

for p=1:n-1

[Y ,j]=max(abs(B(p:n,p)));C=B(p,:);

for k=p+1:n

m=B(k,p)/B(p,p);

B(k,p:n+1)=B(k,p:n+1)-m*B(p,p:n+1)

end

end

b=B(1:n,n+1);A=B(1:n,1:n);X(n)=b(n)/A(n,n);

for q=n-1:-1:1

X(q)=(b(q)-sum(A(q,q+1:n)*X(q+1:n)))/A(q,q);

end

else

disp('Çë×¢Ò⣺ÒòΪRA=RB¡´n£¬ËùÒÔ´Ë·½³ÌÓÐÎÞÇî¶à½â¡£') end

end

键入

A=[1 20 -1 0.001

2 -5 30 -0.1

5 1 -100 -10

2 -100 -1 1];

b=[0;1;0;0];

[RA,RB,n,X]=liezhu(A,b)

结果：

请注意：因为RA=RB=n,所以此方程组有唯一解。

B = 1.0000 20.0000 -1.0000 0.0010 0

0 -45.0000 32.0000 -0.1020 1.0000

5.0000 1.0000 -100.0000 -10.0000 0

2.0000 -100.0000 -1.0000 1.0000 0

B = 1.0000 20.0000 -1.0000 0.0010 0

0 -45.0000 32.0000 -0.1020 1.0000

0 -99.0000 -95.0000 -10.0050 0

2.0000 -100.0000 -1.0000 1.0000 0

B = 1.0000 20.0000 -1.0000 0.0010 0

0 -45.0000 32.0000 -0.1020 1.0000

0 -99.0000 -95.0000 -10.0050 0

0 -140.0000 1.0000 0.9980 0

B = 1.0000 20.0000 -1.0000 0.0010 0

0 -45.0000 32.0000 -0.1020 1.0000

0 0.0000 -165.4000 -9.7806 -2.2000

0 -140.0000 1.0000 0.9980 0