现代数值计算方法实验二.doc

实验二 线性方程组直接实验
一、实验目的
1.运用matlab软件完成线性方程组的直接实验;
2.通过实验,了解列主元和顺序高斯消去法的异同,LU分解的优点及追赶法的优点.
二、实验题目
1.分别用列主元法和顺序高斯消去法求解下面的线性方程组,分析对结果的影响:
.
2.LU分解的优点
实验题目：给矩阵A和向量b:

⑴求A的LU分解,n的值自己确定;
⑵利用A的LU分解求解下列方程组(a)Ax=b, (b)A2x=b, (c)A3x=b.
对方程组(c),若先求LU=A3,再解(LU)x=b有何缺点？
3.追赶法的优点
实验题目：用追赶法分别对n=10,n=100,n=1000解方程组Ax=b,其中

再用LU分解法解此方程组,并对二者进行比较.
三、实验原理
1.⑴列主元法算法原理：
%magauss2.m
function x=magauss2(A,b,flag)
if nargin<3,flag=0;end
n=length(b);
for k=1:(n-1)
%选主元
[ap,p]=max(abs(A(k:n,k))); p=p+k-1;
if p>k
t=A(k,:); A(k,:)=A(p,:); A(p,:)=t;
t=b(k); b(k)=b(p); b(p)=t;
end
%消元
m=A(k+1:n,k)/A(k,k);
A(k+1:n,k+1:n)=A(k+1:n,k+1:n)-m*A(k,k+1:n);
b(k+1:n)=b(k+1:n)-m*b(k);
A(k+1:n,k)=zeros(n-k,1);
if flag~=0, Ab=[A,b], end
end
%回代
x=zeros(n,1);
x(n)=b(n)/A(n,n);
for k=n-1:-1:1
x(k)=(b(k)-A(k,k+1:n)*x(k+1:n))/A(k,k);
end

⑵顺序高斯消去法算法原理：
%magauss.m
function x=magauss(A,b,flag)
if nargin<3,flag=0;end
n=length(b);
%消元
for k=1:(n-1)
m=A(k+1:n,k)/A(k,k);
A(k+1:n,k+1:n)=A(k+1:n,k+1:n)-m*A(k,k+1:n);
b(k+1:n)=b(k+1:n)-m*b(k);
A(k+1:n,k)=zeros(n-k,1);
if flag~=0, Ab=[A,b], end
end
%回代
x=zeros(n,1);
x(n)=b(n)/A(n,n);
for k=n-1:-1:1
x(k)=(b(k)-A(k,k+1:n)*x(k+1:n))/A(k,k);
end

2.LU分解算法原理：
%LU分解
function [l,u]=lufj(A)
[n,m]=size(A);
u=zeros(n,n);l=eye(n,n);
u(1,:)=A(1,:); l(2:n,1)=A(2:n,1)/u(1,1);
for k=2:n
u(k,k:n)=A(k,k:n)-l(k,1:k-1)*u(1:k-1,k:n);
l(k+1:n,k)=(A(k+1:n,k)-l(k+1:n,1:k-1)*u(1:k-1,k))/u(k,k);
end

%解下三角方程组Ly=b
function y=lowt(l,b)
n=length(b);
y=zeros(n,1);
y(1)=b(1);
for k=2:n
y(k)=b(k)-l(k,1:k-1)*y(1:k-1);
end

%解上三角方程组Ux=y
function x=upt(u,y)
n=length(y);
x=zeros(n,1);
x(n)=y(n)/u(n,n);
for k=n-1:-1:1
x(k)=(y(k)-u(k,k+1:n)*x(k+1:n))/u(k,k);
end

malu.m
function [x,l,u]=malu(A,b)
%用途：用LU分解法解方程组Ax=b
%格式：[x,l,u]=malu(A,b) A为系数矩阵，b为右端向量，x返回
% 解向量，l返回下三角矩阵，u返回上三角矩阵
format short
%LU分解
n=length(b);
u=zeros(n,n);l=eye(n,n);
u(1,:)=A(1,:); l(2:n,1)=A(2:n,1)/u(1,1);
for k=2:n
u(k,k:n)=A(k,k:n)-l(k,1:k-1)*u(1:k-1,k:n);
l(k+1:n,k)=(A(k+1:n,k)-l(k+1:n,1:k-1)*u(1:k-1,k))/u(k,k);
end
%解下三角方程组Ly=b
y=zeros(n,1);
y(1)=b(1);
for k=2:n
y(k)=b(k)-l(k,1:k-1)*y(1:k-1);
end
%解上三角方程组Ux=y
x=zeros(n,1);
x(n)=y(n)/u(n,n);
for k=n-1:-1:1
x(k)=(y(k)-u(k,k+1:n)*x(k+1:n))/u(k,k);
end

3. 追赶法的算法原理
%machase.m
function x=machase(a,b,c,d)
%用途：追赶法解三对角方程组Ax=d
%格式：x= machase(a,b,c,d) a为次下对角线元素向量，b主对角元素
% 向量，c为次上对角线元素向量，d为右端向量，x返回解向量
n=length(a);
for k=2:n
b(k)=b(k)-a(k)/b(k-1)*c(k-1);
d(k)=d(k)-a(k)/b(k-1)*d(k-1);
end
x(n)=d(n)/b(n);
for k=n-1:-1:1
x(k)=(d(k)-c(k)*x(k+1))/b(k);
end
四、实验内容
1. 程序：
A=[0.3e-16 59.14 3 1;1 2 1 1;11.2 9 5 2;5.291 -6.13 -1 2]；
b=[51.97;2;1;46.78]；
x1=magauss2(A,b)
x2=magauss(A,b)

2.程序:
n=6;
A=zeros(n,n);
b=zeros(n,1); b(n)=1;
for(i=1:n)
for(j=i:n)
A(i,j)=n-j+i;
end
end
%LU分解
[l,u]=lufj(A)
%解下三角方程组Ly=b
y=lowt(l,b)
%解上三角方程组Ux=y
x1=upt(u,y)
xb1=malu(A,b)
%以下解A^2x=b
%Ly=b
y=lowt(l,b);
%Uz=y
z=upt(u,y);
%Lw=z
w=lowt(l,z);
%Ux=w
x2=upt(u,w)
xb2=malu(A*A,b)
%以下解A^3x=b
%Ly=b
y=lowt(l,b);
%Uz=y
z=upt(u,y);
%Lw=z
w=lowt(l,z);
%Us=w
s=upt(u,w);
%Lt=s
t=lowt(l,s);
%Ux=t
x3=upt(u,t)
xb3=malu(A*A*A,b)

3.程序:
①n=10;
a=-ones(n,1); c=a;
b=4*ones(n,1);
d=5*ones(n,1); d(1)=6;
xb1=machase(a,b,c,d);
x1=xb1'
%用LU分解求解
A=zeros(n,n);
for(i=1:n)
A(i,i)=4;
end
for(i=1:n-1)
A(i,i+1)=-1;
end
for(i=2:n)
A(i,i-1)=-1;
end
x2=malu(A,d)

②n=100;
a=-ones(n,1); c=a;
b=4*ones(n,1);
d=5*ones(n,1); d(1)=6;
xb1=machase(a,b,c,d);
x1=xb1'
%用LU分解求解
A=zeros(n,n);
for(i=1:n)
A(i,i)=4;
end
for(i=1:n-1)
A(i,i+1)=-1;
end
for(i=2:n)
A(i,i-1)=-1;
end
x2=malu(A,d)

③n=1000;
a=-ones(n,1); c=a;
b=4*ones(n,1);
d=5*ones(n,1); d(1)=6;
xb1=machase(a,b,c,d);
x1=xb1'
%用LU分解求解
A=zeros(n,n);
for(i=1:n)
A(i,i)=4;
end
for(i=1:n-1)
A(i,i+1)=-1;
end
for(i=2:n)
A(i,i-1)=-1;
end
x2=malu(A,d)

五、实验结果
1. x1 =
3.7876
1.4680
-15.0620
10.3384
x2 =
NaN
NaN
NaN
NaN

2.⑴
l =

1 0 0 0 0 0
0 1 0 0 0 0
0 0 1 0 0 0
0 0 0 1 0 0
0 0 0 0 1 0
0 0 0 0 0 1

u =

6 5 4 3 2 1
0 6 5 4 3 2
0 0 6 5 4 3
0 0 0 6 5 4
0 0 0 0 6 5
0 0 0 0 0 6
⑵
x1 =

0.0074
0.0063
0.0054
0.0046
-0.1389
0.1667

xb1 =

0.0074
0.0063
0.0054
0.0046
-0.1389
0.1667

x2 =

0.0008
0.0006
0.0005
0.0208
-0.0463
0.0278

xb2 =

0.0008
0.0006
0.0005
0.0208
-0.0463
0.0278

x3 =

0.0001
0.0001
-0.0029
0.0100
-0.0116
0.0046

xb3 =

0.0001
0.0001
-0.0029
0.0100
-0.0116
0.0046

3. ①n=10
x1 =

2.0981
2.3923
2.4711
2.4920
2.4970
2.4960
2.4870
2.4519
2.3205
1.8301

x2 =

2.0981
2.3923
2.4711
2.4920
2.4970
2.4960
2.4870
2.4519
2.3205
1.8301
②n=100;
x1 =

2.0981
2.3923
2.4711
2.4923
2.4979
2.4994
2.4999
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.4999
2.4998
2.4991
2.4965
2.4871
2.4519
2.3205
1.8301

x2 =

2.0981
2.3923
2.4711
2.4923
2.4979
2.4994
2.4999
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.4999
2.4998
2.4991
2.4965
2.4871
2.4519
2.3205
1.8301
③n=1000;

x1 =

2.0981
2.3923
2.4711
2.4923
2.4979
2.4994
2.4999
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000
2.5000