現(xiàn)代數(shù)值計(jì)算（第3版）課件 第1-3章 科學(xué)計(jì)算與MATLAB、線性代數(shù)方程組的直接法、線性代數(shù)方程組的迭代法

第一章AdvancedNumericalComputing科學(xué)計(jì)算與MATLAB現(xiàn)代數(shù)值計(jì)算同濟(jì)大學(xué)數(shù)學(xué)科學(xué)學(xué)院目錄/Contents第一章科學(xué)計(jì)算與MATLAB第一節(jié)科學(xué)計(jì)算的意義第二節(jié)誤差基礎(chǔ)知識(shí)第三節(jié)MATLAB軟件1.1科學(xué)計(jì)算的意義科學(xué)計(jì)算的出現(xiàn)利用現(xiàn)代計(jì)算機(jī)輔助，解決實(shí)際問(wèn)題計(jì)算的挑戰(zhàn)：基因測(cè)序，全球天氣模擬科學(xué)計(jì)算問(wèn)題的主要步驟數(shù)學(xué)建模數(shù)值算法評(píng)價(jià)科學(xué)計(jì)算軟件MATLAB，Mathematica,Maple,Python,1.2誤差基礎(chǔ)知識(shí)實(shí)際問(wèn)題數(shù)學(xué)問(wèn)題(模型誤差)計(jì)算問(wèn)題(截?cái)嗾`差、觀測(cè)誤差)結(jié)果(舍入誤差)1.2.1誤差的來(lái)源實(shí)際問(wèn)題數(shù)學(xué)模型數(shù)值算法結(jié)果編程處理模型誤差觀測(cè)誤差舍入誤差截?cái)嗾`差比較檢驗(yàn)1.2誤差基礎(chǔ)知識(shí)設(shè)有真值,及近似值，稱(chēng)為該近似值的絕對(duì)誤差1.2.2誤差度量稱(chēng)為絕對(duì)誤差限若稱(chēng)為相對(duì)誤差稱(chēng)為相對(duì)誤差限由于真值難以求出，通常也使用后者更加合理1.2誤差基礎(chǔ)知識(shí)十進(jìn)制數(shù)的標(biāo)準(zhǔn)形式(其中),1.2.3有效數(shù)字四舍五入保留位:因此有誤差限稱(chēng)為有效數(shù)字,

稱(chēng)為有效數(shù)問(wèn)題:有效數(shù)字和相對(duì)誤差限有什么關(guān)系?有效數(shù)的誤差限是末位數(shù)單位的一半,其本身就體現(xiàn)了誤差界有效數(shù)末尾不可以隨便增減零1.2誤差基礎(chǔ)知識(shí)1.2.4計(jì)算機(jī)的浮點(diǎn)數(shù)系單精度實(shí)數(shù)由32位二進(jìn)制的浮點(diǎn)數(shù)表示:1位符號(hào),

23位尾數(shù),

8位階數(shù)(本身也有符號(hào))機(jī)器所能表示的數(shù)中,離最近的是和.因此若,則在機(jī)器中記為,即.相對(duì)誤差限最大數(shù),最小數(shù),

上溢,

下溢1.2誤差基礎(chǔ)知識(shí)1.2.5一個(gè)實(shí)例有一艘駁船,寬度為5米,欲駛過(guò)一個(gè)河渠.該河渠有一個(gè)直角彎道,形狀和尺寸如圖所示.試問(wèn),要駛過(guò)這個(gè)河渠,駁船的長(zhǎng)度不能超過(guò)多少米?駁船的長(zhǎng)度有如下關(guān)系極小化問(wèn)題1.2誤差基礎(chǔ)知識(shí)1.2.5一個(gè)實(shí)例這個(gè)過(guò)程中有多少處有誤差？或者求零點(diǎn)問(wèn)題可證,對(duì)任意并可求得1.2誤差基礎(chǔ)知識(shí)1.2.6數(shù)值計(jì)算中應(yīng)注意的幾個(gè)問(wèn)題計(jì)算容易推導(dǎo)出因此,調(diào)用MatLab程序nademo11.2誤差基礎(chǔ)知識(shí)1.2.6數(shù)值計(jì)算中應(yīng)注意的幾個(gè)問(wèn)題計(jì)算,其中,保留四位有效數(shù)字.其它的例子:避免大數(shù)和小數(shù)相加減1.2誤差基礎(chǔ)知識(shí)1.2.6數(shù)值計(jì)算中應(yīng)注意的幾個(gè)問(wèn)題,其中計(jì)算中間步驟簡(jiǎn)化計(jì)算步驟Horner算法或秦九韶算法1.3Matlab軟件1.3.1簡(jiǎn)介全稱(chēng):MatrixLaboratory功能:科學(xué)計(jì)算、符號(hào)計(jì)算、圖形處理等數(shù)據(jù)類(lèi)型:數(shù)、字符串、矩陣、單元型數(shù)據(jù)和結(jié)構(gòu)型數(shù)據(jù)集成界面:命令窗口、命令歷史窗口、當(dāng)前路徑窗口、工作空間變量窗口等提示符>>,換行符...,注釋符

%,默認(rèn)變量ans1.3Matlab軟件1.3.2向量和矩陣的基本運(yùn)算矩陣A=[13;24]向量a=[123456]冒號(hào)a=1:6a:s:b列向量A=[1;2;3]字符串A='hellomatlab'

A='This''smatlab''sworld.'1.3Matlab軟件1.3.2向量和矩陣的基本運(yùn)算常量:在運(yùn)行過(guò)程中不能變化的量科學(xué)記數(shù)法:顯示方式:format(只影響顯示)變量:保存在內(nèi)存(地址),可隨時(shí)變化內(nèi)置變量:i,j,pi,Inf,NaN(NotaNumber)Inf及NaN的運(yùn)算規(guī)律1.3Matlab軟件1.3.2向量和矩陣的基本運(yùn)算矩陣的加(+)、乘(*)、數(shù)乘(*)、冪(^)矩陣的點(diǎn)乘(.*)、點(diǎn)除(./)、點(diǎn)冪(.^)矩陣的左除(X=A\B即求解AX=B)矩陣的右除(X=A/B即求解A=XB)1.3Matlab軟件1.3.2向量和矩陣的基本運(yùn)算>>x=[0pi/6pi/4pi/3pi/2];>>sin(x)ans=

0

0.5000

0.7071

0.8660

1.0000向量功能其他初等函數(shù):三角反三角、指數(shù)對(duì)數(shù)、根號(hào)、絕對(duì)值等等1.3Matlab軟件1.3.2向量和矩陣的基本運(yùn)算>>sqrt([91011])>=pians=

0

1

1>>a=[2300];>>b=[-1010];>>a&bans=

1

0

0

0>>a|bans=

1

1

1

0>>~bans=

0

1

0

11.3Matlab軟件1.3.2向量和矩陣的基本運(yùn)算矩陣運(yùn)算>>A=magic(3)

816>>A(2,1:3)

357>>A(2:end,[1end])492

>>B=[23];>>C=[12;34];>>D=[57]';>>A=[B9;CD]

>>A(A>=4)=0

>>v=1:9;>>v(abs(v-5)<=2)=[]1.3Matlab軟件1.3.3流程控制基本語(yǔ)法ifvalue1,

statement1,elseifvalue2,

statement2,else

statement3end1.3Matlab軟件1.3.3流程控制例如(判別閏年)ifmod(year,400)==0,

fprintf('%disaleapyear.\n',year);elseifmod(year,100)==0,

fprintf('%disnotaleapyear.\n',year);elseifmod(year,4)==0,

fprintf('%disaleapyear.\n',year);else

fprintf('%disnotaleapyear.\n',year);end1.3Matlab軟件1.3.3流程控制基本語(yǔ)法forloopvalue=value,

statement,end和whilevalue,

statement,end1.3Matlab軟件1.3.3流程控制例如(利用計(jì)算圓周率的近似值)>>s=0;>>fork=1:10000,

s=s+1/k^2;

end>>s=sqrt(6*s)1.3Matlab軟件1.3.3流程控制Collatz猜想：輸入一個(gè)正整數(shù)n,如果是偶數(shù)就除以2,是奇數(shù)就乘3加1,如此一直變換,最后會(huì)變成1.n=input('n=

');whilen~=1,

ifmod(n,2)==1,

n=n*3+1;

else

n=n/2;

end

disp(n);end1.3Matlab軟件1.3.3流程控制冒泡排序：把一列數(shù)想象為垂直存放,數(shù)值大的在下方,每輪比較時(shí)從上到下依次比較相鄰的兩個(gè)數(shù),若是上面的數(shù)大,把它們對(duì)調(diào),否則不動(dòng)。直至沒(méi)有對(duì)調(diào)為止。>>done=0;k=1;>>v=input('arowvector:');arowvector:

[1863975024]>>while~done,

done=1;

forp=1:length(v)-k,

ifv(p)>v(p+1),

tmp

=v(p);

v(p)=v(p+1);

v(p+1)=tmp;

%OR

v([pp+1])=v([p+1p]);

done=0;

end

end

k=k+1;

end1.3Matlab軟件1.3.4腳本文件和函數(shù)文件把一系列命令收集在一個(gè)文件里,保存為以.m為后綴的文件.執(zhí)行時(shí)只需要鍵入文件名,不需鍵入后綴.例:>>mysortarowvector:

[1863975024]v=

0

1

2

3

4

5

6

7

8

91.3Matlab軟件1.3.4腳本文件和函數(shù)文件一種封裝的文件,具有特定的頭格式:function[out1,out2,...]=funname(in1,in2,...)函數(shù)名必須和文件名一致與腳本文件的比較例:文件mysort2.m1.3Matlab軟件1.3.4腳本文件和函數(shù)文件函數(shù)頭function[v,s]=mysort3(v)調(diào)用>>d=[53421];>>[r,w]=mysort3(d)傳值方式:在輸出或輸入列表中的位置列表不一樣長(zhǎng)的情形命令global的用法1.3Matlab軟件1.3.4腳本文件和函數(shù)文件變量nargin和nargout的含義用法(例如:根據(jù)輸入計(jì)算面積)functions=zhouchang(a,b,c)ifnargin==1,

s=2*pi*a;elseifnargin==2,

s=2*(a+b);elseifnargin==3,

s=a+b+c;end1.3Matlab軟件1.3.4腳本文件和函數(shù)文件直接或間接地用到了自己例如:Fibonacci數(shù)列定義為functionf=fib(n)ifn>=3,

f=fib(n-1)+fib(n-2);elseifn==1|n==2,

f=1;end1.3Matlab軟件1.3.5幫助系統(tǒng)Help：查看工具箱,函數(shù)可以自己書(shū)寫(xiě)文件的幫助,寫(xiě)在function之后其他查看系統(tǒng)命令用法的工具:doc,lookfor其他幫助命令:which,who等輔助命令:clc,home,clear1.3Matlab軟件1.3.6畫(huà)圖功能>>x=0:0.01:10;>>y=1./(1+x.^2)+sin(x).*exp(x/3);plot(x,y,‘g*-’)畫(huà)函數(shù)的圖像hold命令命令plot中的參數(shù)選項(xiàng)plot(x,y1,'yo--',x,y2,'g*',x,y3,'r+:',x,y4,'bp:’);其他類(lèi)似命令:loglog,semilogx,semilogy1.3Matlab軟件1.3.6畫(huà)圖功能三維線圖>>t=linspace(0,10*pi,2000);>>plot3(sin(t).*t,cos(t).*t,t,'r-');>>view(-17,66)1.3Matlab軟件1.3.6畫(huà)圖功能三維面圖：命令meshgrid>>x=1:4;>>y=5:3:11;>>[X,Y]=meshgrid(x,y)X=

1

2

3

4

1

2

3

4

1

2

3

4Y=

5

5

5

5

8

8

8

8

11

11

11

111.3Matlab軟件1.3.6畫(huà)圖功能三維面圖：命令surf及contour例如:畫(huà)下面函數(shù)的圖像及等高線>>x=linspace(-10,10,200);>>[X,Y]=meshgrid(x);>>Z=exp(-abs(X))+cos(X+Y)+1./(X.^2+Y.^2+1);>>surf(X,Y,Z);>>contour(X,Y,Z,20)1.3Matlab軟件1.3.6畫(huà)圖功能標(biāo)注:坐標(biāo)軸xlabel,曲線legend,圖形標(biāo)題title窗口控制:打開(kāi)figure,關(guān)閉close,清除clf坐標(biāo)軸控制:axis('equal','off',[-1327])或者

axisequal;

axisoff1.3Matlab軟件1.3.7數(shù)據(jù)操作文本方式>>x=0:0.1:1;y=[x;exp(x)];>>fid=fopen('exp.txt','wt');>>fprintf(fid,'%s\n','%---exp.txt---');>>fprintf(fid,'%6.2f%12.8f\n',y);>>fclose(fid);>>fid=fopen('exp.txt');>>s=fscanf(fid,'%c',[117])s=%---exp.txt--->>A=fscanf(fid,'%6f%12f\n',[2inf]);>>A=A'轉(zhuǎn)置關(guān)系1.3Matlab軟件1.3.7數(shù)據(jù)操作二進(jìn)制方式>>loadclown.mat>>whoYourvariablesare:X

caption

map>>image(X);colormap(map)>>savea.matx*A1.3Matlab軟件1.3.7數(shù)據(jù)操作>>A=[120;030;0-16];>>A=sparse(A)A=

(1,1)

1

(1,2)

2

(2,2)

3

(3,2)

-1

(3,3)

6或者1.3Matlab軟件1.3.7數(shù)據(jù)操作>>i=[11233];>>j=[12223];>>s=[123-16];>>A=sparse(i,j,s)A=

(1,1)

1

(1,2)

2

(2,2)

3

(3,2)

-1

(3,3)

6指明階數(shù)>>A=sparse(i,j,s,200,100);找出非零元>>[i,j,s]=find(A);1.3Matlab軟件1.4評(píng)注Matlab參考書(shū)目:[1]MatLab與科學(xué)計(jì)算(第2版),王沫然,電子工業(yè)出版社,2007年8月[2]MatLab工程數(shù)學(xué)應(yīng)用,許波、劉征,清華大學(xué)出版社,2000年4月[3]MatLab數(shù)學(xué)實(shí)驗(yàn),胡良劍、孫曉君,高等教育出版社,2006年6月[4]MatLab高等數(shù)學(xué)實(shí)驗(yàn),章恩棟、馬玉蘭、徐美萍、李雙,電子工業(yè)出版社,2008年11月AdvancedNumericalComputing現(xiàn)代數(shù)值計(jì)算同濟(jì)大學(xué)數(shù)學(xué)科學(xué)學(xué)院學(xué)海無(wú)涯，祝你成功！第二章AdvancedNumericalComputing線性代數(shù)方程組的直接法現(xiàn)代數(shù)值計(jì)算同濟(jì)大學(xué)數(shù)學(xué)科學(xué)學(xué)院目錄/Contents第二章線性代數(shù)方程組的直接法第一節(jié)高斯消去法第二節(jié)矩陣的三角分解第三節(jié)正交矩陣與奇異值分解前言求解大規(guī)模線性方程組:如何利用計(jì)算機(jī)來(lái)快速、穩(wěn)定、有效地求解該問(wèn)題是科學(xué)計(jì)算的核心問(wèn)題之一直接法和迭代法由于浮點(diǎn)運(yùn)算的精度的影響，直接法不可能給出完全精確的計(jì)算解給定

階方陣

和

維向量

,尋找向量

使得：2.1高斯消去法2.1.1基本步驟記線性方程組的分量形式為下面演示用高斯消去法求解上述線性方程組的計(jì)算過(guò)程。計(jì)算過(guò)程分為消去過(guò)程和回代過(guò)程記矩陣

，向量

，它們的元素分別為：一、消去過(guò)程第一步：如果

，用數(shù)

依次乘以方程組的第一行，并加到第

行上去，

其中2.1高斯消去法2.1.1基本步驟第二步：如果

，用數(shù)

依次乘以方程組的第二行，并加到第

行上去，其中2.1高斯消去法2.1.1基本步驟類(lèi)似地，這樣的運(yùn)算過(guò)程一直可作到第

步，結(jié)果轉(zhuǎn)化為一個(gè)上三角形方程組2.1高斯消去法2.1.1基本步驟二、回代過(guò)程如果

，可以逐次回代計(jì)算出線性方程組的解。這就是求解線性方程組的高斯消去法。在沒(méi)有浮點(diǎn)運(yùn)算誤差的情況下，該方法在有限的計(jì)算步驟內(nèi)能夠得到原線性方程組的精確解，是直接法的一種。2.1高斯消去法2.1.1基本步驟【高斯消去法】算法(1)對(duì)

做：

對(duì)

做：

用數(shù)

乘以方程組第

行加到第

行上;

標(biāo)記得到的矩陣及右端向量的上標(biāo) ;(2) 且對(duì)于

做:2.1高斯消去法2.1.1基本步驟消去過(guò)程的第

步，對(duì)矩陣需作

次乘法運(yùn)算及

次除法運(yùn)算，對(duì)右端向量作

次乘法運(yùn)算，在消去過(guò)程總的乘除法運(yùn)算工作量為回代過(guò)程中，計(jì)算每個(gè)

需作

次乘除法運(yùn)算，其工作量為用高斯消去法計(jì)算線性方程組所需要總的乘除法運(yùn)算工作量為2.1高斯消去法2.1.2乘除法的運(yùn)算量高斯消去法能夠順利進(jìn)行到底是有前提條件的，即要求所有的主導(dǎo)元素不等于零。如果某個(gè)主元為零，則高斯消去法中斷?！径x2.1】在計(jì)算中做除數(shù)的元素

，被稱(chēng)為主導(dǎo)元素，簡(jiǎn)稱(chēng)主元。2.1高斯消去法2.1.3選主元策略【例2.1】取

，用高斯消去法計(jì)算下述線性方程組。（假定模型計(jì)算機(jī)具有8位字長(zhǎng)的浮點(diǎn)表示及16位的累加器。）解：首先用高斯消去法，這時(shí)

，對(duì)方程組消元，于是得到2.1高斯消去法2.1.3選主元策略在這個(gè)模型計(jì)算機(jī)上，具體計(jì)算是這樣的：

通過(guò)四舍五入，輸出結(jié)果變成

同理，

的計(jì)算結(jié)果也是

。這樣高斯消去的計(jì)算解為：

實(shí)際上方程組的精確解為

這里未知量

的計(jì)算解的相對(duì)誤差達(dá)到了驚人的

，這就是“大數(shù)吃小

數(shù)”的現(xiàn)象。2.1高斯消去法2.1.3選主元策略由于浮點(diǎn)運(yùn)算誤差的影響，高斯消去法過(guò)程中會(huì)得到錯(cuò)誤的解。為了避免上述不穩(wěn)定的現(xiàn)象，對(duì)一般的線性方程組而言，我們采用選主元的策略。采用列主元素高斯消去法，對(duì)上述例子中方程組進(jìn)行行交換：這時(shí)

，進(jìn)行消元后，可得到在模型計(jì)算機(jī)上，

和

都被算成為

。2.1高斯消去法2.1.3選主元策略若考慮方程此時(shí),選不選主元一樣:因?yàn)檫@種方程的選主元同等變形變得毫無(wú)意義。兩邊同乘了一個(gè)很大的數(shù)2.1高斯消去法2.1.3選主元策略目錄/Contents第一節(jié)高斯消去法第二節(jié)矩陣的三角分解第三節(jié)正交矩陣與奇異值分解第二章線性代數(shù)方程組的直接法2.2矩陣的三角分解常見(jiàn)的矩陣三角分解有：LU分解（杜利脫爾分解，克洛脫分解）LDU分解喬列斯基分解【定義2.1】把一個(gè)

階矩陣分解成結(jié)構(gòu)簡(jiǎn)單的三角形矩陣的乘積稱(chēng)為矩陣的三角分解。高斯消去法的消去過(guò)程與左乘下述矩陣是等價(jià)的：2.2.1LU分解2.2矩陣的三角分解其中易知?jiǎng)t其中2.2.1LU分解2.2矩陣的三角分解這里，是單位下三角矩陣，

是上三角矩陣，這種矩陣分解稱(chēng)為杜利脫爾(Doolittle)分解，或者杜利脫爾三角分解?？寺迕?Crout)分解LDU分解以上三種分解統(tǒng)稱(chēng)為矩陣的三角分解，或者

分解。如果不作特殊說(shuō)明，一般我們所說(shuō)的

分解就是指杜利脫爾三角分解。，這里是下三角矩陣，

是單位上三角矩陣。

，這里

是單位下三角矩陣，

是對(duì)角矩陣，是單位上三角矩陣。2.2.1LU分解2.2矩陣的三角分解矩陣三角分解的存在唯一性【定理2.1】(存在性)利用高斯消去法求解方程組

時(shí)的主元素

的充要條件是

階矩陣

A

的所有順序主子式均不為零。若

A

為

階矩陣，且所有順序主子式均不等于零，則

A

可分解為一個(gè)單位下三角矩陣

L與一個(gè)上三角矩陣

U

的乘積: ，且分解是唯一的?！径ɡ?.2】(唯一性)2.2.1LU分解2.2矩陣的三角分解【算法2.1】杜利脫爾算法(1)對(duì)

做:(2)(3)2.2.1LU分解2.2矩陣的三角分解(1)對(duì)

做:(2)(3)2.2.1LU分解2.2矩陣的三角分解【算法2.2】克洛脫算法【算法2.3】回代(1)對(duì)

做：(2)對(duì)

做：2.2.1LU分解2.2矩陣的三角分解解：【例2.1】利用LU分解求解

：2.2.1LU分解2.2矩陣的三角分解解：【例2.1】利用LU分解求解

：2.2.1LU分解2.2矩陣的三角分解2.2.1LU分解2.2矩陣的三角分解解：【例2.1】利用LU分解求解

：2.2.1LU分解2.2矩陣的三角分解解：【例2.1】利用LU分解求解

：2.2.1LU分解2.2矩陣的三角分解解：【例2.1】利用LU分解求解

：2.2.1LU分解2.2矩陣的三角分解解：【例2.1】利用LU分解求解

：2.2.1LU分解2.2矩陣的三角分解解：【例2.1】利用LU分解求解

：2.2.1LU分解2.2矩陣的三角分解解：【例2.1】利用LU分解求解

：用MATLAB可以計(jì)算矩陣的LU分解，其語(yǔ)法為:>>A=[211;

232;

234];

>>[l,u]=lu(A)

l=

1

0

0

1

1

0

1

1

1

u=

2

1

1

0

2

1

0

0

22.2.1LU分解2.2矩陣的三角分解>>b=[479]';>>y=L\b;

y=

4

3

2>>x=U\y

x=

1

1

12.2.1LU分解2.2矩陣的三角分解【例2.2】已知

，作A的杜利脫爾分解，并求解方程組，其中解：假設(shè)2.2.1LU分解2.2矩陣的三角分解經(jīng)計(jì)算可得對(duì)

進(jìn)行回代可得

；再對(duì)

進(jìn)行回代，則

。用MATLAB可以計(jì)算矩陣的LU分解，其語(yǔ)法為:2.2.1LU分解2.2矩陣的三角分解[L,U]=lu(A)當(dāng)矩陣A為對(duì)稱(chēng)正定時(shí)，它的所有順序主子式都大于零，易知存在唯一的LU分解：由A的對(duì)稱(chēng)性可得

按照分解的唯一性可得

即

2.2.2喬列斯基分解2.2矩陣的三角分解易證這個(gè)三角分解是唯一的，稱(chēng)之為喬列斯基(Cholesky)分解根據(jù)

A

是對(duì)稱(chēng)正定矩陣，有從而

記

則有其中，

G是對(duì)角元均大于零的下三角矩陣2.2.2喬列斯基分解2.2矩陣的三角分解【算法2.4】喬列斯基分解算法(1)對(duì)

(2) (3) 2.2.2喬列斯基分解2.2矩陣的三角分解給定喬列斯基分解，線性方程組

的求解可轉(zhuǎn)化為計(jì)算公式為平方根法2.2.2喬列斯基分解2.2矩陣的三角分解解：設(shè)A的喬列斯基分解

，經(jīng)計(jì)算得【例2.3】利用平方根法求解下述對(duì)稱(chēng)正定方程組2.2.2喬列斯基分解2.2矩陣的三角分解2.2.2喬列斯基分解2.2矩陣的三角分解解：設(shè)

A

的喬列斯基分解

，經(jīng)計(jì)算得【例2.3】利用平方根法求解下述對(duì)稱(chēng)正定方程組2.2.2喬列斯基分解2.2矩陣的三角分解解：設(shè)

A

的喬列斯基分解

，經(jīng)計(jì)算得【例2.3】利用平方根法求解下述對(duì)稱(chēng)正定方程組2.2.2喬列斯基分解2.2矩陣的三角分解解：設(shè)

A的喬列斯基分解

，經(jīng)計(jì)算得【例2.3】利用平方根法求解下述對(duì)稱(chēng)正定方程組2.2.2喬列斯基分解2.2矩陣的三角分解解：設(shè)A的喬列斯基分解

，經(jīng)計(jì)算得【例2.3】利用平方根法求解下述對(duì)稱(chēng)正定方程組2.2.2喬列斯基分解2.2矩陣的三角分解解：設(shè)

A

的喬列斯基分解

，經(jīng)計(jì)算得【例2.3】利用平方根法求解下述對(duì)稱(chēng)正定方程組回代2.2.2喬列斯基分解2.2矩陣的三角分解解：設(shè)

A

的喬列斯基分解

，經(jīng)計(jì)算得【例2.3】利用平方根法求解下述對(duì)稱(chēng)正定方程組【例2.4】求下列矩陣的喬列斯基分解2.2.2喬列斯基分解2.2矩陣的三角分解解：設(shè)

的喬列斯基分解

，經(jīng)計(jì)算得2.2.2喬列斯基分解2.2矩陣的三角分解解：設(shè)

的喬列斯基分解

，經(jīng)計(jì)算得【例2.4】求下列矩陣的喬列斯基分解2.2.2喬列斯基分解2.2矩陣的三角分解解：設(shè)

的喬列斯基分解

，經(jīng)計(jì)算得【例2.4】求下列矩陣的喬列斯基分解2.2.2喬列斯基分解2.2矩陣的三角分解解：設(shè)

的喬列斯基分解

，經(jīng)計(jì)算得【例2.4】求下列矩陣的喬列斯基分解2.2.2喬列斯基分解2.2矩陣的三角分解解：設(shè)

的喬列斯基分解

，經(jīng)計(jì)算得【例2.4】求下列矩陣的喬列斯基分解2.2.2喬列斯基分解2.2矩陣的三角分解解：設(shè)

的喬列斯基分解

，經(jīng)計(jì)算得【例2.4】求下列矩陣的喬列斯基分解利用矩陣的三角分解，很容易導(dǎo)出一些特殊方程組的解法。對(duì)矩陣

作克洛脫分解，得2.2.3追趕法2.2矩陣的三角分解設(shè)

為三對(duì)角矩陣，即(3）(4)(1)設(shè)

(2)對(duì)

算法2.52.2.3追趕法2.2矩陣的三角分解追趕法的Matlab程序tridiagsolver.m如下functionx=tridiagsolver(A,b)[n,n]=size(A);fori=1:nif(i==1) l(i)=a(i,i); y(i)=b(i)/l(i); else l(i)=a(i,i)-a(i,i-1)*u(i-1); y(i)=(b(i)-y(i-1)*a(i,i-1))/l(i); end if(i<n) u(i)=a(i,i+1)/l(i); end2.2.3追趕法2.2矩陣的三角分解endx(n)=y(n)forj=n-1:-1:1 x(j)=y(j)-u(j)*x(j+1);end2.2.3追趕法2.2矩陣的三角分解追趕法的Matlab程序tridiagsolver.m如下【例2.5】用追趕法求解下述三對(duì)角線性方程組解：追的過(guò)程為2.2.3追趕法2.2矩陣的三角分解趕的過(guò)程為因此，原方程組的解為

。2.2.3追趕法2.2矩陣的三角分解設(shè)滿(mǎn)足且

和

分別是分塊下三角矩陣和分塊下三角2.2.3分塊三角分解2.2矩陣的三角分解矩陣經(jīng)計(jì)算即這樣的矩陣分解稱(chēng)為矩陣的分塊三角分解2.2.3分塊三角分解2.2矩陣的三角分解其中，

稱(chēng)為

的舒爾(Schur)補(bǔ)【例2.6】求分塊矩陣的一個(gè)分塊三角分解，其中解：因?yàn)樗裕?.2.3分塊三角分解2.2矩陣的三角分解因此，

分塊矩陣

的一個(gè)分塊三角分解為:2.2.3分塊三角分解2.2矩陣的三角分解目錄/Contents第一節(jié)高斯消去法第二節(jié)矩陣的三角分解第三節(jié)正交矩陣與奇異值分解第二章線性代數(shù)方程組的直接法矩陣除了三角分解以外，還有QR分解為正交矩陣，

為上三角矩陣。奇異值分解和

為正交矩陣，

為對(duì)角矩陣。2.3正交矩陣與奇異值分解2.3.1正交矩陣正交矩陣

有如下性質(zhì)：

的長(zhǎng)度與

的長(zhǎng)度相等矩陣

，滿(mǎn)足

，

我們稱(chēng)這樣的矩陣

為正交矩陣。2.3正交矩陣與奇異值分解單位矩陣任意個(gè)置換矩陣的乘積仍然是正交矩陣置換矩陣：將單位矩陣的任意兩行（列）交換得到的矩陣。

譬如，交換第

行和第

行，得：2.3.1正交矩陣2.3正交矩陣與奇異值分解旋轉(zhuǎn)矩陣（Givens變換）形如的矩陣被稱(chēng)為Givens矩陣或Givens變換，或稱(chēng)（平面）旋轉(zhuǎn)矩陣（或旋轉(zhuǎn)變換）其中

為旋轉(zhuǎn)的角度。2.3.1正交矩陣2.3正交矩陣與奇異值分解平面旋轉(zhuǎn)變換的幾何意義：2.3.1正交矩陣2.3正交矩陣與奇異值分解反射矩陣（Householder變換）設(shè)

，且

，則稱(chēng)為Householder變換，或者Householder矩陣。

2.3.1正交矩陣2.3正交矩陣與奇異值分解

Householder矩陣有如下性質(zhì)：不難驗(yàn)證所以，Householder變換又稱(chēng)鏡面反射變換，Householder矩陣也稱(chēng)初等反射矩陣。(1)

，即

是對(duì)稱(chēng)陣。(2)

，即

是正交陣。(3)設(shè)

，

為過(guò)原點(diǎn)的平面且

，可分解成為

，其中

。2.3.1正交矩陣2.3正交矩陣與奇異值分解鏡面反射變換的幾何意義：2.3.1正交矩陣2.3正交矩陣與奇異值分解則從而其中,其中，

。一個(gè)重要的應(yīng)用是對(duì)

，求Householder矩陣

，使得由

知

，由

的構(gòu)造，有設(shè)

，為了使

計(jì)算時(shí)不損失有效數(shù)位，取2.3.1正交矩陣2.3正交矩陣與奇異值分解則【例3.1】已知

，求Householder矩陣

，使得

，

其中

解：取

，

，

，2.3.1正交矩陣2.3正交矩陣與奇異值分解證：給出構(gòu)造性證明如下：【定理3.1】設(shè)

，則存在正交陣

，使

，其中

為上三角陣。首先，考慮

的第一列

，可找到Householder矩陣

，使得

的除了第

個(gè)元素以外都為零。

同理，找到

使得

的第2列對(duì)角元以下元素為0，而第一列對(duì)角元以下元素與

一樣是0。依次這樣下去，可以得到其中

為上三角矩陣，

為正交陣。定理證畢。2.3.2QR分解2.3正交矩陣與奇異值分解【定理3.2】

，設(shè)

非奇異，則存在正交陣

與上三角陣

，使得

有如下分解：且當(dāng)

的對(duì)角元均為正時(shí)，分解是唯一的。該定理保證了

可分解為

，若

非奇異，則

也非奇異。如果不規(guī)定

的對(duì)角元為正，則分解不是唯一的。2.3.2QR分解2.3正交矩陣與奇異值分解【例3.2】用Householder變換作矩陣

的QR分解則有解：Householder矩陣

，使2.3.2QR分解2.3正交矩陣與奇異值分解和再找

，使

，得且這是一個(gè)下三角矩陣，但對(duì)角元皆為負(fù)數(shù)。2.3.2QR分解2.3正交矩陣與奇異值分解只要令

，則有

是對(duì)角元為正的上三角矩陣，使得

。其中，2.3.2QR分解2.3正交矩陣與奇異值分解

【定理3.3】（奇異值分解）設(shè)

，則存在酉陣

和

，使得：其中

, , 。【定義3.2】設(shè)

，

的

個(gè)特征值的非負(fù)平方根叫作

的奇異值，記為

。2.3.3奇異值分解2.3正交矩陣與奇異值分解奇異值分解與對(duì)稱(chēng)矩陣基于特征向量的對(duì)角化類(lèi)似，但還是有明顯的不同。對(duì)稱(chēng)矩陣特征向量分解的基礎(chǔ)是譜分析，而奇異值分解則是譜分析理論在任意矩陣上的推廣。奇異值分解提供了一些關(guān)于

的信息，例如非零奇異值的數(shù)目（

的階數(shù)）和

的秩相同。一旦秩r確定，那么U的前r列構(gòu)成了

的列向量空間的正交基。奇異值分解非常有用和可靠的分解，但是它比QR分解要花上近十倍的計(jì)算時(shí)間。2.3.3奇異值分解2.3正交矩陣與奇異值分解本章小結(jié)高斯消去法矩陣的三角分解法矩陣的QR分解法直接法相對(duì)來(lái)說(shuō)，工作量小，精度高，但程序復(fù)雜，并且對(duì)于高階線性方程組易于受計(jì)算機(jī)容量的限制，所以它適于求解中小型方程組。對(duì)于高階大型線性方程組，有效的解法是第三章要討論的迭代法。AdvancedNumericalComputing現(xiàn)代數(shù)值計(jì)算同濟(jì)大學(xué)數(shù)學(xué)科學(xué)學(xué)院學(xué)海無(wú)涯，祝你成功！第三章AdvancedNumericalComputing線性代數(shù)方程組的迭代法現(xiàn)代數(shù)值計(jì)算同濟(jì)大學(xué)數(shù)學(xué)科學(xué)學(xué)院目錄/Contents第三章線性方程組的迭代解法第一節(jié)范數(shù)和條件數(shù)第二節(jié)基本迭代法第三節(jié)不定常迭代法迭代方法的優(yōu)點(diǎn)：能夠充分利用系數(shù)矩陣的稀疏性質(zhì)占用內(nèi)存少運(yùn)算方便，計(jì)算程序也較為簡(jiǎn)單前言迭代法分類(lèi)常用定常迭代法雅可比（Jacobi）迭代法高斯-塞德?tīng)?Gauss-Seidel)迭代法超松弛(SOR)迭代法常用不定常迭代法共軛梯度法廣義極小殘量（GMRES）法前言前言基本概念【定義3.1】線性方程組的解是一個(gè)向量，稱(chēng)為解向量【定義3.2】近似解向量與精確解向量之差稱(chēng)為近似解的誤差向量為了估計(jì)誤差的大小，需要引入衡量向量和矩陣大小的度量概念-范數(shù)，滿(mǎn)足：【定義3.4】對(duì)任意

階方陣

，若對(duì)應(yīng)一個(gè)非負(fù)實(shí)數(shù)時(shí)成立；（1）

，等號(hào)當(dāng)且僅當(dāng)（2）對(duì)任意實(shí)數(shù)，

；（3）對(duì)任意兩個(gè)n階方陣

和

，

；（4）對(duì)任意兩個(gè)n階方陣

和

，

；則稱(chēng)

為方陣

的矩陣范數(shù)。為的譜半徑3.1范數(shù)和條件數(shù)3.1.1向量范數(shù)和矩陣范數(shù)【定義3.5】設(shè)

是

階矩陣，則稱(chēng)這里，

為

的特征值（

）對(duì)于

階方陣，定義分別為矩陣的1范數(shù)，2范數(shù)，無(wú)窮范數(shù)和范數(shù)。3.1范數(shù)和條件數(shù)3.1.1向量范數(shù)和矩陣范數(shù)【定義3.6】滿(mǎn)足的向量范數(shù)和矩陣范數(shù)，稱(chēng)為相容的。常用范數(shù)有下列相容關(guān)系：3.1范數(shù)和條件數(shù)3.1.1向量范數(shù)和矩陣范數(shù)考慮線性方程組解是由

和

決定的；數(shù)據(jù)

和

都帶有誤差（測(cè)量誤差，模型誤差，舍入誤差）；3.1范數(shù)和條件數(shù)3.1.2擾動(dòng)分析和條件數(shù)通常

和

的誤差相對(duì)于精確數(shù)據(jù)都是微小的，稱(chēng)之為小的擾動(dòng)，分別記為

和

。一個(gè)很自然的問(wèn)題是：

和

的微小擾動(dòng)將對(duì)計(jì)算線性方程組的解有何影響？考察解為若

，則原方程組變?yōu)榻鉃?.1范數(shù)和條件數(shù)3.1.2擾動(dòng)分析和條件數(shù)有即解的相對(duì)誤差是右端項(xiàng)相對(duì)誤差的10000倍。表明一些線性方程組系數(shù)的微小變化會(huì)引起解的巨大變化。為什么呢？3.1范數(shù)和條件數(shù)3.1.2擾動(dòng)分析和條件數(shù)若

有一個(gè)小擾動(dòng)

，則解

產(chǎn)生一個(gè)擾動(dòng)

，即：

有

兩邊取范數(shù)另得：3.1范數(shù)和條件數(shù)3.1.2擾動(dòng)分析和條件數(shù)若

有一個(gè)小擾動(dòng)

，則解

相應(yīng)產(chǎn)生一個(gè)擾動(dòng)

，即：有兩邊取范數(shù)得得3.1范數(shù)和條件數(shù)3.1.2擾動(dòng)分析和條件數(shù)觀察得，解

的相對(duì)誤差除了與方程組系數(shù)矩陣

和右端

擾動(dòng)有關(guān)外，還與

的大小有關(guān)。條件數(shù)的大小與范數(shù)有關(guān)：【定義3.7】設(shè)

為

階非奇異矩陣，稱(chēng)數(shù)

為線性方程組

的條件數(shù)，或者稱(chēng)為矩陣

的條件數(shù)。3.1范數(shù)和條件數(shù)3.1.3向量范數(shù)與矩陣范數(shù)如果

對(duì)稱(chēng)正定，則和

是最大和最小特征值。條件數(shù)有如下性質(zhì)：3.1范數(shù)和條件數(shù)3.1.3向量范數(shù)與矩陣范數(shù)對(duì)于任意的

階非奇異矩陣

，

成立；對(duì)于任意的正交矩陣

，有

；對(duì)于任意的

階非奇異矩陣

及任意非零常數(shù)

，

成立；對(duì)于任意的

階非奇異矩陣

及任意

階正交矩陣

，有典型的病態(tài)矩陣希爾伯特（Hilbert）矩陣隨著

的增大,條件數(shù)

非常迅速的增加希爾伯特矩陣階數(shù)越高，病態(tài)程度就越為嚴(yán)重。3.1范數(shù)和條件數(shù)3.1.3向量范數(shù)與矩陣范數(shù)目錄/Contents第三章線性方程組的迭代解法第一節(jié)范數(shù)和條件數(shù)第二節(jié)基本迭代法第三節(jié)不定常迭代法假定

有分解為

；其中

是非奇異方陣3.2

基本迭代法考慮有或；其中

且

3.2.1基本知識(shí)【算法】建立迭代公式：即給定初始向量

，按上述公式迭代，可以得到一個(gè)

向量序列

這種方法就是解線性方程組的基本迭代解法若

收斂于確定的向量

，則

，亦即

，

就是方程組

的解3.2

基本迭代法3.2.1基本知識(shí)令其中為對(duì)角矩陣，嚴(yán)格下三角矩陣和嚴(yán)格上三角矩陣通過(guò)不同的構(gòu)造方法，可得到三種基本迭代解法：雅可比迭代法高斯-塞德?tīng)柕⊿OR迭代法3.2

基本迭代法3.2.1基本知識(shí)線性代數(shù)方程組Jacobi迭代:從第

個(gè)方程解出

3.2

基本迭代法3.2.2雅可比（Jacobi）迭代法3.2

基本迭代法3.2.2雅可比（Jacobi）迭代法Jacobi迭代:從第

個(gè)方程解出

即取

和

，得雅可比迭代法:或：迭代矩陣為：【算法3.1】雅可比迭代算法1)選定初值2)對(duì)

計(jì)算:3)如果近似解達(dá)到收斂條件，退出；否則，繼續(xù)第2)步的計(jì)算。3.2

基本迭代法3.2.2雅可比（Jacobi）迭代法function[x,iter]=jacobi(A,b,tol)

D=diag(diag(A));

L=D-tril(A);

U=D-triu(A);

x=zeros(size(b));

foriter=1:500

x=D\(b+L*x+U*x);

error

=norm(b-A*x)/norm(b);

if(error<tol)

break;

end

endMatlab程序3.2

基本迭代法3.2.2雅可比（Jacobi）迭代法Gauss-Seidel迭代:從第

個(gè)方程解出 ,盡量使用新分量3.2

基本迭代法3.2.3高斯-塞德?tīng)?Gauss-Seidel)迭代法如果雅可比迭代法是收斂的，將

的分量投入下一個(gè)迭代方程中使用，這樣可能會(huì)收到更好的效果，即高斯-塞德?tīng)柕?或迭代矩陣為：3.2

基本迭代法3.2.3高斯-塞德?tīng)?Gauss-Seidel)迭代法【算法3.2】高斯-塞德?tīng)柕惴?)如果近似解達(dá)到收斂條件，退出；否則，繼續(xù)第2)步的計(jì)算1)選定初值

2)對(duì)

計(jì)算:3.2

基本迭代法3.2.3高斯-塞德?tīng)?Gauss-Seidel)迭代法function[x,iter]=gs(A,b,tol)

D=diag(diag(A));

L=D-tril(A);

U=D-triu(A);

x=zeros(size(b));

foriter=1:500

x=(D-L)\(b+U*x);

error

=norm(b-A*x)/norm(b);

if(error<tol)

break;

end

endMatlab程序3.2

基本迭代法3.2.3高斯-塞德?tīng)?Gauss-Seidel)迭代法迭代格式為:即：為了加快迭代的收斂速度，等號(hào)右端的第二項(xiàng)前乘以一個(gè)參數(shù)

，得逐次超松弛迭代法，簡(jiǎn)稱(chēng)SOR迭代法，

被稱(chēng)為松弛因子3.2

基本迭代法3.2.4超松弛(SOR)迭代法迭代矩陣為：注意當(dāng)

時(shí)，SOR迭代法其實(shí)就是GS迭代法【算法3.3】SOR迭代算法1)選定初值

,3)如果近似解達(dá)到收斂條件，退出；否則，繼續(xù)第2)步的計(jì)算。2)對(duì)

計(jì)算3.2

基本迭代法3.2.4超松弛(SOR)迭代法function[x,iter]=sor(A,b,omega,tol)

D=diag(diag(A));

L=D-tril(A);

U=D-triu(A);

x=zeros(size(b));

foriter=1:500

x=(D-omega*L)\(b+(1-omega)*D*x+omega*U*x);

error

=norm(b-A*x)/norm(b);

if(error<tol)

break;

end

endMatlab程序3.2

基本迭代法3.2.4超松弛(SOR)迭代法【例3.1】取初值為

，試用雅可比迭代，GS迭代以及 SOR迭代分別計(jì)算線性方程組迭代過(guò)程保留5位有效數(shù)字（精確解

）3.2

基本迭代法3.2.4超松弛(SOR)迭代法建立雅可比迭代格式解：雅可比方法迭代21步就得到了保留5位有效數(shù)字的近似解取

，

對(duì)

計(jì)算可得

3.2

基本迭代法3.2.4超松弛(SOR)迭代法建立GS迭代格式取

，對(duì)

計(jì)算可得GS方法迭代9步就得到了保留5位有效數(shù)字的近似解。對(duì)于這個(gè)例子而言，GS迭代的收斂速度大約是雅可比迭代的2倍。3.2

基本迭代法3.2.4超松弛(SOR)迭代法選取

，建立SOR迭代格式:取

，對(duì)

計(jì)算可得3.2

基本迭代法3.2.4超松弛(SOR)迭代法SOR方法迭代7步就得到了保留5位有效數(shù)字的近似解。對(duì)于這個(gè)例子而言，

時(shí)，SOR迭代法的收斂速度比GS迭代法要快。首先引進(jìn)兩個(gè)重要的概念:【定義3.8】設(shè)

是

階矩陣

，如果存在

階排列矩陣

，

使

，有如下形狀：3.2

基本迭代法3.2.5迭代的收斂性分析和誤差估計(jì)其中，

和

分別為

階和

階方陣

，則稱(chēng)

為可約矩陣，如果不存在這樣的排列陣，則稱(chēng)

為不可約陣。首先引進(jìn)兩個(gè)重要的概念:【定義3.9】設(shè)

是

階矩陣，若

滿(mǎn)足：且其中至少有一個(gè)嚴(yán)格不等式成立，稱(chēng)

是(行)弱對(duì)角占優(yōu)若式中的每一個(gè)不等式都是嚴(yán)格不等號(hào)，則稱(chēng)

是(行)嚴(yán)格對(duì)角占優(yōu)3.2

基本迭代法3.2.5迭代的收斂性分析和誤差估計(jì)【引理3.1】設(shè)

階矩陣

是嚴(yán)格對(duì)角占優(yōu)或不可約弱對(duì)角占優(yōu)矩陣，

則

是非奇異矩陣。證明(一)：假設(shè)

不可逆。則存在非零向量

，使得

不妨設(shè)

且

則

從

的第

個(gè)方程得是嚴(yán)格占優(yōu)的，因此上述不等式不成立，矛盾。3.2

基本迭代法3.2.5迭代的收斂性分析和誤差估計(jì)利用矩陣的Jordan標(biāo)準(zhǔn)型，可以得出下述引理【定理3.10】任一矩陣

的譜半徑均不大于

的任一與某一向量范數(shù)相容的

矩陣范數(shù)，即證明:兩邊取范數(shù)，假設(shè)

是矩陣

的特征對(duì)，則【引理3.2】設(shè)

是任意

階矩陣，則

的

次冪

（當(dāng)

）的

充要條件為譜半徑3.2

基本迭代法3.2.5迭代的收斂性分析和誤差估計(jì)【定理3.11】對(duì)于基本迭代格式，給定初值

，

是真解，有下列收斂

結(jié)果和誤差估計(jì)：(1)迭代格式收斂的充要條件為譜半徑

；(2)如果

，則有如下估計(jì)：3.2

基本迭代法3.2.5迭代的收斂性分析和誤差估計(jì)證明:是任意的，上式等價(jià)于

，即 .3.2

基本迭代法3.2.5迭代的收斂性分析和誤差估計(jì)【定理3.11】對(duì)于基本迭代格式，給定初值

，

是真解，有下列收斂

結(jié)果和誤差估計(jì)：(1)迭代格式收斂的充要條件為譜半徑

；證明:移項(xiàng)即得3.2

基本迭代法3.2.5迭代的收斂性分析和誤差估計(jì)【定理3.11】對(duì)于基本迭代格式，給定初值

，

是真解，有下列收斂

結(jié)果和誤差估計(jì)：(2)如果

，則有如下估計(jì)：【定理3.12】若

是嚴(yán)格對(duì)角占優(yōu)或不可約弱對(duì)角占優(yōu)矩陣，則雅可比迭代和 GS迭代都收斂【定理3.14】SOR迭代收斂的必要條件是

【定理3.13】若

是對(duì)稱(chēng)正定矩陣，則雅可比迭代收斂的充要條件是

也

對(duì)稱(chēng)正定【定理3.15】設(shè)系數(shù)矩陣

對(duì)稱(chēng)正定，則

時(shí)SOR迭代收斂3.2

基本迭代法3.2.5迭代的收斂性分析和誤差估計(jì)目錄/Contents第三章線性方程組的迭代解法第一節(jié)范數(shù)和條件數(shù)第二節(jié)基本迭代法第三節(jié)不定常迭代法本節(jié)將介紹兩類(lèi)最基本的不定常迭代方法：一類(lèi)是求解對(duì)稱(chēng)正定線性方程組的最速下降法和共軛梯度法另一類(lèi)是求解不對(duì)稱(chēng)線性方程組的廣義極小殘量法3.3

不定常迭代法對(duì)于一切

，通過(guò)直接計(jì)算：知函數(shù)

的梯度為：定義

元二次函數(shù)

為：3.3

不定常迭代法3.3.1最速下降法令

為方程組的真解，則有對(duì)于一切

，

，且對(duì)于一切

有

3.3

不定常迭代法3.3.1最速下降法【定理3.16】設(shè)