【正文】 1 第四章 數(shù)值計(jì)算 引言 本章將花較大的篇幅討論若干常見(jiàn)數(shù)值計(jì)算問(wèn)題：線性分析、一元和多元函數(shù)分析、微積分、數(shù)據(jù)分析、以及常微分方程（初值和邊值問(wèn)題）求解等。但與一般數(shù)值計(jì)算教科書不同，本章的討論重點(diǎn)是：如何利用現(xiàn)有的世界頂級(jí)數(shù)值計(jì)算資源 MATLAB。至于數(shù)學(xué)描述，本章將遵循“最低限度自封閉”的原則處理，以最簡(jiǎn)明的方式闡述理論數(shù)學(xué)、數(shù)值數(shù)學(xué)和 MATLAB計(jì)算指令之間的內(nèi)在聯(lián)系及區(qū)別。 對(duì)于那些熟悉其他高級(jí)語(yǔ)言（如 FORTRAN， Pascal， C++）的讀者來(lái)說(shuō)，通過(guò)本章，MATLAB卓越的數(shù)組處理能力、 浩瀚而靈活的 M函數(shù)指令、豐富而友善的圖形顯示指令將使他們體驗(yàn)到解題視野的豁然開朗，感受到擺脫煩瑣編程后的眉眼舒展。 對(duì)于那些經(jīng)過(guò)大學(xué)基本數(shù)學(xué)教程的讀者來(lái)說(shuō)，通過(guò)本章， MATLAB 精良完善的計(jì)算指令，自然易讀的程序?qū)⑹顾麄兏形颉敖坛獭睌?shù)學(xué)的基礎(chǔ)地位和局限性，看到從“理想化”簡(jiǎn)單算例通向科學(xué)研究和工程設(shè)計(jì)實(shí)際問(wèn)題的一條途徑。 對(duì)于那些熟悉 MATLAB基本指令的讀者來(lái)說(shuō)，通過(guò)本章，圍繞基本數(shù)值問(wèn)題展開的內(nèi)容將使他們體會(huì)到各別指令的運(yùn)用場(chǎng)合和內(nèi)在關(guān)系，獲得綜合運(yùn)用不同指令解決具體問(wèn)題的思路和借鑒。 由于 MATLAB的基本運(yùn)算單元是數(shù)組，所以本章內(nèi)容將從矩陣分析、線性代數(shù)的數(shù)值計(jì)算開始。然后再介紹函數(shù)零點(diǎn)、極值的求取，數(shù)值微積分，數(shù)理統(tǒng)計(jì)和分析，擬合和插值， Fourier 分析，和一般常微分方程初值、邊值問(wèn)題。本章的最后討論稀疏矩陣的處理，因?yàn)檫@只有在大型問(wèn)題中，才須特別處理。 從總體上講，本章各節(jié)之間沒(méi)有依從關(guān)系，即讀者沒(méi)有必要從頭到尾系統(tǒng)閱讀本章內(nèi)容。讀者完全可以根據(jù)需要閱讀有關(guān)節(jié)次。除特別說(shuō)明外，每節(jié)中的例題指令是獨(dú)立完整的，因此讀者可以很容易地在自己機(jī)器上實(shí)踐。 MATLAB從 ， 本章內(nèi)容的變化如下： ? MATLAB從 ，其矩陣和特征值計(jì)算指令不再以 LINPACK和 EISPACK 庫(kù)為基礎(chǔ)，而建筑在計(jì)算速度更快、運(yùn)行更可靠的 LAPACK和 ARPACK程序庫(kù)的新基礎(chǔ)上。因此，雖然各種矩陣計(jì)算指令沒(méi)有變化，但計(jì)算結(jié)果卻可能有某些不同。這尤其突出地表現(xiàn)在涉及矩陣分解、特征向量、奇異向量等的計(jì)算結(jié)果上。對(duì)此，用戶不必詫異，因?yàn)闃?gòu)成空間的基向量時(shí)不唯一的，且新版的更可信。本書新版全部算例結(jié)果是在 。 ? 在 ，泛函指令對(duì)被處理函數(shù)的調(diào)用是借助函數(shù)名字符串進(jìn)行的。這種調(diào)用方式在 “過(guò)渡期內(nèi)允許使用但即將被淘汰的調(diào)用方式”；而新的調(diào)用方式是借助“函數(shù)句柄”進(jìn)行的。因此，關(guān)于述泛函指令，本章新版著重講述如何使用“函數(shù)句柄”，同時(shí)兼顧“函數(shù)名字符串”調(diào)用法。 ? MATLAB從 ，提供了一組專門求微分方程“邊值問(wèn)題”數(shù)值解的指令。適應(yīng)這種變化，本章新增第 ，用 2個(gè)算例闡述求解細(xì)節(jié)。 ? quad8已經(jīng)廢止； quadl ； 分指令 triplequad。本章新版對(duì)此作了相應(yīng)的改變。 LU 分解和恰定 方程組的解 2 LU 分解、行列式和逆 恰定方程組的解 【例 】“求逆”法和“左除”法解恰定方程的性能對(duì)比 （ 1） randn(39。state39。,0)。 A=gallery(39。randsvd39。,100,2e13,2)。 x=ones(100,1)。 b=A*x。 cond(A) ans = +013 （ 2） tic xi=inv(A)*b。 ti=toc eri=norm(xxi) rei=norm(A*xib)/norm(b) ti = eri = rei = （ 3） tic。xd=A\b。 td=toc,erd=norm(xxd),red=norm(A*xdb)/norm(b) td = 0 erd = red = 范數(shù)、條件數(shù)和方程解的精度 【例 】 Hilbert矩陣是著名的病態(tài)矩陣。 MATLAB中有專門的 Hilbert矩陣及其準(zhǔn)確逆矩陣的生成函數(shù)。本例將對(duì)方程 bHx? 近 似解和準(zhǔn)確解進(jìn)行比較。 N=[6 8 10 12 14]。 for k=1:length(N) n=N(k)。 H=hilb(n)。 Hi=invhilb(n)。 b=ones(n,1)。 x_approx=H\b。 x_exact=Hi*b。 ndb=norm(H*x_approxb)。nb=norm(b)。 ndx=norm(x_approx x_exact)。nx=norm(x_approx)。 er_actual(k)=ndx/nx。 K=cond(H)。 3 er_approx(k)=K*eps。 er_max(k)=K*ndb/nb。 end disp(39。Hilbert 矩陣階數(shù) 39。),disp(N) format short e disp(39。實(shí)際誤差 er_actual39。),disp(er_actual),disp(39。39。) disp(39。近似的最大可能誤差 er_approx39。),disp(er_approx),disp(39。39。) disp(39。最大可能誤差 er_max39。),disp(er_max),disp(39。39。) Hilbert矩陣階數(shù) 6 8 10 12 14 實(shí)際誤差 er_actual +000 近似的最大可能誤差 er_approx +000 +001 最大可能誤差 er_max +003 +007 +010 矩陣特征值和矩陣函數(shù) 特征值和特征向量的求取 【例 】簡(jiǎn)單實(shí)陣的特征值問(wèn)題。 A=[1,3。2,2/3]。[V,D]=eig(A) V = + D = + 0 0 【例 】本例演示：如矩陣中有元素與截?cái)嗾`差相當(dāng)時(shí)的特性值問(wèn)題。 A=[3 2 2*eps 2 4 1 eps eps/4 eps/2 1 0 1 ]。 [V1,D1]=eig(A)。ER1=A*V1V1*D1 [V2,D2]=eig(A,39。nobalance39。)。ER2=A*V2V2*D2 ER1 = 0 ER2 = * 0 【例 】指令 eig 與 eigs 的比較。 4 rand(39。state39。,1),A=rand(100,100)。 t0=clock。[V,D]=eig(A)。T_full=etime(clock,t0) =1e8。 =0。 t0=clock。[v,d]=eigs(A,1,39。lr39。,options)。 T_part=etime(clock,t0) [Dmr,k]=max(real(diag(D)))。 d,D(1,1) T_full = T_part = d = + ans = + vk1=V(:,k)。 vk1=vk1/norm(vk1)。v=v/norm(v)。 V_err=acos(norm(vk139。*v))*180/pi D_err=abs(D(k,k)d)/abs(d) V_err = D_err = 特征值問(wèn)題的條件數(shù) 【例 】矩陣的代數(shù) 方程條件數(shù)和特征值條件數(shù)。 B=eye(4,4)。B(3,4)=1。B format short e,c_equ=cond(B),c_eig=condeig(B) B = 1 0 0 0 0 1 0 0 0 0 1 1 0 0 0 1 c_equ = +000 Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = . In D:\MATLAB6P1\toolbox\matlab\matfun\ at line 30 c_eig = +000 +000 +015 +015 【例 】對(duì)虧損矩陣進(jìn)行 Jordan分解。 A=gallery(5) [VJ,DJ]=jordan(A)。 [V,D,c_eig]=condeig(A)。c_equ=cond(A)。 DJ,D,c_eig,c_equ A = 9 11 21 63 252 5 70 69 141 421 1684 575 575 1149 3451 13801 3891 3891 7782 23345 93365 1024 1024 2048 6144 24572 DJ = 0 1 0 0 0 0 0 1 0 0 0 0 0 1 0 0 0 0 0 1 0 0 0 0 0 D = Columns 1 through 4 0 0 0 0 + 0 0 0 0 0 0 0 0 + 0 0 0 0 Column 5 0 0 0 0 c_eig = +010 * c_equ = +018 復(fù)數(shù)特征值對(duì)角陣與實(shí)數(shù)塊特征值對(duì)角陣的轉(zhuǎn)化 【例 】把例 中的復(fù)數(shù)特征值對(duì)角陣 D 轉(zhuǎn)換成實(shí)數(shù)塊對(duì)角陣，使VR*DR/VR=A。 A=[1,3。2,2/3]。[V,D]=eig(A)。 [VR,DR]=cdf2rdf(V,D) VR = 0 DR = 矩陣的譜分解和矩陣函數(shù) 【例 】數(shù)組乘方與矩陣乘方的比較。 clear,A=[1 2 3。4 5 6。7 8 9]。 A_Ap=A.^ A_Mp=A^ A_Ap = 6 A_Mp = + + + + + 【例 2】標(biāo)量的數(shù)組乘方和矩陣乘方的比較。（ A取自例 ） pA_A=().^A pA_M=()^A pA_A = pA_M = 【例 】 sin的數(shù)組運(yùn)算和矩陣運(yùn)算比較 。（ A取自例 ） A_sinA=sin(A) A_sinM=funm(A,39。sin39。) A_