解线性方程组的直接解法

解线性方程组的直接解法

一、实验目的及要求

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

求解线性方程组的直接方法有以下几种： 1、利用左除运算符直接求解

线性方程组为Ax?b，则输入x?A\\b即可。 2、列主元的高斯消元法 程序流程图：

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

根据矩阵的秩判断是否有解，若无解停止；否则，顺序进行； 对于p?1:n?1

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

回代求解。（此部分可以参看课本第150页相关算法） 3、利用矩阵的分解求解线性方程组 （1）LU分解

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

[L,U]=lu（A）

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

角。

（2）平方根法

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

R=chol（A）

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

三、研究、解答以下问题

问题1、先将矩阵A进行楚列斯基分解，然后解方程组Ax?b(即利用平方根

法求解线性方程组，直接调用函数)：

1??12?32?6??????323?7?33????， A??b????2?799?6?16?????1?3?619??7?????解答： 程序：

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分解，然后解方程组Ax?b（直接调用函数）：

3?1?276?4?31?30?7?3A??560?13??2165?78?23765162???5?????89??，b???13???15???81???2??5??2?3? ?51?5??71?解答： 程序：

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、利用列主元的高斯消去法，求解下列方程组：??x1?20x2?x3?0.001x4?0??2x1?5x2?30x3?0.1x4?1?5x1?x2?100x 3?10x4?0??2x1?100x2?x3?x4?0解答：

程序：

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('??×￠òa￡oRA~=RB￡??ùò?′?·?3ì×é?T?a?￡') return end

if RA==RB if RA==n

disp('??×￠òa￡oòò?aRA=RB=n,?ùò?′?·?3ì×éóD?¨ò??a?￡') 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('??×￠òa￡oòò?aRA=RB?′n￡??ùò?′?·?3ìóD?T???à?a?￡') 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