数值分析实验一

数值分析 实验一
一、实验目的
熟悉MATLAB 编程。

学习线性方程组数值解法的程序设计算法。

二、实验题目
1.给定线性方程组
（a ）用LU 分解和列主元消去法求解。

输出A=LU 的分解的L 和U ，detA 以及解向量x 。

（b ）将2.099999改为2.1，5.900001改为5.9.用列主元消去法求解。

输出detA 及解向量 x ，并与（a ）的结果比较。

2。

线性方程组Ax=b 的A 及b 为
10 7 8 9 32
7 5 6 5 23
A= 8 6 10 9 b= 33
7 5 9 10 31
则解x=（1,1,1,1）T 用MATLAB 内部函数求detA 及A 的所有特征值和cond （A ）2。

若令
⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡=⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡----15900001.582012151526099999.2310
7104321x x x x
10 7 8.1 7.2
A+δA= 7.08 5.04 6 5
8 5.98 9.89 9
6.99 5 9 9.98
求解（A+δA ）（x+δx ）=b,输出向量δx 和2 从理论结果和实际计算两方面分析线性方程组Ax=b 解的相对误差2 及A 的相对误差2 的关系。

三、实验原理与理论基础
1.本题用LU 分解法及列主元高斯消去法
LU 分解法原理：
解Ax=b 相当于解两个三角方程，即 Ly=b,Ux=y 分解求出x 和y 。

列主元消去法原理：
利用逐次消去未知数的方法，把线性方程组 Ax=b 化为与其等价的三角形线性方程组，求解线性方程组的方法可用回代的方法求。

2.本题通过 计算式：
)
()()(min max 2122T T AA A A A
A A cond λλ==- 来计算A 的谱子条件数。

四、实验内容
1.解：
（a ）LU 分解法：
程序如下：
function [ x ] = zhijiejiefangcheng( A,b )
n=length(A);
%求出LU分解矩阵
for i=1:n
U(1,i)=A(1,i);
end
L(1,1)=A(1,1);
for i=2:n
L(i,1)=A(i,1)/U(1,1);
end
for r=2:n
for i=r:n
p=0;q=0;
for k=1:(r-1)
p=p+L(r,k)*U(k,i);
q=q+L(i,k)*U(k,r);
end
U(r,i)=A(r,i)-p;
if r<n L(i,r)=(A(i,r)-q)/U(r,r);
end
end
end
for i=1:n
for j=1:n
if i==j
L(i,j)=1;
end
end
end
L
U
%计算X的值
y(1)=b(1);
for i=2:n
e=0;
for k=1:(i-1)
e=e+L(i,k)*y(k);
end
y(i)=b(i)-e;
end
X(n)=y(n)/U(n,n);
for i=(n-1):-1:1
f=0;
for k=(i+1):n
f=f+U(i,k)*X(k);
end
X(i)=(y(i)-f)/U(i,i);
end
X
%计算detA的值
detA=det(L*U)
列主元消去法：
程序如下：
function [ x ] = gauss( A,b )
A=[10 -7 0 1;-3 2.099999 6 2;5 -1 5 -1;2 1 0 2];
b=[8;5.900001;5;1];
[n,n]=size(A);
x=zeros(n,1);
Aug=[A,b]; %增广矩阵
for k=1:n-1
[piv,r]=max(abs(Aug(k:n,k))); %找列主元所在子矩阵的行r r=r+k-1; %列主元所在矩阵的行
if r>k
temp=Aug(k,:);
Aug(k,:)=Aug(r,:);
Aug(r,:)=temp;
end
if Aug(k,k)==0,
'error';
end
for p=k+1:n
Aug(p,:)=Aug(p,:)-Aug(k,:)*Aug(p,k)/Aug(k,k);
end
end
A=Aug(:,1:n);b=Aug(:,n+1);
x(n)=b(n)/A(n,n);
for k=n-1:-1:1
x(k)=b(k);
for p=n:-1:k+1
x(k)=x(k)-A(k,p)*x(p);
end
x(k)=x(k)/A(k,k);
end
detA=det(A)
End
（b）改动数字后列主元消去法
程序如下：
function [ x ] = gauss2( A,b )
A=[10 -7 0 1;-3 2.1 6 2;5 -1 5 -1;2 1 0 2];
b=[8;5.1;5;1];
[n,n]=size(A);
x=zeros(n,1);
Aug=[A,b]; %增广矩阵
for k=1:n-1
[piv,r]=max(abs(Aug(k:n,k))); %找列主元所在子矩阵的行r
r=r+k-1; %列主元所在矩阵的行
if r>k
temp=Aug(k,:);
Aug(k,:)=Aug(r,:);
Aug(r,:)=temp;
end
if Aug(k,k)==0,
'error';
end
for p=k+1:n
Aug(p,:)=Aug(p,:)-Aug(k,:)*Aug(p,k)/Aug(k,k);
end
end
A=Aug(:,1:n);b=Aug(:,n+1);
x(n)=b(n)/A(n,n);
for k=n-1:-1:1
x(k)=b(k);
for p=n:-1:k+1
x(k)=x(k)-A(k,p)*x(p);
end
x(k)=x(k)/A(k,k);
end
detA=det(A)
end
2.解：求detA及A的所有特征值和cond（A）2的程序如下：
function [x] = jhls( A )
A=[10 7 8 7;7 5 6 5;8 6 10 9;7 5 9 10];
detA=det(A) %求A的行列式值
[D]=eig(A) %求行列式的特征值
x=cond(A,2) %求cond（A）2
end
求解(A+&A)(x+&x)=b 输出向量2x δδ和x 的程序如下：
>> A=[10 7 8.1 7.2;7.08 5.04 6 5;8 5.98 9.89 9;6.99 5 9 9.98];
>> b=[32 23 33 31]';
>> x=A\b
x =
-9.5863
18.3741
-3.2258
3.5240
>> x0=[1 1 1 1]';
>> x1=x-x0
x1 =
10.5863
17.3741
-4.2258
2.5240
>> A=[10 7 8.1 7.2;7.08 5.04 6 5;8 5.98 9.89 9;6.99 5 9 9.98];
A0=[10 7 8 7;7 5 6 5;8 6 10 9;7 5 9 10];
b=[32 23 33 31]';
x0=[1 1 1 1]';
x=A\b;
x1=x-x0;
norm(x1,2); % ||δx||2/||x||2的值
A1=A-A0; %||δA||2/||A||2的值
>> norm(x1,2)/norm(x0,2)
ans =
10.4661
>> norm(A1,2)/norm(A0,2)
ans =
0.0076
由上式计算结果可知，A相对误差为0.0076，而x的相对误差为10.4661，所以方程组是病态的，A的条件数为2984.1远大于1，当A只有很小的误差就会给结果带来很大的影响。

五、实验结果
1.（a）LU分解法
列主元消去法：
（b）改动数字后计算结果
虽然只对矩阵中两个数做了很小的改动，但是计算结果却有很大差别。

对行列式A值没有太大影响，对解向量x的影响大。

2.解：求得detA及A的所有特征值和cond(A)2的结果为：
由上式计算结果可知，A相对误差为0.0076，而x的相对误差为10.4661，所以方程组是病态的，A的条件数为2984.1远大于1，当A只有很小的误差就会给结果带来很大的影响。

六、实验结果分析与小结
通过本次试验又熟悉了MATLAB软件，来解决一些线性方程问题。

学习了用LU分解法和列主元分解法。

通过对矩阵经一系列运算进行误差分析，分析是否为病态的。

知道了很小的能造成很大的误差，所以计算要认真仔细，避免一些较大错误。