南京邮电大学 数值代数实验

数值代数实验

数值线性代数实验一

一、实验名称：矩阵的LU分解.

二、实验目的：用不选主元的LU分解和列主元LU分解求解线性方程组Ax=b, 并比较这

两种方法.

三、实验内容与要求

（1）用所熟悉的计算机语言将不选主元和列主元LU分解编成通用的子程序，然后用编写的程序求解下面的84阶方程组

将计算结果与方程组的精确解进行比较，并就此谈谈你对Gauss消去法的看法.

（2）写出追赶法求解三对角方程组的过程，并编写程序求该实验中的方程组

Gauss消去法：

用消去法解方程组的基本思想是用逐次消去未知数的方法把原来方程组Ax=b化为与其等价的三角方程组，而求解三角方程组就容易了。换句话说，上述过程就是用行的初等变换将原方程组系数矩阵化为简单形式，从而将求解原方程组的问题转化为求解简单方程组的问题。

利用Gauss消去法对线性方程组Ax=b进行求解。

用MATLAB建立m文件DelGauss.m,程序如下：

function x=DelGauss(a,b)

[n,m]=size(a);

nb=length(b);

det=1;

x=zeros(n,1);

for k=1:n-1

for i=k+1:n

if a(k,k)==0

return

end

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);

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

在matlab中输入如下：

结果如下：

方程组的精确解为x1=x2=…=x84=1.0000，与Gauss消去法求得的解差距很大，所得结果不够准确，计算简单但其消元过程有时不能进行到底而使求解出现解失真的情况。

数值线性代数实验二

一、实验名称：实对称正定矩阵的A的Cholesky分解.

二、实验目的：用平方根法和改进的平方根方法求解线性方程组Ax=b.

三、实验内容与要求

用所熟悉的计算机语言将Cholesky分解和改进的Cholesky分解编成通用的子程序，然后用编写的程序求解对称正定方程组Ax=b，其中

（1）b随机的选取，系数矩阵为100阶矩阵

（2）系数矩阵为40阶Hilbert矩阵，即系数矩阵A的第i行第j列元素为

，向量b的第i个分量为

（3）用实验一的程序求解这两个方程组，并比较所有的计算结果，然后评价各个方法的优劣。

平方根法：

平方根法就是利用对称正定矩阵的三角分解而得到的求解对称正定方程组的一种有效方法。平方根法递推公式可以证明对于对称正定矩阵A,可以唯一地分解成A=LL T,其中L是非奇异下三角形矩阵。

模型二：利用平方根法对线性方程组Ax=b进行求解。

用MATLAB建立m文件pingfg.m,程序如下：

function [x]=pingfg(A,b) %Cholesky分解

[n,n]=size(A);

L=zeros(n,n);%实际上不用为 L 申请空间，使用 A 即可

L(1,1)=sqrt(A(1,1));

for k=2:n

L(k,1)=A(k,1)/L(1,1);

end

for k=2:n-1

L(k,k)=sqrt(A(k,k)-sum(L(k,1:k-1).^2));

for i=k+1:n

L(i,k)=(A(i,k)-sum(L(i,1:k-1).*L(k,1:k-1)))/L(k,k);

end

end

L(n,n)=sqrt(A(n,n)-sum(L(n,1:n-1).^2));%解下三角方程组Ly=b

y=zeros(n,1);

for k=1:n

j=1:k-1;

y(k)=(b(k)-L(k,j)*y(j))/L(k,k);

end%解上三角方程组L'x=y

x=zeros(n,1);

U=L';

for k=n:-1:1

j=k+1:n;

x(k)=(y(k)-U(k,j)*x(j))/U(k,k);

End

模型三：利用改进的平方根法对线性方程组Ax=b进行求解。

用MATLAB建立m文件ave.m,程序如下：

function [x]=ave(A,b,n) %用改进平方根法求解Ax=b

L=zeros(n,n); %L为n*n矩阵

D=diag(n,0); %D为n*n的主对角矩阵

S=L*D;

for i=1:n %L的主对角元素均为1

L(i,i)=1;

end

for i=1:n

for j=1:n %验证A是否为对称正定矩阵

if (eig(A)<=0)|(A(i,j)~=A(j,i)) %A的特征值小于0或A非对称时，输出wrong disp('wrong');break;end

end

end

Hilbert矩阵用MATLAB建立m文件Hil.m,程序如下：

function b=Hil()

for k=1:40

for m=1:40

s=0;

t=s+1/(k+m-1);

s=t;

end

b(k,1)=s;

end

在matlab中输入如下：

输出结果如下：