NA32用矩阵分解法求解线性方程组

数值计算方法兰洲交通大学数理学院1－第二节用矩阵分解法求解线性方程组

一、矩阵的三角分解

1.顺序高斯消元与矩阵的LU分解的等价性

顺序高斯消元的基本思想：将矩阵A的下三角部分消为零，

数值分析数值分析2－数值分析数值分析3－数值分析数值分析4－数值分析数值分析5－

数值分析数值分析6－数值分析数值分析7－数值分析数值分析8－进行到第k步消元时数值分析数值分析9－数值分析数值分析10－数值分析数值分析11－其中数值分析数值分析12－13－14－例3.2矩阵分解法LY=bUx=Yy1=6y2=-4y3=-4x1=-13x2=8x3=2A=LU3/1815－基于高斯消元的LU分解算法A=[2,3,4;3,5,2;4,3,30];b=[6;5;32];[n,r]=size(A);fork=1:n-1fori=k+1:nm=A(i,k)/A(k,k);forj=k+1:nA(i,j)=A(i,j)-m*A(k,j);endA(i,k)=m;endendA=2.003.004.001.500.50-4.002.00-6.00-2.004/1816－fori=2:ns=0;forj=1:i-1s=s+A(i,j)*b(j);endb(i)=b(i)-s;endb(n)=b(n)/A(n,n);fori=n-1:-1:1s=0;forj=i+1:ns=s+A(i,j)*b(j);endb(i)=(b(i)-s)/A(i,i);endb’=-1382triu(A)ans=2.00003.00004.000000.5000-4.000000-2.0000tril(A)+eye(3)-diag(diag(A))ans=1.0000001.50001.000002.0000-6.00001.00005/1817－Doolittle分解若矩阵A有分解：A=LU，其中L为单位下三角阵，U为上三角阵，则称该分解为Doolittle分解，可以证明，当A的各阶顺序主子式均不为零时，Doolittle分解可以实现并且唯一。18－Doolittle分解A的各阶顺序主子式均不为零，即19－Doolittle分解20－1、Doolittle分解L为下三角，U为单位上三角比较第1行：比较第1列：21－Doolittle分解22－Doolittle分解23－例3.5求矩阵的Doolittle分解

12/1824－13/1825－举例（一）解：例：用LU分解求线性方程组Ax=b

的解，其中令A=LU，由LU分解可得回代：解Ly

=b

得：y=[2,8,18,24]T解Ux

=y

得：x

=[-1,1,-1,1]T26－27－28－29－30－31－lupdsv.m%功能：调用列主元三角分解函数[LU,p]=lud(A)%求解线性方程组Ax=b。%解法：PA=LU,Ax=b←→PAx=Pb%LUx=Pb,y=Ux%Ly=f=Pb,f(i)=b(p(i))%输入：方阵A，右端项b（行或列向量均可）%输出：解x（行向量）32－functionx=lupdsv(A,b)n=length(b);[LU,p]=lu(A);y(1)=b(p(1));fori=2:ny(i)=b(p(i))-LU(i,1:i-1)*y(1:i-1)';endx(n)=y(n)/LU(n,n);fori=(n-1):-1:1x(i)=(y(i)-LU(i,i+1:n)*x(i+1:n)')/LU(i,i);end33－34－a=243141264>>[l,up]=lu(a)l=1.00000

00.50001.000001.00001.00001.0000u=2.00004.00003.00000

2.0000-0.50000

0

1.5000p=10001000135－a=264141385>>[l,u,p]=lu(a)l=1.00000

00.33331.000000.66670.50001.0000u=3.00008.0000

5.00000

1.3333

-0.66670

0

1.0000p=00101010036－>>p*aans=385141264>>[l,u,p]=lu(p*a)37－l=1.00000

00.33331.000000.66670.50001.0000u=3.00008.0000

5.00000

1.3333

-0.66670

0

1.0000p=10001000138－39－对称正定矩阵的Cholesky分解AT=A,A的各阶顺序主子式大于零.A的LU分解uii

>0(i=1,···,n)40－AT=A

RTDLT=LDR

RT=L

R=

LT由

ukk>0,记

得其中(2)对称矩阵的LLT分解(1)对称矩阵的

LDLT

分解:A=LDLT在不引起符号混乱时仍记为:A=LLTA=LU

A=LDR

AT=RTDLT41－LDLT分解计算公式d1=a11,l21=a21/d1,······,ln1=an1/d142－A

a21

a21/a11,······,an1=an1/a11(j=2,······,n-1)(i=

j+1,······,n)工作量:n(n–1)(n+4)/343－44－例

.LDLT分解45－例

用分解法求解46－解Ly=b得y=(6,1,-1)T。解LTx=D-1y得x=(2,1,-1)T。

例

用LDLT分解求解47－

三对角方程组求解的追赶法48－49－50－51－52－53－54－55－56－57－58－例4.7用追赶发求解三对角方程组Ax=d，其中

解得由求解公式得59－60－求解Ux=y,x4=0.3333,x3=-0.3333,x2=-1,x1=-1求解Ly=b,y1=1,y2=1.5,y3=1,y4=0.561－

对另一类方程组，在周期样条插值等为题遇到的循环三对角方程组Ax=d，其中我们也可用三对角分解的方法。从矩阵零元素的位置不难验证L和U可写成下面的形式：由此不难得到L和U的元素的计算公式，这里不在介绍。62－解三对角矩阵线性方程组的追赶法程序框图63－64－ccc

用追赶法求例3.4implicitreal*8(a-h,o-z)

parameter(n=4)dimensiona(n),b(n),c(n),d(n)dataa/0,-1,-1,-1/,b/4*2/datac/-1,-1,-1,0/,d/1,0,0,1/ccc

追的过程

c(1)=c(1)/b(1)d(1)=d(1)/b(1)do10i=2,n

b(i)=b(i)-a(i)*c(i-1)

c(i)=c(i)/b(i)10d(i)=(d(i)-a(i)*d(i-1))/b(i)ccc

赶的过程

do20i=1,n-120d(n-i)=d(n-i)-c(n-i)*d(n-i+1)ccc

结果输出

write(*,*)'三对角方程组的解为’'do40i=1,nwrite(*,30)'x(',i,')=',d(i)30format(1x,a2,