1、“數(shù)值計(jì)算方法”上機(jī)實(shí)驗(yàn)指導(dǎo)書實(shí)驗(yàn)一 誤差分析實(shí)驗(yàn)1.1（病態(tài)問題）實(shí)驗(yàn)?zāi)康模核惴ㄓ小皟?yōu)”與“劣”之分，問題也有“好”與“壞”之別。對(duì)數(shù)值方法的研究而言，所謂壞問題就是問題本身對(duì)擾動(dòng)敏感者，反之屬于好問題。通過本實(shí)驗(yàn)可獲得一個(gè)初步體會(huì)。數(shù)值分析的大部分研究課題中，如線性代數(shù)方程組、矩陣特征值問題、非線性方程及方程組等都存在病態(tài)的問題。病態(tài)問題要通過研究和構(gòu)造特殊的算法來解決，當(dāng)然一般要付出一些代價(jià)（如耗用更多的機(jī)器時(shí)間、占用更多的存儲(chǔ)空間等）。問題提出：考慮一個(gè)高次的代數(shù)多項(xiàng)式顯然該多項(xiàng)式的全部根為1,2,20共計(jì)20個(gè)，且每個(gè)根都是單重的。現(xiàn)考慮該多項(xiàng)式的一個(gè)擾動(dòng)其中是一個(gè)非常小的數(shù)。這相當(dāng)

2、于是對(duì)（1.1）中的系數(shù)作一個(gè)小的擾動(dòng)。我們希望比較（1.1）和（1.2）根的差別，從而分析方程（1.1）的解對(duì)擾動(dòng)的敏感性。實(shí)驗(yàn)內(nèi)容：為了實(shí)現(xiàn)方便，我們先介紹兩個(gè)MATLAB函數(shù)：“roots”和“poly”。其中若變量a存儲(chǔ)n+1維的向量，則該函數(shù)的輸出u為一個(gè)n維的向量。設(shè)a的元素依次為，則輸出u的各分量是多項(xiàng)式方程的全部根；而函數(shù) 的輸出b是一個(gè)n+1維向量，它是以n維向量v的各分量為根的多項(xiàng)式的系數(shù)?？梢姟皉oots”和“poly”是兩個(gè)互逆的運(yùn)算函數(shù)。上述簡單的MATLAB程序便得到（1.2）的全部根，程序中的“ess”即是（1.2）中的。實(shí)驗(yàn)要求：（1） 選擇充分小的ess，反復(fù)

3、進(jìn)行上述實(shí)驗(yàn)，記錄結(jié)果的變化并分析它們。如果擾動(dòng)項(xiàng)的系數(shù)很小，我們自然感覺（1.1）和（1.2）的解應(yīng)當(dāng)相差很小。計(jì)算中你有什么出乎意料的發(fā)現(xiàn)？表明有些解關(guān)于如此的擾動(dòng)敏感性如何？（2） 將方程（1.2）中的擾動(dòng)項(xiàng)改成或其它形式，實(shí)驗(yàn)中又有怎樣的現(xiàn)象出現(xiàn)？（3） （選作部分）請(qǐng)從理論上分析產(chǎn)生這一問題的根源。注意我們可以將方程（1.2）寫成展開的形式， 同時(shí)將方程的解x看成是系數(shù)的函數(shù)，考察方程的某個(gè)解關(guān)于的擾動(dòng)是否敏感，與研究它關(guān)于的導(dǎo)數(shù)的大小有何關(guān)系？為什么？你發(fā)現(xiàn)了什么現(xiàn)象，哪些根關(guān)于的變化更敏感？思考題一：（上述實(shí)驗(yàn)的改進(jìn)）在上述實(shí)驗(yàn)中我們會(huì)發(fā)現(xiàn)用roots函數(shù)求解多項(xiàng)式方程的精度不高

4、，為此你可以考慮用符號(hào)函數(shù)solve來提高解的精確度，這需要用到將多項(xiàng)式轉(zhuǎn)換為符號(hào)多項(xiàng)式的函數(shù)poly2sym,函數(shù)的具體使用方法可參考MATLAB的幫助。思考題二：（二進(jìn)制產(chǎn)生的誤差）用MATLAB計(jì)算。結(jié)果居然有誤差！因?yàn)閺氖M(jìn)制數(shù)角度分析，這一計(jì)算應(yīng)該是準(zhǔn)確的。實(shí)驗(yàn)反映了計(jì)算機(jī)內(nèi)部的二進(jìn)制本質(zhì)。思考題三：（一個(gè)簡單公式中產(chǎn)生巨大舍入誤差的例子）可以用下列式子計(jì)算自然對(duì)數(shù)的底數(shù) 這個(gè)極限表明隨著n的增加，計(jì)算e值的精度是不確定的?，F(xiàn)編程計(jì)算與exp(1)值的差。n大到什么程度的時(shí)候誤差最大？你能解釋其中的原因嗎？相關(guān)MATLAB函數(shù)提示：poly(a) 求給定的根向量a生成其對(duì)應(yīng)的多項(xiàng)式系

5、數(shù)（降序）向量roots(p) 求解以向量p為系數(shù)的多項(xiàng)式（降序）的所有根poly2sym(p) 將多項(xiàng)式向量p表示成為符號(hào)多項(xiàng)式（降序）sym(arg) 將數(shù)字、字符串或表達(dá)式arg轉(zhuǎn)換為符號(hào)對(duì)象syms arg1 arg2 argk 將字符arg1,arg2,argk定義為基本符號(hào)對(duì)象solve(eq1) 求符號(hào)多項(xiàng)式方程eq1的符號(hào)解實(shí)驗(yàn)二 插值法實(shí)驗(yàn)2.1（多項(xiàng)式插值的振蕩現(xiàn)象）問題提出：考慮一個(gè)固定的區(qū)間上用插值逼近一個(gè)函數(shù)。顯然拉格朗日插值中使用的節(jié)點(diǎn)越多，插值多項(xiàng)式的次數(shù)就越高。我們自然關(guān)心插值多項(xiàng)式的次數(shù)增加時(shí)，是否也更加靠近被逼近的函數(shù)。龍格（Runge）給出一個(gè)例子是極著名

6、并富有啟發(fā)性的。設(shè)區(qū)間-1,1上函數(shù) 實(shí)驗(yàn)內(nèi)容：考慮區(qū)間-1,1的一個(gè)等距劃分，分點(diǎn)為 則拉格朗日插值多項(xiàng)式為 其中的是n次拉格朗日插值基函數(shù)。實(shí)驗(yàn)要求：（1） 選擇不斷增大的分點(diǎn)數(shù)目n=2,3.，畫出原函數(shù)f(x)及插值多項(xiàng)式函數(shù)在-1，1上的圖像，比較并分析實(shí)驗(yàn)結(jié)果。（2）選擇其他的函數(shù)，例如定義在區(qū)間-5，5上的函數(shù)重復(fù)上述的實(shí)驗(yàn)看其結(jié)果如何。（3）區(qū)間a,b上切比雪夫點(diǎn)的定義為 以為插值節(jié)點(diǎn)構(gòu)造上述各函數(shù)的拉格朗日插值多項(xiàng)式，比較其結(jié)果,試分析原因。實(shí)驗(yàn)2.2（樣條插值的收斂性）問題提出：多項(xiàng)式插值是不收斂的，即插值的節(jié)點(diǎn)多，效果不一定就好。對(duì)樣條函數(shù)插值又如何呢？理論上證明樣條插值的

7、收斂性是比較困難的，但通過本實(shí)驗(yàn)可以驗(yàn)證這一理論結(jié)果。實(shí)驗(yàn)內(nèi)容：請(qǐng)按一定的規(guī)則分別選擇等距或者非等距的插值節(jié)點(diǎn)，并不斷增加插值節(jié)點(diǎn)的個(gè)數(shù)?？紤]實(shí)驗(yàn)2.1中的函數(shù)或選擇其他你有興趣的函數(shù)。實(shí)驗(yàn)要求：（1）隨節(jié)點(diǎn)個(gè)數(shù)增加，比較被逼近函數(shù)和樣條插值函數(shù)誤差的變化情況。分析所得結(jié)果并與拉格朗日多項(xiàng)式插值比較(可以用MATLAB的函數(shù)“spline”作此函數(shù)的三次樣條插值，取n=10、20，分別畫出插值函數(shù)及原函數(shù)的圖形)。（2）樣條插值的思想是早產(chǎn)生于工業(yè)部門。作為工業(yè)應(yīng)用的例子考慮如下問題：某汽車制造商用三次樣條插值設(shè)計(jì)車門的曲線，其中一段的數(shù)據(jù)如下：xk012345678910yk0.00.791

8、.532.192.713.033.272.893.063.193.29yk0.80.2要求：i.自己編程計(jì)算（用三彎矩、三轉(zhuǎn)角方程均可） ii.主函數(shù)myspline(x,y,邊界類型，邊界值，xi ) 其中：x 節(jié)點(diǎn) y 節(jié)點(diǎn)上的函數(shù)值 xi 未知節(jié)點(diǎn) 返回：S(xi) iii.三對(duì)角方程組用追趕法求解(書P160)。實(shí)驗(yàn)2.3：(一維插值的應(yīng)用畫圖)畫你自己的手的形狀，在MATLAB中輸入figure(position,get(0,screensize)axes(position,0 0 1 1)x,y=ginput;將你的手掌張開放在計(jì)算機(jī)屏幕上，然后使用計(jì)算機(jī)鼠標(biāo)選取一系列點(diǎn)勾勒出手的

9、輪廓，按回車鍵結(jié)束ginput過程，這樣就獲得了一系列你的手掌外形數(shù)據(jù)點(diǎn)。也可以這樣獲得數(shù)據(jù)點(diǎn)，先把手放在一張白紙上，并用筆畫出它的輪廓線，然后將紙貼在計(jì)算機(jī)屏幕上，透過紙能看到平面上的鼠標(biāo)，并通過ginput記錄下輪廓上的點(diǎn)。將和坐標(biāo)值看作是兩個(gè)獨(dú)立變量的函數(shù)，獨(dú)立變量的取值為從1到記錄的點(diǎn)的數(shù)目。利用MATLAB的插值函數(shù)進(jìn)行插值，并畫出你的手掌外形輪廓。思考題：（二維插值）1. 在一丘陵地帶測(cè)量高程，x和y方向每隔100米測(cè)一個(gè)點(diǎn)，得高程數(shù)據(jù)如下。試用MATLAB的二維插值函數(shù)“interp2”進(jìn)行插值，并由此找出最高點(diǎn)和該點(diǎn)的高程。yx1002003004001006366976244

10、782006987126304783006806745984124006626265523342. 畫出山區(qū)地貌圖。利用MATLAB的peaks函數(shù)生成某山區(qū)的一些地點(diǎn)及其高度三維數(shù)據(jù)（單位：m）。命令格式：x,y,z=peaks(n)，生成的n階矩陣x,y,z為測(cè)量的山區(qū)地點(diǎn)三維數(shù)據(jù)（n=30）。根據(jù)peaks函數(shù)生成的數(shù)據(jù)，利用Matlab二維插值畫出該山區(qū)的地貌圖和等值線圖(提示函數(shù)：interp2、meshgrid、plot3等)。相關(guān)MATLAB函數(shù)提示：plot(x,y) 作出以數(shù)據(jù)(x(i),y(i)為節(jié)點(diǎn)的折線圖，其中x,y為同長度的向量subplot(m,n,k) 將圖形窗口

11、分為m*n個(gè)子圖，并指向第k幅圖yi=interp1(x,y,xi) 根據(jù)數(shù)據(jù)(x,y)給出在xi的分段線性插值結(jié)果yipp=spline(x,y) 返回樣條插值的分段多項(xiàng)式(pp)形式結(jié)構(gòu)pp=csape(x,y,邊界類型,邊界值) 生成各種邊界條件的三次樣條插值yi=ppval(pp,xi) pp樣條在xi的函數(shù)值ZI=interp2(x,y,z,xi,yi) x,xi為行向量，y,yi為列向量，z為矩陣的雙線性二維插值ZI=interp2(,spline) 使用二元三次樣條插值ZI=griddata(x,y,z,xi,yi) x,y,z均為向量（不必單調(diào)）表示數(shù)據(jù)，xi,yi為網(wǎng)格向量的

12、三角形線性插值（不規(guī)則數(shù)據(jù)的二維插值）實(shí)驗(yàn)三 函數(shù)逼近與曲線擬合實(shí)驗(yàn)3.1（曲線逼近方法的比較）問題提出：曲線的擬合和插值，是逼近函數(shù)的基本方法，每種方法具有各自的特點(diǎn)和特定的適用范圍，實(shí)際工作中合理選擇方法是重要的。實(shí)驗(yàn)內(nèi)容：考慮實(shí)驗(yàn)2.1中的著名問題。下面的MATLAB程序給出了該函數(shù)的二次和三次擬合多項(xiàng)式。x=-1:0.2:1;y=1/(1+25*x.*x);xx=-1:0.02:1;p2=polyfit(x,y,2);yy=polyval(p2,xx);plot(x,y,o,xx,yy);xlabel(x);ylabel(y);hold on;p3=polyfit(x,y,3);yy=

13、polyval(p3,xx);plot(x,y,o,xx,yy);hold off;適當(dāng)修改上述MATLAB程序，也可以擬合其他你有興趣的函數(shù)。實(shí)驗(yàn)要求：（1）將擬合的結(jié)果與拉格朗日插值及樣條插值的結(jié)果比較。（2）歸納總結(jié)數(shù)值實(shí)驗(yàn)結(jié)果，試定性地說明函數(shù)逼近各種方法的適用范圍，及實(shí)際應(yīng)用中選擇方法應(yīng)注意的問題。實(shí)驗(yàn)3.2：（最小二乘擬合的經(jīng)驗(yàn)公式和模型）1.（已知經(jīng)驗(yàn)公式）：某類疾病發(fā)病率為和年齡段(每五年為一段，例如05歲為第一段，610歲為第二段)之間有形如的經(jīng)驗(yàn)關(guān)系，觀測(cè)得到的數(shù)據(jù)表如下1234567890.8982.383.071.842.021.942.222.774.02101112

14、131415161718194.765.466.5310.916.522.535.750.661.681.8實(shí)驗(yàn)要求：（1）用最小二乘法確定模型中的參數(shù)和(提示函數(shù)：lsqcurvefit，lsqnonlin)。（2）利用MATLAB畫出離散數(shù)據(jù)及擬合函數(shù)圖形。（3）利用MATLAB畫出離散點(diǎn)處的誤差圖，并計(jì)算相應(yīng)的均方誤差。2（最小二乘擬合模型未知）某年美國轎車價(jià)格的調(diào)查資料如表，其中表示轎車的使用年數(shù)，表示相應(yīng)的平均價(jià)格，實(shí)驗(yàn)要求：試分析用什么形式的曲線來擬合表中的數(shù)據(jù)，并預(yù)測(cè)使用4.5年后轎車的平均價(jià)格大致為多少？1 2 3 4 5 6 7 8 9 102615 1943 1494 10

15、87 765 538 484 290 226 204實(shí)驗(yàn)3.3（研究最佳平方逼近多項(xiàng)式的收斂性質(zhì)）實(shí)驗(yàn)內(nèi)容和要求：取函數(shù)，在-1，1上以勒讓德多項(xiàng)式為基函數(shù)，對(duì)于構(gòu)造最佳平方逼近多項(xiàng)式，令，將的曲線畫在一個(gè)圖上。令，畫出的曲線。做出之間的最小二乘曲線，能否提出關(guān)于收斂性的猜測(cè)。思考題一：（病態(tài)）考慮將0,130等分節(jié)點(diǎn)，用多項(xiàng)式生成數(shù)據(jù)，再用polyfit求其3次、5次、10次、15次擬合多項(xiàng)式，并分析誤差產(chǎn)生的原因。思考題二：調(diào)研Matlab擬合工具箱的使用。相關(guān)MATLAB函數(shù)提示：p=polyfit(x,y,k) 用k次多項(xiàng)式擬合向量數(shù)據(jù)(x,y)，返回多項(xiàng)式的降冪系數(shù)，當(dāng)k=n-1時(shí)，

16、實(shí)現(xiàn)多項(xiàng)式插值，這里n是向量維數(shù)實(shí)驗(yàn)四 數(shù)值積分與數(shù)值微分實(shí)驗(yàn)4.1 （高斯數(shù)值積分方法用于積分方程求解）問題提出：線性的積分方程的數(shù)值求解，可以被轉(zhuǎn)化為線性代數(shù)方程組的求解問題。而線性代數(shù)方程組所含未知數(shù)的個(gè)數(shù)，與用來離散積分的數(shù)值方法的節(jié)點(diǎn)個(gè)數(shù)相同。在節(jié)點(diǎn)數(shù)相同的前提下，高斯數(shù)值積分方法有較高的代數(shù)精度，用它通常會(huì)得到較好的結(jié)果。實(shí)驗(yàn)內(nèi)容：求解第二類Fredholm積分方程首先將積分區(qū)間a,b等分成n份，在每個(gè)子區(qū)間上離散方程中的積分就得到線性代數(shù)方程組。實(shí)驗(yàn)要求：分別使用如下方法，離散積分方程中的積分1.復(fù)化梯形方法；2.復(fù)化辛甫森方法；3.復(fù)化高斯方法。求解如下的積分方程。1），方程的

17、準(zhǔn)確解為；2），方程的準(zhǔn)確解為；比較各算法的計(jì)算量和誤差以分析它們的優(yōu)劣。實(shí)驗(yàn)4.2（高維積分?jǐn)?shù)值計(jì)算的蒙特卡羅方法）問題提出：高維空間中的積分，如果維數(shù)不很高且積分區(qū)域是規(guī)則的或者能等價(jià)地寫成多重積分的形式，可以用一元函數(shù)積分的數(shù)值方法來計(jì)算高維空間的積分。蒙特卡羅方法對(duì)計(jì)算復(fù)雜區(qū)域甚至不連通的區(qū)域上的積分并沒有特殊的困難。實(shí)驗(yàn)內(nèi)容：對(duì)于一般的區(qū)域，計(jì)算其測(cè)度（只要理解為平面上的面積或空間中的體積）的一般方法是：先找一個(gè)規(guī)則的區(qū)域A包含，且A的測(cè)度是已知的。生成區(qū)域A中m個(gè)均勻分布的隨機(jī)點(diǎn)，如果其中有n個(gè)落在區(qū)域中，則區(qū)域的測(cè)度m()為n/m。函數(shù)f(x)在區(qū)域上的積分可以近似為：區(qū)域的測(cè)度

18、與函數(shù)f(x)在中n個(gè)隨機(jī)點(diǎn)上平均值的乘積，即實(shí)驗(yàn)要求：假設(shè)冰琪淋的下部為一錐體而上面為一半球，考慮冰琪淋體積問題：計(jì)算錐面上方和球面內(nèi)部區(qū)域的體積。如果使用球面坐標(biāo)，該區(qū)域可以表示為如下的積分：用蒙特卡羅方法可以計(jì)算該積分。另一方面，顯然這樣的冰琪淋可以裝在如下立方體的盒子里而該立方體的體積為8。只要生成這個(gè)盒子里均勻分布的隨機(jī)點(diǎn)，落入冰琪淋錐點(diǎn)的個(gè)數(shù)與總點(diǎn)數(shù)之比再乘以8就是冰琪淋錐的體積。比較兩種方法所得到的結(jié)果。類似的辦法可以計(jì)算復(fù)雜區(qū)域的測(cè)度（面積或體積）。試求由下列關(guān)系所界定區(qū)域的測(cè)度：相關(guān)MATLAB函數(shù)提示：diff(x) 如果x是向量，返回向量x的差分；如果x是矩陣，則按各列作

19、差分diff(x,k) k階差分q=polyder(p) 求得由向量p表示的多項(xiàng)式導(dǎo)函數(shù)的向量表示qFx=gradient(F,x) 返回向量F表示的一元函數(shù)沿x方向的導(dǎo)函數(shù)F(x)，其中x是與F同維數(shù)的向量z=trapz(x,y) x表示積分區(qū)間的離散化向量；y是與x同維數(shù)的向量，表示被積函數(shù)；z返回積分的近似值z(mì)=guad(fun,a,b,tol) 自適應(yīng)步長Simpson積分法求得Fun在區(qū)間a,b上的定積分，F(xiàn)un為M文件函數(shù)句柄，tol為積分精度z=dblquad(fun,a,b,c,d,tol,method) 求得二元函數(shù)Fun(x,y)的重積分z=triplequad(fun,a

20、,b,c,d,e,f,tol,method) 求得三元函數(shù)Fun(x,y,z)的重積分實(shí)驗(yàn)五 解線性方程組的直接方法實(shí)驗(yàn)5.1 （主元的選取與算法的穩(wěn)定性）問題提出：Gauss消去法是我們?cè)诰€性代數(shù)中已經(jīng)熟悉的。但由于計(jì)算機(jī)的數(shù)值運(yùn)算是在一個(gè)有限的浮點(diǎn)數(shù)集合上進(jìn)行的，如何才能確保Gauss消去法作為數(shù)值算法的穩(wěn)定性呢？Gauss消去法從理論算法到數(shù)值算法，其關(guān)鍵是主元的選擇。主元的選擇從數(shù)學(xué)理論上看起來平凡，它卻是數(shù)值分析中十分典型的問題。實(shí)驗(yàn)內(nèi)容：考慮線性方程組 編制一個(gè)能自動(dòng)選取主元，又能手動(dòng)選取主元的求解線性方程組的Gauss消去過程。實(shí)驗(yàn)要求：（1）取矩陣，則方程有解。取n=10計(jì)算矩

21、陣的條件數(shù)。讓程序自動(dòng)選取主元(順序消元)，結(jié)果如何？（2）現(xiàn)選擇程序中手動(dòng)選取主元的功能。每步消去過程總選取按模最小或按模盡可能小的元素作為主元，觀察并記錄計(jì)算結(jié)果。若每步消去過程總選取按模最大的元素作為主元，結(jié)果又如何？分析實(shí)驗(yàn)的結(jié)果。（3）取矩陣階數(shù)n=20或者更大，重復(fù)上述實(shí)驗(yàn)過程，觀察記錄并分析不同的問題及消去過程中選擇不同的主元時(shí)計(jì)算結(jié)果的差異，說明主元素的選取在消去過程中的作用。（4）選取其他你感興趣的問題或者隨機(jī)生成矩陣，計(jì)算其條件數(shù)。重復(fù)上述實(shí)驗(yàn)，觀察記錄并分析實(shí)驗(yàn)結(jié)果。實(shí)驗(yàn)5.2（線性代數(shù)方程組的性態(tài)與條件數(shù)的估計(jì)）問題提出：理論上，線性代數(shù)方程組的攝動(dòng)滿足 矩陣的條件數(shù)確

22、實(shí)是對(duì)矩陣病態(tài)性的刻畫，但在實(shí)際應(yīng)用中直接計(jì)算它顯然不現(xiàn)實(shí)，因?yàn)橛?jì)算通常要比求解方程還困難。實(shí)驗(yàn)內(nèi)容：MATLAB中提供有函數(shù)“condest”可以用來估計(jì)矩陣的條件數(shù)，它給出的是按1-范數(shù)的條件數(shù)。首先構(gòu)造非奇異矩陣A和右端，使得方程是可以精確求解的。再人為地引進(jìn)系數(shù)矩陣和右端的攝動(dòng)，使得充分小。實(shí)驗(yàn)要求：（1）假設(shè)方程Ax=b的解為x，求解方程，以1-范數(shù)，給出的計(jì)算結(jié)果。（2）選擇一系列維數(shù)遞增的矩陣（可以是隨機(jī)生成的），比較函數(shù)“condest”所需機(jī)器時(shí)間的差別.考慮若干逆是已知的矩陣，借助函數(shù)“eig”很容易給出cond2(A)的數(shù)值。將它與函數(shù)“cond(A,2)”所得到的結(jié)果進(jìn)

23、行比較。（3）利用“condest”給出矩陣A條件數(shù)的估計(jì)，針對(duì)（1）中的結(jié)果給出的理論估計(jì)，并將它與（1）給出的計(jì)算結(jié)果進(jìn)行比較，分析所得結(jié)果。注意，如果給出了cond(A)和的估計(jì)，馬上就可以給出的估計(jì)。（4）估計(jì)著名的Hilbert矩陣的條件數(shù)。思考題一：（Vadermonde矩陣）設(shè) ，其中，（1）對(duì)n=2,5,8，計(jì)算A的條件數(shù)；隨n增大，矩陣性態(tài)如何變化？（2）對(duì)n=5，解方程組Ax=b；設(shè)A的最后一個(gè)元素有擾動(dòng)10-4，再求解Ax=b（3）計(jì)算（2）擾動(dòng)相對(duì)誤差與解的相對(duì)偏差，分析它們與條件數(shù)的關(guān)系。（4）你能由此解釋為什么不用插值函數(shù)存在定理直接求插值函數(shù)而要用拉格朗日或牛頓插

24、值法的原因嗎？相關(guān)MATLAB函數(shù)提示：zeros(m,n) 生成m行，n列的零矩陣ones(m,n) 生成m行，n列的元素全為1的矩陣eye(n) 生成n階單位矩陣rand(m,n) 生成m行,n列(0,1)上均勻分布的隨機(jī)矩陣diag(x) 返回由向量x的元素構(gòu)成的對(duì)角矩陣tril(A) 提取矩陣A的下三角部分生成下三角矩陣triu(A) 提取矩陣A的上三角部分生成上三角矩陣rank(A) 返回矩陣A的秩det(A) 返回方陣A的行列式inv(A) 返回可逆方陣A的逆矩陣V，D=eig(A) 返回方陣A的特征值和特征向量norm(A,p) 矩陣或向量的p范數(shù)cond(A,p) 矩陣的條件數(shù)

25、L，U，P=lu(A) 選列主元LU分解R=chol(X) 平方根分解Hi=hilb(n) 生成n階Hilbert矩陣實(shí)驗(yàn)六 解線性方程組的迭代法實(shí)驗(yàn)6.1（病態(tài)的線性方程組的求解）問題提出：理論的分析表明，求解病態(tài)的線性方程組是困難的。實(shí)際情況是否如此，會(huì)出現(xiàn)怎樣的現(xiàn)象呢？實(shí)驗(yàn)內(nèi)容：考慮方程組Hx=b的求解，其中系數(shù)矩陣H為Hilbert矩陣， 這是一個(gè)著名的病態(tài)問題。通過首先給定解（例如取為各個(gè)分量均為1）再計(jì)算出右端b的辦法給出確定的問題。實(shí)驗(yàn)要求：（1）選擇問題的維數(shù)為6，分別用Gauss消去法、列主元Gauss消去法、J迭代法、GS迭代法和SOR迭代法求解方程組，其各自的結(jié)果如何？將

26、計(jì)算結(jié)果與問題的解比較，結(jié)論如何？（2）逐步增大問題的維數(shù)（至少到100），仍然用上述的方法來解它們，計(jì)算的結(jié)果如何？計(jì)算的結(jié)果說明了什么？（3）討論病態(tài)問題求解的算法實(shí)驗(yàn)6.2 書上P211計(jì)算實(shí)習(xí)題2，其中N=100，第二小問改為用Jacobi迭代、G-S迭代、紅黑排序的G-S迭代求解，并比較他們之間的收斂速度；進(jìn)一步，用BSOR迭代求解，試找出最優(yōu)松弛因子。實(shí)驗(yàn)七 非線性方程求根實(shí)驗(yàn)7.1（迭代法、初始值與收斂性）實(shí)驗(yàn)?zāi)康模撼醪秸J(rèn)識(shí)非線性問題的迭代法與線性問題迭代法的差別，探討迭代法及初始值與迭代收斂性的關(guān)系。問題提出：迭代法是求解非線性方程的基本思想方法，與線性方程的情況一樣，其構(gòu)造方

27、法可以有多種多樣，但關(guān)鍵是怎樣才能使迭代收斂且有較快的收斂速度。實(shí)驗(yàn)內(nèi)容：考慮一個(gè)簡單的代數(shù)方程 針對(duì)上述方程，可以構(gòu)造多種迭代法，如 在實(shí)軸上取初始值x0，請(qǐng)分別用迭代（7.1）-（7.3）作實(shí)驗(yàn)，記錄各算法的迭代過程。實(shí)驗(yàn)要求：（1）取定某個(gè)初始值，分別計(jì)算（7.1）-（7.3）迭代結(jié)果，它們的收斂性如何？重復(fù)選取不同的初始值，反復(fù)實(shí)驗(yàn)。請(qǐng)自選設(shè)計(jì)一種比較形象的記錄方式（如利用MATLAB的圖形功能），分析三種迭代法的收斂性與初值選取的關(guān)系。（2）對(duì)三個(gè)迭代法中的某個(gè)，取不同的初始值進(jìn)行迭代，結(jié)果如何？試分析迭代法對(duì)不同的初值是否有差異？（3）線性方程組迭代法的收斂性是不依賴初始值選取的。

28、比較線性與非線性問題迭代的差異，有何結(jié)論和問題。相關(guān)MATLAB函數(shù)提示：x=fzero(fun,x0) 返回一元函數(shù)fun的一個(gè)零點(diǎn)，其中fun為函數(shù)句柄，x0為標(biāo)量時(shí)，返回在x0附近的零點(diǎn)；x0為向量a,b時(shí)，返回函數(shù)在a,b中的零點(diǎn)x,f,h=fsolve(fun,x0) 返回一元或多元函數(shù)x0附近fun的一個(gè)零點(diǎn)，其中fun為函數(shù)句柄，x0為迭代初值；f返回fun在x的函數(shù)值，應(yīng)該接近0； h返回值如果大于0，說明計(jì)算結(jié)果可靠，否則不可靠實(shí)驗(yàn)八 常微分方程初值問題數(shù)值解法實(shí)驗(yàn)8.1（Lorenz問題與混沌）問題提出：考慮著名的Lorenz方程 其中s,r,b為變化區(qū)域有一定限制的實(shí)參數(shù)

29、。該方程形式簡單，表面上看并無驚人之處，但由該方程提示出的許多現(xiàn)象，促使“混沌”成為數(shù)學(xué)研究的嶄新領(lǐng)域，在實(shí)際應(yīng)用中也產(chǎn)生了巨大的影響。實(shí)驗(yàn)內(nèi)容：先取定初值y0=(0,0,0),參數(shù)s=10,r=28,b=8/3,用MATLAB的數(shù)值求常微分方程函數(shù)ods45編程對(duì)（8.1）進(jìn)行求解實(shí)驗(yàn)要求：（1）對(duì)目前取定的參數(shù)值s,r和b，選取不同的初值y0進(jìn)行運(yùn)算，觀察計(jì)算的結(jié)果有什么特點(diǎn)？解的曲線是否有界？解的曲線是不是周期的或趨于某個(gè)固定點(diǎn)？（2）在問題允許的范圍內(nèi)適當(dāng)改變其中的參數(shù)值s,r,b，再選取不同的初始值y0進(jìn)行運(yùn)算，觀察并記錄計(jì)算的結(jié)果有什么特點(diǎn)？是否發(fā)現(xiàn)什么不同的現(xiàn)象？思考題一：考慮如

30、下常微分方程 先固定其中參數(shù)：r=0.9,b=3,a=1,k0=0,m=1,初始值為原點(diǎn)，t的范圍為0-500分別取不同的k1,如1，1.5，2，2.5，3，3.5等，繪圖觀察相位圖（x,x）的變化情況，能找到其中的規(guī)律嗎？試試改變其它參數(shù)呢？相關(guān)MATLAB函數(shù)提示：t,y=ode45(odefun,tspan,y0) odefun表示f(t,y)的函數(shù)句柄，t是標(biāo)量，y是標(biāo)量或向量；tspan是二維向量t0,tf，表示自變量初值t0和終值tf；y0表示初值向量，若無輸出參數(shù)，則作出圖形附錄MATLAB簡介這里介紹MATLAB一些入門知識(shí)，包括MATLAB桌面和窗口，MATLAB命令格式、數(shù)

31、據(jù)格式、數(shù)據(jù)文件和變量管理，MATLAB的數(shù)組和矩陣運(yùn)算，MATLAB的字符串、元胞和結(jié)構(gòu)等數(shù)據(jù)類型，MATLAB的程序設(shè)計(jì)方法，MATLAB作圖方法在線幫助的使用和程序文件和目錄的管理等。表一 MATLAB的基本命令主題詞含義主題詞含義format設(shè)置數(shù)據(jù)顯示格式feval函數(shù)求值who顯示變量名input提示輸入whos顯示變量信息disp輸出clear清除內(nèi)存變量tic啟動(dòng)秒表save保存工作變量到文件toc時(shí)間讀數(shù)（秒）load從文件裝載變量help幫助linspace區(qū)間等分lookfor查找length獲取數(shù)組長度type列程序清單size矩陣大小which查找文件目錄max最大值

32、double雙精度min最小值str2num字符串轉(zhuǎn)化為數(shù)值sum求和num2str數(shù)值轉(zhuǎn)化為字符串find條件檢索一、MATLAB桌面啟動(dòng)MATLAB后，就進(jìn)入MATLAB的桌面，圖1為MATLAB6.1的默認(rèn)（Default）桌面。第一行為菜單欄，第二行為工具欄，下面是三個(gè)最常用的窗口。右邊最大的是命令窗口（Command Window），左上方前臺(tái)為發(fā)行說明書窗口（Launch pad），后臺(tái)為工作空間（Workspace），左下方為命令歷史（Command History）后臺(tái)為當(dāng)前目錄（Current Directory）。1.窗口（1）命令窗口該窗口是進(jìn)行MATLAB操作最主要的窗

33、口。窗口中“”為命令輸入提示符，其后輸入運(yùn)算命令，按回車鍵就可執(zhí)行運(yùn)算，并顯示運(yùn)算結(jié)果.。圖1（2）發(fā)行說明書窗口發(fā)行說明書窗口是MATLAB所特有的，用來說明用戶所擁有的Mathworks公司產(chǎn)品的工具包、演示以及幫助信息。（3）工作空間在默認(rèn)桌面，位于左上方窗口前臺(tái)，列出內(nèi)存中MATLAB工作空間的所有變量的變量名、尺寸、字節(jié)數(shù)。用鼠標(biāo)選中變量，擊右鍵可以打開、保存、刪除、繪圖等操作。（4）當(dāng)前目錄在默認(rèn)桌面，位于左下方窗口后臺(tái)，用鼠標(biāo)點(diǎn)擊可以切換到前臺(tái)。該窗口列出當(dāng)前目錄的程序文件（.m）和數(shù)據(jù)文件（.mat）等。用鼠標(biāo)選中文件，擊右鍵可以進(jìn)行打開、運(yùn)行、刪除等操作。（5）命令歷史（Co

34、mmand History）該窗口列出在命令窗口執(zhí)行過的MATLAB命令行的歷史記錄。用鼠標(biāo)選中命令行，擊右鍵可以進(jìn)行復(fù)制、執(zhí)行（Evaluate Selection）、刪除等操作。除上述窗口外，MATLAB常用窗口還有編程器窗口、圖形窗口等。二、數(shù)據(jù)和變量1.表達(dá)式在命令窗口作一些簡單的計(jì)算，就如同使用一個(gè)功能強(qiáng)大的計(jì)算器，使用變量無須預(yù)先定義類型。例如，設(shè)球半徑為r=2，求球的體積。r=2 %表達(dá)式將2賦予變量rr= %系統(tǒng)返回r的值2 v=4/3*pi*r3 %pi為內(nèi)置常量，乘方用表示v= 33.5103幾個(gè)表達(dá)式可以寫在一行，用分號(hào)（；）或逗號(hào)（，）分割，用分號(hào)（；）使該表達(dá)式運(yùn)算結(jié)

35、果不顯示，而逗號(hào)（，）則顯示結(jié)果。也可以將一個(gè)長表達(dá)式分在幾行上寫，用三點(diǎn)（）續(xù)行。若需要修改已執(zhí)行過的命令行，可以在命令歷史中找到該命令行復(fù)制，再粘貼至命令窗口修改。也可以直接使用鍵盤調(diào)出已執(zhí)行過的命令行修改。2.數(shù)據(jù)顯示格式MATLAB默認(rèn)的數(shù)據(jù)顯示格式為短格式（short）：當(dāng)結(jié)果為整數(shù)，就作為整數(shù)顯示；當(dāng)結(jié)果是實(shí)數(shù)，以小數(shù)點(diǎn)后四位的長度顯示。若結(jié)果的有效數(shù)字超出一定范圍，以科學(xué)計(jì)數(shù)法顯示（如3.2000e-006表示）。數(shù)據(jù)顯示格式可使用命令Format改變。例如： format long;v %長格式，16位v = 33.112 format short;v %短格式v = 33.5

36、103 format rational;v %有理格式，近似分?jǐn)?shù)v = 6501/194 3.復(fù)數(shù)MATLAB中復(fù)數(shù)可以如同實(shí)數(shù)一樣，直接輸入和計(jì)算。例如： a=1+2i;b=5-4*i;c=a/bc = -0.0732 + 0.3415i4.預(yù)定義變量MATLAB有一些預(yù)定義變量（表1），啟動(dòng)時(shí)就已賦值，可以直接使用，如前我們使用的圓周率pi和虛數(shù)單位i.表1 常用預(yù)定義變量變量名說 明i或j虛數(shù)單位pi圓周率3.14159eps浮點(diǎn)數(shù)識(shí)別精度2(-52)=realmin最小正實(shí)數(shù)realmax最大正實(shí)數(shù)inf無窮大NaN沒有意義的數(shù)預(yù)定義變量在工作空間觀察不到。如果預(yù)定義變量被用戶重新賦值

37、，則原來的功能暫不能使用。當(dāng)這些用戶變量被清除（clear）或MATLAB重新啟動(dòng)后，這些功能得以恢復(fù)。5.用戶變量MATLAB變量名總以字母開頭，以字母、數(shù)字或下劃線組成，區(qū)分大小寫，有效字符長度為63個(gè)。如A,a,a1,a_b都是合法的，且a與A表示不同變量。在Command Window中使用的變量一旦被賦值，就會(huì)攜帶這個(gè)值存在于工作空間，直到被清除或被賦予新的值。ans是系統(tǒng)一個(gè)特別的變量名。若一個(gè)表達(dá)式運(yùn)算結(jié)果沒有賦予任何變量，系統(tǒng)自動(dòng)用ans存放答案。例如： A=5+4i;b=5-4*i;B=1;A*b %沒有定義A*b的輸出變量ans = 41%ans來接受計(jì)算結(jié)果，注意這是大寫

38、A與小寫b的乘積，盡管我們可以使用工作空間來查詢和清除變量，但使用下列命令方式更快捷： whos%查詢Workspace中的變量列表 Name Size Bytes Class A 1x1 16 double array (complex) B 1x1 8 double array a 1x1 16 double array (complex) ans 1x1 8 double array b 1x1 16 double array (complex) c 1x1 16 double array (complex)Grand total is 6 elements using 80 bytes

39、 A %查詢變量A的值A(chǔ) = 5.0000 + 4.0000i clear A%清除變量A A%再查詢A的值，已經(jīng)不存在了? Undefined function or variable A. clear%清除Workspace中所有變量 whos%Workspace中已沒有任何變量了 三、數(shù)組和矩陣運(yùn)算MATLAB基本數(shù)據(jù)單元是無需指定維數(shù)的數(shù)組。數(shù)組運(yùn)算是MATLAB最鮮明的特點(diǎn)，一方面可以使得計(jì)算程序簡明易讀，另一方面可以提高計(jì)算速度。1.數(shù)組的輸入最常用的數(shù)組是雙精度數(shù)值數(shù)組（double array）。一維數(shù)組相當(dāng)于向量，二維數(shù)組相當(dāng)于矩陣，一維數(shù)組可以視為二維數(shù)組的特例。二維數(shù)組的

40、第一維稱為“行”，第二維稱為“列”。MATLAB數(shù)組無需預(yù)先定義維數(shù)。直接輸入數(shù)組的元素，用中括號(hào)（）表示一個(gè)數(shù)組，同行元素間用空格或逗號(hào)分隔，不同行間用分號(hào)或回車分隔，例如： clear;a=1,2,3;4,5,6;7,8,9a = 1 2 3 4 5 6 7 8 9或 a=1 2 3%這種方式特別適用于大型矩陣4 5 67 8 9a = 1 2 3 4 5 6 7 8 9對(duì)于等差數(shù)列構(gòu)造的一維數(shù)組，可用冒號(hào)運(yùn)算生成，也可用函數(shù)linspace生成。 b=0:3:10%初值：增量：終值b = 0 3 6 9 b=0:10%增量為1可省略b = 0 1 2 3 4 5 6 7 8 9 10 b

41、=10:-3:0%遞減b = 10 7 4 1 b=linspace(0,10,4)%將區(qū)間0，10等分為4-1=3份b = 0 3.3333 6.6667 10.0000 length(b)%查詢b的長度ans = 4 b(3)%查詢b的第三個(gè)元素ans = 6.6667 b(1,end)%查詢b的首和尾元素ans = 0 10二維數(shù)組元素雙下標(biāo)編址按通常方式，單下標(biāo)編址按列排序。 size(a)%查詢數(shù)組a的尺寸ans = 3 3 a(3,2),a(6)ans = 8ans = 8 c=a(1 3,2 3)%提取a的第一、第三行和第二、第三列（分塊矩陣）c = 2 3 8 9 d=a(2,

42、:)%提取a的第二行d = 4 5 6 a(:)%將a所有元素按單下標(biāo)順序排為列向量ans = 1 4 7 2 5 8 3 6 9一些特殊的二維數(shù)組可以用函數(shù)產(chǎn)生，例如： a=zeros(2,4)%生成2行4列零矩陣a = 0 0 0 0 0 0 0 0 b=ones(1,4)%生成1行4列1矩陣b = 1 1 1 1 c=a;b%拼接c = 0 0 0 0 0 0 0 0 1 1 1 1 c(2,1)=100%修改部分元素c = 0 0 0 0 100 0 0 0 1 1 1 1 reshape(c,2,6)%按2行6列重排矩陣元素ans = 0 1 0 0 1 0 100 0 1 0 0

43、1注意：數(shù)組下標(biāo)對(duì)應(yīng)矩陣的行和列，編址一律從1開始，不能用0.矩陣輸入也可用“l(fā)oad”命令從外部數(shù)據(jù)文件導(dǎo)入2.數(shù)組運(yùn)算數(shù)組運(yùn)算是指數(shù)組對(duì)應(yīng)元素之間的運(yùn)算，也稱點(diǎn)運(yùn)算。矩陣的乘法、乘方和除法有特殊的數(shù)學(xué)含義，并不是數(shù)組對(duì)應(yīng)元素的運(yùn)算，所數(shù)組乘法、乘方和除法的運(yùn)算符前特別加了一個(gè)點(diǎn)。特別要區(qū)分?jǐn)?shù)組運(yùn)算在乘法、乘方和除法上的意義和表示上與矩陣運(yùn)算的不同。表2 數(shù)組運(yùn)算符運(yùn)算符號(hào)說明數(shù)組加與減A+B與A-B對(duì)應(yīng)元素之間加減數(shù)乘數(shù)組k*A或A*kk乘A的每個(gè)元素?cái)?shù)與數(shù)組加減k+A或k-Ak加（減）A的每個(gè)元素?cái)?shù)組乘數(shù)組A.*B點(diǎn)運(yùn)算只有點(diǎn)乘、點(diǎn)乘方、點(diǎn)除三個(gè)，表示對(duì)應(yīng)元素之間的運(yùn)算；（.*）是一個(gè)整

44、體，點(diǎn)(.)不能漏掉，（.）和（*）之間也不能有空格數(shù)組乘方A.k,k.A數(shù)除以數(shù)組k./A數(shù)組除法左除A.B,右除B./A clear;A=1 -1;0 2;B=0 1;1 -1; A.*B%注意不是A*Bans = 0 -1 0 -2 A.B,A./BWarning: Divide by zero.ans = 0 -1.0000 Inf -0.5000Warning: Divide by zero.ans = Inf -1 0 -2 A.2ans = 1 1 0 4 1./AWarning: Divide by zero.ans = 1.0000 -1.0000 Inf 0.50003.矩

45、陣運(yùn)算矩陣是一個(gè)二維數(shù)組，所以矩陣的加、減、數(shù)乘等運(yùn)算與數(shù)組運(yùn)算是一致的。但是有兩點(diǎn)需要注意：（1）對(duì)于乘法、乘方和除法等三種運(yùn)算，矩陣運(yùn)算與數(shù)組運(yùn)算的運(yùn)算符及含義不同：矩陣運(yùn)算按線性變換定義，使用通常符號(hào)；數(shù)組運(yùn)算按對(duì)應(yīng)元素運(yùn)算定義，使用點(diǎn)運(yùn)算符；（2）數(shù)與矩陣加減、矩陣除法在數(shù)學(xué)上是沒有意義的，在MATLAB中為簡便起見，定義了這兩類運(yùn)算，其含義見表3.表3 矩陣運(yùn)算符運(yùn)算符號(hào)說明轉(zhuǎn)置A加與減A+B與A-B同數(shù)組運(yùn)算數(shù)乘矩陣k*A或A*k同數(shù)組運(yùn)算矩陣乘法A*B矩陣乘方Ak數(shù)與矩陣加減k+A與k-Ak+A等價(jià)于k*ones(size(A)+A矩陣除法左除AB，右除B/A它們分別為矩陣方程A

46、X=B和XA=B的解 A=1 2;3 4;B=4 3;2 1; 100+Aans = 101 102 103 104 A*B,A.*B%注意矩陣運(yùn)算和數(shù)組運(yùn)算的區(qū)別ans = 8 5 20 13ans = 4 6 6 4 AB,B/A,A.B,B./A%注意矩陣運(yùn)算和數(shù)組運(yùn)算的區(qū)別ans = -6.0000 -5.0000 5.0000 4.0000ans = -3.5000 2.5000 -2.5000 1.5000ans = 4.0000 1.5000 0.6667 0.2500ans = 4.0000 1.5000 0.6667 0.25004.數(shù)學(xué)函數(shù)數(shù)組的數(shù)學(xué)函數(shù)也是按每個(gè)元素的運(yùn)算

47、，使用通常的函數(shù)符號(hào)，常用數(shù)學(xué)函數(shù)見表4表4 數(shù)學(xué)函數(shù)函數(shù)意義函數(shù)意義sin正弦fix向0取整cos余弦mod模余tan正切rem除法余數(shù)cot余切abs絕對(duì)值（模）asin反正弦exp指數(shù)函數(shù)acos反余弦log自然對(duì)數(shù)sqrt開方log10以10為底的對(duì)數(shù) A=4 -1;3 2; B=exp(A)B = 54.5982 0.3679 20.0855 7.3891 C=fix(B)C = 54 0 20 7 D=sin(C)D = -0.5588 0 0.9129 0.6570 E=log(D)Warning: Log of zero.E = -0.5820 + 3.1416i -Inf -

48、0.0911 -0.4201 5.關(guān)系與邏輯運(yùn)算MATLAB的關(guān)系運(yùn)算和邏輯運(yùn)算符都是對(duì)于元素的操作，其結(jié)果是特殊的邏輯數(shù)組（logical array）表5，“真”用1表示，“假”用0表示，而邏輯運(yùn)算中，所有非零元素作為1（真）處理。表5 關(guān)系運(yùn)算和邏輯運(yùn)算運(yùn)算符含義運(yùn)算符含義小于&與大于非=大于等于all= =等于any=不等于 A=-2:4,B=4:-1:-2A = -2 -1 0 1 2 3 4B = 4 3 2 1 0 -1 -2 ABans = 0 0 0 0 1 1 1 A=Bans = 0 0 0 1 0 0 0 A&B %邏輯運(yùn)算中，所有非零元素作為1（真）處理ans = 1 1 0 1