数值计算报告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，所以此方程组有无穷多解.'
升级会员 