数值分析 第四章 基于MATLAB的科学计算—解线性方程组的迭代法
数值分析综述报告

淮阴工学院《数值分析》考试──基于Matlab的方法综合应用报告班级:金融1121姓名:姚婷婷学号:1124104129成绩:数理学院2014年6月7日《数值分析》课程综述报告前言:数值分析也称计算方法,它与计算工具的发展密切相关。
数值分析是一门为科学计算提供必需的理论基础和有效、实用方法的数学课程,它的任务是研究求解各类数学问题的数值方法和有关的理论。
正文:第一章 近似计算与误差分析1、产生误差的原因:①模型误差;②观测误差;③截断误差;④舍入误差。
2、四则运算的误差: ①加减法运算()()()****x y x y δδδ±=+②乘法运算()()()***************xy x y xy xy xy x y x y y y x x x y x y y x δδδ-=-+-≤-+-⇒=+③ 除法运算:()()()()()**********************2**x x xy x y y y yy xy x y x y x yyy x x yy y x yy x y y x x y y δδδ--=-+-=-+-=+⎛⎫⇒≈⎪⎝⎭3、科学表示法、有效数字、近似值的精度 任何一个实数都可以表示成如下的形式:其中:是正整数,是整数,如果是数的近似值并且则称该近似值具有位有效数字(significant digit )。
此时,该近似值的相对误差为另一方面,若已知()()*1111021nr x a δ-≤+那么,()()***1112110.10211102r m n n m n x x x x a a a a δ----≤⨯=+≤即:*x 至少有n 位有效数字。
例: 3.141592653589793...π= 取其近似值如下: x*=3.14 x * =3.14159 x*=3.1415 x*=3.141**213100.314110.0016...0.005101022x x π--=⨯-=<=⨯=⨯**516100.314159110.0000026...0.000005101022x x π--=⨯-=<=⨯=⨯**314100.31415110.000092...0.0001101022x x π--=⨯-=<<⨯=⨯**213100.3141110.00059...0.001101022x x π--=⨯-=<<⨯=⨯第二章 线性方程组在科学计算中,问题的本身就是求解线性方程组,许多问题的求解需要最后也归结为线性方程组的求解,所以线性方程组的求解是科学计算中最常见的问题。
数值分析实验报告--解线性方程组的迭代法及其并行算法

disp('请注意:高斯-塞德尔迭代的结果没有达 到给定的精度,并且迭代次数已经超过最大迭 代次数max1,方程组的精确解jX和迭代向量X 如下: ') X=X';jX=jX' end end X=X';D,U,L,jX=jX'
高斯-塞德尔的输入为:
A=[10 2 3;2 10 1;3 1 10]; b=[1;1;2]; X0=[0 0 0]'; X=gsdddy(A,b,X0,inf, 0.001,100) A=[10 2 3;2 10 1;3 1 10]; 请注意:因为对角矩阵 D 非奇异,所以此方程组有解.
0.0301 0.0758 0.1834
8.心得体会:
这已经是第三次实验了, 或多或少我已经对 MATLAB 有了更多的了 解与深入的学习。通过这次实验我了解了雅可比迭代法和高斯- 塞德尔迭代法的基本思想,虽然我们不能熟练编出程序,但还是 能看明白的。运行起来也比较容易,让我跟好的了解迭代法的多 样性,使平常手算的题能得到很好的验证。通过这次实验让我对 MATLAB 又有了更深一层的认识,使我对这门课兴趣也更加浓厚。
运行雅可比迭代程序输入: A=[10
b=[1;1;2];X0=[0 0 0]'; X=jacdd(A,b,X0,inf,0.001,100)
2 3;2 10 1;3 1 10];
结果为:
k= 1 X=
0.1000 k= 2 X= 0.0200 k= 3 X= 0.0400 k= 4 X= 0.0276 k= 5 X= 0.0314 k= 6 X= 0.0294 k= 7 X= 0.0301 k= 8 X= 0.0297
6、 设计思想:先化简,把对角线的项提到左边,其它项
matlab 解线性方程组的迭代法

迭代过程本质上就是计算极限的过程,一般不能 得到精确解。
迭代法的优点是程序简单,适合于大型方程组求 解,但缺点是要判断迭代是否收敛和收敛速度问题 。 1. 雅可比(Jacobi(1804-1851))迭代法(简单迭代法) 2. 赛得尔 (Seidel (1821 - 1896))迭代法
2、简单迭代法
while(norm(x-x1)>eps) x1=x; x=(I-A)*x1+b; n = n + 1; if(n>=M) disp('Warning: 迭代次数太多,现
在退出!'); return;
end end
例:求解方程组
clear all; A =[ 1.0170 -0.0092 0.0095;
遗传算法是一种基于自然选择的用于求解有约束和无约束 最优问题的方法。遗传算法反复修改包含若干个体的种群 。遗传算法在每一步中,随机从当前种群中选择若干个个 体作为父辈,并用它们产生下一代子辈。在若干代之后, 种群就朝着最优解“进化”。我们可以利用遗传算法去解决 各种最优化问题,包括目标函数是不连续、不可微、随机 或者高度非线性的问题。
若不满足收敛条件,适当调整方程次序或作一 定的线性组合,就可能满足收敛条件。
5、MATLAB的线性方程组求解函数 2
格式
solve('eqn1','eqn2',...,'eqnN','var1,var2,...,varN')
基于Matlab的解线性方程组的几种迭代法的实现及比较

基于Matlab的解线性方程组的几种迭代法的实现及比较线性方程组的解法有很多种,其中一类常用的方法是迭代法。
迭代法根据一个初值逐步逼近方程组的解,在每一次迭代中利用现有的信息产生新的近似值,并不断地修正。
下面介绍基于Matlab的三种迭代法:雅可比迭代法、高斯-赛德尔迭代法和超松弛迭代法,并进行比较。
1. 雅可比迭代法雅可比迭代法是迭代法中最简单的一种方法。
对于线性方程组Ax=b,雅可比迭代法的迭代公式为:x_{i+1}(j)=1/a_{jj}(b_j-\\sum_{k=1,k\eq j}^n a_{jk}x_i(k))其中,i表示迭代次数,j表示未知数的下标,x_i表示第i次迭代的近似解,a_{jk}表示系数矩阵A的第j行第k列元素,b_j 表示方程组的常数项第j项。
在Matlab中,可以使用以下代码实现雅可比迭代:function [x,flag]=jacobi(A,b,X0,tol,kmax)n=length(b);x=X0;for k=1:kmaxfor i=1:nx(i)=(b(i)-A(i,:)*x+A(i,i)*x(i))/A(i,i);endif norm(A*x-b)<tolflag=1;returnendendflag=0;return其中,参数A为系数矩阵,b为常数项列向量,X0为初值列向量,tol为迭代误差容许值(默认为1e-6),kmax为最大迭代次数(默认为1000)。
函数返回值x为近似解列向量,flag表示是否满足容许误差要求。
2. 高斯-赛德尔迭代法高斯-赛德尔迭代法是雅可比迭代法的改进。
其基本思想是,每次迭代时,利用已经求出的新解中的信息来更新其他未知数的值。
迭代公式为:x_{i+1}(j)=(1/a_{jj})(b_j-\\sum_{k=1}^{j-1}a_{jk}x_{i+1}(k)-\\sum_{k=j+1}^n a_{jk}x_i(k))与雅可比迭代法相比,高斯-赛德尔迭代法的每一次迭代都利用了前面已求得的近似解,因此可以更快地收敛。
matlab超松弛迭代法求方程组

一、介绍MATLAB(Matrix Laboratory)是一种用于数值计算和数据可视化的专业软件。
在MATLAB中,超松弛迭代法是解决线性方程组的一种有效算法。
本文将介绍MATLAB中超松弛迭代法的基本原理和实现方法,并给出一个具体的例子进行演示。
二、超松弛迭代法的基本原理超松弛迭代法是一种逐步迭代的算法,用于求解线性方程组。
它的基本原理是通过不断迭代更新方程组的解,直到达到满足精度要求的解。
超松弛迭代法的公式如下:X(k+1) = (1-w)X(k) + w*(D-L)⁻¹*(b+U*X(k))其中,X(k)代表第k次迭代的解向量,X(k+1)代表第k+1次迭代的解向量,D、L和U分别代表方程组的对角线元素、下三角元素和上三角元素构成的矩阵,b代表方程组的右端向量,w代表松弛因子。
超松弛迭代法的关键在于选择合适的松弛因子w,一般情况下,可以通过试验选取一个合适的值。
在MATLAB中,可以使用sor函数来实现超松弛迭代法。
三、MATLAB中超松弛迭代法的实现方法在MATLAB中,可以通过调用sor函数来实现超松弛迭代法。
sor 函数的语法格式如下:[X,flag,relres,iter,resvec] = sor(A,b,w,tol,maxit)其中,A代表线性方程组的系数矩阵,b代表右端向量,w代表松弛因子,tol代表迭代的精度要求,maxit代表最大迭代次数,X代表迭代求解得到的解向量,flag代表迭代的结果标志,relres代表相对残差的大小,iter代表迭代次数,resvec代表迭代过程中的残差向量。
以下是一个使用sor函数求解线性方程组的示例:A = [4 -1 0 -1 0 0; -1 4 -1 0 -1 0; 0 -1 4 0 0 -1; -1 0 0 4 -1 0; 0 -1 0 -1 4 -1; 0 0 -1 0 -1 4];b = [1; 0; -1; 0; 1; 0];w = 1.25;tol = 1e-6;maxit = 100;[X,flag,relres,iter,resvec] = sor(A,b,w,tol,maxit);通过调用sor函数,可以得到方程组的解向量X,迭代的结果标志flag,相对残余resrel和迭代次数iter。
MATLAB代码解线性方程组的迭代法

MATLAB代码解线性方程组的迭代法解线性方程组的迭代法1.rs里查森迭代法求线性方程组Ax=b的解function[x,n]=rs(A,b,x0,eps,M)if(nargin==3)eps=1.0e-6;%eps表示迭代精度M=10000;%M表示迭代步数的限制值elseif(nargin==4) M=10000;endI=eye(size(A));n=0;x=x0;tol=1;%迭代过程while(tol>eps)x=(I-A)*x0+b;n=n+1;%n为最终求出解时的迭代步数tol=norm(x-x0); x0=x;if(n>=M)disp('Warning:迭代次数太多,可能不收敛!'); return;endend2.crs里查森参数迭代法求线性方程组Ax=b的解function[x,n]=crs(A,b,x0,w,eps,M)if(nargin==4)eps=1.0e-6;%eps表示迭代精度M=10000;%M表示迭代步数的限制值elseif(nargin==5)M=10000;endI=eye(size(A));n=0;x=x0;tol=1;%迭代过程while(tol>eps)x=(I-w*A)*x0+w*b;n=n+1;%n为最终求出解时的迭代步数tol=norm(x-x0); x0=x;if(n>=M)disp('Warning:迭代次数太多,可能不收敛!'); return;endend3.grs里查森迭代法求线性方程组Ax=b的解function[x,n]=grs(A,b,x0,W,eps,M)if(nargin==4)eps=1.0e-6;%eps表示迭代精度M=10000;%M表示迭代步数的限制值elseif(nargin==5)M=10000;endI=eye(size(A));n=0;x=x0;tol=1;%前后两次迭代结果误差%迭代过程while(tol>eps)x=(I-W*A)*x0+W*b;%迭代公式n=n+1;%n为最终求出解时的迭代步数tol=norm(x-x0); x0=x;if(n>=M)disp('Warning:迭代次数太多,可能不收敛!'); return;endend4.jacobi雅可比迭代法求线性方程组Ax=b的解function[x,n]=jacobi(A,b,x0,eps,varargin)if nargin==3eps=1.0e-6;M=200;elseif nargin<3errorreturnelseif nargin==5M=varargin{1};endD=diag(diag(A));%求A的对角矩阵L=-tril(A,-1);%求A的下三角阵U=-triu(A,1);%求A的上三角阵B=D\(L+U);f=D\b;x=B*x0+f;n=1;%迭代次数while norm(x-x0)>=epsx0=x;x=B*x0+f;n=n+1;if(n>=M)disp('Warning:迭代次数太多,可能不收敛!');return;endend5.gauseidel高斯-赛德尔迭代法求线性方程组Ax=b的解function[x,n]=gauseidel(A,b,x0,eps,M)if nargin==3eps=1.0e-6;M=200;elseif nargin==4M=200;elseif nargin<3errorreturn;endD=diag(diag(A));%求A的对角矩阵L=-tril(A,-1);%求A的下三角阵U=-triu(A,1);%求A的上三角阵G=(D-L)\U;f=(D-L)\b;x=G*x0+f;n=1;%迭代次数while norm(x-x0)>=epsx0=x;x=G*x0+f;n=n+1;if(n>=M)disp('Warning:迭代次数太多,可能不收敛!');return;endend6.SOR超松弛迭代法求线性方程组Ax=b的解function[x,n]=SOR(A,b,x0,w,eps,M)if nargin==4eps=1.0e-6;M=200;elseif nargin<4errorreturnelseif nargin==5M=200;endif(w<=0||w>=2)error;return;endD=diag(diag(A));%求A的对角矩阵L=-tril(A,-1);%求A的下三角阵U=-triu(A,1);%求A的上三角阵B=inv(D-L*w)*((1-w)*D+w*U);f=w*inv((D-L*w))*b;x=B*x0+f;n=1;%迭代次数while norm(x-x0)>=epsx0=x;x=B*x0+f;n=n+1;if(n>=M)disp('Warning:迭代次数太多,可能不收敛!');return;endend7.SSOR对称逐次超松弛迭代法求线性方程组Ax=b的解function[x,n]=SSOR(A,b,x0,w,eps,M)if nargin==4eps=1.0e-6;M=200;elseif nargin<4errorreturnelseif nargin==5M=200;endif(w<=0||w>=2)error;return;endD=diag(diag(A));%求A的对角矩阵L=-tril(A,-1);%求A的下三角阵U=-triu(A,1);%求A的上三角阵B1=inv(D-L*w)*((1-w)*D+w*U);B2=inv(D-U*w)*((1-w)*D+w*L);f1=w*inv((D-L*w))*b;f2=w*inv((D-U*w))*b;x12=B1*x0+f1;x=B2*x12+f2;n=1;%迭代次数while norm(x-x0)>=epsx0=x;x12=B1*x0+f1;x=B2*x12+f2;n=n+1;if(n>=M)disp('Warning:迭代次数太多,可能不收敛!');return;endend8.JOR雅可比超松弛迭代法求线性方程组Ax=b的解function[x,n]=JOR(A,b,x0,w,eps,M)if nargin==4eps=1.0e-6;M=10000;elseif nargin==5M=10000;endif(w<=0||w>=2)%收敛条件要求error;return;endD=diag(diag(A));%求A的对角矩阵B=w*inv(D);%迭代过程x=x0;n=0;%迭代次数tol=1;%迭代过程while tol>=epsx=x0-B*(A*x0-b);n=n+1;tol=norm(x-x0);x0=x;if(n>=M)disp('Warning:迭代次数太多,可能不收敛!'); return;endend9.twostep两步迭代法求线性方程组Ax=b的解function[x,n]=twostep(A,b,x0,eps,varargin) if nargin==3eps=1.0e-6;M=200;elseif nargin<3errorreturnelseif nargin==5M=varargin{1};endD=diag(diag(A));%求A的对角矩阵L=-tril(A,-1);%求A的下三角阵U=-triu(A,1);%求A的上三角阵B1=(D-L)\U;B2=(D-U)\L;f1=(D-L)\b;f2=(D-U)\b;x12=B1*x0+f1;x=B2*x12+f2;n=1;%迭代次数while norm(x-x0)>=epsx0=x;x12=B1*x0+f1;x=B2*x12+f2;n=n+1;if(n>=M)disp('Warning:迭代次数太多,可能不收敛!');return;endend10.fastdown最速下降法求线性方程组Ax=b的解function[x,n]=fastdown(A,b,x0,eps)if(nargin==3)eps=1.0e-6;endx=x0;n=0;tol=1;while(tol>eps)%以下过程可参考算法流程r=b-A*x0;d=dot(r,r)/dot(A*r,r);x=x0+d*r;tol=norm(x-x0);x0=x;n=n+1;end11.conjgrad共轭梯度法求线性方程组Ax=b的解function[x,n]=conjgrad(A,b,x0)r1=b-A*x0;p=r1;n=0;for i=1:rank(A)%以下过程可参考算法流程if(dot(p,A*p)< 1.0e-50)%循环结束条件break;endalpha=dot(r1,r1)/dot(p,A*p);x=x0+alpha*p;r2=r1-alpha*A*p;if(r2< 1.0e-50)%循环结束条件break;endbelta=dot(r2,r2)/dot(r1,r1);p=r2+belta*p;n=n+1;end12.preconjgrad预处理共轭梯度法求线性方程组Ax=b的解function[x,n]=preconjgrad(A,b,x0,M,eps)if nargin==4eps=1.0e-6;endr1=b-A*x0;iM=inv(M);z1=iM*r1;p=z1;n=0;tol=1;while tol>=epsalpha=dot(r1,z1)/dot(p,A*p);x=x0+alpha*p;r2=r1-alpha*A*p;z2=iM*r2;belta=dot(r2,z2)/dot(r1,z1);p=z2+belta*p;n=n+1;tol=norm(x-x0);x0=x;%更新迭代值r1=r2;z1=z2;end13.BJ块雅克比迭代法求线性方程组Ax=b的解function[x,N]=BJ(A,b,x0,d,eps,M)if nargin==4eps=1.0e-6;M=10000;elseif nargin<4errorreturnelseif nargin==5M=10000;%参数的默认值endNS=size(A);n=NS(1,1);if(sum(d)~=n)disp('分块错误!');return;endbnum=length(d);bs=ones(bnum,1);for i=1:(bnum-1)bs(i+1,1)=sum(d(1:i))+1;%获得对角线上每个分块矩阵元素索引的起始值endDB=zeros(n,n);for i=1:bnumendb=bs(i,1)+d(i,1)-1;DB(bs(i,1):endb,bs(i,1):endb)=A(bs(i,1):endb,bs(i,1):endb);%求A的对角分块矩阵endfor i=1:bnumendb=bs(i,1)+d(i,1)-1;invDB(bs(i,1):endb,bs(i,1):endb)=inv(DB(bs(i,1):endb,bs(i,1):e ndb));%求A的对角分块矩阵的逆矩阵endN=0;tol=1;while tol>=epsx=invDB*(DB-A)*x0+invDB*b;%由于LB+DB=DB-AN=N+1;%迭代步数tol=norm(x-x0);%前后两步迭代结果的误差x0=x;if(N>=M)disp('Warning:迭代次数太多,可能不收敛!');return;endend14.BGS块高斯-赛德尔迭代法求线性方程组Ax=b的解function[x,N]=BGS(A,b,x0,d,eps,M)if nargin==4eps=1.0e-6;M=10000;elseif nargin<4errorreturnelseif nargin==5M=10000;endNS=size(A);n=NS(1,1);bnum=length(d);bs=ones(bnum,1);for i=1:(bnum-1)bs(i+1,1)=sum(d(1:i))+1;%获得对角线上每个分块矩阵元素索引的起始值endDB=zeros(n,n);for i=1:bnumendb=bs(i,1)+d(i,1)-1;DB(bs(i,1):endb,bs(i,1):endb)=A(bs(i,1):endb,bs(i,1):endb); %求A的对角分块矩阵endLB=-tril(A-DB);%求A的下三角分块阵UB=-triu(A-DB);%求A的上三角分块阵N=0;tol=1;while tol>=epsinvDL=inv(DB-LB);x=invDL*UB*x0+invDL*b;%块迭代公式N=N+1;tol=norm(x-x0);x0=x;if(N>=M)disp('Warning:迭代次数太多,可能不收敛!');return;end15.BSOR块逐次超松弛迭代法求线性方程组Ax=b的解function[x,N]=BSOR(A,b,x0,d,w,eps,M)if nargin==5eps=1.0e-6;M=10000;elseif nargin<5errorreturnelseif nargin==6M=10000;%参数默认值endNS=size(A);n=NS(1,1);bnum=length(d);bs=ones(bnum,1);for i=1:(bnum-1)bs(i+1,1)=sum(d(1:i))+1;%获得对角线上每个分块矩阵元素索引的起始值endDB=zeros(n,n);for i=1:bnumendb=bs(i,1)+d(i,1)-1;DB(bs(i,1):endb,bs(i,1):endb)=A(bs(i,1):endb,bs(i,1):endb); %求A的对角矩阵endLB=-tril(A-DB);%求A的下三角阵UB=-triu(A-DB);%求A的上三角阵N=0;iw=1-w;while tol>=epsinvDL=inv(DB-w*LB);x=invDL*(iw*DB+w*UB)*x0+w*invDL*b;%块迭代公式N=N+1;tol=norm(x-x0);x0=x;if(N>=M)disp('Warning:迭代次数太多,可能不收敛!'); return;endend。
基于matlab的线性方程组迭代法(实验报告)

基于matlab 的线性方程组迭代法实验题目:实验要求:(1)分别试用 Jacobi 和Gauss-Seidel 迭代法计算,要求达到的精度为:(1)()510k k x x +-∞->(2)观测得到的迭代序列是否收敛?若收敛,记录迭代次数并分析计算结果。
实验流程一、迭代法简介 1、 Jacobi 迭代法对于方程组Ax b =有A 非奇异情况下且0ij a ≠时,A 分裂为A D L U =--,可得到:0x B x f =+,其中1110(),B I D A D L U f D b ---=-=+=,得到雅克比迭代法:(0)(1)()0()k k x xB x f +⎧⎪⎨=+⎪⎩初始向量 2、 Gauss-Seidel 迭代法(0)(1)()()k k x x Gx f +⎧⎪⎨=+⎪⎩初始向量 其中11(),()G D L U f D L b --=-=-。
其迭代法优点为只需一组存储单元。
3、 超松弛迭代法(SOR)Gauss-Seidel 迭代法的一种加速方法,ω松弛因子。
(0)(1)()(1)(1))()(1)k k k k k x x Gx f x x x ωω+++⎧⎪⎪=+⎨⎪=+-⎪⎩(初始向量 其中11(),()G D L U f D L b --=-=-。
二、迭代法的matlab 程序1、 Jacobi 迭代法Jacobi.mfunction [y,n]= Jacobi( A,b,x0,e )%JACOBI ÇëÔÚ´Ë´¦ÊäÈ뺯Êý¸ÅÒªif(nargin<4)e=1e-5;endD=diag(diag(A));I=eye(size(A));B=I-D\A;f=D\b;y=x0+2*e;n=0;while norm(y-x0,inf)>ey=x0;x0=B*y+f;n=n+1;endnend2、Gauss-Seidel迭代法GaussSeidel.mfunction [y,n]= GaussSeidel( A,b,x0,e ) %GS ÇëÔÚ´Ë´¦ÊäÈ뺯Êý¸ÅÒªif(nargin<4)e=1e-5;endD=diag(diag(A));I=eye(size(A));L=D-tril(A);U=D-triu(A);f=(D-L)\b;G=(D-L)\U;y=x0+2*e;n=0;while norm(y-x0,inf)>ey=x0;x0=G*y+f;n=n+1;endnend3、超松弛迭代法(SOR) SOR.mfunction [y,n]= SOR( A,b,w,x0,e )%SORÇëÔÚ´Ë´¦ÊäÈ뺯Êý¸ÅÒªif(nargin<5)e=1e-5;endD=diag(diag(A));I=eye(size(A));L=D-tril(A);U=D-triu(A);f=(D-L)\b;G=(D-L)\U;y=x0+2*e;n=0;while norm(y-x0,inf)>ex0=y;x1=G*x0+f;y=(1-w)*x0+w*x1;n=n+1;endnend4、变量初始化creatMatrix.mclear;clc;a=diag(3*ones(1,20));b=diag(-0.5*ones(1,19),1);c=diag(-0.25*ones(1,18),2);A=a+b+b'+c+c';%ϵÊý¾ØÕób=ones(20,1)*7/4;b(1)=9/4;b(20)=9/4;x0=zeros(20,1);A,b,x0,w=1.5建立A数组以及初始化b,松弛因子w,迭代初值x05、程序运行和结果记录solve.mclc;tic,s1=Jacobi(A,b,x0),toctic,s2=GaussSeidel(A,b,x0),toctic,s3=SOR(A,b,w,x0),toc三、计算结果运行程序得到几种方法的计算结果。
数值分析MATLAB科学计算—线性方程组

科学计算—理论、方法及其基于MATLAB的程序实现与分析 三、 解线性方程组(线性矩阵方程)解线性方程组是科学计算中最常见的问题。
所说的“最常见”有两方面的含义:1) 问题的本身是求解线性方程组;2) 许多问题的求解需要或归结为线性方程组的求解。
关于线性方程组B A x B Ax 1-=⇒=(1)其求解方法有两类:1) 直接法:高斯消去法(Gaussian Elimination ); 2) 间接法:各种迭代法(Iteration )。
1、高斯消去法1) 引例考虑如下(梯形)线性方程组:()⎪⎩⎪⎨⎧==+==+-=⇒⎪⎪⎪⎭⎫ ⎝⎛=⎪⎪⎪⎭⎫ ⎝⎛⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡--⇔⎪⎩⎪⎨⎧==-=+-5.0141315.3221122004301211214322332321321332321x x x x x x x x x x x x x x x 高斯消去法的求解思路:把一般的线性方程组(1)化成(上或下)梯形的形式。
2)高斯消去法——示例考虑如下线性方程组:⎪⎪⎪⎭⎫ ⎝⎛-=⎪⎪⎪⎭⎫ ⎝⎛⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡---⇔⎪⎩⎪⎨⎧=++-=-+-=+-306015129101.2001.221113060129501.2001.221321321321321x x x x x x x x x x x x 1) 第一个方程的两端乘12加到第二个方程的两端,第一个方程的两端乘-1加到第三个方程的两端,得⎪⎪⎪⎭⎫⎝⎛-=⎪⎪⎪⎭⎫ ⎝⎛⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡--3060031110001.0001.00111321x x x2) 第二个方程的两端乘001.010-加到第三个方程的两端,得 ⎪⎪⎪⎭⎫⎝⎛-=⎪⎪⎪⎭⎫ ⎝⎛⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡--60600311010001.0001.00111321x x x3) 从上述方程组的第三个方程依此求解,得()⎪⎩⎪⎨⎧==+-==+-=600300001.03100024011332321x x x x x x 3)高斯消去法的不足及其改进——高斯(全、列)主元素消去法在上例中,由于建模、计算等原因,系数2.001而产生0.0005的误差,实际求解的方程组为⎪⎪⎪⎭⎫⎝⎛-=⎪⎪⎪⎭⎫ ⎝⎛⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡---306015129101.20005.22111321x x x ⎪⎩⎪⎨⎧===⇒70.4509.30142.2565321x x x注:数值稳定的算法高斯列主元素消去法就是在消元的每一步选取(列)主元素—一列中绝对值最大的元取做主元素,高斯列主元素消去法是数值稳定的方法。
- 1、下载文档前请自行甄别文档内容的完整性,平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益,请联系客服反馈,我们会尽快为您处理(人工客服工作时间:9:00-18:30)。
科学计算—理论、方法及其基于MATLAB 的程序实现与分析三、 解线性方程组的迭代法(Iteration )线性方程组的理论求解公式bA x 1-= (1)在应用于实际问题的计算时,通常面临两方面的问题 1、计算过程复杂, 2、不能保证算法的稳定性;此外,当初始数据(可能)存在误差时,按公式(1)即使求出了“精确解”意义也不大,因此,对于存在初始数据误差、特别是大型的线性方程组求解,需要寻求能达到精度要求的、操作和计算过程相对简单的求解方法。
下面将要介绍的迭代法就属于这类方法。
迭代法求解线性方程组的基本思想是1) 不追求“一下子”得到方程组的解,而是在逐步逼近方程组的精确解的迭代过程中获得满足精度要求的近似解,这一点与直接法不同;2) 通过对问题的转化,避免(困难的)矩阵求逆运算。
用迭代法求解线性方程组,首先要把线性方程组写成等价的形式f Mx x b Ax +=⇔= (2)式(2)的右端称为迭代格式,由迭代格式(2)确定如下的迭代算法:,2,101=⎩⎨⎧∈∀+=+k R x fMx x nk k (3)对于给定的线性方程组,可以写成不同的(无穷多)迭代格式,有意义的(可用的)迭代格式应具有收敛性―生成的解向量序列{}x n 收敛于方程组的解;而好的迭代法应具有较高的收敛速度。
关于迭代法收敛性的两个判别条件:a 、充分必要条件是:矩阵M 的谱半径(){}1,,2,1max <==n i M iiλρb 、充分条件是:矩阵M 的某个算子范数M <1。
设x 是方程组(2)的解,{}m x 是迭代法(3)生成的任一序列,因为f Mx x +=,f Mx x m m +=+1所以()()()0221x x Mx x Mx x M x x mm m m -==-=-=--- (4)设11--=⇔=TJTM J MTT ,其中矩阵J 是矩阵M 的Jordan 标准型,那么容易验证1-=TTJ M m m,并且()[]()()()1lim ,2,10lim lim 0lim lim 010<⇔=⇔==⇔-⇔=-⇔=+∞→+∞→-+∞→+∞→+∞→M Mn i x x TJT x x Mx x mm mi n mm mm m m ρρλ (5)此外,因为()02211x x Mx x M x x M x x M x x m m m m m -≤≤-≤-≤-=---- (6)所以x x x x MM m m m m mm =⇔=-⇒=⇒<+∞→+∞→+∞→lim 0lim0lim1(7)注:迭代格式(2)所确定的迭代法收敛与否,完全由系数矩阵M 决定,而与常数项f 无关.常用的迭代法1、Jacobian 迭代法:UL D a a a a a a a a a a a a a a a a a a A n n n nn n nn nn n n nn --=⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡----⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡----⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡=⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡=--0000000000000111211212211212222111211()()⎩⎨⎧=+=⇒++=⇒⎩⎨⎧--==----+b D f U L D M b D x U L D x U L D A b Ax m m 11111 (8)例1 解下面方程组(精确解为T x )1,1,1(*=).⎪⎩⎪⎨⎧=++-=+-=++.14103,53102,14310321321321x x x x x x x x x解 1) 改写成等价形式⎪⎪⎪⎩⎪⎪⎪⎨⎧--=++=----=--=).314(101),325(101)325(101),314(10121321312321x x x x x x x x x x x 2) 构造迭代公式,即为雅可比迭代公式⎪⎪⎪⎩⎪⎪⎪⎨⎧=--=++=--=+++.,2,1,0),314(101),325(101),314(101)(2)(1)1(3)(2)(1)1(2)(3)(2)1(1 k x x x x x x x x x k k k k k k k k k 3) 取初始向量T x )0,0,0()0(=,即,0)0(3)0(2)0(1===x x x 代入上式,求出4.11014,5.0105,4.11014)1(3)1(2)1(1=====x x x .再代回公式中,求出11.1)4.15.0314(101)2(1=-⨯-=x , 2.1)5.034.125(101)2(2=⨯+⨯+=x ,11.1)2.1311.114(101)2(3=⨯--=x .依次迭代,计算结果如表4-1.Example intera_j.mitera_j2、Gauss-Seidel 迭代法:()()()()⎩⎨⎧-=-=⇒-+-=⇒⎩⎨⎧--==----+b L D f U L D M b L D Ux L D x U L D A b Ax n n 11111 (9)根据GS 迭代法(9),可进一步得到()()()bD Ux Lx Dx b Ux x L D bL D Ux L D x n n n n n n n 11111111)(-+-++--+++=⇔+=-⇔-+-= (10)即⎪⎪⎭⎪⎪⎬⎫⎪⎪⎩⎪⎪⎨⎧+⎪⎪⎪⎪⎪⎭⎫⎝⎛⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡---+⎪⎪⎪⎪⎪⎭⎫ ⎝⎛⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡---=⎪⎪⎪⎪⎪⎭⎫⎝⎛-+++--+++b x x x a a a x x x a a a D x x x Nnn n n n n N n n n nn n N n n n2111121211111211121110000000000(11)式(11)表明:Gauss-Seidel 迭代法在计算第k 个迭代值k n x 1+时,及时地利用了在此步迭代中得到的新的迭代值:in x1+,1,,2,1-=k i ,由于第1+n 步的迭代值通常比第n 步的迭代值更接近方程组的精确解,所以,在Jacobian 迭代法和Gauss-Seidel 迭代法都收敛的情况下,Gauss-Seidel 迭代法的收敛速度比Jacobian 迭代法的收敛速度快。
例2 解下面方程组(与例1相同,精确解为T x )1,1,1(*=).⎪⎩⎪⎨⎧=++-=+-=++.14103,53102,14310321321321x x x x x x x x x 解 1) 原方程组改为等价方程组⎪⎪⎪⎩⎪⎪⎪⎨⎧--=++=--=).314(101),325(101),314(101213212321x x x x x x x x x . 2) 构造迭代公式,即为高斯-赛德尔迭代公式⎪⎪⎪⎩⎪⎪⎪⎨⎧=--=++=--=++++++.,2,1,0),314(101),325(101),314(101)1(2)1(1)1(3)(2)1(1)1(2)(3)(2)1(1k x x x x x x x x x k k k k k k k k k . 3) 取初始向量T x )0,0,0()0(=,即,0)0(3)0(2)0(1===x x x 代入上式,求出4.1)00314(101)1(1=-⨯-=x , 78.0)034.125(101)1(2=⨯+⨯+=x , 026.1)78.034.114(101)1(3=⨯--=x .迭代计算下去,得表4-2.Example itera_gs.mitera_gs对Gauss-Seidel 迭代法做进一步的研究,式(9)还可以写成如下的形式:()()()()()n n Jn SG n n n n n n n n n n n n n n x x L Dx x x x L Db D x U L D x b D Lx Ux Lx Lx D x b Ux x L D bL D Ux L D x -+=⇔-+++=⇔+-++=⇔+=-⇔-+-=+-++---+-+-++--+11_11111111111111)()( (12)即Gauss-Seidel 迭代法是在Jacobian 迭代法的基础上增加了一个修正项:()n n x x L D-+-11并且这个修正项起到了加速的作用,受这一事实的启发,为获得更快的收敛速度,我们对Gauss-Seidel 迭代法再做进一步的研究()()()()()()b x D U Lx Dx x b x D U Lx Dx Dx bUx x L D bL D Ux L D x n n n n n n n n n n n n +-++=⇔+-++=⇔+=-⇔-+-=+-++++--+111111111 (13)式(13)表明Gauss-Seidel 迭代法的第1+n 步迭代值是在第n 步迭代值的基础上增加了一个修正项:()()b x D U Lx Dn n +-++-11(14)并且这个修正项使第1+n 步迭代值更接近方程组的精确解,受这一事实的启发,我们希望通过用一个适当的因子乘修正项(14)的办法达到获得更快收敛速度的目的:()()()()()()[]()()[]()bwLD w x wU Dw wLD x wbx wU D w x wL D b x D U Lx w Dx Dx b x D U Lx wD x x n n n n n n n n n n n n 11111111111--+++++-+-++--=⇔++-=-⇔+-++=⇔+-++= (15)从而得到3、SOR (Successive Over Relaxation Method ) 迭代法:()()[]()x D w L w D w U x w D w L b n n+--=--++-1111 (16)其中w 称为松弛因子,由于()()[]()()()⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡------⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡=+--=----nn n n n nn nn n a w wa a w wa wa a w a wawa a wa a wUDw wL D M 101010001122112111112221111(17) ()()()()∏∏∏====-=-=⇔nk knnk kknk kkn M w aa w M 11111det λ (18)()()1120111>∃⇔>⇔>∨<<⇔>-⇔∏=MMw w wk nk k nλλ (19) 所以,为保证SOR 迭代法的收敛性,要求20<<w 。
例3 用超松驰迭代法解下面方程组,取松驰因子4.1=ω.⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡=⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡------0101211210012100124321x x x x .解 方程组的精确解为:T x )8.0,6.1,4.1,2.1(=.取初值T x )1,1,1,1()0(=,用高斯-赛德尔迭代10次,得T x )798216.0,5964336.1,3955917.1,1966324.1()10(=.利用SOR 方法,构造迭代公式⎪⎪⎪⎪⎩⎪⎪⎪⎪⎨⎧-+=+-++=+-+=+-+=+++++++).2(2),21(2),2(2),21(2)(4)1(3)(4)1(4)(4)(3)1(2)(3)1(3)(3)(2)1(1)(2)1(2)(2)(1)(1)1(1k k k k k k k k k k k k k k k k k k x x x x x x x x x x x x x x x x x x ωωωω与高斯-赛德尔方法相同,初值为T x )1,1,1,1()0(=.迭代计算结果列于表4-3.表4-3 SOR 迭代的数值结果k)(1k x)(2k x)(3k x)(1k x1 2 3 4 51 11.343 1.195451.203472 11.49 1.4753 1.402361.4028735 1.7 1.616 1.65095 1.60198151.59390590.79 0.8152 0.829585 0.789553 0.7999129由表4-3可知,用超松驰迭代法只迭代了5次,结果与G-S 法迭代10次的结果大体相同,可见,SOR 方法的松驰因子起到了加速收敛的重要作用.Example itera_sor.mitera_sor关于迭代法收敛性的两个判别条件:a 、充分必要条件是:矩阵M 的谱半径(){}1,,2,1max <==n i M iiλρb 、充分条件是:矩阵M的某个算子范数M <1。