NA32用矩陣分解法求解線性方程組

數(shù)值計算方法蘭洲交通大學(xué)數(shù)理學(xué)院1－第二節(jié)用矩陣分解法求解線性方程組

一、矩陣的三角分解

1.順序高斯消元與矩陣的LU分解的等價性

順序高斯消元的基本思想：將矩陣A的下三角部分消為零，

數(shù)值分析數(shù)值分析2－數(shù)值分析數(shù)值分析3－數(shù)值分析數(shù)值分析4－數(shù)值分析數(shù)值分析5－

數(shù)值分析數(shù)值分析6－數(shù)值分析數(shù)值分析7－數(shù)值分析數(shù)值分析8－進行到第k步消元時數(shù)值分析數(shù)值分析9－數(shù)值分析數(shù)值分析10－數(shù)值分析數(shù)值分析11－其中數(shù)值分析數(shù)值分析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分解，可以證明，當(dāng)A的各階順序主子式均不為零時，Doolittle分解可以實現(xiàn)并且唯一。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(yù)=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%功能：調(diào)用列主元三角分解函數(shù)[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用追趕發(fā)求解三對角方程組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－

對另一類方程組，在周期樣條插值等為題遇到的循環(huán)三對角方程組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

結(jié)果輸出

write(*,*)'三對角方程組的解為’'do40i=1,nwrite(*,30)'x(',i,')=',d(i)30format(1x,a2,