第7节 迭代法

1、第三章第三章 线性方程组求解的数值方法线性方程组求解的数值方法第四节第四节 迭代法迭代法迭代法（iterative method ） 迭代法是从迭代法是从初始猜测值初始猜测值开始，通过开始，通过依次逼近依次逼近获得问题解的获得问题解的方法方法。 应用：应用： 线性方程组求解、非线性方程求解、非线性方程组求解、特征值线性方程组求解、非线性方程求解、非线性方程组求解、特征值求解、最优化方法、分形求解、最优化方法、分形 数字信号处理：维纳滤波、卡拉曼滤波数字信号处理：维纳滤波、卡拉曼滤波 典型方法：典型方法： 数值计算：数值计算：Jacobi法、法、Gauss-Seidel法、法、Newton迭代法

2、、最快迭代法、最快下降法（下降法（Steepest Descent）分形 （fractal） 迭代法能从迭代法能从简单的规则简单的规则出发，通过迭代得到出发，通过迭代得到复杂的结果复杂的结果。最直观的展示就是最直观的展示就是“分形图形分形图形” 分形是一种粗糙的或分裂的几何形状，其每一部分都是整分形是一种粗糙的或分裂的几何形状，其每一部分都是整体的缩小比例的复制（自相似性）体的缩小比例的复制（自相似性）NOVA纪录片：纪录片： hunting the hidden dimension “寻找隐藏的维度寻找隐藏的维度” 分形应用：分形应用： 计算机虚拟现实技术；计算机虚拟现实技术； 天线设计；天

3、线设计； 医学；医学； 分形 fractal科赫雪花（科赫雪花（Koch snowflake）：）：1、初始为正三角形；、初始为正三角形；2、将每条边等分三段；、将每条边等分三段；3、以中段为底向外做等边三角形、以中段为底向外做等边三角形4、移除被三等分的边的中段。、移除被三等分的边的中段。5、如此迭代。、如此迭代。分形图形分形图形分形图形分形图形分形图形分形图形clear allclose allclcshgclf resetset(gcf, color, white, menubar, none, numbertitle, off, name, ?)x = -0.1; 0.5h = plo

4、t(x(1), x(2), .)mycolor = 0, 2/3,0;set(h, markersize, 10, color, mycolor, erasemode, none);axis(-3, 3, 0, 10);axis off% stop = uicontrol(style, toggle, string, stop, background, white)% drawnowp = 0.85, 0.92, 0.99, 1.00;A1 = 0.85, 0.04; -0.04, 0.85;b1 = 0; 1.6;A2 = 0.20, -0.26; 0.23, 0.22;b2 = 0; 1.

5、6;A3 = -0.15, 0.28; 0.26, 0.24;b3 = 0; 0.44;A4 = 0, 0; 0, 0.16;for iii = 1 : 10000 r = rand; if r p(1) x = A1 * x + b1; elseif r p(2) x = A2 * x + b2; elseif r 2 ()() +-()kkkkkkkkkkkkxKkKxxxxxxxxxxxxxxxf x，：，收敛则有:使得时，满足 即：收敛 又有：，则： ( )( )( )( )*()- lim()-0()0 ()kkkkkf xxf xxf xxf xx ， x*称为函数 f（x）的不动

6、点不动点（fixed point）证明：证明：迭代函数的构造 例：用迭代法求解方程 2 x = 1(1)( )1(1)( )2(1)( )3(1)( )4( )11( )(1)/3(1)/3( )2(1 1.5 )2(1 1.5)( )(1 10 )/12(1 10)/12kkkkkkkkxxxxxxxxxxxxxxxx ( )0( )( )0( )f xf xxxa f xa f xxx 步骤二、构造迭代函数： 下列方程（组）等价。 步骤一、写出等价方程：11 102112(1 1.5 )312xxxxxxxxx ( )0( )f xxx010203000.20.40.60.81f(x) =

7、 1 - x0102030-1-0.500.511.52x 108f(x) = 2 * (1 - 1.5 x)010203000.10.20.30.40.5f(x) = (1 + 10 x) / 12010203000.10.20.30.40.5f(x) = (1 + x) / 3振荡序列振荡序列收敛序列收敛序列发散序列发散序列收敛序列收敛序列 迭代法的关键问题： 迭代函数收不收敛迭代函数收不收敛 收敛速度快不快收敛速度快不快迭代法程序clear allclose allclcformat long% 迭代法 求解2x=1x0 = 0 % 迭代初值N = 30 % 迭代次数% 方法1：f(x)

8、 = 1 - xx1(1) = x0;for iii = 2 : N x1(iii) = 1 - x1(iii - 1);endfigure;plot(x1, o-)title(f(x) = 1 - x)% 方法2：f(x) = (1 + x) / 3x2(1) = x0;for iii = 2 : N x2(iii) = (1 + x2(iii - 1) / 3;endfigure;plot(x2, o-)title(f(x) = (1 + x) / 2)% 方法3：f(x) = 2 * (1 - 1.5 x)x3(1) = x0;for iii = 2 : N x3(iii) = (1 -

9、 2 * x3(iii - 1);endfigure;plot(x3, o-)title(f(x) = 2 * (1 - 1.5 x)% 方法4：f(x) = (1 + 10 x) / 12x4(1) = x0;for iii = 2 : N x4(iii) = (1 + 10 * x4(iii - 1) / 12;endfigure;plot(x4, o-)title(f(x) = (1 + 10 x) / 12)迭代法求解线性方程组 迭代公式的构建：迭代公式的构建： 将方程将方程Ax=b改写为：改写为：x=Mx+c，M称为迭代矩阵。称为迭代矩阵。迭代法求解线性方程组 迭代公式一：最简单的构

10、造迭代公式一：最简单的构造0()AxbAxbAI xbx(1)( )()kkxAI xbMAI 迭代公式二：最有效的构造迭代公式二：最有效的构造111, =A AxA bxxA bM00迭代法求解线性方程组 迭代公式三（Jacobi 迭代）：11221331111112211211222222211233221112 nnnnnnnnnnnnnngxb xb xb xa x a xa xbga x a xa xbxb xb xb xa x a xa xb112233 nnnnngxb x b xb x11000.000nnaDa定义对角矩阵：)11(1)( () kkMIJacD Aobixb

11、xggDM：迭代公式为迭代法求解线性方程组迭代次数迭代次数x的第一分量的第一分量x的第二分量的第二分量00.00.510.60.220.360.6830.7440.48840.59040.795250.83620.672360.73790.868970.89510.7903 例：利用例：利用Jacobi迭代求解方程组迭代求解方程组54,(1,1) .45TAb100.8(),0.80JMDLU 1(0.2,0.2)TJgD b 迭代法求解线性方程组 迭代公式三（迭代公式三（Gauss-Seidel迭代）：迭代）：一般认为新近似解要比老近似解更接近真实解，将已计算出的一般认为新近似解要比老近似解

12、更接近真实解，将已计算出的x（k+1）分量替换分量替换Jacobi 迭代公式中迭代公式中x（k）相应分量即可得到相应分量即可得到Gauss-Seidel迭代。迭代。迭代法求解线性方程组 例：利用例：利用Gauss-Seidel迭代求解方程组迭代求解方程组541,.451Ab 100.8(),0.80JMDLU 1(0.2,0.2)TJgD b 步骤步骤1、构造、构造Jacobi迭代公式迭代公式 步骤步骤2、选择初值、选择初值(0)0,0.5x 步骤步骤3、利用利用Jacobi迭代公式计算迭代公式计算一次迭代的第一分量：一次迭代的第一分量：(1)1000.80.20.60.5x迭代法求解线性方程

13、组 步骤步骤4、将步骤、将步骤3得到的得到的一次迭代的第一分量替换初值的第一一次迭代的第一分量替换初值的第一分量，计算一次迭代的第二分量：分量，计算一次迭代的第二分量：(1)20.800.20.680.50.6x 步骤步骤5、如果第三分量存在，如果第三分量存在，利用利用一次迭代的第一、二分量一次迭代的第一、二分量计算第三分量，直到计算出所有迭代向量分量；计算第三分量，直到计算出所有迭代向量分量； 步骤步骤6、重复步骤、重复步骤3-5，进行迭代，进行迭代；迭代法求解线性方程组迭代次数迭代次数x的第一分量的第一分量x的第二分量的第二分量00.00.510.60.6820.7440.79530.83

14、620.868940.89510.916150.93280.946360.95700.965670.97250.9780JacobiGauss与相比，只需一组工作单元存放近似解。迭代法求解线性方程组 Gauss-Seidel迭代矩阵：(1)( )( )( )112213311(1)( )( )22332(12)211(1(1)112 kkkknnkkkknnnnnkkgxb xb xb xgxb xxbbxxxb(1)(1)2,11 kkn nnnb xbxg1212,000U000000nnbbb21120000000 00nnLbbb(1)(1)( )(1)( )(1)1( )1 ()()

15、()kkkkkkkxLxUxgIL xUxgxILUxILg11()JacobixMxgLU xgIL xUxgxILUxILg迭代法求解线性方程组 Gauss-Seidel迭代矩阵：12131212323132112111111(1)1( 000000000 , () ()nnnnnnnnkkaaaaaaLaaUaaaaLD LUD UILD DD LDDLxILUx 如果用矩阵A来表示，记则由)(1)1( )111() () ()(-)kkxDLUxDILgMDLUGaussdeblLSei。式中矩阵为迭代法的迭代矩阵迭代法求解线性方程组 迭代公式四（超松弛迭代 / SOR）：( )(1)

16、(1)( )(1)(1( )( )(1)( )(1)()( )1( )- (1) ()kkGkkGkkkkkkkkGGGkkSORGauss SeidelxxGauss SexxidelxxxSORxxxxxxxxxx 是 迭代法的一种加速法。 假设 :已知，为 迭代法结果，定义 :，则,迭代法的表达式为：11-1Gauss Seidel其中， 称为松弛因子。当时称为低松弛；是迭代；时称为超松弛法。 相当于在高斯迭代法的基相当于在高斯迭代法的基础上再向前（础上再向前（11）或向后或向后（11）移动一段距离。移动一段距离。迭代法求解线性方程组 超松弛迭代 / SOR迭代矩阵：(1)1(1)1(

17、)1(1)( )1(1)1( )1( )(1)( )( )1(1)1( )1( )( )1(1)1( ) = =() () (1)kkkkkkkkkkkkkkkkkxD LxD UxD bxxxD LxD U xD bxxxxxD LxD UxD bxxD LxD UxD根据高斯迭代法的矩阵表示：1-1-11-1(1)1)111(-1-(-)() (1)( () (1) ()kkbID LID LDLxDLDU xMDLDUgDLDLbb因为,故（）与存在， 有：松弛法的迭代矩阵为:Code 12clear allclose allclcshgclf resetset(gcf, color,

18、white, menubar, none, numbertitle, off, name, ?)x = -0.0; 0.5h = plot(x(1), x(2), .)mycolor = 1, 0, 0;set(h, markersize, 20, color, mycolor, erasemode, none);axis(0.8, 1.05, 0.8, 1.05);axis ongrid on% stop = uicontrol(style, toggle, string, stop, background, white)% drawnowA = 5, -4.0; -4.0, 5;b = 1

19、; 1;xdir = A b% % 方法 1% M1 = A + eye(2);% g1 = -1 * b;% % for iii = 1 : 50% x = -1 * M1 * x + g1;% set(h, xdata, x(1), ydata, x(2)% drawnow% e1(iii) = norm(x - xdir);% pause(0.1)% end% % figure% % plot(e1)% Jacobix = -0.0; 0.5Mja = A / 5 - eye(2);gja = b / 5;for iii = 1 : 30 x = -1 * Mja * x + gja;

20、set(h, xdata, x(1), ydata, x(2) drawnow eja(iii) = norm(x - xdir); pause(0.5)end% figure% plot(eja)% Gauss-seidelx = -0.0; 0.5set(h, xdata, x(1), ydata, x(2) mycolor = 0,0.8,0set(h, markersize, 18, color, mycolor, erasemode, none);x = -0.0; 0.5D = A(1, 1), 0; 0, A(2, 2)L = -1 * 0, 0; A(2, 1), 0U = -1 * 0, A(1, 2); 0, 0Mgs = inv(D - L) * Uggs = inv(D - L) * bfor iii = 1 : 15 x = Mgs * x + ggs; set(h, xdata, x(1), y