计算方法作业第六章.docx
- 文档编号:3981907
- 上传时间:2022-11-26
- 格式:DOCX
- 页数:16
- 大小:30.86KB
计算方法作业第六章.docx
《计算方法作业第六章.docx》由会员分享,可在线阅读,更多相关《计算方法作业第六章.docx(16页珍藏版)》请在冰豆网上搜索。
计算方法作业第六章
1.考虑两个线性方程组,其系数矩阵如下
问题的真解均取为
,线性方程组的右端项用这个真解计算出来。
相应的问题分别称为问题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);
fori=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);
fori=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)解上三角方程组源程序:
%解上三角方程组
functionx=SolveUpTriangle(A,b)
N=size(A);
n=N
(1);
x(n)=b(n)/A(n,n);
fork=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);
fori=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方法解线性方程组
functionx=Cholesky(A,b);
N=size(A);
n=N
(1);
L(1,1)=sqrt(A(1,1));
%L为下三角矩阵
forj=2:
n
L(j,1)=A(1,j)./L(1,1)
end
fork=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的转置
forj=1:
n
for k=1:
n
LT(j,k)=L(k,j);
end
end
y=SolveDownTriangle(L,b);
x=SolveUpTriangle(LT,y);
解上三角方程组源程序:
%解上三角方程组
functionx=SolveUpTriangle(A,b)
N=size(A);
n=N
(1);
x(n)=b(n)/A(n,n);
fork=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
解下三角方程组源程序:
functionx=SolveDownTriangle(A,b)
N=size(A);
n=N
(1);
x
(1)=b
(1)/A(1,1);
fork=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:
计算程序:
n=100;
A=2*eye(n);
fori=1:
n-1
A(i+1,i)=-1;
A(i,i+1)=-1;
end
b=0;
b
(1)=1;
b(100)=1;
x3=GJCholesky(A,b);
(GJCholesky.m):
Cholesky法源程序
%用改进Cholesky法求解线性方程组
%记LD为U,U为下三角矩阵,LT为单位上三角矩阵
functionx=GJCholesky(A,b);
N=size(A);
n=N
(1);
fork=1:
n
U(k,1)=A(k,1);
end
fori=2:
n
L(k,1)=A(1,k)./U(1,1);
end
fork=2:
n
for j=k:
n
s=0;
for i=1:
k-1
s=s+U(j,i).*L(j,i);
end
U(j,k)=A(j,k)-s;
end
for j=k+1:
n
sum=0;
for i=1:
k-1
sum=sum+U(k,i).*L(j,i);
end
L(j,k)=(A(k,j)-sum)./U(k,k);
end
end
fork=1:
n
D(k,k)=U(k,k);
L(k,k)=1;
LT(k,k)=1;
for j=k+1:
n
LT(k,j)=L(j,k);
end
end
y=SolveDownTriangle(L,b);
fori=1:
n
z(i)=y(i)./D(i,i);
end
x=SolveUpTriangle(LT,z);
解上三角方程组源程序:
%解上三角方程组
functionx=SolveUpTriangle(A,b)
N=size(A);
n=N
(1);
x(n)=b(n)/A(n,n);
fork=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
解下三角方程组源程序:
functionx=SolveDownTriangle(A,b)
N=size(A);
n=N
(1);
x
(1)=b
(1)/A(1,1);
fork=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
结果:
n=8
x1=[0,0,0,…..0,0,0.9891]T;x2=[0,0,0,…..0,0,0.9889]T;
n=10
x1=[0,0,0,…..0,0,0.9889]T;x2=[0,0,0,…..0,0,0.9885]T;
n=12
x1=[0,0,0,…..0,0,0.9886]T;x2=[0,0,0,…..0,0,0.9885]T;
n=14
x1=[0,0,0,…..0,0,0.9884]T;x2=[0,0,0,…..0,0,0.9884]T;
n=16
x1=[0,0,0,…..0,0,0.9881]T;x2=[0,0,0,…..0,0,0.9882]T;
n=18
x1=[0,0,0,…..0,0,0.9878]T;x2=[0,0,0,…..0,0,0.9879]T;
n=20
x1=[0,0,0,…..0,0,0.9875]T;x2=[0,0,0,…..0,0,0.9875]T;
追赶法:
计算程序:
n=100;
A=2*eye(n);
fori=1:
n-1
A(i+1,i)=-1;
A(i,i+1)=-1;
end
b=0;
b
(1)=1;
b(100)=1;
x4=zhuigan(A,b);
追赶法源程序
%用追赶法求解线性方程组
functionM=zhuigan(A,g)
n=length(A);
L=eye(n);
U=zeros(n);
fori=1:
n-1
U(i,i+1)=A(i,i+1);
end
U(1,1)=A(1,1);
fori=2:
n
L(i,i-1)=A(i,i-1)/U(i-1,i-1);
U(i,i)=A(i,i)-L(i,i-1)*A(i-1,i);
end
Y
(1)=g
(1);
fori=2:
n
Y(i)=g(i)-L(i,i-1)*Y(i-1);
end
M(n)=Y(n)/U(n,n);
fori=n-1:
-1:
1
M(i)=(Y(i)-A(i,i+1)*M(i+1))/U(i,i);
end
(2)
Gauss消去:
计算程序:
n=6;
fori=1:
n
s=0;
for j=1:
n
A(i,j)=1/(i+j-1);
s=s+A(i,j);
end
b(i)=s;
end
[x1,XA]=GaussJordanXQ(A,b);
Gauss消元法源程序
%用Gauss消元法解线性方程组
function[x,XA]=GaussJordanXQ(A,b)
N=size(A);
n=N
(1);
fori=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;
%消元后的系数矩阵
解上三角方程组源程序:
%解上三角方程组
functionx=SolveUpTriangle(A,b)
N=size(A);
n=N
(1);
x(n)=b(n)/A(n,n);
fork=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
列主元Gauss消去:
计算程序:
n=6;
fori=1:
n
s=0;
for j=1:
n
A(i,j)=1/(i+j-1);
s=s+A(i,j);
end
b(i)=s;
end
[x2,XA]=GaussLineMainXQ(A,b);
列主元Gauss消去法源程序
%列主元Gauss消去法解线性方程组
functionx=GaussLinemainXQ(A,b);
N=size(A);
n=N
(1);
index=0;
fori=1:
n-1
main=max(abs(A(1:
n,i)));
%选取列主元
for k=i:
n
if (abs(A(k,i)==main)
index=k;
break;
end
%保存主元所在行的指标
end
temp=A(i,1:
n);
A(i,1:
n)=A(index,1:
n);
A(index,1:
n)=temp;
btemp=b(index);
b(index)=b(i);
b(i)=btemp;
%交换主元所在行
for j=(i+1):
n
if (A(i,i)==0)
disp('对角线元素为0’);
return;
end
l=A(j.i);
m=A(i,i);
A(j,1:
n)=A(j,1:
n)-A(i,1:
n).*l/m;
b(j)=b(j)-b(i).*l/m;
end
end
x=SolveUpTriangle(A,b);
解上三角方程组源程序:
%解上三角方程组
functionx=SolveUpTriangle(A,b)
N=size(A);
n=N
(1);
x(n)=b(n)/A(n,n);
fork=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
2.考虑矩阵
,其中
。
利用Gauss-Jordan方法求解逆矩阵,并验证逆矩阵与A的乘积是否为单位矩阵
Gauss-Jordan方法程序:
functionx=GaussJordan(A)
[n,n]=size(A);
whiledet(A)==0
break
end
s=1;P=zeros(1,n);
whiles<=n
max=abs(A(s,s));
big=0;
form=s:
n
ifmax max=abs(A(m,s)); q=m; big=1; elsecontinue end end ifbig==1 P(s)=q; D=A(s,: ); A(s,: )=A(q,: ); A(q,: )=D; end fori=1: n ifi~=s forj=1: n ifj~=s A(i,j)=A(i,j)-A(s,j)*A(i,s)/A(s,s); elsecontinue end end elsecontinue end A(i,s)=-A(i,s)/A(s,s); end fork=1: n ifk A(s,k)=A(s,k)/A(s,s); end ifk==s A(s,k)=1/A(s,s); end ifk>s A(s,k)=A(s,k)*A(s,s); end end s=s+1; end fori=n: (-1): 1 ifP(i)~=0 p=P(i); d=A(: i); A(: i)=A(: p); A(: p)=d; elsecontinue end end A 调用程序为: formatrat n=6; a=zeros([n,n]); fori=1: n forj=1: n a(i,j)=i-j; End end b=eye([n,n])+a B=GaussJordan(b) B*b 结果为: B= 51/106-39/106-23/106-7/1069/10625/106 -41/10675/106-21/106-11/106-1/1069/106 -27/106-23/10687/106-15/106-11/106-7/106 -13/106-15/106-17/10687/106-21/106-23/106 1/106-7/106-15/106-23/10675/106-39/106 15/1061/106-13/106-27/106-41/10651/106 B*b ans= 1.00000.00000-0.0000-0.0000-0.0000 -0.00001.0000-0.0000-0.00000.0000-0.0000 -0.0000-0.00001.00000.0000-0.00000.0000 -0.0000-0.0000-0.00001.0000-0.0000-0.0000 -0.0000-0.0000-0.0000-0.00001.00000.0000 0.00000.00000.00000.00000.00001.0000 是单位矩阵
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 计算方法 作业 第六