数值计算报告2Word下载.docx
- 文档编号:21119221
- 上传时间:2023-01-27
- 格式:DOCX
- 页数:14
- 大小:61.09KB
数值计算报告2Word下载.docx
《数值计算报告2Word下载.docx》由会员分享,可在线阅读,更多相关《数值计算报告2Word下载.docx(14页珍藏版)》请在冰豆网上搜索。
close(10)
2理论简述
本文采用了三种方式进行matlab仿真——gauss消元法、Jacob迭代、SOR超松弛迭代,经过实验可知道Jacob迭代是失败的,并不收敛。
2.1gauss消元法
对于线性方程组Ax=b,经过n-1次行初等变换,将B=(A,b)变成上三角阵。
对上式进行回带,从而化简成为单位I阵,从而得到解形式如下所示。
其中,
。
该方法对线性方法是最原始的最准确的代数解法。
2.2Jacobi迭代法
如果假设如下两系数矩阵,那么可以得到Jacobi迭代公式
式子中A为方程系数矩阵,b为右端常数项。
2.3SOR迭代公式
ω为可选择的松弛因子,那么得到如下迭代公式。
在松弛迭代公式里有向前和向后的超松弛迭代方法,这里就不详细阐述。
3仿真实验
3.1数据提取程序
通过matlab程序提取txt文档里将其保存在A,b中,从而提取Ax=b的方程,程序如下所示。
%%初始化
clearall;
clc;
%%数据提取
A=textread('
);
%A=textread('
test.txt'
N=A(1,1);
fori=1:
N
forj=1:
N+1
B(i,j)=A((i-1)*(N+1)+j+1,1);
end
end
%清除不必要数据节约内存
C=B(:
1:
N);
f=B(:
N+1);
clearABij;
得到的矩阵C和f
3.2Gauss消元法
对于一个3x3的矩阵Ax=b进行验证计算,首先我们已知该方程的矩阵A,b,以及解x;
理论解:
用gauss解出结果如所示
图1方程解分布图
通过计算得到解如图2所示,程序如gaus_program.m所示。
图2线性Ax=b对应解
3.3Jacobi迭代法
取eps<
0.01为精确解条件,通过迭代后发现该种方法并不收敛,随着迭代的次数不断的增加误差也越来越大,最后超出内存,程序代码如Jacob_iteration.m所示。
3.4SOR—迭代法
迭代程序如SOR_iteration.m和SSOR_iteration.m所示,其中,松弛因子ω取不同的值会得到有不同的迭代精度和迭代时间,得到ω=0.5情况下验证矩阵解和收敛趋势分别如图3和图4所示。
图3ω=0.5时方程的解
图4ω=0.5时的线性方程的收敛趋势
对于matrix.txt文档里的数据根本是不收敛的,松弛迭代方法在w=0.5,1.5都不能收敛,其收敛曲线ω=0.5迭代误差曲线如图5所示。
图5不收敛解趋势
4结论分析
从结果上看
1.Guass消元法能很好的求出结果,不存在收敛与否的问题;
2.SOR和SSOR在针对解线性方程上有一定的局限,不能保证所有的矩阵都能收敛,本矩阵就是一个例子。
5代码附录
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%SSOR%%%%%%%%%%%%%%%%%
%%清除数据
%%读取数据
%提取系数矩阵和常数矩阵Cx=f
A=C;
B=[Af];
b=f;
%%计算上三角、下三角、对角矩阵
L=-tril(A,-1);
D=triu(A,0)-triu(A,1);
U=-triu(A,1);
%松弛因子w
w=0.5;
%矩阵计算
D1=inv(D);
E=D1*L;
F=D1*U;
I=eye(N);
E1=inv(I-w*E);
F1=inv(I-w*F);
Lw=E1*((1-w)*I+w*F);
Uw=F1*((1-w)*I+w*E);
Sw=Lw*Uw;
kw=w*(2-w)*E1*F1*D1*b;
%初值选取
Y=10*ones(N,1);
%%迭代程序
eps=100;
kk=1;
whileeps>
0.001
X=Sw*Y+kw;
eps_temp=sum(abs(X-Y))
eps(kk)=eps_temp;
kk=kk+1;
Y=X;
%%画出解的图,便于点数多后查看
figure
(1);
plot(X);
gridon;
xlabel('
X(n)中的n'
figure
(2);
plot(eps);
ylabel('
X的数值解'
w=1.5;
P=inv(D-w*L);
Gw=P*[(1-w)*D+w*U];
fw=w*P*b;
%X的初值
X=Gw*Y+fw;
eps(kk)=sum(abs(X-Y))
%%%清除数据
%clearall;
%clc;
%%%读取数据
%%A=textread('
%N=A(1,1);
%fori=1:
%forj=1:
%B(i,j)=A((i-1)*(N+1)+j+1,1);
%end
%%提取系数矩阵和常数矩阵Cx=f
%C=B(:
%f=B(:
%clearABij;
%A=C;
%B=[Af];
%b=f;
%系数矩阵
Bj=inv(D)*(L+U);
fj=inv(D)*b;
Y=0.01*ones(N,1);
X=Bj*Y+fj;
eps=sum(abs(X-Y))
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%高斯解法计算
clearCf;
n=length(b);
RA=rank(A);
RB=rank(B);
judge=RB-RA;
ifjudge>
0,
disp('
因为AB秩不同,所以此方程组无解.'
return
ifRA==RB
ifRA==n
因为AB秩相同且为N,所以此方程组有唯一解.'
%RA==n
X=zeros(n,1);
C=zeros(1,n+1);
forp=1:
n-1
fork=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);
b=B(1:
n,n+1);
A=B(1:
n,1:
n);
X(n)=b(n)/A(n,n);
forq=n-1:
-1:
1
X(q)=(b(q)-sum(A(q,q+1:
n)*X(q+1:
n)))/A(q,q);
else
因为R(A)=R(B)<
n,所以此方程组有无穷多解.'
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 数值 计算 报告