基于MATLAB的数学实验——线性代数中的若干问题(二)2.1引言用迭代法(Iteration)寻求问题的答案,是逐步逼近思想的具体体现,其基本思想是:不追求“一下子”得到(非)线性方程(组)的解,而是在逐步逼近方程组的精确解的迭代过程中获得满足精度要求的近似解,这一点与直接法不同;通过对问题的转化,避免(困难的)矩阵求逆运算。在不能或者很难一下子求得结果的情况下,使用迭代法就是一种自然的选择,因此了解并学会使用迭代法是十分有益的。在本节你将了解到1.什么是迭代;2.什么样的迭代是可用的;3.线性迭代、非线性迭代以及二者的差异;2.2(非)线性方程的场合数列是我们熟知的数学概念和对象,下面就是一例222,,,(1)这个数列产生于这样的过程:从2x开始,计算函数xxf在x2点处的值得到数列的第一项,然后将计算结果再代入函数的表达式中得到数列的第二项,数列(1)重复这种操作(迭代)的结果,其中2x称为迭代的初值。迭代公式(定义了一种迭代算法)为:nnxx1。由此可见,任何函数(包括一元函数和多元函数)在其定义域内都对应一种迭代,线性函数、非线性函数确定的迭代分别称为线性迭代和非线性迭代。面对这无穷多的迭代算法,在实际使用时如何选择呢?下面我们通过一个具体的代数方程求根的例子来说明选择的依据和准则,例用迭代法求代数方程根的例子:让我们来求如下方程的根fxxxxx3223030.sin.(2)首先要确定适当的包含根的区间,这可以依据闭区间上连续函数的介值定理来确定,例如,f1110sin,f222090sin.,所以方程(2)至少有一个实根属于区间]21[,,图1表明区间]21[,中只含有一个根,显然方程(2)的根不易直接求得。-1-0.500.511.52-2.5-2-1.5-1-0.500.51y=f(x)图1下面,我们采用迭代法求方程(2)位于区间]21[,中的根,为此构造迭代算法如下(原因稍后加以说明):nnnnnxxxxgx3.2sin3.01,,,21n(3)在区间]21[,中任取一个迭代初值0x,如取初值20x。执行下面的程序:EqutIteration.m:x=[];x(1)=2;n=100;fori=1:100;x(i+1)=(0.3+x(i)*sin(x(i)))/(x(i)*(2.3-x(i)));ifabs(x(i+1)-x(i))1e-15;breakelseN=i+1;endendNxplot([1:length(x)],x)title('ITERATIONFORSOLVINGEQUATION')gtext('CurveofIteration')共迭代42次,得到结果如下:xn2.00003.53100.23940.72330.68280.66180.65190.64730.64530.64440.64410.64390.64380.64380.64380.64380.64380.64380.64380.64380.64380.64380.64380.64380.64380.64380.64380.64380.64380.64380.64380.64380.64380.64380.64380.64380.64380.64380.64380.64380.64380.64380.6438可以看出n=13以后,在取4位有效数字时xn不再变化,下图也表明了这一事实,事实上(能够证明)用迭代所求方程的根,误差不超过9.91514e-020,这是相当精确的。05101520253035404500.511.522.533.54ITERATIONFORSOLVINGEQUATION图2以上我们已经看到,按迭代算法(3)生成的数列不仅收敛而且收敛于所求方程的根,这种收敛性是对迭代算法的基本要求,下面的定理给出了迭代算法收敛的一个充分条件:定理设方程0xf在][ba,上存在唯一解,xgx是方程的等价形式,如果1、xg在][ba,上连续可微;2、对任何][bax,,][baxg,;3、1Lxg,则对任何][0bax,,由迭代算法xgxnn1,(4)生成的序列nx收敛于方程0xf在][ba,上的唯一解。问题与实验1:用迭代算法(3)能否求出方程(2)的所有解?不妨用程序EqutIteration.m先做实验,再说明原因。问题与实验2:就方程(2)或选择别的例子构造不同的收敛迭代法,通过实验和比较,你可能会有所感悟。2.3(非)线性方程组的场合用迭代法求解线性方程组,首先要把线性方程组写成等价的形式fMxxbAx(5)式(5)的右端称为迭代格式,由迭代格式(5)确定如下的迭代算法:,2,101mRxfMxxnmm(6)对于给定的线性方程组,可以写成不同的(无穷多)迭代格式,有意义的(可用的)迭代格式应具有收敛性―生成的解向量序列xn收敛于方程组的解;而好的迭代法应具有较高的收敛速度。关于迭代法收敛性的两个判别条件:a、充分必要条件是:矩阵M的谱半径1,,2,1maxniMiib、充分条件是:矩阵M的某个算子范数M1。设x是方程组(5)的解,mx是迭代法(6)生成的任一序列,因为fMxx,fMxxmm1所以0221xxMxxMxxMxxmmmm(7)设11TJTMJMTT,其中矩阵J是矩阵M的Jordan标准型,那么容易验证1TTJMmm,并且10lim,2,10limlim0limlim010MMnixxTJTxxMxxmmminmmmmmm(8)此外,因为02211xxMxxMxxMxxMxxmmmmm(9)所以xxxxMMmmmmmmlim0lim0lim1(10)注:迭代格式(6)所确定的迭代法收敛与否,完全由系数矩阵M决定,而与常数项f无关.例用迭代法求线性方程组根的例子:让我们来求5阶线性方程组:其中A=[4.39990.6686-1.60410.5287-1.0106;0.69005.19080.25730.21930.6145;0.8156-1.20254.0565-0.92190.5077;0.7119-0.01981.41513.17071.6924;1.2902-0.1567-0.8051-0.05925.5913]b=[-0.4326;-1.6656;0.1253;0.2877;-1.1465]构造迭代法(Jacobi迭代法)如下:xMxfnn1(5)其中Ddiadaaa112255,,,,MDDA1=[0-0.15200.3646-0.12020.2297-0.13290-0.0496-0.0422-0.1184-0.20110.296400.2273-0.1252-0.22450.0062-0.44630-0.5338-0.23080.02800.14400.01060]fDb1=[-0.0983-0.32090.03090.0907-0.2051]使用程序ITERA_J:A=[4.39990.6686-1.60410.5287-1.0106;0.69005.19080.25730.21930.6145;0.8156-1.20254.0565-0.92190.5077;0.7120-0.01981.41513.17071.6924;1.2902-0.1567-0.8051-0.05925.5913];b=[-0.4326;-1.6656;0.1253;0.2877;-1.1465]D=diag(diag(A));LU=D-A;M=D\LU;f=D\b;x=[];z=[];x(:,1)=eye(5,1);N=2000;fori=1:N;ifnorm(A*x(:,i)-b)1e-010;m=i;breakelsex(:,i+1)=M*x(:,i)+f;z=x(:,i+1);endendme=norm(A*z-b)plot([1:length(x)],x)title('JACOBIITERATIONOFLINEAREQUATIONS')gtext('x1')gtext('x2')gtext('x3')gtext('x4')gtext('x5')迭代42次,误差向量eAxb427.0162e-011,并且看出n23以后,在取4位有效数字时xn不再变化,迭代过程如图3所示。051015202530354045-0.500.51JACOBIITERATIONOFLINEAREQUATIONSx1x2x3x4x5图3问题与实验3:一元线性迭代的收敛性条件怎样表述?问题与实验4:在本例中,12M,这时迭代序列是收敛的,就本例或选择别的例子,按12M和12M构造不同的迭代法,通过实验和比较,并给出你对实验结果的解释(如关于收敛性、收敛速度等),当然这需要你首先知道矩阵范数的概念,并且对它有比较好的理解。1.4不动点、k循环与混沌我们已经知道任何函数(包括一元函数和多元函数)在其定义域内都对应一种迭代,对于非线性函数迭代来说,不动点是下面所定义的k—循环的特例(1—循环):定义:设由迭代算法nnxgx1,,,10n,生成的序列nx具有如下性质:Nk(自然数集合),0x][ba,,000xN,当00xNn时,knnxx,innxx,ki(6)那么称knnnxxx,,,1是由迭代算法nnxgx1生成的k—循环。引例关于迭代xaxann12(7)实验与讨论。容易求得迭代(7)有两个不动点:0*x和12*ax,图4给出了a2及迭代初值分别取0100.x和900.x时,迭代收敛于12*ax的情形:使用程序(ITERA1:)a=2;x=[];y=[];x(1)=input('pleaseinputthefirsevalueofiterationx(1)=:');%x(1)=0.01y(1)=input('pleaseinputthefirsevalueofiterationy(1)=:');%y(1)=0.9n=input('pleaseinputthetotalnumberofiterationn=:');%n=150fori=1:n-1;x(i+1)=a-(x(i)-sqrt(a))^2;y(i+1)=a-(y(i)-sqrt(a))^2;endsubplot(2,1,1)plot(1:n,x,'bo')title('ITERATIONFORSOLVINGEQUATION')subplot(2,1,2)plot(1:n,y,'r*')title('ITERATIONFORSOLVINGEQUATION')05010015000.511.52ITERATIONFORSOLVINGEQUATION0501001500.511.52ITERATIONFORSOLVINGEQUATION图4图5给出了a3时,出现4—循环的情形:使用程序(ITERA2:)a=3;x=[];x(1)=input('pleaseinputthefirsevalueofiterationx(1)=:');%2.5n=input('pleaseinputthetotalnumberofiterationn=:');%n=200fori=1:n-1;x(i+1)=a-(x(i)-sqrt(a))^2;endplot(1:n,x,'r*')title('ITERATIONFORSOLVINGEQUATION')0204060801001201401601802001.41.61.822.22.42.62.83ITERATI