数值线性代数第二版徐树方高立张平文上机习题第一章实验报告(供参考)

上机习题

1.先用你所熟悉的的计算机语言将不选主元和列主元Gauss 消去法编写成通用的子程序；然后用你编写的程序求解84阶方程组；最后将你的计算结果与方程的精确解进行比较，并就此谈谈你对Gauss 消去法的看法。

Sol ：

（1）先用matlab 将不选主元和列主元Gauss 消去法编写成通用的子程序，得到P U L ,,： 不选主元Gauss 消去法：[])(,A GaussLA U L =得到U L ,满足LU A =

列主元Gauss 消去法：[])(,,A GaussCol P U L =得到P U L ,,满足LU PA =

（2）用前代法解()Pb or b Ly =，得y

用回代法解y Ux =，得x

求解程序为()P U L b A Gauss x ,,,,=（P 可缺省，缺省时默认为单位矩阵）

（3）计算脚本为ex1_1

代码

%算法

function [L,U]=GaussLA(A)

n=length(A);

for k=1:n-1

A(k+1:n,k)=A(k+1:n,k)/A(k,k);

A(k+1:n,k+1:n)=A(k+1:n,k+1:n)-A(k+1:n,k)*A(k,k+1:n);

end

U=triu(A);

L=tril(A);

L=L-diag(diag(L))+diag(ones(1,n));

%算法

function [L,U,P]=GaussCol(A)

n=length(A);

for k=1:n-1

[s,t]=max(abs(A(k:n,k)));

p=t+k-1;

temp=A(k,1:n);

A(k,1:n)=A(p,1:n);

A(p,1:n)=temp;

u(k)=p;

if A(k,k)~=0

A(k+1:n,k)=A(k+1:n,k)/A(k,k);

A(k+1:n,k+1:n)=A(k+1:n,k+1:n)-A(k+1:n,k)*A(k,k+1:n); else

break;

end

end

L=tril(A);U=triu(A);L=L-diag(diag(L))+diag(ones(1,n));

P=eye(n);

for i=1:n-1

temp=P(i,:);

P(i,:)=P(u(i),:);

P(u(i),:)=temp;

end

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

function x=Gauss(A,b,L,U,P)

if nargin<5

P=eye(length(A));

end

n=length(A);

b=P*b;

for j=1:n-1

b(j)=b(j)/L(j,j);

b(j+1:n)=b(j+1:n)-b(j)*L(j+1:n,j);

end

b(n)=b(n)/L(n,n);

y=b;

for j=n:-1:2

y(j)=y(j)/U(j,j);

y(1:j-1)=y(1:j-1)-y(j)*U(1:j-1,j);

end

y(1)=y(1)/U(1,1);

x=y;

end

ex1_1

clc;clear;

%第一题

A=6*eye(84)+diag(8*ones(1,83),-1)+diag(ones(1,83),1);

b=[7;15*ones(82,1);14];

%不选主元Gauss 消去法

[L,U]=GaussLA(A);

x1_1=Gauss(A,b,L,U);

%列主元Gauss 消去法

[L,U,P]=GaussCol(A);

x1_2=Gauss(A,b,L,U,P);

%解的比较

subplot(1,3,1);plot(1:84,x1_1,'o-');title('Gauss');

subplot(1,3,2);plot(1:84,x1_2,'.-');title('PGauss');

subplot(1,3,3);plot(1:84,ones(1,84),'*-');title('精确解');

结果为（其中Gauss 表示不选主元的Gauss 消去法，PGauss 表示列主元Gauss

消去法，精确解为[]'⨯8411,,1 ）：

由图，显然列主元消去法与精确解更为接近，不选主元的Gauss 消去法误差比列主元消去法大，且不如列主元消去法稳定。

Gauss 消去法重点在于A 的分解过程，无论A 如何分解，后面两步的运算过程不变。

2.先用你所熟悉的的计算机语言将平方根法和改进的平方根法编写成通用的子程序；然后用你编写的程序求解对称正定方程组Ax=b 。

Sol ：

（1）先用matlab 将平方根法和改进的平方根法编写成通用的子程序，得到L ，(D)： 平方根法：L=Cholesky(A)

改进的平方根法：[L,D]=LDLt(A)

（2）求解得b Ly =

求解得y x DL or y x L T

T ==

求解程序为x=Gauss(A,b,L,U,P)（T

T DL U or L U == ，P 此时缺省，缺省时默认为单位矩阵）

（3）计算脚本为ex1_2

代码

%算法

function L=Cholesky(A)

n=length(A);

for k=1:n

A(k,k)=sqrt(A(k,k));

A(k+1:n,k)=A(k+1:n,k)/A(k,k);

for j=k+1:n

A(j:n,j)=A(j:n,j)-A(j:n,k)*A(j,k);

end

end

L=tril(A);

end

%计算LDL ‘分解：改进的平方根法

function [L,D]=LDLt(A)

n=length(A);

for j=1:n

for i=1:n

v(i,1)=A(j,i)*A(i,i);

end

A(j,j)=A(j,j)-A(j,1:j-1)*v(1:j-1,1);

A(j+1:n,j)=(A(j+1:n,j)-A(j+1:n,1:j-1)*v(1:j-1,1))/A(j,j);