整理LU分解法列主元高斯法Jacobi迭代法GaussSeidel法的原理及Matlab程序.docx
- 文档编号:14792213
- 上传时间:2023-06-27
- 格式:DOCX
- 页数:13
- 大小:118.72KB
整理LU分解法列主元高斯法Jacobi迭代法GaussSeidel法的原理及Matlab程序.docx
《整理LU分解法列主元高斯法Jacobi迭代法GaussSeidel法的原理及Matlab程序.docx》由会员分享,可在线阅读,更多相关《整理LU分解法列主元高斯法Jacobi迭代法GaussSeidel法的原理及Matlab程序.docx(13页珍藏版)》请在冰点文库上搜索。
整理LU分解法列主元高斯法Jacobi迭代法GaussSeidel法的原理及Matlab程序
精品文档
、实验目的及题目
1.1实验目的:
(1)学会用高斯列主元消去法,LU分解法,Jacobi迭代法和Gauss-Seide迭代法解线性
方程组。
(2)学会用Matlab编写各种方法求解线性方程组的程序。
1.2实验题目:
用列主元消去法解方程组:
iX”2+3x4=412为+X2-X3+&=1
12X|—x?
-X3+3x^=-3[-为+2X2+3X3-X4=4
1.
2.
用LU分解法解方程组
广48
-24
0-
12
2、
—24
24
12
,b=
4
0
6
20
2|
-2
<~6
6
2
16
、乂)
Ax=b,其中
A=
3.
分别用Jacobi迭代法和
Gauss-Seide迭代法求解方程组:
10%—X2+2x3=-11
8X2-X3+3x4=-11
2x^x2+10x3=6
一為+3x2—X3+11%=25
二、实验原理、程序框图、程序代码等
2.1实验原理
2.1.1高斯列主元消去法的原理
Gauss消去法的基本思想是一次用前面的方程消去后面的未知数,从而将方程组化为等
价形式:
》1必+Gx2+III+gXn=g
Ib22X^|(+b2nX^g2
bnnXn=gn
这个过程就是消元,然后再回代就好了。
具体过程如下:
对于k=1,2,川,n-1,若akJ工0,依次计算
精品文档
mik=aik)/akk)
ai(k—aj(k)-mikak:
)
b(5=b(k)-mkbkk),i,j=k+1,川,n
然后将其回代得到:
Xn尽)/^
n
}xk=(bkk)-2akjk)Xj)/akk),k=n—1,n—2,川,1
jM
以上是高斯消去。
但是高斯消去法在消元的过程中有可能会出现akk^0的情况,这时消元就无法进行了,
即使主元数akk^O,但是很小时,其做除数,也会导致其他元素数量级的严重增长和舍入误差的扩散。
因此,为了减少误差,每次消元选取系数矩阵的某列中绝对值最大的元素作为主元素。
然后换行使之变到主元位置上,再进行销元计算。
即高斯列主元消去法。
2.1.2直接三角分解法(LU分解)的原理
先将矩阵A直接分解为A=LU则求解方程组的问题就等价于求解两个三角形方程组。
直接利用矩阵乘法,得到矩阵的三角分解计算公式为:
|u1i=a1i,i二1,2,川,n
山=ai1/U11,i=2,ill,n
「kA
|Ukj=akj-送lkmUmj,,j=k,k+1,川,n
ml,k=2,3川In
lik=(ak—送limUmk)/Ukk,i=k+1,k+2」||,n且kHn
mzt
由上面的式子得到矩阵A的LU分解后,求解Ux=y的计算公式为
i-4
yi=D
IVi=bi—送lijyj,i=2,3,川n
[y
Xn=Vn/Unn
*n
Xi=(yi-2UijXj)/Uii,i=n—1,n—2,川,1j止十
以上为LU分解法。
精品文档
2.1.3Jacobi迭代法和Gauss-Seidel迭代法的原理
(1)Jcaobi迭代
设线性方程组
Ax=b
的系数矩阵A可逆且主对角元素a11,a22,...,ann均不为零,令
D=diag(an,a22,...,ann
并将A分解成
A=(A-D)中D
从而
(1)可写成
X=BiX+fi
-1
11
其中B=1-DA,f1=Db
以Bi为迭代矩阵的迭代法(公式)
称为雅可比Jacobi)迭代法,其分量形式为
(kdt)1r_J(k)
X=—Lb-壬aX八i^iiJj
a••J
u"J#
「=1,2,...n,k=01,2,...
其中x(J=(xFW),..人(0斤为初始向量.
(2)Gauss-Seidel迭代
由雅可比迭代公式可知,在迭代的每一步计算过程中是用X(k)的全部分量来计算x(k+的
(k+)(k+)(k+)
所有分量,显然在计算第i个分量Xi时,已经计算出的最新分量X1,...,Xi~没有被利用。
把矩阵A分解成
A=D-L-U
其中D二diag(a11,a22,...,ann\~L"U分别为A的主对角元除外的下三角和上三角部分,
精品文档
于是,方程组
(1)便可以写成
(D-L风=Ux+b
X=B2X+f2
其中
以B2为迭代矩阵构成的迭代法(公式)
称为高斯一塞德尔迭代法,用分量表示的形式为
functionGauss(A,b)%A为系数矩阵,b为右端项矩阵
[m,n]=size(A);
n=length(b);
fork=1:
n-1
[pt,p]=max(abs(A(k:
n,k)));%找出列中绝对值最大的数
p=p+k-1;
ifp>k
t=A(k,:
);A(k,:
)=A(p,:
);A(p,:
)=t;%交换行使之变到主元位置上
t=b(k);b(k)=b(p);b(p)=t;
A(k+1:
n,k+1:
n)=A(k+1:
n,k+1:
n)-m*A(k,k+1:
n);
ifflag~=0
精品文档
精品文档
Ab=[A,b];
end
end
fork=n-1:
-1:
1
x(k)=(b(k)-A(k,k+1:
n)*x(k+1:
n))/A(k,k);
end
fork=1:
n
fprintf('x[%d]=%f\n',k,x(k));
end
L(2:
n,1)=A(2:
n,1)/U(1,1);
fork=2:
n
U(k,k:
n)=A(k,k:
n)-L(k,1:
k-1)*U(1:
k-1,k:
n);
L(k+1:
n,k)=(A(k+1:
n,k)-L(k+1:
n,1:
k-1)*U(1:
k-1,k))/U(k,k);
end
%输出L矩阵
%输出U矩阵
y=zeros(n,1);%开始解方程组Ux=y
y
(1)=b
(1);
fork=2:
n
y(k)=b(k)-L(k,1:
k-1)*y(1:
k-1);
endx=zeros(n,1);精品文档
精品文档x(n)=y(n)/U(n,n);fork=n-1:
-1:
1
x(k)=(y(k)-U(k,k+1:
n)*x(k+1:
n))/U(k,k);endfork=1:
n
fprintf('x[%d]=%f\n',k,x(k));
end
2.2.3Jacobi迭代法的程序
[m,n]=size(A);
temp=1;
x=zeros(m,1);
k=0;
whileabs(max(x)-temp)>eps
temp=max(abs(x));
k=k+1;
%记录循环次数
%雅克比迭代公式
x=-inv(D)*(L+U)*x+inv(D)*b;
endfork=1:
n
fprintf('x[%d]=%f\n',k,x(k));
end
精品文档
2.2.4Gauss-Seidel迭代程序
functionGauss_Seidel(A,b,eps)%A为系数矩阵,b为后端项矩阵,epe为精度
[m,n]=size(A);
temp=1;
x=zeros(m,1);
k=0;
whileabs(max(x)-temp)>eps
temp=max(abs(x));
endfork=1:
n
fprintf('x[%d]=%f\n',k,x(k));
end
三、实验过程中需要记录的实验数据表格
3.1第一题(高斯列主元消去)的数据>>A=[1103;21-11;3-1-13;-123-1];>>b=[4;1;-3;4];
>>Gauss(A,b)
x[1]=-1.333333
x[2]=2.333333
x[3]=-0.333333
x[4]=1.0000003.2第二题(LU分解法)数据
>>A=[48-240-12;-24241212;06202;-66216];
>>b=[4;4;-2;-2];
>>LU(A,b)
精品文档
1.0000
-0.5000
0
-0.1250
0
1.0000
0.5000
0.2500
0
0
1.0000
-0.0714
0
0
0
1.0000
-24.0000
12.0000
0
0
-12.0000
6.0000
-1.0000
12.9286
x[1]=0.521179
x[2]=1.005525
x[3]=-0.375691
x[4]=-0.259669
>>A=[10-120;08-13;2-1100;-13-111];b=[-11;-11;6;25];
Jacobi(A,b,0.00005)
x[1]=-1.467396
x[2]=-2.358678
x[3]=0.657604
x[4]=2.8423973.4第三题用Gauss_Seide迭代的数据>>A=[10-120;08-13;2-1100;-13-111];
>>b=[-11;-11;6;25];
>>Gauss_Seidel(A,b,0.00005)
x[1]=-1.467357
x[2]=-2.358740
x[3]=0.657597
x[4]=2.842405
四、实验中存在的问题及解决方案
4.1存在的问题
(1)第一题中在matlab中输入》Gauss(A,b)(数据省略)得到m=4n=4?
?
?
Undefinedfunctionorvariable"Ab".Errorin==>Gaussat8[ap,p]=max(abs(Ab(k:
n,k)));没有得到想要的结
精品文档
果。
2)第二题中在matlab中输入>>y=LU(A,b)得到y=4.00006.0000-5.0000-3.3571不是方程组的解。
(3)第三题中在用高斯赛德尔方法时在matlab中输入>>Gauss-Seidel(A,b,eps结果程序报
错?
?
?
Errorusing==>GaussToomanyoutputargumentS得不至U想要的结果。
4.2解决方案
(1)针对第一题中由于程序的第二行漏了一个分号导致输出了m和n的值,第8行中将
Ab改为A问题就解决了。
(2)由于程序后面出现了矩阵y故输出的事矩阵y的值,但是我们要的事X的值,故只需
要将y改成X,或者直接把y去掉就解决了问题。
3)在function文件中命名不能出现“-”应该将其改为下划线“_”,所以将M文件名
Gauss-Seidel(A,b,eps”)改成“Gauss_Seidel(A,b,eps”)就解决问题了。
五、心得体会
本次试验涉及到了用高斯列主元消去法,LU分解法,Jacobi迭代法以及Gauss-seide迭
代法等四种方法。
需要对这些方法的原理都要掌握才能写出程序,由于理论知识的欠缺,我
花了很大一部分时间在看懂实验的原理上,看懂了实验原理之后就开始根据原理编写程序,程序中还是出现了很多的低级错误导致调试很久才能运行。
通过这次试验使我深刻的体会到理论知识的重要性,没有理论知识的支撑是写不出程序精品文档
精品文档
来的。
写程序时还会犯很多低级的错误,以后一定要加强理论知识的学习,减少编程时低级错误的产生。
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 整理 LU 解法 列主元高斯法 Jacobi 迭代法 GaussSeidel 原理 Matlab 程序
链接地址:https://www.bingdoc.com/p-14792213.html