Matlab70使用說明-數(shù)值計(jì)算部分

1、代號9- 06A模塊計(jì)算機(jī)輔助設(shè)計(jì)層次基礎(chǔ)型Mat lab 7.0使用說明（數(shù)值計(jì)算部分）華中科技大學(xué)國家機(jī)械基礎(chǔ)課程教學(xué)基地2010年9月i第一部分基本操作1§1.Matlab 的使用11-1.直接輸入命令11-2.用M文件開發(fā)程序 12M文件程序的主要語句和主要函數(shù) 22-1.Matlab的數(shù)字特征22-2.主要語句32-3.常用函數(shù)42-4.幾個(gè)常用的命令5矩陣的有關(guān)計(jì)算53-1.矩陣的輸入53-2.矩陣/向量的運(yùn)算63-3.矩陣的范數(shù)63-4.向量的范數(shù)73-5.矩陣的條件數(shù)73-6.矩陣的特征值和特征向量 8§4.Matlab 繪圖94-1.繪圖的基本命令94-2

2、.圖形的交互編輯 11第二部分?jǐn)?shù)值計(jì)算 12§1 .方程求根121-1.牛頓迭代法121-2.圖解法確定迭代的初始點(diǎn) 13§2 性方程組132-1.迭代法的收斂性 132-2.線性方程組白病態(tài)問題142-3.求解線性方程組15§3 插值和擬合163-1.Lagrange 插值163-2.代數(shù)多項(xiàng)式插值173-3.插值誤差 173-4.分段線性插值183-5.數(shù)據(jù)的曲線擬合 18§4 .數(shù)值積分204-1.復(fù)合梯形求積公式204-2.復(fù)合 Simpson 求積公式20§5 .常微分方程的數(shù)值解法 215-1.Euler 方法215-2.改進(jìn)的Eu

3、ler方法235-3.四階龍格-庫塔方法24習(xí)題27一、方程求根 27二、線性方程組 27三、插值與擬合 28四、數(shù)值積分29五、常微分方程30計(jì)算方法實(shí)驗(yàn)報(bào)告 31一、方程求根 31二、線性方程組31三、插值與擬合 32四、數(shù)值積分32五、常微分方程33III第一部分基本操作§1.Matlab的使用Matlab的使用方法有兩種：（1）在 Matlab的命令窗口（ Matlab Command Windows）中直接輸入命 令，即可得到結(jié)果；（2）在Matlab的編輯窗口（ Matlab Editor）內(nèi)編寫M文件，然后在命令窗口執(zhí)行該文件，得到所需的結(jié)果。1-1.直接輸入命令在命令

4、窗口中，直接輸入命令。例如：鍵入x=3+5顯示x=8若鍵入x=3+5 ;不能直接顯示結(jié)果，還須鍵入x方可顯示x=81-2.用M文件開發(fā)程序1 .設(shè)置當(dāng)前目錄M文件分為腳本 M文件（相當(dāng)于主程序）和函數(shù) M文件（相當(dāng)于子程序或函數(shù)）。子程序或函數(shù) 調(diào)用時(shí)，須在當(dāng)前目錄內(nèi)操作，故建議將用戶創(chuàng)建的子目錄設(shè)置為當(dāng)前目錄。這樣，所有程序及命令的 操作都在用戶子目錄內(nèi)，較為方便。以下操作只須設(shè)置一次，以后再進(jìn)入Matlab系統(tǒng)時(shí)，當(dāng)前目錄即為已設(shè)定的目錄。設(shè)置當(dāng)前目錄步驟為： 將鼠標(biāo)移至Matlab圖標(biāo)，按右鍵彈出下拉菜單； 點(diǎn)擊屬性：彈出Matlab屬性窗口； 將該窗口內(nèi) 起始位置”欄中的路徑更改為用戶

5、創(chuàng)建的子目錄路徑; 點(diǎn)擊應(yīng)用：最后點(diǎn)擊確定2 .如何編寫程序?qū)τ贛atlab命令窗口的上方菜單條，點(diǎn)擊File New一 M-file ,則彈出Matlab的編輯窗口。在編輯窗口內(nèi)鍵入M文件，最后點(diǎn)擊該窗口上方菜單條中的File-Save保存文件。3 .例示一一 221）計(jì)算函數(shù) f （x,為）=2x+7& - 3x& - 5x + 3& - 12編輯窗口中輸入該函數(shù) M文彳functiony=demof_1(x1,x2) y=2*x1A2+7*x2A2-3*x1*x2-5*x1+3*x2-12;輸入完畢后存盤(默認(rèn)文件名為demof_1.m)在命令窗口中調(diào)用該函數(shù)鍵入

6、y=demof_1(2,3)顯示y=40鍵入y=demof_1(3,-5)顯示 y=196Bia2)計(jì)算一組數(shù)據(jù)x i(i = 1,2,?, ne S =1- 和S2 = 產(chǎn)nn 在編輯窗口中鍵入該函數(shù)的 M文件function s1,s2=demof_2(x)n=length(x);s1=sum(x)/n;s2=sqrt(sum(x.A2)/n);這是一個(gè)具有多個(gè)返回值的函數(shù)，調(diào)用語句的左邊須為向量。輸入完畢后保存文件(默認(rèn)文件名為demof_2.m)。 在命令窗口中調(diào)用函數(shù)輸入數(shù)組x=1,2,3,4,5輸入調(diào)用t1,t2=demof_2(x)顯示t1=3t2=3.3166或者，通過主程序調(diào)

7、用該函數(shù)i) 首先，在編輯窗口輸入腳本M文件如下：x=input('x= =');s1,s2=demof_2(x);fprintf('s1=%12.5f s2=%12.5f ' ,s1,s2);ii) 保存該文件，備下次使用。iii)運(yùn)行該腳本M文件，即在命令窗口點(diǎn)擊File一Open,調(diào)用該文件(該文件在編輯窗口)；點(diǎn)擊Debug一Run,在命令窗口顯示待輸入的數(shù)值變量提示x=輸入一組數(shù)據(jù)，如 2,-3,0,9,13,55則顯示結(jié)果s1=12.66667s2=23.409402 M文件程序的主要語句和主要函數(shù)2-1.Matlab的數(shù)字特征1 .數(shù)字類型Matl

8、ab中所有變量均為雙精度，整型和實(shí)型之間沒有差別。幾個(gè)特殊的數(shù)字，如pi表示兀的值、inf表示8、eps表示計(jì)算機(jī)精度 2.2204 X10-16等。2 .算術(shù)運(yùn)算符加減乘除乘方+-*/A3 .邏輯比較符大于小于大于等于小于等于andorif語句中的不等于if語句中的等于><>=<=&|=4 .數(shù)組變量（1）數(shù)組的表示 一維數(shù)組（等同于向量），例如：x=1,3,-4,7,21y=3：7（等同于 y=3,4,5,6,7）z=3:0.5:6（等同于 z=3,3.5,4,4.5,5,5.5,6）v='g','k','m'

9、二維數(shù)組（等同于矩陣），例如一個(gè)3X2數(shù)組：m=0.1,0.2,0.3;0.7,0.8,0.9（2 ）數(shù)組的元素對應(yīng)于上述數(shù)組，x（2）=3, y（3）=5,m（2,1）=0.7等等。二維數(shù)組的整行或整列可以一個(gè)冒號表示，例如：m（1,:）=0.1 0.2 0.3,m（:,2）=0.20.8（3 ）數(shù)組的運(yùn)算兩數(shù)組的加（+）、減（-）運(yùn)算符與向量的運(yùn)算相同，而乘（.*）、除（./）、乘方（人）運(yùn)算符不同。例如,對于數(shù)組a和b相加Z=a+ba和b的對應(yīng)兀素相加相減Z=a-ba和b的對應(yīng)兀素相減相乘Z=a.*ba和b的對應(yīng)兀素相乘相除Z=a./ba和b的對應(yīng)兀素相除乘方Z=a.Al.3a的所白兀素

10、的1.3次方2-2.主要語句1.if-end 語句每個(gè)if語句必須以一個(gè)end結(jié)束，二者一一對應(yīng)。例如: if n= =2,price=17; end if n %price=15;else price=12; end if a>5,tap=10;elseif a<5,tap=-10;（如有必要，可多次使用elseif）else tap=1; end2. for-end 語句 for x=1:0.5:9 y=xA2-5*x-3;end上述語句一次計(jì)算 x=1, 1.5, 2,，9時(shí)y的值。若改為for x=-2,0,15,則依次計(jì)算 x=-2,0,1 5時(shí) 的值。 for x=0:

11、0.1:10 y=sin（x）; if y<0,y=0;（循環(huán)中可以插入if-end語句）endend3. while-end 語句 i=1;c=0;x=-8,0,12,33,-6,5,-7; while i<=length(x)if x(i)<0,c=c+1; end i=i+1;endfprintf('C=%d',c);上述語句為統(tǒng)計(jì)數(shù)組x中小于0的元素個(gè)數(shù)。 while 1 ?if x>xlimit,break; end ?endwhile 1-end將可以無限循環(huán)下去，當(dāng)條件 x>xlimit得到滿足時(shí)，通過 break語句中止循環(huán)。4.

12、輸入、輸出語句輸入(input)語句舉例如下：z=input('input z=');輸入一個(gè)數(shù)zp=input('Your name=','s');輸入一個(gè)字符或字符串輸出(fprintf)與C語言的輸出語句類似，如fprintf('Volume=%12.5fn',vol);除120 格式，還可以輸出 12.5e 格式、%d 格式、%s格式等。2-3.常用函數(shù)1 .常用教學(xué)函數(shù) 三角 函數(shù)： sin(x)、cos(x)、tan(x)、asin(x)、acos(x)、atan(x) 初等函數(shù)：abs(x)、sqrt(x)、roun

13、d(x)(四舍五入取整)、fix(x)(去尾數(shù)取整)、mod(x,y)(取余數(shù))、 exp(x)(以e為底的指數(shù))、log(x)(以e為底的對數(shù))、log10(x)(以10為底的對數(shù))2 .常用功能函數(shù)sum(x)求向量/數(shù)組元素之和 max(x)求向量/數(shù)組的最大值 min(x)求向量/數(shù)組的最小值rand(n)生成隨機(jī)數(shù)，當(dāng)n=1時(shí)返回一個(gè)隨機(jī)數(shù)；n>1時(shí)，返回nXn隨機(jī)矩陣 feval(f_name,x)計(jì)算以x為參數(shù)，名為f_name的函數(shù)，例如,s_name='sin'貝Uy=feval(s_name,x)等價(jià)于 y=sin(x)eval(s)執(zhí)行字符串s所代表

14、的命令例如s='x=cos(pi/3)'則eval(s)等價(jià)于執(zhí)行x=cos(pi/3)2-4.幾個(gè)常用的命令1. cd 顯示或改變當(dāng)前目錄cd 顯示當(dāng)前目錄 cd c:matlabrllwork將此目錄設(shè)置為當(dāng)前目錄2. dir 列出當(dāng)前目錄或列出指定目錄中的文件dir顯示當(dāng)前目錄中的內(nèi)容 dir c:matlabrllbin顯示此目錄中的內(nèi)容3. disp顯示變量內(nèi)容4. type列出指定文件的全部內(nèi)容5. clear清除內(nèi)容中的變量和函數(shù)6. home 清屏并將光標(biāo)移至左上角7. exit 或 quit 退出 Matlab做矩陣的有關(guān)計(jì)算3-1.矩陣的輸入1 .在命令窗口

15、內(nèi)直接輸入?1 2 3?對于a=?4 5 6?,可采用如下輸入歹8 9?a=1,2,3;4,5,6;7,8,9或 a=1 2 3;4 5 6;7 8 9或 a=1 2 34 5 67 8 92 .創(chuàng)建M文件對于大型矩陣直接輸入不方便，可以采用創(chuàng)建 M文件的方式，即在編輯窗口中輸入該矩陣(輸入 格式同3-1的1),然后保存在當(dāng)前目錄，每次進(jìn)入 Matlab時(shí)該矩陣已自動(dòng)調(diào)入內(nèi)容。3 .生成特殊矩陣的命令 zeros生成0陣?yán)纾?a=zeros(3)生成 3X3 的 0 陣 a=zeros(3,2)生成 3 >2 的 0 陣 b=zeros(size(a)生成與矩陣a維數(shù)相同的0陣格式同z

16、eros格式同zeros ones生成1陣eye生成單位陣4 .向量的輸入行向量x=1,2,3,4,5列向量x=1,2,3,4,5'3-2.矩陣/向量的運(yùn)算1 .加減c=a+bc=a-bc=a+3（矩陣亦可與數(shù)量運(yùn)算）2 .乘c=a*bc=3*a3 .求逆c=inv（a） （a須為非奇異方陣）4 .除 矩陣的左除c=ab（等效于a-1*b） 矩陣的右除c=a/b（等效于b*a-1）5 .行列式的值b=det（a） （a須為方陣）6 .轉(zhuǎn)置c=a'對于矢邱車A=?3O?6即3-3.矩陣的范數(shù)111.無窮范數(shù)即矩陣元素的絕對值按行相加的最大值7H 義： A = max1110010

17、可j=1輸入命令p=norm(A,inf)顯示p=192.1- 范數(shù)定義:IAn=max Ta11勺叼臺(tái)F即矩陣元素絕對值按列相加的最大值輸入命令p=norm(A,1)顯示p=223.2- 范數(shù)定義：IA| 2二(ATA),入max 6A)表示為ATA的最大特征值輸入命令p=norm(A,2)或 p=norm(A)顯示p=15.86873-4.向量的范數(shù)對于向量V=1 3 6 -71 .無窮范數(shù)定義：VII= max v1輸入命令 p=norm(V,inf)即 max(abs(V)顯示p=7輸入命令p=norm(V,-inf)即min(abs(V)顯示p=11.1- 范數(shù)n定義：V 1 = 1

18、> i=1輸入命令p=norm(V,1)顯示p=173.2- 范數(shù)定義：V2輸入命令p=norm(V,2)或 p=norm(V)顯示p=9.74683-5.矩陣的條件數(shù)定義：A為非奇異矩陣，cond(A)r = A1 口 A r為矩陣A的條件數(shù)。其中r= 8,1,2 ；分別稱為無窮 范數(shù)條彳數(shù)、1-范數(shù)條件數(shù)、2-范數(shù)條件數(shù)。對于矢邱車A =? ? ? 09. 4-171 .無窮范數(shù)條件數(shù)輸入命令p=cond(A,inf)顯示p=29.52632.1-范數(shù)條件數(shù)輸入命令p=cond(A,1)顯示p=30.31583.2-范數(shù)條件數(shù)輸入命令p=cond(A) 或 p=cond(A,2)顯示

19、p=19.05753-6.矩陣的特征值和特征向量定義：設(shè)A是nxn的矩陣，滿足關(guān)系式AXi = AX，則(i=1,2,n)為A的特征值，向量 X為人對應(yīng)的特征向量。?3 2?對于矢I陣A= ?1 2?1 .特征值輸入命令d=eig(A)顯示d= 41(即 1=4、 2=1)2 .特征向量和特征值輸入命令v,d=eig(A)顯示v= 0.8944-0.70710.44720.7071d= 4001?0.8944?以上結(jié)果即為1=4、其對應(yīng)的特征向量為 ?；1W4472? ,?- 0.7071?2=1、其對應(yīng)的特征向量為 ?。? 0.7071 ?§4.Matlab 繪圖4-1.繪圖的基本

20、命令l.plot畫出點(diǎn)數(shù)據(jù)集合的曲線1) 基本格式plot(x,y) x,y是點(diǎn)的坐標(biāo)數(shù)組例如 x=0:0.2:10y=x.A3-12*x.A2-7*x+12Plot(x,y)2) 擴(kuò)展格式 plot(x,y,'linewidth',3)設(shè)置線的粗細(xì) plot(x,y,'r')設(shè)置線的顏色3時(shí)，線變細(xì)，反之亦然，此數(shù)的默認(rèn)值是0.5。'linewidth'和3控制所繪線的粗細(xì)，當(dāng)用較小的數(shù)代替r表示畫出紅色的線，其余顏色設(shè)置見下表:紅黃紫青綠蘭白里 八、rymcgbwk plot(x,y,'+')設(shè)置繪出線的標(biāo)記上述命令表示用標(biāo)記

21、+繪出線，設(shè)置標(biāo)記見下表：星形占八、十字形叉形菱形圓形五角形正方形三角形*.+XDOPSA顏色與標(biāo)記可以組合使用，例如plot(x,y,' g')2 .grid 畫網(wǎng)格格式 grid on繪圖形加上網(wǎng)格grid off去掉網(wǎng)格3 .hold保持圖形格式 hold on保持圖形hold off取消此功能使用plot繪出一條曲線后，再繪出另一條曲線之前，需用此命令將原曲線保持下來。4 .xlabel和ylabel給坐標(biāo)軸加標(biāo)注給x和y坐標(biāo)軸加上標(biāo)注例如 xlabel('x')給x軸加上標(biāo)注xylabel('y=sin(x)')給 y 軸加上標(biāo)注 y=

22、sin(x)5 .title 給圖形加上標(biāo)題例如 title('pressure Rario')6 .clf和cla清除圖形窗口中的全部內(nèi)容7 .subplot繪制圖形陣列調(diào)用格式為 subplot(m,n,k)在圖形窗口內(nèi)繪制 mxn個(gè)圖形陣列，k是圖形窗口的序號。見下述例2。8.示例例1 輸入下面的M文件并運(yùn)行之clf; x=0:0.2:10;10.6y1=sin(x);y2=exp(-0.4*x);y3=y1.*y2;plot(x,y1,'r');hold onplot(x,y2,'g');plot(x,y3,'linewidth&

23、#39;,3);grid onxlabel('X');ylabel('Y1=sin(X) Y2=exp(-0.4*X)Y3=Y1*Y2');例2輸入下面的M文件并運(yùn)行之-oa6 4 2 D 2 4 6 0.OI.Q 。-O.-OI1*>=£>-(各 FUKMm UkiV 二clf; t=0:0.1:30;subplot(2,2,1);plot(t,sin(t);title('SubPlot No:1');SubPlot Nd:1xlabel('t');ylabel('Y=sin(t)');su

24、bplot(2,2,2);plot(t,t.*sin(t);title('SubPlot No:2');xlabel('t');ylabel('Y=t*sin(t)');subplot(2,2,3);plot(t,t.*sin(t)A2);title('SubPlot No:3');xlabel('t');ylabel('Y=t*sin(t)A2');subplot(2,2,4);plot(t,t.A2.*sin(t).A2);title('SubPlot No:4');xlabel

25、('t');ylabel('Y=tA2*sin(t)A2');垂直疊放的兩個(gè)圖形可通過如下方式繪制:subplot(2,1,1);plot(;subplot(2,1,2);plot(；20口-i J410102030tSuhPhrt Nd:330r-4C 01000Figure 14-2.圖形的交互編輯圖4-2-1所示為圖形窗口左上角局部。1.改變曲線的屬性用鼠標(biāo)點(diǎn)擊菜單欄中按鈕移動(dòng)鼠標(biāo)使其指向所要編輯曲線上的任意一點(diǎn)。 按下鼠 標(biāo)右鍵，曲線變成連續(xù)的黑點(diǎn)，并彈出下拉菜File Edit View (ns：art Tools Desktop口方0昌|恃I默RO&

26、#174;SubPlot No.1Window He 制400.5-0.6) |單。點(diǎn)擊菜單中各盤命令,可以改變曲線的 顏色、線寬、標(biāo)記和曲線類型等屬性。圖.'42T.i30o o o 2 4 a三產(chǎn)192.其他的編輯命令通過編輯窗口內(nèi)的Tools菜單項(xiàng)亦可完成若干編輯功能。讀者可以自行摸索。ii第二部分?jǐn)?shù)值計(jì)算§1.方程求根1-1.牛頓迭代法x- 0.01 一例 1求解萬程(0.01x+1)sin x-2+1- - 0.0096 = 0在 xo=4 附近的一個(gè)根。1 .編寫牛頓迭代法的函數(shù)M文件function x=Newt_n(f_name,x0)x=x0;xb=x+1;

27、k=0;n=50;del_x=0.01;while abs(x-xb)>0.000001k=k+1;xb=x;if k>=n break;end; y=feval(f_name,x);y_driv=(feval(f_name,x+del_x)-y)/del_x;x=xb-y/y_driv;fprintf('k=%3.0f,x=%12.5e,y=%12.5e,yd=%12.5en',k,x,y,y_driv);end;fprintf('n Final answer=%12.6en',x);用默認(rèn)文件名 Newt_n.m存盤。2 .編寫待求方程的函數(shù)M文

28、件function y=eqn_1(x)y=(0.01*x+1)*sin(x)-(x-0.01)*(xA2+1)A(-1)-0.0096;用默認(rèn)文件名eqn_1.m存盤。3 .求解步驟在命令窗口中鍵入 Newt_n('eqn_1',4)顯示： k= 1,x=2.36795e+000,y=-1.03138e+000,yd=-6.31954e-001k= 2,x=2.92631e+000,y=3.48815e-001,yd=-6.24716e-001k= 6,x=2.82170e+000,y=-9.55313e-009,yd=-8.88819e-001Final answer=2.

29、821702e+000 注意：Newt_n.m和eqn_1.m應(yīng)保存在當(dāng)前目錄。1-2.圖解法確定迭代的初始點(diǎn)x- 0.01萬程(0.0僅+1)sin x- 0.0096 = 0具有多重根，若要求出某一區(qū)間(如xC|0,10|)的所x +1有根，可用圖解法確定其初始點(diǎn)，然后再調(diào)用Newt_n函數(shù)解之，具體步驟為：1 .繪方程曲線該方程的解也就是? y產(chǎn)(0.0僅+1) sin xx- 0.010.0096的交點(diǎn)，故應(yīng)繪出曲線y1和y2,觀測兩條曲線交點(diǎn)的坐16標(biāo)值。輸入M文件如下:clf;hold on; x=0:0.01:10;y1=(0.01*x+1).*sin(x);y2=(x-0.01

30、)./(x.A2+1)-0.0096;plot(x,y1,'r');plot(x,y2,'g');觀測圖形，兩條曲線在 xe0,10有3個(gè)交點(diǎn)，交點(diǎn)x坐標(biāo)值約為3、6.5和9.5,可將這些點(diǎn)作為牛 頓迭代法的初始點(diǎn)。2 .求解鍵入：Newt_n('eqn_1',3)顯示：ans=2.8217鍵入：Newt_n('eqn_1',6.5)顯示：ans=6.4351 鍵入：Newt_n('eqn_1',9.5)顯示：ans=9.31892線性方程組2-1.迭代法的收斂性?10x -為-2x = 7.2 ?例 2對于? -

31、x+-10x+; 523x=34=28.3建立 Jocobi 迭代格式為：1)?x1k+1 ? ?0 ? k+1 ?為?=夕.19 k+1解法一:.2k -0.1 0.2?x1 ?0.72?00.2?>2 k?+ ?0.83?,判斷其收斂性。0.20 ?xk ? ?0.42?0 0.1 0.2? ? 輸入迭代矩陣 B= ?0.100.2?如 0.2 0 ?調(diào)用求矩陣特征值命令(eig命令)輸入：eig(B)顯示：ans= - 0.10000.3372-0.2372矩陣B的譜半徑均為該矩陣的最大特征值，故B)= 0.3372<1,該迭代格式收斂。2) 解法二：調(diào)用求矩陣范數(shù)的命令(n

32、orm命令)若矩陣B的任意一種算子的范數(shù)小于1,則迭代收斂。故可計(jì)算B |、|闿1或B|2中的任意一種，只要小于1,即可判斷出迭代收斂。輸入：norm(B,1)顯示：ans=0.4000即|B | =0.4000<1 ,該迭代格式收斂。亦可調(diào)用norm(B,2)或norm(B,inf)來判斷。例3對于例2的線性方程組，若建立迭代格式為?x1k+1? ?010- 2?x1 ? ?- 8.3? 9o? ,? 9 o?2k+1 ?= ?- 105 ?為 k?+ ? 4.2?,其收斂性如何？?x：+1 ? ?5 -0.50 ?x3? ?-3.6?2-2.線性方程組的病態(tài)問題?11 ?x ? ?2

33、?例4分析方程組? ?= ? ?的病態(tài)性。.夕 1.000?%? ?2?輸入該方程的系數(shù)矩陣A=?11 ?夕 1.0001?cond命令) 調(diào)用求矩陣條件數(shù)命令(輸入：cond(A,1)顯示： ans=4.0004e+004即cond(A)= A-1 |1? A |1 =4.0004 X104>>1 ,該方程組是病態(tài)方程組。亦可調(diào)用 cond(A)或cond(A,inf)命令來判斷。2-3.求解線性方程組1 .調(diào)用x=Ab命令求解方程組對于線性方程組 AX=b （其中A為n3的矩陣、x和b均為n維列向量），有X=A -1b,該式與Matlab 的矩陣左除運(yùn)算等效，即 X=Ab 。?

34、3x + 2為=-1 例5求解? X - & =1?3輸入系數(shù)矩陣A= ?12 ?- 1?和吊數(shù)歹U向重b= ? ?1 ?-1?輸入：X=Ab 顯示：X=0.2000-0.80002 .利用LU分解求解方程組1) LU分解? 2x + x - 3%=2一一 ?對于-x + 3x + 2x = 0?1 -3?2? 3x +&-敘=1?2? ? 輸入系數(shù)矩陣A=? 1 32 ?和常數(shù)項(xiàng)列向量b= ?0? ?3 1 -3?1?調(diào)用LU分解命令（lu命令）輸入：l,u,p=lu(A)顯示：1=1.000000-0.33331.000000.66670.10001.0000u=3.000

35、01.0000-3.000003.33331.000000-1.1000?123I2p=0 01其中l(wèi)為單位下三角陣100（對角線元素均為1）, u為上三角；p為置換矩陣，它記錄了選主元高斯消去法中行交換的信息。它們之間的關(guān)系為p*A=l*u2）利用LU分解求解線性方程組 對于AX=b ,用p左乘方程兩邊，得pAX=pb21uX=yly=pb通過對原方程組系數(shù)矩陣 操作步驟為：在命令窗口內(nèi)輸入插值節(jié)點(diǎn)坐標(biāo)x=2,4,6,8,10y=5,7,9,11,13xi=2.5,5,7.3,9.1x、y以及待插值的數(shù)據(jù)xi如下:(a)(b)A的LU分解，將求解AX=b的問題轉(zhuǎn)化為求解方程組(a)和方程組(

36、b)。其D求方程組(b)中的未知變量y輸入：y=l(p*b) 顯示：y=1.00000.3333 1.3000D求方程組(a)中的未知變量X輸入：X=uy顯示：X=-1.00000.4545-1.1818亦可將上述 、的命令合并為 X=出(l(p*b)§3.插值和擬合3-1.Lagrange 插值1 .編寫Lagrange函數(shù)M文件function fi=lagrange(x,y,xi)fi=zeros(size(xi);npl=length(y);for i=1:nplz=ones(size(xi);for j=1:nplif i=j,z=z.*(xi-x(j)/(x(i)-x(j

37、);end;end;fi=fi+z*y(i);end注意： 插值節(jié)點(diǎn)坐標(biāo)x、y為數(shù)組，待插值的xi可以是一個(gè)標(biāo)量，亦可以是數(shù)組。Lagrange函數(shù)M文件必須保存在當(dāng)前目錄內(nèi)。2 .調(diào)用Lagrange插值函數(shù)例6已知插值節(jié)點(diǎn)xk246810yk .、._ _ 5 _一 7,一、9 1113調(diào)用Lagrange函數(shù)輸入：yi=lagrange(x,y,xi)顯示：yi=5.5000 8.000010.3000 12.10003-2.代數(shù)多項(xiàng)式插值1 .代數(shù)多項(xiàng)式的表示方法已知(n+1)個(gè)節(jié)點(diǎn)信息 保/N=0,1,n可以構(gòu)造一個(gè)n次代數(shù)插值多項(xiàng)式n ,n-1Pn = a X + an- 1X+?

38、+&x+a)如，watiabx5上峪用桃示的數(shù)舞麗,4喬素為多項(xiàng)式的系數(shù)，并且從左至右按降哥排列。例2 .代數(shù)多項(xiàng)式的插值計(jì)算分成二步進(jìn)行：通過polyfit 命令，又已知(n+1)個(gè)節(jié)點(diǎn)唯一確定一個(gè)n次代數(shù)插值多項(xiàng)式；利用這個(gè)代數(shù)插值多項(xiàng)式計(jì)算待插點(diǎn)處的函數(shù)值。1)構(gòu)造n次代數(shù)插值多項(xiàng)式 調(diào)用格式polyfit(x,y,n)其中數(shù)組x,y是已知節(jié)點(diǎn)坐標(biāo)，n是代數(shù)插值多項(xiàng)式的階數(shù)。利用該命令可求得n次代數(shù)多項(xiàng)式p(x) = c x+cn-1x n-1+?+cx+G 的系數(shù)。例如，輸入節(jié)點(diǎn)數(shù)據(jù)：x=1.1,2.3,3.9,5.1y=3.887,4.276,4.651,2.117 調(diào)用命令

39、：a=polyfit(x,y,length(x)-1)顯示：a=-0.20151.4385-2.74775.4370其中l(wèi)ength(x)-1=4-1=3,故所得的多項(xiàng)式是：y=-0.2015x3+1.438或-2.7477x+5.43702)代數(shù)多項(xiàng)式的插值計(jì)算調(diào)用格式 polyval(a,xi)其中a為插值多項(xiàng)式的系數(shù)數(shù)值，xi為待插點(diǎn)數(shù)組。若要計(jì)算xi=1.8,1.95,2.93,4.8 時(shí)的函數(shù)值，調(diào) 用該命令即可。輸入：xi=1.8,1.95,2.93,4.8yi=polyval(a,xi)顯示：yi=3.97714.05524.66853.1127若調(diào)用Lagrange插值，可得到

40、與上面相同的結(jié)果，說明了插值多項(xiàng)式存在的唯一性。3 - 3.插值誤差以下用作圖方法觀察插值的誤差。x例7已知曲線f(>) = e.5- 2sin x,采用3個(gè)節(jié)點(diǎn)x1.12.35.1yf(i.i)f(2.3)f(5.1)構(gòu)造插值多項(xiàng)式P2(x)作為f(x)的近似，通過作圖的方法觀察在x41.1,5.1的誤差。運(yùn)行以下腳本M文件clear,clf;hold on;x=1.1,2.3,5.1;y=exp(x/1.5)-2*sin(x);plot(x,y,'o');n=length(x)-1;coeff=polyfit(x,y,n);xp=1.1:0.05:5.1;yp=pol

41、yval(coeff,xp);plot(xp,yp,'r');yh=exp(xp/1.5)-2*sin(xp);plot(xp,yh,'k');xlabel('X');ylabel('f(x) P(x),data points:。');x圖中黑線為曲線f(R = e，5 - 2sin x,紅線為由3個(gè)節(jié)點(diǎn)成的2次插值多項(xiàng)式，圓圈表示插值節(jié)點(diǎn)。3-4.分段線性插值調(diào)用格式 yi=interpl(x,y,xi,'linear')其中數(shù)組x和y為已知節(jié)點(diǎn)坐標(biāo)，xi是待插值的標(biāo)量或數(shù)組，對應(yīng)的 yi值通過線性插值算出。注意

42、，調(diào)用格式中的數(shù)組全部用列向量表示。例8已知數(shù)據(jù)xk135.281013yk46389.1-4求xi=4,5.5；11.7時(shí)線性插值函數(shù)之值。在命令窗口內(nèi)輸入列數(shù)組x、y和xix=1,3,5.2,8,10,13' y=4,6,3,8,9.1,-4' xi=4,5.5,11.7'調(diào)用分段線性插值命令yi=interp1(x,y,xi,'linear') 顯示：yi=4.63643.53571.67673-5.數(shù)據(jù)的曲線擬合例9已知一組數(shù)據(jù)，作擬合曲線i012345xi0.10.40.50.70.70.9yi0.610.920.991.521.472.031

43、）方法一：求解正規(guī)方程 將上述一組數(shù)據(jù)標(biāo)在坐標(biāo)紙上，觀察到各點(diǎn)在一條直線 附近，可設(shè)擬合曲線為y=co+cix其正規(guī)方程為：?G，?0 ) (%?)夜?_ ?0, f)?1?)?,?)? ? ?(? , f)?，其中?0 (X)=1,?1 (x)= x55(?0, ?o)=刀(x)?0(x)= B?1 = 6 i=0=055(?0, ?i)= (?1, ?0)= X? (x)?i(x)=匯1?x = 3.3 ""i=0i=055?, ?)= ?(x )?(x)=x2 = 2.21112 1 i 1 i 2 i% f)= Z?：(x)?y=歸&=7.54 =0i=05

44、5(?1, f)= X? (x)?y= »? = 4.844 i=0=0在Matlab命令窗口中輸入： a=6,3.3;3.3,2.21b=7.54,4.844'c=ab顯示： c=0.28621.7646即所要求的擬合曲線為y=0.2862+1.7646x2)方法二：調(diào)用 polyfit(x,y,n)調(diào)用polyfit(x,y,n),可求得擬合代數(shù)多項(xiàng)式4為=c義+cn.1x *+?+Qx+G的系數(shù)，其中x,y為節(jié)點(diǎn)坐標(biāo)數(shù)組，n為擬合多項(xiàng)式的階數(shù)。對于例9中的節(jié)點(diǎn)數(shù)據(jù)，求擬合直線方程，操作如下：輸入： x=0.1,0.4,0.5,0.7,0.7,0.9,y=0.61,0.9

45、2,0.99,1.52,1.47,2.03c=polyfit(x,y,1)顯示： c=1.76460.2862即所要求的擬合曲線為y=1.7646x+0.28624數(shù)值積分4-1.復(fù)合梯形求積公式復(fù)合梯形求積公式為24I = f(x)dx-h?f(a+2 n-1 f(x ) + f(b)? Z2 k1)2)3)其中h=b-?k=1=a+k k=1,2,?n-1。n輸入復(fù)合梯形公式的函數(shù)function t=trapez限件ia(a,b,n,f_name)h=(b-a)/n;fa=feval(f_name,a);fb=feval(f_name,b);t1=0;for k=1:n-1 x(k)=a

46、+k*h;t1=t1+feval(f_name,x(k);end;t=h*(fa+2*t1+fb)輸入待積分的函數(shù)M文件例10中的被積函數(shù)為可勾=4y，其M文件如下，并以funjts.m文件名存盤。1 + x2function y=func_ts(x)y=4/(1+xA2);調(diào)用復(fù)合梯形公式若取 n=8,則輸入：t=trapezia(0,1,8,'func_ts')顯示：t=3.13904-2.復(fù)合Simpson求積公式復(fù)合Simpson求積公式為b h?1I = /(Rdx ?2nf(a + f(切+。2 f%1) + f %)?k=1其中 h= b- a,x n=a+ kh

47、rk= 1,2?,2n。21)輸入復(fù)合Simpson公式的函數(shù) M文件 function s=simpson(a,b,n,f_name) h=(b-a)/n;for k=1:2*nx(k)=a+k*h/2;end;fa=feval(f_name,a); fb=feval(f_name,b); s1=0;s2=0;for k=1:ns1=s1+feval(f_name,x(2*k-1);s2=s2+feval(f_name,x(2*k);end;s=h*(0.5*(fa-fb)+2*s1+s2)/3;2) 輸入待積分的函數(shù)M文件此處仍以例10中的被積函數(shù)為例，且以 fun_ts.m文件名存盤。3

48、) 調(diào)用復(fù)合Simpson公式若取 n=4,則輸入：s=simpson(0,1,4,'func_ts') 顯示：s=3.1416通過對復(fù)合梯形公式和復(fù)合Simpson公式的比較，對于同一個(gè)被積函數(shù)，前者劃分為8等分、計(jì)算結(jié)果為3.1390;后者劃分為 4等分、計(jì)算結(jié)果為3.141。兩種算法計(jì)算被積函數(shù)值的次數(shù)相近(復(fù)合梯形公式計(jì)算被積函數(shù)值9次、復(fù)合Simpson公式10次)，但復(fù)合Simpson公式的精度(代數(shù)精度3)高于復(fù)合梯形公式(代數(shù)精度1)。§5,常微分方程的數(shù)值解法5-1.Euler 方法1)計(jì)算公式,、r一？y = f(x y1,對于初值問題?, Eul

49、er公式y(tǒng)+1 = y + h?f(x, y)。以下通過仞題，說明Euler方法的?y(x) = %使用。例11質(zhì)量m=70kg的人在t=0時(shí)刻從飛機(jī)上跳出，假設(shè)跳傘者垂直下降，且初始時(shí)刻垂直速度速度v(t1)=0,空氣阻力F=cv2 (N), c=0.27kg/m,取步長h=0.1。試求速度，并繪出 00t 020s的解的圖 形。2)分析由牛頓第二定律，mFv= - F + ng,式中重力加速度g=9.8m/s2。那么，建立初值問題為： dt22?V= f (t, V = - + g(5-1)?m? v(t) = 0根據(jù)例題要求，tea,t=0,20,將求解區(qū)間離散如下：11切131用I|1

50、I0=Q1=20其中等分?jǐn)?shù)n= -b-ah對于所求的(5-1)式，采用Euler公式，可得下述的計(jì)算格式:? = M+=(0 1)hi = 2,3,?,n+1i=2,3?,n+1? v = 0?V= v-i + h* f(tiNi)3) 編寫f(t,v)函數(shù)程序function f=descent(t,v) c=0.27,m=70,g=9.8;f=-c*vA2/m+g;4) 編寫該問題的Euler法求解程序clear,clc;a=0,b=20,h=0.1;n=(b-a)/h;t(1)=0；v(1)=0;for i=2:n+1t(i)=t(1)+(i-1)*h;v(i)=v(i-1)+h*des

51、cent(t(i-1),v(i-1);end;for i=1:n+1fprintf('Time=%f,t(i);fprintf(' Velocity=%fn',v(i);end; plot(t,v);xlabel('Time (s)');ylabel('Velocity (m/s)');5) 運(yùn)行求解例11的結(jié)果如下：Time=0.000000Time=2.000000Time=4.000000Time=6.000000Time=8.000000Time=10.000000Time=12.000000Time=14.000000Time=16.000000Time=18.000000Time=20.000000Velocity=0.000000Velocity=18.730327Velocity=32.989678Velocity=41.671830Velocity=46.250031Velocity=48.480593Velocity=49.525261Velocity=50.005441Velocity=50.224250Velocity=50.323563Velocity=50.3685585-2.改進(jìn)的Euler