求解线性方程组的直接解法.doc

求解线性方程组的直接解法5.2LU分解 Gauss消去法实现了LU分解顺序消元结束时的上三角矩阵U和所用的乘数,严格下三角矩阵。将下三角矩阵的对角元改成1,记为L,则有A=LU,这事实是一般的,我们不难从消去的第k个元素时的矩阵k行及k列元素的历史得到这一点.因为从消元的历史有ukj=akj-mk1u1j- mk2u2j - mk,k-1uk-1,j,j=k,k+1,nmik=(aik-mi1u1k- mi2u2k -mi,k-1uk-1,k)/ukki=k+1,k+2,n于是akj=mk1u1j+mk2u2j+mk,k-1uk-1,j+ukj,j=k,k+1,naik=mi1u1k+mi2u2k+mi,k-1uk-1,k+mikukki=k+1,k+2,n从前面两个式子我们可以直接计算L和U(见下段).将矩阵分解为单位下三角矩阵和上三角矩阵之积称为矩阵的LU分解.顺序消元实现了LU分解,同时还求出了g, Lg=b的解. 直接LU分解上段我们得到(lij=mij)ukj=akj-lk1u1j-lk2u2j - lk,k-1uk-1,j,j=k,k+1,nlik=(aik-li1u1k-li2u2k -li,k-1uk-1,k)/ukki=k+1,k+2,n二式也可从A=LU的n2个等式解出.下面以n=3为例说明.u11=a11u12=a12u13=a13l21=a21/u11u22=a22-l21u12u23=a23-l21u13l31=a31/u11l32=(a32-l31u12)/ u22u33=a33-l31u13- l32u23从表上看到每个元素由所在位置的元素减去同行L左边诸元素与上方U诸元素对应乘积,只不过算L的元素时还要除以同列对角元.这一规律很容易记住.可写成算法(L和U可存放于A):fork=1:n-1forj=k:nukj=akj-lk1u1j-lk2u2j - lk,k-1uk-1,jendfori=k+1:nlik=(aik-li1u1k-li2u2k -li,k-1uk-1,k)/ukkendend这一算法也叫Gauss消去法的紧凑格式,可一次算得L,U的元素,不需逐步计算存储.考察上面的表格会发现还可安排其它计算次序,只要在这一次序下每个元素左边的L的元素与上方的U的元素已计算在先。例如, 逐行自左而右的次序, 逐列自上而下的次序, 易知g的计算规律同U.利用LU分解解Ax=b分三步:1分解A=LU2解Lg=b求g3解Ux=y求x例3. 用直接LU分解法解解 用分解公式计算得 求解 其它分解我们用顺序消元和直接分解两种方法实现了LU分解.还有更一般的三角分解,比如,下三角矩阵和单位上三角矩阵之积,又如单位下三角矩阵,对角矩阵,单位上三角矩阵之积,等等.下面给出第二种分解形式的算法LDR分解法。A=LDR,L是单位下三角矩阵,D是对角矩阵,R是单位上三角矩阵.逐列计算(逐列作LU分解,再用U的对角元素除各行),结果存入A。forj=1:nfori=2:jaij=aij-ai1a1j-ai2a2j -ai,i-1ai-1,jendfori=j+1:naij=(aij-ai1a1j-ai2a2j -ai,j-1aj-1,j)/ajjendfori=1:j-1aij= aij/aiiendend 列主元素的LU分解对照顺序消元和LU分解,列主元素法也可得列主元素的LU分解:PA=LUP是行交换结果的排列阵,L和U同前.例4. 列主元素法解方程组并写出系数矩阵的LU分解.括号内是乘数,k=2时2,3行交换.因而有直接作列主元素LU分解,因为在k步要先选主元素,所以作如下改变:fork=1:n-1fori=k:naik=aik-li1u1k-li2u2k -li,k-1uk-1,kend找p:p行k行ik=pforj=k+1:nukj=akj-lk1u1j-lk2u2j - lk,k-1uk-1,jendfori=k+1:nlik= aik/ukkendend可将lik存于aik,ukj存于akj.二 实验部分本章实验内容：实验题目：Gauss消元法,追赶法,范数。实验内容：编制用Gauss消元法求解线性方程组Ax=f的程序。 编制用追赶法求解线性方程组Ax=f的程序。 编制向量和矩阵的范数程序。实验目的：了解Gauss消元法原理及实现条件,熟练掌握Gauss消元法解方程组的算法,并能计算行列式的值。 掌握追赶法，能利用追赶法求解线性方程组。理解向量和矩阵范数定义,性质并掌握其计算方法.编程要求：利用Gauss消元法,追赶法解线性方程组。分析误差。计算算法:Gauss消元法:1. 消元过程设,对计算回代过程追赶法：1.分解Ax=f: ( )2.解Lg=f,求g: ()3解Ux=g,求x: （）范数：常用向量范数有：（令x=( x1,x2,xn)）1-范数： x1=x1+x2+xn2-范数： x2=(x12+x22+xn2)1/2-范数： x=max(x1,x2,xn)常用的三种向量范数导出的矩阵范数是：1-范数:A1= maxAx1/x1=1=2-范数:A2=maxAx2/x2=1=,1是ATA的最大特征值.-范数:A=maxAx/x=1=实验例题：用Gauss消元法解方程组 实验例题: 用追赶法解三对角方程组Ax=b,其中 实验例题：设 ， 计算A的行列范数.程序:Gauss消元法function x=Gauss(A,b) %A是线性方程组的系数矩阵,b为自由项.n=length(A)for k=1:n-1 m(k+1:n,k)=A(k+1:n,k)/A(k,k); A(k+1:n,k+1:n)=A(k+1:n,k+1:n)-m(k+1:n,k)*A(k,k+1:n); b(k+1:n)=b(k+1:n)-m(k+1:n,k)*b(k);endx=zeros(n,1);x(n)=b(n)/A(n,n);for k=n-1:-1:1 x(k)=(b(k)-A(k,k+1:n)*x(k+1:n)/A(k,k);endx=x;disp(sprintf(k x(k);for i=0:n disp(sprintf(%d %f ,i,x(i+1);end数值结果：x=Gauss(A,b)n =3k x(k)0 1.111111 1 0.777778 2 2.555556程序:追赶法function x,y,beta=zhuiganfa(a,b,c,f)%a,b,c是三对角阵的对角线上的元素,f是自由项.n=length(b);beta(1)=c(1)/b(1);for i=2:n beta(i)=c(i)/(b(i)-a(i)*beta(i-1);endy(1)=f(1)/b(1);for i=2:n y(i)=(f(i)-a(i)*y(i-1)/(b(i)-a(i)*beta(i-1);endx(n)=y(n);for i=n-1:-1:1 x(i)=y(i)-beta(i)*x(i+1);enddisp(sprintf(k x(k) y(k) beta(k);for i=0:n disp(sprintf(%d %f ,i,x(i+1),y(i+1),beta(i+1);end数值结果:a=0 -1 -1 -1 -1;b=2 2 2 2 2;c=-1 -1 -1 -1 0;f=1 0 0 0 0;x,y,beta=zhuiganfa(a,b,c,f)k x(k) y(k) beta(k)0 0.833333 5.000000e-001 -0.500000 1 0.666667 3.333333e-001 -0.666667 2 0.500000 2.500000e-001 -0.750000 3 0.333333 2.000000e-001 -0.800000 4 0.166667 1.666667e-001 0.000000程序： 1.列范数:function fan=lie(A)%A为已知矩阵n=length(A)for j=1:n x(j)=0 for i=1:n x(j)=x(j)+abs(A(i,j); endendfan=max(x)disp(sprintf(n x(n);for i=0:ndisp(sprintf( %d %f,i,x(i+1);end数值结果: fan=lie(A) fan =0.8000n x(n)0 0.7000001 0.8000002.行范数:function fan=hang(A)%A为已知矩阵n=length(A)for i=1:n x(i)=0 for j=1:n x(i)=x(i)+abs(A(i,j); endendfan=max(x)disp(sprintf(n x(n);for i=0:n disp(sprintf( %d %f,i,x(i+1);end数值结果: fan=hang(A) fan =1.1000 n x(n)0 1.1000001 0.400000总结:从代数上看,直接分解法和Gauss消去法本质上一样,但如果我们采用”双精度累加”计算,那么直接三角分解法的精度要比Gauss消去法为高. 求线性方程组的直接法,其算式繁杂,给人以枯燥沉闷的感觉.为了改善教学效果,本章着重介绍了三对角方程组的追赶法.三对角方程组以及其拓广形式的带状方程组有着广泛的实际应用.追赶法是解三对角线方程组(对角元占优势)的有效方法,它具有计算量少,方法简单,算法稳定等优点,具有鲜明的对称美.复习思考题 1. 用消去法解线性方程组为什么最好选主元？怎样的方程组可以不用选主元？ 2.