整理LU分解法列主元高斯法Jacobi迭代法GaussSeidel法的原理与Matlab程序文档格式.docx
《整理LU分解法列主元高斯法Jacobi迭代法GaussSeidel法的原理与Matlab程序文档格式.docx》由会员分享,可在线阅读,更多相关《整理LU分解法列主元高斯法Jacobi迭代法GaussSeidel法的原理与Matlab程序文档格式.docx(12页珍藏版)》请在冰豆网上搜索。
然后换行使之变到主元位置上,再进行销元计算。
即高斯列主元消去法。
2.1.2直接三角分解法(LU分解)的原理
先将矩阵A直接分解为
则求解方程组的问题就等价于求解两个三角形方程组。
直接利用矩阵乘法,得到矩阵的三角分解计算公式为:
由上面的式子得到矩阵A的LU分解后,求解Ux=y的计算公式为
以上为LU分解法。
2.1.3Jacobi迭代法和Gauss-Seidel迭代法的原理
(1)Jcaobi迭代
设线性方程组
(1)
的系数矩阵A可逆且主对角元素
均不为零,令
并将A分解成
(2)
从而
(1)可写成
令
.(3)
以
为迭代矩阵的迭代法(公式)
(4)
称为雅可比(Jacobi)迭代法,其分量形式为
(5)
为初始向量.
(2)Gauss-Seidel迭代
由雅可比迭代公式可知,在迭代的每一步计算过程中是用
的全部分量来计算
的所有分量,显然在计算第i个分量
时,已经计算出的最新分量
没有被利用。
把矩阵A分解成
(6)其中
分别为
的主对角元除外的下三角和上三角部分,于是,方程组
(1)便可以写成
即
(7)
为迭代矩阵构成的迭代法(公式)
(8)
称为高斯—塞德尔迭代法,用分量表示的形式为
2.2程序代码
2.2.1高斯列主元的代码
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;
end
m=A(k+1:
n,k)/A(k,k);
%开始消元
A(k+1:
n,k+1:
n)=A(k+1:
n)-m*A(k,k+1:
n);
b(k+1:
n)=b(k+1:
n)-m*b(k);
n,k)=zeros(n-k,1);
ifflag~=0
Ab=[A,b];
end
x=zeros(n,1);
%开始回代
x(n)=b(n)/A(n,n);
fork=n-1:
-1:
1
x(k)=(b(k)-A(k,k+1:
n)*x(k+1:
n))/A(k,k);
fork=1:
n
fprintf('
x[%d]=%f\n'
k,x(k));
2.2.2LU分解法的程序
functionLU(A,b)%A为系数矩阵,b为右端项矩阵
%初始化矩阵A,b,L和U
L=eye(n,n);
U=zeros(n,n);
U(1,1:
n)=A(1,1:
%开始进行LU分解
L(2:
n,1)=A(2:
n,1)/U(1,1);
fork=2:
U(k,k:
n)=A(k,k:
n)-L(k,1:
k-1)*U(1:
k-1,k:
L(k+1:
n,k)=(A(k+1:
n,k)-L(k+1:
n,1:
k-1,k))/U(k,k);
L%输出L矩阵
U%输出U矩阵
y=zeros(n,1);
%开始解方程组Ux=y
y
(1)=b
(1);
y(k)=b(k)-L(k,1:
k-1)*y(1:
k-1);
x=zeros(n,1);
x(n)=y(n)/U(n,n);
fork=n-1:
x(k)=(y(k)-U(k,k+1:
n))/U(k,k);
2.2.3Jacobi迭代法的程序
functionJacobi(A,b,eps)%A为系数矩阵,b为后端项矩阵,epe为精度
D=diag(diag(A));
%求矩阵D
L=tril(A)-D;
%求矩阵L
U=triu(A)-D;
%求矩阵U
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;
%雅克比迭代公式
end
2.2.4Gauss-Seidel迭代程序
functionGauss_Seidel(A,b,eps)%A为系数矩阵,b为后端项矩阵,epe为精度
L=D-tril(A);
U=D-triu(A);
x=inv(D-L)*U*x+inv(D-L)*b;
%Gauss_Seidel的迭代公式
三、实验过程中需要记录的实验数据表格
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.000000
3.2第二题(LU分解法)数据
A=[48-240-12;
-24241212;
06202;
-66216];
4;
-2;
-2];
LU(A,b)
L=
1.0000000
-0.50001.000000
00.50001.00000
-0.12500.2500-0.07141.0000
U=
48.0000-24.00000-12.0000
012.000012.00006.0000
0014.0000-1.0000
00012.9286
x[1]=0.521179
x[2]=1.005525
x[3]=-0.375691
x[4]=-0.259669
3.3第三题Jacobi迭代法的数据
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.842397
3.4第三题用Gauss_Seidel迭代的数据
b=[-11;
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=4Undefinedfunctionorvariable"
Ab"
.Errorin==>
Gaussat8[ap,p]=max(abs(Ab(k:
没有得到想要的结果。
(2)第二题中在matlab中输入>
y=LU(A,b)得到y=4.00006.0000-5.0000-3.3571不是方程组的解。
(3)第三题中在用高斯赛德尔方法时在matlab中输入>
Gauss-Seidel(A,b,eps)结果程序报错Errorusing==>
GaussToomanyoutputarguments.得不到想要的结果。
(6)列出选定的评价方法,并作简单介绍。
4.2解决方案
(5)为保障评价对象建成或实施后能安全运行,应从评价对象的总图布置、功能分布、工艺流程、设施、设备、装置等方面提出安全技术对策措施;
从评价对象的组织机构设置、人员管理、物料管理、应急救援管理等方面提出安全管理对策措施;
从保证评价对象安全运行的需要提出其他安全对策措施。
对策措施的建议应有针对性、技术可行性和经济合理性,可分为应采纳和宜采纳两种类型。
(1)针对第一题中由于程序的第二行漏了一个分号导致输出了m和n的值,第8行中将Ab改为A问题就解决了。
(2)由于程序后面出现了矩阵y故输出的事矩阵y的值,但是我们要的事x的值,故只需要将y改成x,或者直接把y去掉就解决了问题。
(1)内涵资产定价法(3)在function文件中命名不能出现“-”应该将其改为下划线“_”,所以将M文件名“Gauss-Seidel(A,b,eps)”改成“Gauss_Seidel(A,b,eps)”就解决问题了。
(一)环境影响评价的概念
1.环境影响评价依据的环境标准体系
3)选择价值。
选择价值(OV)又称期权价值。
我们在利用环境资源的时候,并不希望它的功能很快消耗殆尽,也许会设想未来该资源的使用价值会更大。
五、心得体会
本次试验涉及到了用高斯列主元消去法,LU分解法,Jacobi迭代法以及Gauss-seidel迭代法等四种方法。
需要对这些方法的原理都要掌握才能写出程序,由于理论知识的欠缺,我花了很大一部分时间在看懂实验的原理上,看懂了实验原理之后就开始根据原理编写程序,程序中还是出现了很多的低级错误导致调试很久才能运行。
4.建设项目环境影响评价文件的分级审批通过这次试验使我深刻的体会到理论知识的重要性,没有理论知识的支撑是写不出程序来的。
写程序时还会犯很多低级的错误,以后一定要加强理论知识的学习,减少编程时低级错误的产生。