计算方法作业第六章

1．考虑两个线性方程组，其系数矩阵如下

1211

11...23211111...1212341,121

11

11...3452........................

...

(12111)

1 (1)

2

21n n A A n n n n n ⎡

⎤⎢⎥⎢

⎥-⎡⎤⎢⎥⎢⎥--⎢⎥

+⎢⎥

⎢⎥⎢⎥==--⎢⎥⎢⎥⎢⎥+⎢⎥⎢⎥

⎢⎥-⎣⎦

⎢

⎥⎢⎥⎢⎥++-⎣⎦

问题的真解均取为[1,1,1,1,...1]T

x =，线性方程组的右端项用这个真解计算出来。相应的问题分别称为问题I 和问题II 。请进行如下数值实验：

（1） 对问题I 分别用Gauss 消元法，Cholesky 方法，修改的LDLT 算法，追赶法四

种方法求解，其中n=100；

（2） 对问题II 分别用Gauss 消去法，列主元Gauss 消去法，不做行交换的列主元

Gauss 消去法求解，其中n=6；

（3） 不断增加问题II 的矩阵阶数n=6,8,10,…,20，重复（2）的工作，看看会有什么

问题发生？解释其原因。

（1） Gauss ： 计算程序： n=100; A=2*eye(n); for i=1:n-1 A(i+1,i)=-1; A(i,i+1)=-1; end b=0; b(1)=1; b(100)=1;

[x,XA]=GaussJordanXQ(A,b); Gauss 消元法源程序：

%用Gauss 消元法解线性方程组 function [x,XA]=GaussJordanXQ(A,b) N = size(A); n = N(1); for i=1:(n-1)

for j=(i+1):n

if(A(i,i)==0)

disp('对角元素为0！');

%防止对角元素为0

return;

end

l = A(j,i);

m = A(i,i);

A(j,1:n)=A(j,1:n)-l*A(i,1:n)/m;

%消元方程

b(j)=b(j)-l*b(i)/m;

end

end

x=SolveUpTriangle(A,b);

%通用的求上三角系数矩阵线性方程组的函数XA = A;

%消元后的系数矩阵

（SolveUpTriangle.m）解上三角方程组源程序：%解上三角方程组

function x=SolveUpTriangle(A,b)

N = size(A);

n = N(1);

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

for k=n-1:1

s=0;

for i=k+1:n

s=s+A(k,i)*x(i);

end

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

end

结果：

x1=[0,0,0,…..0,0,1]T

x2=[0,0,0,…..0,0,0.3820]T

x3=[0,0,0,…..0,0,0.9900]T

x4=[1,1,1,…..1,1,1]T

Cholesky：

计算程序：

n=100;

A=2*eye(n);

for i=1:n-1

A(i+1,i)=-1;

A(i,i+1)=-1;

end

b=0;

b(1)=1;

b(100)=1;

x2=Cholesky(A,b); (Cholesky.m):Cholesky法源程序%用Cholesky方法解线性方程组function x=Cholesky(A,b);

N=size(A);

n=N(1);

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

%L为下三角矩阵

for j=2:n

L(j,1)=A(1,j)./L(1,1)

end

for k=2:n

s=0;

for i=1:k-1

s=s+L(k,i);

end

L(k,k)=sqrt(A(k,k)-s);

for j=(k+1):n

sum=0;

for i=1:k-1

sum=sum+L(j,i)*L(k,i);

end

L(j,k)=(A(j,k)-sum)./L(k,k); end

end

%定义L的转置

for j=1:n

for k=1:n

LT(j,k)=L(k,j);

end

end

y=SolveDownTriangle(L,b);

x=SolveUpTriangle(LT,y);

解上三角方程组源程序：

%解上三角方程组

function x=SolveUpTriangle(A,b)

N = size(A);

n = N(1);

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

for k=n-1:1

s=0;

for i=k+1:n

s=s+A(k,i)*x(i);

end

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

End

解下三角方程组源程序：function x=SolveDownTriangle(A,b) N = size(A);

n = N(1);

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

for k=2:n

s=0;

for i=1:k-1

s=s+A(k,i)*x(i);

end

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

end

结果：

x1=[0,0,0,…..0,0,0.9894]T;

x2=[0,0,0,…..0,0,0.9899]T;

LTDT：

计算程序：