《计算方法》课内实验报告
学生姓名: 及 学 号： 学 院: 班 级: 课程名称： 实验题目: 指导教师 姓名及职称：
张学阳 理学院 数学101 计算方法
解线性方程组的直接方法和迭代法
宋云飞 讲 师 朱秀丽 讲 师 尚宝欣 讲 师
1009300132
2012年12月10日
目 录
一、实验题目........................................................................................ - 1 - 二、实验目的........................................................................................ - 1 - 三、实验内容........................................................................................ - 1 - 四、实验结果........................................................................................ - 2 - 五、实验体会或遇到问题 ................................................................... - 8 -
一、实验题目
解线性方程组的直接方法和迭代法
二、实验目的
1．熟悉Matlab编写及运行数值计算程序的方法。
2．进一步理解求解线性方程组的直接方法和迭代法基础理论。 3．进一步掌握应用不同的方法求解线性方程组的收敛速度及误差分析。
三、实验内容
1．用LU分解及列主元高斯消去法解线性方程组
?701??x1??8?10???32.09999692??x2??5.9000?01??????. ??15?1??x3??5?5???????102??x4??1?2?输出Ax?b中系数A?LU分解法的矩阵L及U，解向量x及detA；列主元法的行交换次序，解向量x及detA；比较两种方法所得的结果. 2．线性方程组Ax?b的A及b为
?10?7 A???8??7787??32??23?565??,b???,
6109??33????5910??31?则解x?(1,1,1,1)T.用Matlab内部函数求detA及A的所有特征值和cond(A)2.若令
78.17.2??10?7.085.0465??， A??A??5.989.899??8??99.98??6.995求解(A??A)(x??x)?b，输出向量?x和?x2，从理论结果和实际计算两方面分析线性方程组Ax?b解的相对误差?x2x2及A的相对误差?A2A2的关系。
- 1 -
3．给出线性方程组Hnx?b，其中系数矩阵Hn为希尔伯特矩阵： Hn??hij??Rn?n，hij?T1,i,j?1,2,?,n.
i?j?1假设x*??1,1,?,1??Rn,b?Hnx*，若取n?6,8，分别用雅可比迭代法及SOR迭代法(??1,1.25,1.5)求解，比较计算结果。
四、实验结果
1.文件程序：
function hl=zhjLU(A,b)
[n n] =size(A); RA=rank(A); if RA~=n
disp('请注意：因为A的n阶行列式hl等于零，所以A不能进行LU分解.A的秩RA如下:'), RA,hl=det(A); return end
if RA==n for p=1:n
h(p)=det(A(1:p, 1:p)); end
hl=h(1:n); for i=1:n if h(1,i)==0
disp('请注意：因为A的r阶主子式等于零，所以A不能进行LU分解.A的秩RA和各阶顺序主子式值hl依次如下：'), hl;RA return end end
if h(1,i)~=0
disp('请注意：因为A的各阶主子式都不等于零，所以A能进行LU分解.A的秩RA和各阶顺序主子式值hl依次如下：') for j=1:n
U(1,j)=A(1,j); end
for k=1:n for i=2:n
for j=2:n
L(1,1)=1;L(i,i)=1; if i>j
L(1,1)=1;L(2,1)=A(2,1)/U(1,1); L(i,1)=A(i,1)/U(1,1);
- 2 -
L(i,k)=(A(i,k)- L(i,1:k-1)*U(1:k-1,k))/U(k,k); else
U(k,j)=A(k,j)-L(k,1:k-1)*U(1:k-1,j); end end end end
hl;RA,U,L end end
文件程序：
function [RA,RB,n,X]=liezhu(A,b) B=[A b]; n=length(b); RA=rank(A); RB=rank(B);zhica=RB-RA; if zhica>0,
disp('请注意：因为RA~=RB，所以此方程组无解.') return end
if RA==RB if RA==n
disp('请注意：因为RA=RB=n，所以此方程组有唯一解.') X=zeros(n,1); C=zeros(1,n+1); for p= 1:n-1
[Y,j]=max(abs(B(p:n,p))); C=B(p,:); B(p,:)= B(j+p-1,:); B(j+p-1,:)=C; for k=p+1:n
m= B(k,p)/ B(p,p);
B(k,p:n+1)= B(k,p:n+1)-m* B(p,p:n+1); end end
b=B(1:n,n+1);A=B(1:n,1:n); X(n)=b(n)/A(n,n); for q=n-1:-1:1
X(q)=(b(q)-sum(A(q,q+1:n)*X(q+1:n)))/A(q,q); end else
disp('请注意：因为RA=RB 窗口输入输出 >> A=[10 -7 0 1;-3 2.099999 6 2;5 -1 5 -1;2 1 0 2]; >> b=[8 5.9000001 5 1]'; >> zhjLU(A,b) 请注意：因为A的各阶主子式都不等于零，所以A能进行LU分解.A的秩RA和各阶顺序主子式值hl依次如下： - 3 - RA =4 %U矩阵 U = 10.0000 -7.0000 0 1.0000 0 2.1000 6.0000 2.3000 0 0 -2.1429 -4.2381 0 -0.0000 0 12.7333 %L矩阵 L = 1.0000 0 0 0 -0.3000 1.0000 0 0 0.5000 1.1905 1.0000 -0.0000 0.2000 1.1429 3.2000 1.0000 %矩阵顺序主子试 ans = 10.0000 -0.0000 -150.0001 -762.0001 >> y=inv(L)*b >> x=inv(U)*y %解向量 x = -0.2749 -1.3299 1.2969 1.4398 %行列式 >> det(L)*det(U) ans = -573.0100 >> [RA,RB,n,X]=liezhu(A,b) 请注意：因为RA=RB=n，所以此方程组有唯一解. A = 10.0000 -7.0000 0 1.0000 0 2.5000 5.0000 -1.5000 0 0 6.0000 2.3000 0 0 0 5.0800 RA =4 RB = 4 n =4 %解向量 X = 0.0000 -1.0000 1.0000 - 4 - 1.0000 %行列式 >> det(A) ans = -762.0001 2． >> A=[10 7 8 7;7 5 6 5;8 6 10 9;7 5 9 10]; >> b=[32 23 33 31]'; >> cond(A,2)%求cond(A)2. ans =2.9841e+003 >> det(A)%求行列式 ans =1.0000 >> eig(A)%求A的特征值 ans =0.0102 0.8431 3.8581 30.2887 >> A1=[10 7 8.1 7.2;7.08 5.04 6 5;8 5.98 9.89 9;6.99 5 9 9.98]; >> x=[1 1 1 1]'; >> I=[1 0 0 0;0 1 0 0;0 0 1 0;0 0 0 1]; >> x1=-inv(I+inv(A)*(A1-A))*inv(A)*(A1-A)*x%求?x x1 = -10.5863 17.3741 -4.2258 2.5240 >> norm(x1,2)%求?x2 ans = 20.9322 >> norm(x1,2)/norm(x,2)%相对误差?x2x2 ans = 5.2931 >> norm(A1-A,2)/norm(A,2)%A的相对误差?A2A2 ans =0.0076 3. 雅可比迭代文件一 function X=jacdd(A,b,X0,P,wucha,max1) [n m]=size(A); for j=1:m a(j)=sum(abs(A(:,j)))-2*(abs(A(j,j))); end for i=1:n if a(i)>=0 disp('请注意：系数矩阵A不是严格对角占优的，此雅可比迭代不一定收敛') %return end end - 5 - if a(i)<0 disp('请注意：系数矩阵A是严格对角占优的，此方程组有唯一解，且雅可比迭代收敛 ') end for k=1:max1 k for j=1:m X(j)=(b(j)-A(j,[1:j-1,j+1:m])*X0([1: j-1,j+1:m]))/A(j,j); end X,djwcX=norm(X'-X0,P); xdwcX=djwcX/(norm(X',P)+eps); X0=X';X1=A\\b; if (djwcX disp('请注意：雅可比迭代收敛,此方程组的精确解jX和近似解X如下：') return end end if (djwcX>wucha)&(xdwcX>wucha) disp('请注意：雅可比迭代次数已经超过最大迭代次数max1 ') end a,X=X;jX=X1', 文件二 n=8; x0=[]; for i=1:n for j=i:n A(i,j)=1/(i+j-1); A(j,i)=A(i,j); end x0=[x0 1]; end x0=x0'; b=A*x0; jacdd(A,b,x0,2,0.001,100) 运行结果 n=6 ans = 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 n=8 ans = Columns 1 through 6 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 Columns 7 through 8 1.0000 1.0000 超松弛迭代文件一 function X=jacdd(A,b,X0,P,wucha,max1,w) [n m]=size(A); for j=1:m a(j)=sum(abs(A(:,j)))-2*(abs(A(j,j))); - 6 - end for i=1:n if a(i)>=0 disp('请注意：系数矩阵A不是严格对角占优的，此雅可比迭代不一定收敛') %return end end if a(i)<0 disp('请注意：系数矩阵A是严格对角占优的，此方程组有唯一解，且雅可比迭代收敛 ') end for k=1:max1 k for j=1:m X(j)=(b(j)-A(j,[1:j-1,j+1:m])*X0([1: j-1,j+1:m]))/A(j,j); end X,djwcX=norm(X'-X0,P); xdwcX=djwcX/(norm(X',P)+eps); X0=X';X1=A\\b; if (djwcX disp('请注意：雅可比迭代收敛,此方程组的精确解jX和近似解X如下：') return end end if (djwcX>wucha)&(xdwcX>wucha) disp('请注意：雅可比迭代次数已经超过最大迭代次数max1 ') end a,X=X;jX=X1', 文件二 n=8; w=1; x0=[]; for i=1:n for j=i:n A(i,j)=1/(i+j-1); A(j,i)=A(i,j); end x0=[x0 1]; end x0=x0'; b=A*x0; jacdd(A,b,x0,2,0.001,100,w) 运行结果 n=6,w=1 ans =1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 w=1.25 - 7 - ans = 1.0e+072 * 0.8986 1.9393 2.5849 3.0372 3.3749 3.6378 (不收敛) w=1.5 ans =1.0e+080 * 1.2736 2.7485 3.6635 4.3046 4.7832 5.1558 (不收敛) n=8,w=1 ans = Columns 1 through 6 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 Columns 7 through 8 1.0000 1.0000 w=1.25 ans = 1.0e+087 * Columns 1 through 6 0.3650 0.8261 1.1329 1.3580 1.5322 1.6716 Columns 7 through 8 1.7860 1.8818 (不收敛) w=1.5 ans =1.0e+095 * Columns 1 through 6 0.5139 1.1629 1.5948 1.9118 2.1569 2.3531 Columns 7 through 8 2.5142 2.6490 （不收敛） 五、实验体会或遇到问题 通过这次课内试验熟悉了Matlab编写及运行数值计算程序的方法，进一步理解了求解线性方程组的直接方法和迭代法基础理论，掌握了应用不同的方法求解线性方程组的收敛速度及误差分析，感觉收获很大。 - 8 -