5 线性方程组的直接解法

实验五 线性方程组的直接解法
一、实验内容
1、用列主元素法求解方程组
15
123459.170.31059.43
146.785.291
6.3112111.295221211x x x x -⎡⎤⨯⎡⎤⎡⎤⎢⎥⎢⎥⎢⎥--⎢⎥⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦⎣⎦⎣⎦ 并计算误差b-Ax
，分析结果的好坏；
2、 用改进Cholesky 方法求对称正定阵线性方程组 1234248.72171013.741090.7x x x -⎡⎤⎡⎤⎡⎤⎢⎥⎢⎥⎢⎥-=⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥-⎣⎦⎣⎦⎣
⎦
并计算误差b-Ax ，分析结果的好坏；
3、 用追赶法解方程组 123421006132010121000351x x x x -⎡⎤⎡⎤⎡⎤⎢⎥⎢⎥⎢⎥--⎢⎥⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥--⎢⎥⎢⎥⎢⎥-⎣⎦⎣⎦
⎣⎦ 二、要求
1、 对上述三个方程组分别利用Gauss 列主元消去法；Cholesky 方法；追赶法求解（选择其一）；
2、 应用结构程序设计编出通用程序；
3、 比较计算结果，分析数值解误差的原因；
三、目的和意义
1、通过该课题的实验，体会模块化结构程序设计方法的优点；
2、运用所学的计算方法，解决各类线性方程组的直接算法；
3、提高分析和解决问题的能力，做到学以致用；
4、 通过三对角形线性方程组的解法，体会稀疏线性方程组解法的特点。

四、实验学时：2学时
五、实验步骤：
1．进入matlab 开发环境；
2．根据实验内容和要求编写程序；
3．调试程序；
4．运行程序；
5．撰写报告,讨论分析实验结果．
六、程序
1、Gauss列主元素消去法
function x=Gauss_pivot(A,b)
%用Gauss列主元素法求解线性方程组Ax=b
%x是未知向量
n=length(b);
x=zeros(n,1);
c=zeros(1,n);
d1=0;
%消元计算
for i=1:n-1
max=abs(A(i,i));
m=i;
for j=i+1:n
if max<abs(A(j,i));
max=abs(A(j,i));
m=j;
end
end
if m~=i
for k=i:n
c(k)=A(i,k);
A(i,k)=A(m,k);
A(m,k)=c(k);
end
d1=b(i);
b(i)=b(m);
b(m)=d1;
end
for k=i+1:n
for j=i+1:n
A(k,j)=A(k,j)-A(i,j)*A(k,i)/A(i,i);
end
b(k)=b(k)-b(i)*A(k,i)/A(i,i);
A(k,i)=0;
end
end
%回代求解
x(n)=b(n)/A(n,n);
for i=n-1:-1:1
sum=0;
for j=i+1:n
sum=sum+A(i,j)*x(j);
end
x(i)=(b(i)-sum)/A(i,i);
end
2、Cholesky方法
function x=Chol_decompose(A,b)
%用cholesky求解线性方程组Ax=b
%A是对称矩阵
%L是单位下三角阵
%D是对角阵
%对角矩阵A进行三角分解：A=LDL'
N=length(A);
L=zeros(N,N);D=zeros(1,N);
for i=1:N
L(i,i)=1;
end
D(1)=A(1,1);
for i=2:N
for j=1:i-1
if j==1
L(i,j)=A(i,j)/D(j);
else
sum1=0;
for k=1:j-1
sum1=sum1+L(i,k)*D(k)*L(j,k);
end
L(i,j)=(A(i,j)-sum1)/D(j);
end
end
sum2=0;
for k=1:i-1
sum2=sum2+L(i,k)^2*D(k);
end
D(i)=A(i,i)-sum2;
end
%分别求线性方程组Ly=b;L'x=y/D
y=zeros(1,N);
y(1)=b(1);
for i=2:N
sumi=0;
for k=1:i-1
sumi=sumi+L(i,k)*y(k);
end
y(i)=b(i)-sumi;
end
x=zeros(1,N);
x(N)=y(N)/D(N);
for i=N-1:-1:1
sumi=0;
for k=i+1:N
sumi=sumi+L(k,i)*x(k);
end
x(i)=y(i)/D(i)-sumi;
end
3、追赶法
function x=threedia(a,b,c,f)
%求解线性方程组Ax=f，其中A是三对角阵
%a是矩阵A的下对角线元素a(1)=0
%b是矩阵A的对角线元素
%c是矩阵A的上对角线元素c(N)=0
%f是方程组的右端向量
N=length(f);
x=zeros(1,N);y=zeros(1,N);
d=zeros(1,N);u=zeros(1,N);
%预处理
d(1)=b(1);
for i=1:N-1
u(i)=c(i)/d(i);
d(i+1)=b(i+1)-a(i+1)*u(i); % u在这里表示就是公式中的 end
%追的过程
y(1)=f(1)/d(1);
for i=2:N
y(i)=(f(i)-a(i)*y(i-1))/d(i);
end
%赶的过程
x(N)=y(N);
for i=N-1:-1:1
x(i)=y(i)-u(i)*x(i+1);
end
x
七、数值实验
1、第1题运行程序
A=[0.3*^10-15,59.14,3,1;5.291,-6.13,-1,2;11.2,9,5,2;1,2,1,1];
b=[59.17,46.78,1,2];
x=Gauss_pivot(A,b);
b-A*x
2、第2题运行程序
A=[4,-2,4;-2,17,10;4,10,9];
b=[8.7,13.7,-0.7];
x=Chol_decompose(A,b);
3、第3题运行程序
a=[0,-1,-1,-3];
b=[2,3,2,5];
c=[-1,-2,-1,0];
f=[6,1,0,1];
x=threedia(a,b,c,f)
八、思考题
根据高斯列主元素消去法写出高斯消去法的程序，用高斯消去法解算第1题，数值结果与高斯列主元消去法进行比较优缺点。