ch9MATLAB在数学中的应用

1、19.1 多项式与插值多项式与插值9.1.1 插值问题与插值多项式插值问题与插值多项式若有函数若有函数y=f(x)，给定一个区间，给定一个区间a,b上的上的n+1个点个点x0,x1,xn，且有，且有ax0 x1xnb,要求用一个简要求用一个简单的便于计算的函数单的便于计算的函数p(x)近似代替近似代替f(x)，使，使p(xi)=f(xi)=yi,则则p(x)就称为就称为f(x)的插值多项式。的插值多项式。第第9章章 MATLAB在数学中的应用2假设假设p(x)=a0 xn+a1xn-1+an-1x+an 为所求的插为所求的插值多项式，问题的关键是求出值多项式，问题的关键是求出a0,a1,a2,

2、an，由，由于有于有p(xi)=f(xi)=yi，于是可以解一个，于是可以解一个n+1个未知个未知数的线性方程组。数的线性方程组。但直接求解较复杂，下面将用插值法。但直接求解较复杂，下面将用插值法。39.1.2 Lagrange插值插值一、线性插值一、线性插值已知两点（已知两点（x0,y0）及）及(x1,y1)，通过此两点的线性插值多项式为，通过此两点的线性插值多项式为(可由两点公式求得可由两点公式求得)L1(x)= y0+ y1 显然显然L1(x0)=f(x0)=y0,L1(x1)=f(x1)=y1,满足插值条件，所以满足插值条件，所以 L1(x)就是线性插值。就是线性插值。若记若记l0(x

3、)= ,l1(x)= ,则称则称l0(x)和和l1(x)为为x0与与x1的线性插的线性插值基函数。值基函数。101xxxx010 xxxx101xxxx010 xxxx于是有于是有: L1(x)=l0(x)y0+l1(x)y14 二、二次插值二、二次插值当当n=2时，已给时，已给3点点(x0,y0)，(x1,y1)，(x2,y2),依相同原理依相同原理l0(x)=l1(x)=l2(x)=120102()()()()xxxxxxxx021012()()()()xxxxxxxx012021()()()()xxxxxxxx5于是得到二次插值多项式于是得到二次插值多项式L2(x),可表示为可表示为:L

4、2(x)=l0(x)y0+l1(x)y1+l2(x)y2三、n次Lagrange插值已知n+1个点(xi,f(xi)(i=0,1,n)的插值多项式Ln(x),应满足的条件为Ln(xi)=f(xi)=yi。用插值基函数方法可得n次Lagrange插值多项式为Ln(x)= ，其中li(x)=niiixfxl0)()().()().().()().()110110niiiiiiniixxxxxxxxxxxxxxxx6li(x)称为关于x0,x1,xn的n次插值基函数，它满足条件li(xj)=1(当i=j时)或0(当i不等于j)。定理9-1n次Lagrange插值多项式Ln(x)是存在且唯一的。证明假

5、设，n次Lagrange插值多项式Ln(x)为a0 xn+ a1xn - 1+ + an - 1x + an, 已 知 n + 1 个 点(xi,f(xi)(i=0,1,n)的插值多项式Ln(x),应满足的条件为Ln(xi)=f(xi)=yi。将这n+1个已知条件带入，则得到一个有n+1个方程n+1个未知数的线性方程组，由线性代数的知识可知道，该方程组有唯一解。故Ln(x)是存在且唯一的。7定理9-2设f(x)Cn+1a,b(表示f(x)在a,b上连续且有n+1阶导数)，且节点ax0 x1xnb，则插值多项式Ln(x)与f(x)的误差Rn(x)=f(x)-Ln(x)=，其中，ab, n+1(x

6、)=(x-x0)(x-x1)(x-xn)证明:由插值条件可知Rn(xi)=0(i=0,1,n)，故对任何xa,b有Rn(x)=k(x)(x-x0)(x-x1)(x-xn)=k(x) n+1(x)其中，k(x)是依赖x的待定函数，将xa,b看做区间a,b上任一固定点，作函数(t)=f(t)-Ln(t)-k(x)(t-x0)(t-x1)(t-xn)显然(xi)=0，(i=0,1,n),且(x)=0,它表示(t)在a,b上有n+2个零点x0,x1,xn,x，由rolle定理可知(t)在a,b上至少有n+1个零点。反复应用rolle定理，可知(n+1)在a,b上至少有一个零点a,b，使)()!1()(

7、1)1(xnfnn8(n+1)()=f(n+1)()-(n+1)!k(x)=0即k(x)=,得到所需结论。例9-1 已知sin0.32=0.314567,sin0.34=0.333487,sin0.36=0.352274，用线性插值及二次插值计算sin0.3367的近似值并估计误差。解 由题意y=f(x)=sin(x)，x0=0.32,y0=0.314567,x1=0.34,y1=0.333487,x2=0.36,y2=0.352274，线性插值为(1)( )(1)!nfn101xxxx010 xxxxL1(x)=y0+y19=（x-0.34）/(0.32-0.34)*y0+(x-0.32)/

8、(0.34-0.32)y1将x=0.3367代入，有Sin0.3367=L1(0.3367)=0.330365二次插值为L2(x)=l0(x)y0+l1(x)y1+l2(x)y2误差为|R1(x)|=f(x)/2*(x-x0)(x-x1)sinx1/2*0.0167*0.00330.92X10-5将x=0.3367代入，有Sin0.3367=L2(0.3367)=0.330374误差为|R2(x)|=f(x)/6*(x-x0)(x-x1)(x-x2)0.204X10-610四、n次Lagrange插值算法自定义函数lagan.m如下(c为插值多项式的系数)functionc,l=lagan(x

9、,y)w=length(x);n=w-1;l=zeros(w,w);fork=1:n+1v=1;forj=1:n+1 ifk=jv=conv(v,poly(x(j)/(x(k)-x(j)endendl(k,:)=vendc=y*l;11MATLAB命令窗口调用如下（以例9-1为例），得到的结果如图9-1所示。图9-1lagrange插值示范从图可知，求得sin0.3367=0.3304。129.1.3 Newton插值插值1.Newton插值利用Lagrange插值基函数求Lagrange插值多项式时，基函数li(x)的计算比较复杂，若增加节点数目时，所有基函数要重新计算。于是，我们可以给出另

10、一种便于计算的插值多项式Nn(x)，它可以表示为Nn(x)=a0+a1(x-x0)+a2(x-x0(x-x1)+an(x-x0)(x-xn-1)，其中ai(i=0,n)为待定常数。所有的ai可以根据插值条件Nn(xi)=f(xi)=yi求得。13例如，若x=x0时，a0=f(x0)=y0;若x=x1时，a1=为了求出a2,a3,an，先引进以下定义。定义9-1记fxm=f(xm)为f的零阶均差，零阶均差的差商记为fx0,xm=称为函数关于点x0,xm的一阶均差。一般地，记k-1阶均差的差商为fx0,xk-1,xk=称为f关于点x0,x1,.xk-1,xk的k阶均差。0101)()(xxxfxf

11、 00f xmf xxm x110 1, 2,., 0, 1,.,kkkkf x xxxf xxxxx14 根据均差定义，把x看成a,b上一点，可得f(x)=f(x0)+fx,x0(x-x0)fx,x0=fx0,x1+fx,x0,x1(x-x1)fx,x0,xn-1=fx0,x1,xn+fx,x0,xn(x-xn)将后一式代入前一式，就会得到f(x)=f(x0)+fx0,x1(x-x0)+fx0,x1,x2(x-x0)(x-x1)+fx0,x1,xn(x-x0)(x-xn-1)+fx,x0,xn=Nn(x)+Rn(x)15其中Nn(x)=f(x0)+fx0,x1(x-x0)+fx0,x1,x2

12、(x-x0)(x-x1)+fx0,x1,xn(x-x0)(x-xn-1)为Newton插值，Rn(x)=fx,x0,xn为插值的误差。162.Newton插值举例 例 9 - 2 已 知 s i n 0 . 3 2 = 0 . 3 1 4 5 6 7 ,sin0.33=0.324043,sin0.34=0.333487,sin0.35=0.342898,sin0.36=0.352274,求四次Newton插值多项式，并求出sin0.3367的近似值。解依题意有x0=0.32,x1=0.33,x2=0.34,x3=0.35,x4=0.36,f(x0)=0.314567,f(x1)=0.32404

13、3,f(x2)=0.333487,f(x3)=0.342898,f(x4)=0.352274。于是可得f(x)的均差表如表9-1所示。17表9-1f(x)的均差表于是，可得四次Newton插值多项式为N4(x)=0.314567+0.947600(x-x0)-0.160004(x-x0)(x-x1)-0.166667(x-x0)(x-x1)(x-x2)-4.166667(x-x0)(x-x1)(x-x2)(x-x3)将x=0.3367代入得到N4(0.3367)=0.330374。0.320.3145670.330.3240430.9476000.340.3334870.944399-0.16

14、00040.350.3428980.941100-0.164999-0.1666670.360.3522740.937600-0.175000-0.333333-4.166667183.Newton插值算法实现自定义函数如下（c表示多项式系数，d表示均差表）functionC,D=lagan1(x,y)n=length(x);D=zeros(n,n);D(:,1)=y;forj=2:nfork=j:nD(k,j)=(D(k,j-1)-D(k-1,j-1)/(x(k)-x(k-j+1);endendC=D(n,n);fork=(n-1):-1:1C=conv(C,poly(x(k);m=leng

15、th(C);C(m)=C(m)+D(k,k);end19在Matlab命令窗口运行结果如图9-2所示。从图可知，求得sin0.3367=0.3304。图9-2Newton插值示范209.2 数值积分与数值微分数值积分与数值微分9.2.1 数值积分数值积分在工程应用中我们经常会遇到定积分的计算，设F(x)是f(x)在a,b上的原函数，则有=F(b)-F(a)。但是，在实际问题中，常常会遇到如下三方面的问题：第一，f(x)不是连续函数，甚至也不是解析函数，而是一些离散点。第二，f(x)的原函数不能用初等函数表示。第三，f(x)的原函数很难求得。( )baf x dx( )baf x dx21所以在

16、应用中，需要构造一种积分方法，避免求原函数的计算，又能在误差允许范围内求出积分，这就是数值积分所要解决的问题。数值积分的基本思想是用被积函数f(x)在积分区间a,b上某些点处的函数值的线性组合来近似代替定积分，即有求积公式)()()(0fExfAdxxfniiiba其中xia,b称为求积节点，Ai称为求积系数，它与求积节点xi有关，与f(x)的具体表达形式无关，E(f)称为余项（误差）。221.梯形求积用梯形面积近似代替曲边梯形面积求得定积分（f(x)用线性插值代替）+E(f)其中E(f)=-)()(2)(bfafabdxxfba 3)(121fab 232.抛物型求积（1/3simpson求

17、积）f(x)用二次插值代替，得到抛物型求积。+E(f)其中E(f)=-( )( ( )4 ()( )62babaabf x dxf aff b5()2880baf243.Newton-Cotes求积将 区 间 a , b 分 成 n 等 分 ， 取 等 距 节 点 ：ax0 x1xnb，其中xi=，i=1,n+1则Newton-Cotes求积公式可表示为+E(f)(1)baian11( )( )nbiiaiIf x dxA f x25其中，h=当n=1,2,3时，Newton-Cotes求积分别变成了梯形求积，抛物型求积（1/3simpson求积）和3/8simpson求积。Newton-Co

18、tes求积精确度比梯形求积和抛物型求积要高，误差要小，但是求出所有的Ai比较麻烦。0( 1)(1).(2)().()(1)!(1)!nnihAt ttititn dtini nab 264.复合梯形求积若将区间a,b分成n等分，对每一个小区间用梯形求积，则得到复合梯形求积。+E(f)其中h=(b-a)/n，xi=a+i-1)/h，fi=f(xi)，i=1,2,n+1，E(f)=-121( )(2.2)2bnnahIf x dxffff212bah f275.复合抛物型求积若将区间a,b分成n等分（n为偶数），对每二个小区间用抛物型求积，则得到复合抛物型求积（复合1/3simpson求积）。+E

19、(f)其中，h=(b-a)/n，分fi=f(a+(i-1)*h)，E(f)=-(b-a)123411( )(424.24)3bnnnahIf x dxfffffff4180hf286.求积方法举例例9-3分别用梯形求积、1/3simpson求积、n=6的复合梯形求积、n=6的复合1/3simpson求积求定积分，并与真值比较。解:由高等数学知识=1/3+1+1=2.3333用梯形求积=(1-0)/2*(1+4)=2.5120(21)xxdx用1/3simpson求积=1/6*(1+9+4)=2.3333用n=6的复合梯形求积= 1 / 1 2 * ( f1+ 2 f2+ 2 f3+ 2 f4+

20、 2 f5+ 2 f6+ f7)=1/12*(1+49/18+64/18+9/2+50/9+121/18+4)=2.3380120(21)xxdx120(21)xxdx120(21)xxdx120(21)xxdx29用n=6的复合1/3simpson求积= 1 / 1 8 * ( f1+ 4 f2+ 2 f3+ 4 f4+ 2 f5+ 4 f6+ f7)=1/18*(1+49/9+64/18+9+50/9+121/9+4)=2.33337.求积Matlab算法实现用算法求出例9-3的结果120(21)xxdx30自定义函数fff.m如下functionz=fff(x)z=x.2+2*x+1;复

21、合梯形求积函数如下：functionI=fhtx(a,b,n)h=(b-a)/n;I=0;fori=1:n+1x(i)=a+(i-1)*h;I=I+2*fff(x(i);endI=I-fff(x(1)-fff(x(n+1);I=h*I/2;31在命令窗口中调用如下:I=fhtx(0,1,6)则得到结果I=2.3380若用I=fhtx(0,1,1)则得到结果I=2.5相当于梯形求积。若要对其它函数求定积分，只要将函数fff的内容改掉即可，若要变换区间，将函数调用中第1、2个参数改掉即可。32复合1/3simpson求积函数如下：functionI=fhsimpson(a,b,n)h=(b-a)/

22、n;I=0;fori=1:n+1x(i)=a+(i-1)*h;ifi=1,I=I+fff(x(i);elseifi=n+1,I=I+fff(x(i);elseifmod(i,2)=0,I=I+4*fff(x(i);elseI=I+2*fff(x(i);endendI=h*I/3;33在命令窗口中调用如下:I=fhsimpson(0,1,6)则得到结果I=2.3333若用I=fhsimpson(0,1,2)则得到结果I=2.3333相当于1/3simpson求积。349.2.2 数值微分数值微分 将f(x)在x0处进行Taylor展开f(x)=f(x0)+(x-x0)f(x0)+若令h=xk+1

23、-xk,fk=f(xk)导数定义为2( )0000()()().().2!nnxxxxfxfxn()0()()limkkkfxhf xhf xh35则得到两点公式的数值微分上面的O(h)来源于的。若希望提高精度，可设法消除这一项1( )kkkfffO hh1()()()().2!kkkkf xf xhfxfxhh()2!kfx36将上面的两个式子相加，则得到三点的数值微分若希望进一步提高精度，则可得到五点的数值微分1( )1()()()().().2!nnkkkkf xf xhhfxfxhn211()2kkkfffO hh421121(88)()12kkkkkfffffO hh37同理，对于二

24、阶导数，也可以用三点公式和五点公式。三点公式：五点公式：21122()kkkkffffO hh4211221(163016)()12kkkkkkffffffO hh389.3 非线性方程的求根非线性方程的求根9.3.1 概述概述数学物理中的很多问题经常归结为解方程f(x)=0如果有x*使得f(x*)=0，则称x*为方程f(x)=0的根或函数f(x)的零点，如果函数f(x)能写成如下形式： f(x)=(x-x*)mg(x)39其中g(x*)0，m为正整数，当m=1时，称x*为f(x)=0的单根或f(x)的单零点；当m2时，称x*为f(x)=0的m重根或f(x)的m重零点。如果f(x)为超越函数，

25、则称f(x)=0为超越方程；如果f(x)为n次多项式，则称f(x)=0为n次代数方程。一般来说，对于次数较高的代数方程，它的根不能用方程系数的解析式表示。而对于一般的超越方程更没有求根的公式可套。因此，研究非线性方程的数值解法非常有必要。非线性方程的求根通常分为两个步骤：一是对根的搜索，分析方程存在多少个根，找出每个根所在区间；二是根的精确化，求得根的足够精确的近似值。401.根的搜索(1)图解法例如，方程3x-cosx-1=0，等价于3x-1=cosx,在同一坐标系中画出y=3x-1和y=cosx的图。由图可知两曲线交点的横坐标即为方程的根x*1/3,1。(2)近似方程法例如方程3x-cos

26、x+0.01sinx-1=0的根，接近于方程3x-cosx-1=0的根。(3)解析法根据函数的连续性、介值定理以及单调性等去寻找有根区间和有唯一根的区间。41(4)定步长搜索法在某一区间上以适当的步长h,去考察函数值f(xi)(xi=x0+ih,i=0,1,2,)的符号，当f(x)连续且f(xi-1).f(xi)0时，则区间xi-1,xi为有根区间，又若在此区间内f(x)不变号，则在此区间有唯一根。2.二分法二分法是以连续函数的介值定理为基础。考虑方程f(x)=0。设函数f(x)Ca,b且f(a)f(b)0，则方程在a,b内至少存在一个根。42二分法的基本思想是：用对分区间的方法根据分点处函数

27、f(x)的符号逐步将有根区间缩小，使在足够小的区间内，方程有且仅有一根。为方便起见，记a0=a,b0=b。用中点x0=(a0+b0)/2将区间a0,b0分成2个小区间a0,x0和x0,b0。计算f(x0)。若f(x0)=0，则x0为f(x)=0的根，计算结束；否则若f(a0)f(x0)0，令a1=a0,b1=x0;若f(b0)f(x0)0，令a1=x0,b1=b0,不管那种情况，都有f(a1)f(b1)0,break,endmax1=-1+round(log(b-a)-log(delta)/log(2);fork=1:max1c=(a+b)/2;yc=feval(f,c);ifyc=0 a=c

28、;b=c;elseifyb*yc044b=c; yb=yc;elsea=c;ya=yc;endifb-adelta,break,endendc=(a+b)/2;err=abs(b-a);yc=feval(f,c);45例9-4求f(x)=x3-x-1=0在1，1.5内的根，要求精确度达到10-4。可以先定义函数f如下:functiony=f(x)y=x3-x-1;再调用函数i,j,k=bisect(f,1,1.5,1e-4)，得到结果如下：I=1.3248,j=2.4414*10-4,k=4.7404*10-4i为求得的根，j为最后一个区间的长度，k为所求点的函数值（接近0）。469.3.2

29、简单迭代法简单迭代法1.迭代公式的构造若有方程f(x)=0在区间a,b内有1个根x*。可以将方程f(x)=0改写成等价的形式x=(x)，取x0a,b，利用递推公式xk+1=(xk)可以求得序列x1,x2,xn。如果迭代方法是收敛的，则xn就是所求的根。例9-5用迭代法求f(x)=x3-x-1=0在x0=1.5附近的根。47解若将方程改为x=x3-1,构造迭代公式,k=0,1,由 x0= 1 . 5 , 通 过 迭 代 公 式 可 得 到x1=2.375,x2=12.39,可知，xk显然不收敛。若将方程改为x=(x+1)1/3,构造迭代公式Xk+1=(xk+1)1/3,k=0,1,由 x0= 1

30、 . 5 , 通 过 迭 代 公 式 可 得 到x1=1.35721,x2=1.33086,x3=1.32588,x4=1.32494,x5=1.32476,x6=1.32473,x7=1.32472,x8=1.32472。可知，方程的根为x*=1.32472。482.迭代公式的收敛性定理9-3设Ca,b,如果对xa,b有a(x)b，且存在常数L(0,1)使|(x)-(y)|L|x-y|,x,ya,b，则在区间a,b上存在惟一点x*,使生成的迭代序列xk对任何x0a,b收敛于x*，并有误差估计 *10|1kkLxxxxL49证明先证明x*的存在性，记f(x)=x-(x)，由定理有f(a)=a-

31、(a)0及f(b)=b-(b)0，若有一等号成立，则f(a)=0或f(b)=0，即x*存在。再证明惟一性，设Ca,b都是的不动点，且不相等，则由定理得与假设矛盾，所以，即不动点x*是惟一的。*12,xx*12121212| | ()()| |xxxxL xxxx*12xx50下面证明收敛性因为0L1,故|xk-x*|=0，即xk=x*利用定理有|xk+p-xk|=|xk+p-xk+p-1+xk+p-1-xk+p-1+xk+1-xk|xk+p-xk+p-1|+|xk+p-1-xk+p-2|+|xk+1-xk|(Lp-1+Lp-2+L+1)|xk+1-xk|=*110| | ()()| .|kkk

32、kxxxxL xxLxxlimklimk1101|11pkkkLLxxxxLL51推论若C1a,b（表示在a,b上一阶导数连续 ） ， 则 定 理 9 - 3 中 的 条 件 可 以 改为，则定理10.3中结论成立。例9-6构造不同迭代法求x2-3=0的根x*=解(1)(k=0,1,),(x)=，不满足收敛条件。0|( )|1maxx bxL 313kkxx3x*23( ),()1xxx 52(2)(k=0,1,)，迭代收敛。（k=0,1,），211(3)4kkkxxx2*113( )(3),( )1,()( 3)11422xxxxxx 113()2kkkxxx*21313( )(),( )(

33、1),()( 3)0122xxxxxx迭代收敛。533.迭代法的算法实现以例9-6为例来说明(1)的m文件如下functiony1=g1(x,n)fori=1:nz=3/x;x=zendy1=x在Matlab命令窗口中调用为X=2Y1=g1(x,8)则得到的结果为2，1.5，2，1.5，2，1.5，2，1.5，不收敛。54(2)的m文件如下functiony2=g2(x,n)fori=1:nz=x-1/4*(x2-3);x=zendy2=x在Matlab命令窗口中调用为X=2Y2=g2(x,8)则得到的结果为2，1.75，1.7344，1.7324，1.7321，1.7321，1.7321，1

34、.7321，收敛。得到的根为1.7321。55(3)的m文件如下functiony3=g3(x,n)fori=1:nz=1/2*(x+3/x);x=zendy3=x在Matlab命令窗口中调用为X=2Y3=g3(x,8)则得到的结果为2，1.75，1.7321，1.7321，1.7321，1.7321，1.7321，收敛。得到的根为1.7321。569.3.3 Newton法法1.Newton法对于方程f(x)=0，要求出它的根x*，如果已知它的一个近似值xk，利用Taylor展开式求出f(x)在xk附近的线性近似，即有忽略余项后，可以得到近似公式f(x)f(xk)+f(xk)(x-xk)=0

35、，解出方程，得到，记为xk+1，于是得到Newton法迭代公式2( )( )()()()()2!kkkkff xf xfxxxxx()()kkkf xxxfx1()()kkkkf xxxfx572.Newton法的算法实现functionp1=newton(f,df,p0,epsilon,max1)%f表示函数f(x),df表示函数f(x)，p0表示迭代的初值，epsilon表示迭代达到的精度%max表示迭代的最大次数fork=1:max1 p1=p0-feval(f,p0)/feval(df,p0)p0=p1;y=feval(f,p0);if(abs(y)maxmax=a(i,k);row=

36、i;endendtemp=a(k,:);a(k,:)=a(row,:);a(row,:)=temp;t=b(k);b(k)=b(row);b(row)=t;69b(k)=b(k)/a(k,k);forj=n:-1:ka(k,j)=a(k,j)/a(k,k);endfori=k+1:nb(i)=b(i)-a(i,k)*b(k);forj=n:-1:ka(i,j)=a(i,j)-a(i,k)*a(k,j);endendend70 x(n)=b(n);fork=n-1:-1:1t=0;forj=k+1:nt=t+a(k,j)*x(j);endx(k)=b(k)-t;end71例如，给定如下方程组在M

37、atlab命令窗口中调用函数，得到的结果如图9-3所示。图9-3列主元guass消去法求解723.无回带过程的列主元Gauss消去法若在消元的过程中，将系数矩阵消元成一个单位矩阵，则方程组的解可以直接得到。73算法描述如下functionx=c84(a,b,n)fork=1:nrow=k;max=a(k,k);fori=k+1:nifa(i,k)maxmax=a(i,k);row=i;endendtemp=a(k,:);a(k,:)=a(row,:);a(row,:)=temp;t=b(k);b(k)=b(row);b(row)=t;b(k)=b(k)/a(k,k);74forj=n:-1:k

38、a(k,j)=a(k,j)/a(k,k);endfori=k+1:nb(i)=b(i)-a(i,k)*b(k);forj=n:-1:ka(i,j)=a(i,j)-a(i,k)*a(k,j);endendend75fori=n:-1:2fork=i-1:-1:1b(k)=b(k)-a(k,i)*b(i);a(k,i)=0;endendx=b;x例如，给定如下方程组76在Matlab命令窗口中调用如图9-4所示。图9-4无回代的列主元Gauss消去法774.LU分解(1)LU分解原理若将系数矩阵A分解成单位下三角矩阵L和上三角矩阵U，解方程组Ax=b，可变成LUx=b，这时，可以分解成解方程组Ly

39、=b和Ux=y。可利用如下公式求出L和U矩阵。78111111,1,2,.,iiiiaualinu11,2,3,., ,1,.,rririrkkikual urn ir rn11()/,2,3,., ,1,.,riririkkrrrklal uurn ir rn79再利用Ly=b求出y最后利用Ux=y求出x11,1,2,.,kkkkrrrybl y kn11(),1,.,1nkkkrrr kkkxyu xkn nu LU分解也称杜利特尔(Doolitte)分解。80(2)LU分解的算法实现functionx=LU1(a,b,n)fori=1:nu(1,i)=a(1,i);l(i,1)=a(i,

40、1)/u(1,1);endforr=2:nfori=r:ns=0;fork=1:r-1s=s+l(r,k)*u(k,i);end81u(r,i)=a(r,i)-s;s=0;fork=1:r-1s=s+l(i,k)*u(k,r);endl(i,r)=(a(i,r)-s)/u(r,r);endendfork=1:ns=0;forr=1:k-1s=s+l(k,r)*y(r);end82y(k)=b(k)-s;endfork=n:-1:1s=0;forr=k+1:ns=s+u(k,r)*x(r);endx(k)=(y(k)-s)/u(k,k);endy=y;x=xa;839.4.2 解线性方程组的迭代

41、法解线性方程组的迭代法1.Jacobi迭代法设有n阶线性代数方程组84若系数矩阵非奇异，且aii0(i=1,2,n)，上面方程组可改成85于是可以得到迭代公式86写成向量形式为X(k+1)=BX(k)+F其中，B=I-D-1A,F=D-1b,D为矩阵A的对角线元素所组成的对角矩阵，该迭代称为Jacobi迭代。算法描述如下functionz=jacobi(a,b,n,m,x)%a为系数矩阵，b为常量项，n为方程的数目，m为最大的迭代次数,x为迭代的初始值87fork=1:mfori=1:ns=0;forj=1:nifj=i,s=s+a(i,j)*x(j);endendy(i)=(b(i)-s)/

42、a(i,i);endx=y;endz=x88例如，给定如下方程组在Matlab命令窗口中调用如图9-5所示。图9-5jacobi迭代892.Gauss-Seidel迭代在Jacobi迭代中，如果迭代是收敛的，x应该比x更接近于原方程的解x（i=1,2,n），因此，在迭代过程中及时地用已经求出的x来代替x(i=1,2,n-1)，可以使迭代的效果更佳。于是可以将Jacobi迭代改成如下的形式：(1)ki(1)ki( )ki( )ki*i90该迭代称为Gauss-Seidel迭代，可以简写成下面的形式nixaxabaxnijkjijijkjijiiiki,.,2 , 1),(11)(11)1()1(

43、91前 面 介 绍 的 J a c o b i 迭 代 写 成 矩 阵 形 式 为x(k+1)=Bx(k)+F，Gauss-Seidel迭代也可以写成矩阵形式x(k+1)=(I-L)-1Ux(k)+(I-L)-1F，其中B=L+U，L为下三角矩阵，U是上三角矩阵。Gauss-Seidel迭代算法描述如下：functionz=G_S(a,b,n,m,x)%a为系数矩阵，b为常量项，n为方程的数目，m为最多的迭代次数，x为迭代的初值92fork=1:mfori=1:ns=0;forj=1:i-1s=s+a(i,j)*y(j);ends1=0;forj=i+1:ns1=s1+a(i,j)*x(j);endy(i)=(b(i)-s-s1)/a(i,i);endx=y;endz=x93例如，给定如下方程组在Matlab命令窗口中调用如图9-6所示。图9-6Gauss-Seidel迭代943.逐次超松驰迭代在Gauss-Seidel迭代中，将写成将迭代进行加速（i=1,2,n）得到nixaxabaxnijkjijijkjijiiiki,.,2 , 1),(11)(11)1()1(nixaxabanijkjij