matlab數(shù)值運(yùn)算課件

1、第四章第四章 數(shù)值計(jì)算功能數(shù)值計(jì)算功能4.1 多項(xiàng)式計(jì)算多項(xiàng)式計(jì)算4.2 數(shù)值插值和曲線擬合數(shù)值插值和曲線擬合4.3 函數(shù)極值函數(shù)極值4.4 非線性方程問題求解非線性方程問題求解4.5 常微分方程的數(shù)值求解常微分方程的數(shù)值求解4.6 數(shù)值微分與積分?jǐn)?shù)值微分與積分4.7 概率統(tǒng)計(jì)概率統(tǒng)計(jì)4.1 多項(xiàng)式運(yùn)算多項(xiàng)式運(yùn)算 1多項(xiàng)式表示法多項(xiàng)式表示法 2多項(xiàng)式運(yùn)算函數(shù)多項(xiàng)式運(yùn)算函數(shù) 1多項(xiàng)式表示法多項(xiàng)式表示法 (1) 直接法直接法:多項(xiàng)式多項(xiàng)式p 多項(xiàng)式系數(shù)用行向量行向量表示，按降冪排列降冪排列；p 缺某冪次項(xiàng)，其冪次項(xiàng)系數(shù)為零冪次項(xiàng)系數(shù)為零P=a0,a1,.an-1,an符串形式符串形式,x多項(xiàng)式變量

2、y=poly2str(P, x)(2) 字符串表示法：字符串表示法：y= poly2sym(P, x) (3) 符號(hào)多項(xiàng)式表示法：符號(hào)多項(xiàng)式表示法：完整形式完整形式p=1,3,0,4,5y=poly2str(p,x)y = x4 + 3 x3 + 4 x + 5syms xy=poly2sym(p, x)y=x4 + 3*x3 + 4*x + 5 p=1,3,0,4,5p = 1 3 0 4 5 y=poly2str(p,x)y = x4 + 3 x3 + 4 x + 5 rp=roots(p)rp = -3.2346 0.5594 + 1.1980i 0.5594 - 1.1980i -0.

3、8843例例2.多項(xiàng)式的運(yùn)算函數(shù)多項(xiàng)式的運(yùn)算函數(shù)多項(xiàng)式的值；根和微分；擬合曲線；部分分式多項(xiàng)式的值；根和微分；擬合曲線；部分分式polyvalpolyvalpolyvalmcovolution n次多項(xiàng)式次多項(xiàng)式n個(gè)根，或個(gè)根，或?qū)嵏鶎?shí)根，或若干對，或若干對共軛復(fù)根共軛復(fù)根。p P：多項(xiàng)式系數(shù)向量；：多項(xiàng)式系數(shù)向量；p r：n個(gè)根個(gè)根的的。 r=roots(P)r = -7.6998 -0.5572 + 1.1335i -0.5572 - 1.1335i 0.8141 例例求多項(xiàng)式求多項(xiàng)式x4+8x3+3x2+4x-10的根的根P=1,8,3,4,-10;r=roots(P)多項(xiàng)式根的r向量：

4、 創(chuàng)建多項(xiàng)式創(chuàng)建多項(xiàng)式P=poly(r)PP = 1.0000 8.0000 3.0000 4.0000 -10.0000避免浮點(diǎn)顯示 p=poly(r)p = 1 8 3 4 -10多項(xiàng)式的向量系數(shù)：，則求多項(xiàng)式在該點(diǎn)，則求多項(xiàng)式在該點(diǎn) 的值；的值；每元素每元素x(i i,j)多項(xiàng)式求值。多項(xiàng)式求值。Y=a0*X.n+a1*X.(n-1)an+1 p=1,3,0,2; poly2str(p,x)ans = x3 + 3 x2 + 2 a=1,2;3,4a =1 2 3 4 y=polyval(p,a)y =6 22 56 114 y=polyval(p,2)y = 22例例例例 Y = po

5、lyvalm(P,X)矩陣多項(xiàng)式求值矩陣多項(xiàng)式求值Y=a0*Xn+a1*Xn-1an+1p X必為方陣，作必為方陣，作自變量自變量代入多項(xiàng)式求值；代入多項(xiàng)式求值；p 結(jié)果階數(shù)還是保持方陣階數(shù)。結(jié)果階數(shù)還是保持方陣階數(shù)。X設(shè)A為方陣，P代表多項(xiàng)式x3-5x2-2：pp=polyval(P,A)的含義pp=A.*A.*A -5*A. *A-2Pm=polyvalm(P,A)的含義：Pm=A*A*A-5*A*A-2 p=1,-5,0,-2; a=2,4;1,0a = 2 4 1 0 pp=polyval(p,a)pp = -14 -18 -6 -2 pm=polyvalm(p,a)pm =-18 -

6、8 -2 -14 例例例例14p 相同次數(shù)多項(xiàng)式，加減其系數(shù)向量相同次數(shù)多項(xiàng)式，加減其系數(shù)向量, , p 不同次數(shù)不同次數(shù), ,低次多項(xiàng)式中低次多項(xiàng)式中。 p2 = x3 - 0 x2+2x + 1 p1 + p2 = 2x3 - x2 + 2x + 4 2, -1, 0, 3 2, 1 0, 0, 2, 1 2, -1, 2, 4 p1 = 2x3 - x2 + 3例例例例多項(xiàng)式x4+8x3-10與2x2-x+3的乘積。 p1=1,8,0,0,-10; p2=2,-1,3; p12=conv(p1,p2)p12 = Columns 1 through 5 2 15 -5 24 -20 Col

7、umns 6 through 7 10 -30 例例例例 h = 3,2,1,-2,1,0,-4,0,3; x = 1,-2,3,-4,3,2,1; y = conv(h,x); n = 0:14; stem(n,y); %桿圖 xlabel(Time index n); %標(biāo)坐標(biāo)軸 ylabel(Amplitude); title(Output Obtained by Convolution)%圖形標(biāo)題grid;卷積是兩個(gè)變量在某范圍內(nèi)相乘后求和的結(jié)果。卷積的變量是序列x(n)和h(n)，y=conv(x,h) 來實(shí)現(xiàn)卷級(jí)的，輸出的來實(shí)現(xiàn)卷級(jí)的，輸出的結(jié)果個(gè)數(shù)等于結(jié)果個(gè)數(shù)等于x x的長度與的

8、長度與h h的長度之和減去的長度之和減去1 1。卷積公式：卷積公式：z(n)=x(n)z(n)=x(n)* *y(n)= x(m)y(n-m)dm.y(n)= x(m)y(n-m)dm.例例p Q：多項(xiàng)式P1除以P2的商式，p r：P1除以P2的余式，p Q和r仍是多項(xiàng)式系數(shù)向量。積這兩類問題都稱作解卷探、石油勘探等問題。地震信號(hào)處理、地質(zhì)勘求多項(xiàng)式求多項(xiàng)式x4+8x3-10除以除以2x2-x+3的結(jié)果。的結(jié)果。 q,r=deconv(p12,p2)q = 1 8 0 0 -10 r = Columns 1 through 5 0 0 0 0 0 Columns 6 through 7 0 0

9、 q,r=deconv(p2,p12)q = 0 r = 2 -1 3 例例 多項(xiàng)式乘積PQ的導(dǎo)函數(shù)：多項(xiàng)式P的導(dǎo)函數(shù)：導(dǎo)函數(shù)的分子存入p，分母存入q。多項(xiàng)式除法P/Q的導(dǎo)函數(shù)：例：求有理分式的導(dǎo)數(shù)。命令如下：P=1;Q=1,0,5;p,q=polyder(P,Q)例例例例21I=polyint(2,-1,0,3,5); 例：例： p(x) = 2x3 - x2 + 3求求 常數(shù)項(xiàng)常數(shù)項(xiàng) 5( ) dp xx 不定積分，常數(shù)項(xiàng)取不定積分，常數(shù)項(xiàng)取 c不定積分，常數(shù)項(xiàng)不定積分，常數(shù)項(xiàng)取取 零零例例曲線擬合和數(shù)據(jù)插值曲線擬合和數(shù)據(jù)插值4.2.1 曲線擬合曲線擬合 4.2.2 數(shù)據(jù)插值數(shù)據(jù)插值 p

10、最小二乘法(誤差平方和最小) ；p x、y已知數(shù)據(jù)的橫、縱坐標(biāo)；p n為多項(xiàng)式次數(shù);N次次p S結(jié)構(gòu)體數(shù)組(struct)，估計(jì)預(yù)測誤差，含R，df和normr。p R：先輸入x構(gòu)建范德蒙矩陣V，后QR分解，得上三角矩陣。p df：自由度，df=length(y)-(n+1)。df0時(shí)，超定方程組求解，擬合點(diǎn)數(shù)比未知數(shù)(p(1)p(n+1)多。p normr：標(biāo)準(zhǔn)偏差、殘差范數(shù)，normr=norm(y-V*p)， 此處的p為求解之后的數(shù)值。 年份年份1900191019201930194019501960197019801990人口人口75.9991.97105.71123.20131.66

11、150.69179.32203.21226.50249.63 對不同年份人口數(shù)據(jù)分別進(jìn)行1次、2次和9次多項(xiàng)式擬合（polyfit），用poly2str表示多項(xiàng)式完整形式法； 1次、2次和9次多項(xiàng)式估計(jì)2000年人口(polyval)，結(jié)合美國2000年人口普查截至2000年4月1日美國人口2.81421906億數(shù)據(jù) 繪制1900 2000年間的時(shí)間-人口數(shù)(polyval)曲線，要求用plot不同線型（LineSpec）繪制原始數(shù)據(jù)點(diǎn)、擬合的1次、2次和9次多項(xiàng)式曲線，標(biāo)注圖例 ，說明是否階數(shù)越到高越好 美國從1900年1990年每隔10年進(jìn)行人口普查的數(shù)據(jù)如下表所示：（百萬）例例線型線型說

12、明說明數(shù)據(jù)點(diǎn)數(shù)據(jù)點(diǎn)標(biāo)記符標(biāo)記符說明說明顏色顏色說明說明LineSpec例如r-.*、-.r*、*-.r等形式是等效的 plot(x,y,-gh) clearn=1900:10:1990; r=7599,9197,10571,12320,13166,15069,17932,20321,22650,24963; nrf1=polyfit(n,r,1) %1次多項(xiàng)式擬合nrfs1=poly2str(nrf1,n) %1次多項(xiàng)式完整字符串表達(dá)式字符串表達(dá)式nrf2=polyfit(n,r,2) %2次多項(xiàng)式擬合nrf9=polyfit(n,r,9) %9次多項(xiàng)式擬合nrf10=polyfit(n,r,

13、10) %10次多項(xiàng)式擬合nrf1_2000=polyval(nrf1,2000) %一次擬合得到的2000年人口數(shù)nrf2_2000=polyval(nrf2,2000) %2次擬合得到的2000年人口數(shù)nrf9_2000=polyval(nrf9,2000) %9次擬合得到的2000年人口數(shù)n20=1900:4:2000; % 1900年到2000年的線性等分?jǐn)?shù)組nrfv1=polyval(nrf1,n20); % 從1900年到2000年間1次擬合人口數(shù)nrfv2=polyval(nrf2,n20); % 從1900年到2000年間2次擬合人口數(shù) nrfv9=polyval(nrf9,n

14、20); % 從1900年到2000年間9次擬合人口數(shù) plot(n,r,or, n20,nrfv1,- p4=polyfit(x,y,11) y4=polyval(p4,xx);例例hold on;plot(x,y, *); plot(xx,y1,r-);plot(xx,y2,g-);plot(xx,y3,b-)多項(xiàng)式次數(shù)要適當(dāng)，多項(xiàng)式次數(shù)要適當(dāng)，過低誤差大，過高波過低誤差大，過高波動(dòng)明顯動(dòng)明顯已知：10個(gè)實(shí)驗(yàn)數(shù)據(jù)，分別采用2次、5次、9次和11次擬合，選出最佳擬合次數(shù)3. 3. 函數(shù)擬合函數(shù)擬合 clear x=0:5; y=0.2097,0.3523,0.4339,0.5236,0.75

15、90,0.8998; ly=log(y); p=polyfit(x,ly,1);b=p(1);la=p(2);a=exp(la);Xx=linspace(0,5,30);Yy=a*exp(b*Xx);plot(x,y,ro,Xx,Yy)例例*eb xyaln( )ln( )*yab xly 221/()1/yaxbxcdyyaxbxc例例 插值構(gòu)造的函數(shù)必須通過已知數(shù)據(jù)點(diǎn)； 擬合則不要求,只要均方差最小。 都需根據(jù)已知數(shù)據(jù)構(gòu)造函數(shù)； 可使用得到函數(shù)計(jì)算未知點(diǎn)的函數(shù)值。4.2.2 數(shù)據(jù)插值數(shù)據(jù)插值p構(gòu)造函數(shù)近似表達(dá)式的方法。p常用多項(xiàng)式作插值函數(shù),稱多項(xiàng)式插值。p多項(xiàng)式插值基本思想：高次多項(xiàng)式或

16、分段的低次多項(xiàng)式為被插值函數(shù)f(x)的近似解析表達(dá)式。1. 插值法插值法拉格朗日插值法 牛頓插值法 埃爾米特插值法分段插值法樣條插值法等 函數(shù)函數(shù)y=f(x) 一維插值原理：一維插值原理： 調(diào)用格式調(diào)用格式： yi = interp1(x,y,xi,method)p 已知X,Y兩個(gè)等長向量，采樣點(diǎn)和樣本值；p Xi是向量或標(biāo)量，欲插值點(diǎn)；p Yi是一個(gè)與Xi等長的插值結(jié)果。采用差值法估計(jì)美國的人口數(shù)量(1)2000年人口及 1900 1996年每隔4的人口數(shù)(interp1)；(2)將上題原始數(shù)據(jù)點(diǎn)、2階擬合曲線和插值曲線繪制在同一圖，標(biāo)注坐標(biāo)軸為時(shí)間和人口（xlabel、ylabel）、圖形

17、標(biāo)題（title）為插值和擬合和圖例legend年份年份1900191019201930194019501960197019801990人口人口75.9991.97105.71123.20131.66150.69179.32203.21226.50249.63例例 clear n=1900:10:1990; r=7599,9197,10571,12320,13166,15069,17932,20321,22650,24963;n4=1900:4:1996; % 1900 1996年每隔4的線性等分?jǐn)?shù)組nrf2=polyfit(n,r,2) %2次多項(xiàng)式擬合nrfv2=polyval(nrf2,

18、n4)nr4=interp1(n,r,n4,spline) % 插值1900 1996年每隔4的人口數(shù)nr2000=interp1(n,r,2000) % 線性插值2000年人口數(shù)nr2000=interp1(n,r,2000,spline) % 插值2000年人口數(shù)plot(n,r,or,n4,nrfv2,-*k, n4,nr4,-dg) %繪圖legend(原始數(shù)據(jù),2階多項(xiàng)式擬合曲線,插值曲線) %圖例 例例 插值點(diǎn)樣本值落在兩相鄰數(shù)據(jù)點(diǎn)的連線上. (缺省). 兩點(diǎn)間插值點(diǎn)對應(yīng)值就是離兩點(diǎn)最近那點(diǎn)值。 已知數(shù)據(jù)求3次多項(xiàng)式，用多項(xiàng)式求插值。 每分段內(nèi)構(gòu)造一3次多項(xiàng)式，插值函數(shù)除滿足插值條

19、件和節(jié)點(diǎn)處光滑條件，再按照樣條函數(shù)插值。yy = spline(x,Y,xx)clearx=0:1:10;y=sin(x);xx=0:0.2:10;yy=interp1(x,y,xx)subplot(1,4,1)plot(x,y,-*,xx,yy,dr);subplot(1,4,2);y2=interp1(x,y,xx,nearest); plot(x,y,-*,xx,y2,dr);subplot(1,4,3);y3=interp1(x,y,xx,cubic );plot(x,y,-*,xx,y3,dr)subplot(1,4,4);y4=interp1(x,y,xx,spline);plot

20、(x,y,-*,xx,y4,dr)X X1 1取值范圍不能超出取值范圍不能超出X X給給定范圍，否則會(huì)給出定范圍，否則會(huì)給出“NaN”“NaN”錯(cuò)誤。錯(cuò)誤。例例38.02522138.07526738.12531838.17548738.22581538.275150238.325316238.375620738.425841438.475815638.525597538.575347338.625160138.67571038.72538638.77532538.82524738.87523238.92520738.97518239.025157例例3838.238.438.638.8393

21、9.239.40100020003000400050006000700080009000plot(xrd1(:,1),xrd1(:,2) xx=linspace(38,39,80);yy=interp1(xrd1(:,1),xrd1(:,2),xx,spline); plot(xx,yy,-*g) x=xy(:,1)x = 38.0250 38.0750 38.1250 38.1750 38.2250 38.2750 38.3250 38.3750 38.4250 38.4750 38.5250 38.5750 38.6250 38.6750 38.7250 38.7750 38.8250 3

22、8.8750 38.9250 38.9750 39.0250 y=xy(:,2)y = 221 267 318 487 815 1502 3162 6207 8414 8156 5975 3473 1601 710 386 325 247 232 207 182 1573838.138.238.338.438.538.638.738.838.9390100020003000400050006000700080009000 yys=spline(xrd1(:,1),xrd1(:,2),xx) 函數(shù)函數(shù)z=f(x,y) 二維插值的原理：二維插值的原理：p 向量的采樣點(diǎn)向量的采樣點(diǎn)X,Y，與采樣點(diǎn)對

23、應(yīng)函數(shù)值，與采樣點(diǎn)對應(yīng)函數(shù)值Z；p 向量或標(biāo)量的欲插值點(diǎn)向量或標(biāo)量的欲插值點(diǎn)Xi,Yi，插值結(jié)果，插值結(jié)果Zi。p Xi,Yi取值不超取值不超X,Y給定范圍，否則給定范圍，否則“NaN”錯(cuò)誤。錯(cuò)誤。 zi = interp2(x,y,z,xi,yi,method)例例 運(yùn)行結(jié)果如下圖所示。p 向量的采樣點(diǎn)向量的采樣點(diǎn)X,Y，與采樣點(diǎn)對應(yīng)函數(shù)值，與采樣點(diǎn)對應(yīng)函數(shù)值Z；p 向量或標(biāo)量的欲插值點(diǎn)向量或標(biāo)量的欲插值點(diǎn)Xi,Yi，插值結(jié)果，插值結(jié)果Zi。griddata 算法p Xi,Yi取值不超取值不超X,Y給定范圍，否則給定范圍，否則“NaN”錯(cuò)誤。錯(cuò)誤。 xi,yi,zi = griddata(x

24、,y,z,xi,yi,method)輸入?yún)⒘浚╔I,YI）通常是規(guī)則的格點(diǎn) clear p=rand(30,3) X=p(:,1)Y=p(:,2)Z=p(:,3)Xi,Yi=meshgrid(linspace(min(X),max(X),100)Xi,Yi,Zi=griddata(X,Y,Z,Xi,Yi,)mesh(Xi,Yi,Zi)Xi,Yi,Zi=griddata(X,Y,Z,Xi,Yi,v4)mesh(Xi,Yi,Zi) 隨機(jī)生成30*3矩陣?yán)?.3 函數(shù)極值函數(shù)極值調(diào)用格式一樣調(diào)用格式一樣最小值的解最小值的解或零點(diǎn)零點(diǎn)根根x、該點(diǎn)的函數(shù)值、該點(diǎn)的函數(shù)值fval、程序退出標(biāo)志、程序退出

25、標(biāo)志exitflag輸出結(jié)果輸出結(jié)果output、待求根函數(shù)fun、初始值x0、左右邊界x1,x2x,fval,exitflag,output=)迭代解法迭代解法：最小值的x解或零點(diǎn)的根(自變量值)；：最小值或零點(diǎn)時(shí)函數(shù)值f(x)；待求最小值或根函數(shù)f(x)：函數(shù)名(少用)；字符串表達(dá)式，匿名對象、函數(shù)句柄、內(nèi)聯(lián)函數(shù)；程序退出標(biāo)志, 1收斂到解x；選定輸出結(jié)果；搜索初始值；左右邊界 output = 迭代次數(shù) iterations: 24函數(shù)評價(jià)次數(shù) funcCount: 48算法 algorithm: bisection, interpolation對分插值退出信息 message: Zer

26、o found in the interval -28, 190.51l 配置優(yōu)化參數(shù), optimset函數(shù)定義參數(shù);l 最優(yōu)化工具箱的(20多個(gè))選項(xiàng)設(shè)定;l 將option顯示出來;l 改變其中某個(gè)選項(xiàng);選項(xiàng)函數(shù)調(diào)用時(shí)中間結(jié)果顯示方式l off不顯示;l iter每步都顯示;l final只顯示最終結(jié)果。optimset(display,off)options = optimset(param1,value1.)求區(qū)間求區(qū)間x1 x2內(nèi)函數(shù)內(nèi)函數(shù)f(x)最小值時(shí)最小值時(shí)p X：返回最小值的解;p fval=F(x)，即最小值；p fun：用于定義需求解函數(shù);p x1，x2：范圍；p op

27、tion：的選項(xiàng)設(shè)定x = 2.5148f = -0.4993e = 1o = iterations: 13 funcCount: 14 algorithm: golden section search, parabolic interpolation message: 優(yōu)化已終止: 當(dāng)前的 x 滿足使用 1.000000e-04 的 OPTIONS.Tol.function f=mmfun(x)f=exp(-0.1*x)*(sin(x)2)-0.5*(x+0.1)*sin(x)endx,f,e,o=fminbnd(mmfun,-10,10,optimset(display,iter)0.12

28、sin( )0.5 (0.1) sin( )xyexxx求在（-10，10） 最小值例例x,f,e,o=fminbnd(mmfun,6,10,optimset(display,iter)x = 8.0236f = -3.5680e = 1o = iterations: 9 funcCount: 10 algorithm: golden section search, parabolic interpolation message: 優(yōu)化已終止: 當(dāng)前的 x 滿足使用 1.000000e-04 的 OPTIONS.Tol.黃金分割搜索，拋物線插值 x=-10:0.05:10; y=exp(-0.

29、1*x).*(sin(x).2)-0.5*(x+0.1).*sin(x);plot(x,y) plot(x,y)最小最小最小最小0.12sin( )0.5 (0.1) sin( )xyexxx求在（6，10）最小值-10-8-6-4-20246810-4-3-2-101234 x=-10:0.05:10; y=exp(-0.1*x).*(sin(x).2)-0.5*(x+0.1).*sin(x);plot(x,y) plot(x,y)x,fval=fminbnd(exp(-0.1*x)*(sin(x)2)-0.5*(x+0.1)*sin(x),-10,10)x = 2.514797840754

30、235fval = -0.499312445280039最小最小x,fval=fminbnd(exp(-0.1*x)*(sin(x)2)-0.5*(x+0.1)*sin(x),6,10)x = 8.0236；fval =-3.5680例例x,fval=fminbnd(-exp(-0.1*x)*(sin(x)2)+0.5*(x+0.1)*sin(x),-8,-2)x = -4.830203748934343fval = -3.947274022275747-10-8-6-4-20246810-4-3-2-101234最大值求解：最大值求解：-f(x)在區(qū)間(a,b)上的最小值 x=-10:0.0

31、5:10; y=-exp(-0.1*x).*(sin(x).2)+0.5*(x+0.1).*sin(x);plot(x,y) plot(x,y)最大最大例例p在初始x0附近尋找局部最小值；p使用options選項(xiàng)來指定優(yōu)化器的參數(shù)；著名著名Rosenbrocks Banana 測試函數(shù)測試函數(shù),理論極小值理論極小值x=1,y=1.求求極小值點(diǎn)極小值點(diǎn)x = 1 1fval = 0exitflag =1iterations: 24 funcCount: 49 algorithm: Nelder-Mead simplex direct search message: 優(yōu)化已終止: 當(dāng)前的 x 滿足

32、使用 1.000000e-04 的 OPTIONS.TolX 的終止條件，.例例x,fval,exitflag,output=fminsearch(fxy,1;1)function f=fxy(x)f=100*(x(2)-x(1).2).2+(1-x(1).2endx,fval,exitflag,output=fminsearch(100*(x(2)-x(1).2).2+(1-x(1).2,1;1) 代數(shù)方程代數(shù)方程微分方程微分方程線性方程線性方程非線性方程非線性方程常微分方程常微分方程偏微分方程偏微分方程4.4 非線性方程數(shù)值求解非線性方程數(shù)值求解4.4.1 單變量非線性方程求解單變量非線性

33、方程求解 ()4.4.2 非線性方程組的求解非線性方程組的求解()迭代解法迭代解法(fzero)x,fval,exitflag,output= fzero(fun, x1,x2, options)x,fval,exitflag,output=fzero(fun, x0, options)2.調(diào)用格式：調(diào)用格式：函數(shù)是否穿越橫軸決定零點(diǎn)數(shù)值解，求方程函數(shù)是否穿越橫軸決定零點(diǎn)數(shù)值解，求方程f(x)=0根根 1.解方程原理：解方程原理：3. 求根函數(shù)求根函數(shù)fun【f(x)】的的調(diào)用調(diào)用求f(x)=x-10 x+2=0在x0=0.5附近的根。例例最優(yōu)化的結(jié)構(gòu)體最優(yōu)化的結(jié)構(gòu)體output，最優(yōu)化取下列

34、域：，最優(yōu)化取下列域： algorithm 使用算法使用算法 funcCount 函數(shù)個(gè)數(shù)評估函數(shù)個(gè)數(shù)評估 interval iterations 發(fā)現(xiàn)區(qū)間的迭代次數(shù)發(fā)現(xiàn)區(qū)間的迭代次數(shù) iterations 發(fā)現(xiàn)零值點(diǎn)迭代次數(shù)發(fā)現(xiàn)零值點(diǎn)迭代次數(shù) message 退出信息退出信息 x,f,e,o=fzero(fun,0.5,optimset(display,iter)x =0.3758；f =0；e = 1o = interval iterations: 8 iterations: 5 funcCount: 21 algorithm: bisection, interpolation二分法 me

35、ssage: 在區(qū)間 0.34, 0.613137 中發(fā)現(xiàn)零建立函數(shù)文件fun.m。function fv=fun(x) fv=x-10.x+2;x=-4:0.05:1;f=x-10.x+2;plot(x,f)x,y=ginput(1)ans = -2.0012x,fval=fzero(fun,-1)x = -1.989761447718557fval = 0 (1) 建立函數(shù)文件建立函數(shù)文件fun.m。function fv=fun(x)fv=x-10.x+2;(2) 調(diào)用調(diào)用fzero函數(shù)求根。函數(shù)求根。 x,fval,e,output=fzero(fun,0.5)x = 0.3758 f

36、val=0 x,y=fzero(fun,-3,0.9)x = -1.9898y = 0例例x,f,e,o=fzero(fun,8,optimset(display,iter) x = NaNfval = NaNexitflag =-3正在退出 fzero: 終止搜索包含符號(hào)變化的區(qū)間 因?yàn)樵谒阉髌陂g遇到 NaN 或 Inf 函數(shù)值。請檢查函數(shù)或使用其他起始值重試。o = intervaliter: 22 iterations: 0 funcCount: 44 algorithm: bisection, interpolation message: 正在退出 fzero: 終止搜索包含符號(hào)變化的

37、區(qū)間 因?yàn)樵谒阉髌陂g遇到 NaN 或 Inf 函數(shù)值。例例無根為Nan,（Function handle) x,y=fzero(fun,0.9)x = 0.3758y = 2.2204e-016 x,fval=fzero(x)x-10 x+2,-1,6)x =0.3758fval = 2.2204e-016class(f)function_handle 待求最小值或根函數(shù)f(x)：函數(shù)名(少用)；字符串表達(dá)式，匿名對象、函數(shù)句柄、內(nèi)聯(lián)函數(shù)； x,y,e,o=fzero(x-10.x+2,-1,6,optimset(Display,iter)x =-1.9898y =0e =1o = inter

38、val iterations: 12 iterations: 6 funcCount: 31 algorithm: bisection, interpolation message: Zero found in the interval 0.28, -2.28例例例例inline本質(zhì)和本質(zhì)和function一樣，它直接內(nèi)嵌在命令一樣，它直接內(nèi)嵌在命令行里，不用單獨(dú)定義行里，不用單獨(dú)定義function。函數(shù)表達(dá)式函數(shù)表達(dá)式inlineinline(expression(expression) ) f=inline(x-10.x+2); f(3) ans = -995f=inline(x-10.

39、x+2,x)Z=fzero(f,0.5)或或z=fzero(inline(x-10.x+2),0.5)hd = funfa =(x)x-10 x+2fs =x-10 x+2fi=inline(x-10 x+2) Name Size Bytes Class fa 1x1 16 function_handle fi 1x1 832 inline fs 1x8 16 char hd 1x1 16 function_handle class(f)ans =inline例例 roots(1,0,-2,-5)ans = 2.0946 -1.0473 + 1.1359i -1.0473 - 1.1359i

40、x,fval=fzero(x3-2*x-5,3)x = 2.0946fval = -8.8818e-016例例4.4.2 非線性方程組求解非線性方程組求解x向量；向量；F(x) 函數(shù)向量函數(shù)向量。p x：返回的向量向量解;p fval=F(x)：函數(shù)值向量；函數(shù)值向量；p Jacobian：解x處的Jacobian陣；p fun：用于定義需求解函數(shù);p x0：初始估計(jì)值，列向量列向量；p option：的選項(xiàng)設(shè)定方程組的標(biāo)準(zhǔn)形式：方程組的標(biāo)準(zhǔn)形式：F(x) = 0 x,fval,exitflag,output,jacobian=fsolve(fun,x0,options) function F

41、 = myfun(x)F = 2*x(1)-x(2)-exp(-x(1); -x(1)+2*x(2)-exp(-x(2);endx,fval = fsolve(myfun, -5; -5) 求方程組解求方程組解例例Function F = myfun (x)。 % x為自變量所構(gòu)成的數(shù)組。為自變量所構(gòu)成的數(shù)組。F=表達(dá)式表達(dá)式1；表達(dá)式；表達(dá)式2；表達(dá)式表達(dá)式m 122 12012 20 xxxxexxe在在x1=-5,x2=-5解解 (1) 建立函數(shù)文件建立函數(shù)文件myfun.m function y=myfun(x)y(1)=x(1)-0.6*sin(x(1)-0.3*cos(x(2);y

42、(2)=x(2)-0.6*cos(x(1)+0.3*sin(x(2);1 0.6 sin( 1) 0.3 cos( 2)02 0.6 cos( 1) 0.3 sin( 2)0 xxxxxxx,fval=fsolve(myfun,0.5,0.5,optimset(Display,iter)在在(0.5,0.5) 附近解。附近解。例例(2) 初值初值x0=0.5,y0=0.5下，調(diào)用下，調(diào)用fsolve求方程的根。求方程的根。 x=fsolve(myfun,0.5,0.5,optimset(display,iter)x=0.6354 0.3734將求得的解代回原方程，檢驗(yàn)結(jié)果是否正確，將求得的解代

43、回原方程，檢驗(yàn)結(jié)果是否正確， 可見得到了較高精度的結(jié)果?？梢姷玫搅溯^高精度的結(jié)果。q=myfun(0.6354,0.3734)q = 1.0e-004 * -0.2744 -0.6564例例4.5 常微分方程常微分方程微分方程：y=f(t,y), t0ttn初始條件：y( t0 )= y01. 一階常微分方程一階常微分方程 ： 方程的初值問題數(shù)值解：p 求解y(t)在節(jié)點(diǎn)t0,t1 ,t2,t3.tn, 處近似值y0, y1 , y2，y3. yn；p 采用等距節(jié)點(diǎn)tn=t0+nh，n=0,1,.n,h是步長是步長 微分方程描述：時(shí)間列向量;對應(yīng)t的數(shù)值解(列向量);顯式方程y=f(t,y),

44、或含混合矩陣方程M(t,y), 時(shí)間t范圍，形式t0,tf, 或 初始條件的值,用odeset設(shè)置可選參數(shù) 龍格庫塔法龍格庫塔法 t,y=ode23(odefun,tspan,y0, options) t,y=ode45(odefun,tspan,y0, options)一階一階常微分方程數(shù)值解 設(shè)有初值問題，y=(y2-t-2)/4/(t+1); 0t1;y0(t=0)=2 試求其數(shù)值解，并與精確解(y(t)=sqrt(t+1)+1)比較。(1) 建立函數(shù)文件funt.m。function yp=funt(t,y)yp=(y2-t-2)/4/(t+1);(2) 求解微分方程。t0=0;tf=

45、1;y0=2;t,y=ode23(funt,t0,tf,y0); %求數(shù)值解y1=sqrt(t+1)+1; %求精確解tyy1 y為數(shù)值解，y1為精確值，顯然兩者近似。t=linspace(0,1,11)例例求解著名的求解著名的Rossler微分方程組微分方程組a=b=0.2,c=5.7,x(0)=y(0)=z(0)=0function ddx=rossler(t,x,flag,a,b,c ) ddx=-x(2)-x(3);x(1)+a*x(2);b+(x(1)-c)*x(3);end a=0.2;b=0.2;c=5.7; x0=0,0,0;t,y=ode45(rossler,0,100,x0

46、,a,b,c) plot(t,y) figure plot3(y(:,1),y(:,2),y(:,3)0102030405060708090100-15-10-50510152025-10-5051015-20-100100510152025( )( )( )( )( )( )( ) ( ) ( )x ty tz ty tx tay tz tbx tc z t 例例例例2. 二階常微分方程二階常微分方程 ：x=f(t,x,x)求解振蕩器的經(jīng)典求解振蕩器的經(jīng)典Ver der pol的微分方程的微分方程 令y(1)=x, y(2)=d x/dt2222(1) y(2)(1x )x= (1y(1)

47、)y(2)y(1)(2)dydxdtdtdxdyd xdtdtdt一階微分方程組一階微分方程組:初始條件初始條件：x(t=0)=1，y(1)(t=0)=1x(t=0)=0, y(2)(t=0)=00510152025303540-3-2-10123global MUMU=1t,y=ode23(verderpol1,0,40,1;0);plot(t,y)function yy=verderpol1(t,y)global MUyy=y(2); MU*(1-y(1).2).*y(2)-y(1);0510152025303540-6-4-20246MU=32. 高于高于2階常微分方程階常微分方程的的求

48、解求解基本過程基本過程 （1）規(guī)律、定律、公式）規(guī)律、定律、公式的的描述描述形式：形式：( )( ,., )0nF y y yyt(1)000101( ),( ),.( )nny tyy tyyty初始條件初始條件：微分方程微分方程：（2）(1)12,.,nnyy yyyy高階方程高階方程(組組)轉(zhuǎn)換一階：轉(zhuǎn)換一階：),(),(),(2121ytfytfytfyyyynn110002010)()()()(nnyyytytytyty一階微分方程組： 初始條件： 列向量列向量 (2)(1)(1)( )( )(2)=z(1)(3)z(2)-nz(1)(2)kknnkyyyzzzzzzzzzzyyyy

49、yfzzz(k1)=(k-2)(k)=(k-1)(n)( -1)（，. (n),t）=(n).(k+1)=(k) (0)y=y=z(1)f(z(1)(2)zzz，. (n),t）= (n)( )(1)( )(1)=kkkkyyyzzyzz(k+1)(k)(k+1) (k)n( -1)f,.,nyy y yyt （）（3）根據(jù)(1)與(2)，編寫導(dǎo)數(shù)M函數(shù)文件；（4）將M文件與初始條件傳遞給求解器Solver;（5）運(yùn)行后得到ODE的、在指定時(shí)間區(qū)間解列向量y（其中包含y及不同階的導(dǎo)數(shù)）。t,y=solver(odefun,tspan,y0)solversolver求解器求解器solver奇異矩

50、陣奇異矩陣 奇異矩陣奇異矩陣 函數(shù)指令含 義函 數(shù)含 義求解器S o lverode23普通2-3階法解ODEodefile包含ODE的文件ode23s低階法解剛性O(shè)DE選項(xiàng)odeset創(chuàng) 建 、 更 改Solver選項(xiàng)ode23t解適度剛性O(shè)DEodeget讀取Solver的設(shè)置值ode23tb低階法解剛性O(shè)DE輸出odeplotODE的時(shí)間序列圖ode45普通4-5階法解ODEodephas2ODE的二維相平面圖ode15s變階法解剛性O(shè)DEodephas3ODE的三維相平面圖ode113普通變階法解ODEodeprint在命令窗口輸出結(jié)果求 解 器SolverODE類型特點(diǎn)說明ode45

51、非剛性一步算法；4，5階Runge-Kutta方程；累計(jì)截?cái)嗾`差達(dá)(x)3大部分場合的首選算法ode23非剛性一步算法；2，3階Runge-Kutta方程；累計(jì)截?cái)嗾`差達(dá)(x)3使用于精度較低的情形ode113非剛性多步法；Adams算法；高低精度均可到10-310-6計(jì)算時(shí)間比ode45短ode23t適度剛性采用梯形算法適度剛性情形ode15s剛性多步法；Gears反向數(shù)值微分；精度中等若ode45失效時(shí)，可嘗試使用ode23s剛性一步法；2階Rosebrock算法；低精度當(dāng)精度較低時(shí)，計(jì)算時(shí)間比ode15s短ode23tb 剛性梯形算法；低精度當(dāng)精度較低時(shí)，計(jì)算時(shí)間比ode15s短6. 常

52、微分方程組解法器參數(shù)常微分方程組解法器參數(shù)4.6 數(shù)值積分和微分?jǐn)?shù)值積分和微分4.6.1 數(shù)值積分?jǐn)?shù)值積分 4.6.2 數(shù)值微分?jǐn)?shù)值微分 將區(qū)間將區(qū)間a,b等分等分n個(gè)子區(qū)間個(gè)子區(qū)間xi,xi+1，i=1n，x1=a，xn+1=b。 求定積分就求和問題。求定積分就求和問題。 數(shù)值積分方法：數(shù)值積分方法： 簡單的梯形法簡單的梯形法trapz 辛普生辛普生(Simpson)法法 牛頓柯特斯牛頓柯特斯(Newton-Cotes)法法 1.數(shù)值定積分基本原理數(shù)值定積分基本原理 向量自變量X和應(yīng)變量Y(1) 梯形法數(shù)值積分梯形法數(shù)值積分用 trapz函數(shù)計(jì)算定積分。 X=1:0.01:2.5;Y=exp

53、(-X); %生成函數(shù)關(guān)系數(shù)據(jù)向量trapz(X,Y)ans = 0.285796824163932. 數(shù)值積分的實(shí)現(xiàn)方法數(shù)值積分的實(shí)現(xiàn)方法z=trapz(Y)z =28.5797Z = trapz(Y) Z=trapz(X,Y) q,n=quad(exp(-x),1,2.5)q = 0.2858n = 13例例 ：被積函數(shù)； ：被積函數(shù)調(diào)用次數(shù)；：定積分下限和上限；：控制積分精度，缺省1e-6；q,n=quad(fun,a,b,tol,trace) q=quad(fun,a,b,tol,trace) 非非0展現(xiàn)積分過程；展現(xiàn)積分過程；0不展現(xiàn)，缺省時(shí)不展現(xiàn)，缺省時(shí)trace=0；返回參數(shù)返回

54、參數(shù)I即定積分值即定積分值是否展現(xiàn)積分過程是否展現(xiàn)積分過程P為壓力kpa，V為體積，m3；n為摩爾數(shù)kmol；R為理想氣體常數(shù)8.314kpam3/kmolK;T為溫度。氣缸內(nèi)1mol氣體，溫度為300k，氣體在整個(gè)過程恒溫，體積由V1=1m3膨脹到V2=5m3(/)ln( 2/ 1)WPdVnRT V dVnRTVVPV nRT求解氣缸活塞做功n=input(n=);T=input(T=);R=8.314;pp=(v)n*R*T./v;W=quad(pp,1,5)例例 (1) 建立被積函數(shù)建立被積函數(shù)fesin.m。 function f=fesin(x) f=exp(-0.5*x).*si

55、n(x+pi/6); (2) 調(diào)用數(shù)值積分函數(shù)調(diào)用數(shù)值積分函數(shù)quad S,n=quad(fesin,0,3*pi) S = 0.9008 n =77求定積分求定積分quad(-0.5*x).*sin(x+pi/6),0,3*pi)q,n=quad(sqrt(1+cos(x).2),0,pi/2,1.0e-6,1)例例s=quad(0.2+sin(x),0,pi/2,1e-8); x=0:pi/10:pi/2; y=0.2+sin(x); s1=sum(y); ss=s1*pi/10; st=trapz(x,y)format short disp(quad積分積分,blanks(4),sum求

56、積分求積分,blanks(4),trapz求積分求積分)disp(s,ss,st)hold onplot(x,y,linewidth,2) stairs(x,y,-r*) stem(x,y,-gh)quad積分 sum求積分 trapz求積分 1.3142 1.3578 1.3059例例 參數(shù)含義和參數(shù)含義和quad相似；相似； tol缺省值缺省值10-6； 調(diào)用調(diào)用步數(shù)明顯小于步數(shù)明顯小于quad，高效率高效率求值；求值； 積分精度更高。積分精度更高。 求定積分求定積分(1) 被積函數(shù)文件被積函數(shù)文件fx.m。function f=fx(x)f=x.*sin(x)./(1+cos(x).*c

57、os(x);(2) 調(diào)用函數(shù)調(diào)用函數(shù)quad8求定積分。求定積分。I=quad8(fx,0,pi)I = 2.4674例例分別用分別用quad和和quadl求定積分的近似值，并在相同的求定積分的近似值，并在相同的積分精度下，比較函數(shù)的調(diào)用次數(shù)。積分精度下，比較函數(shù)的調(diào)用次數(shù)。調(diào)用函數(shù)quad求定積分：format long;fx=inline(exp(-x);I,n=quad(fx,1,2.5,1e-10)I = 0.28579444254766n = 65調(diào)用函數(shù)quad8求定積分：format long;fx=inline(exp(-x);I,n=quad8(fx,1,2.5,1e-10)

58、I = 0.28579444254754n = 33例例 在在 xmin,xmax ymin,ymax 區(qū)域二重定積分；區(qū)域二重定積分； 參數(shù)參數(shù)toltol，tracetrace的用法與函數(shù)的用法與函數(shù)quadquad相同；相同； toltol或或methodmethod可以忽略?？梢院雎?。q = dblquad(fun,xmin,xmax,ymin,ymax,q = dblquad(fun,xmin,xmax,ymin,ymax,tol,methodtol,method) ) (1) 建立函數(shù)文件建立函數(shù)文件fxy.m：function f=fxy(x,y)global ki;ki=ki+

59、1; %ki用于統(tǒng)計(jì)被積函數(shù)的調(diào)用次數(shù)f=exp(-x.2/2).*sin(x.2+y);(2) 調(diào)用調(diào)用dblquad函數(shù)函數(shù)global ki;ki=0;I=dblquad(fxy,-2,2,-1,1)kiI = 1.57449318974494ki = 1038 計(jì)算二重積分計(jì)算二重積分f=inline(exp(-x.2/2).*sin(x.2+y),x,y);dblquad(f,-2,2,-1,1)ans =1.5745dblquad(f,-2,2,-1,1 )ans=1.5745dblquad(exp(-x.2/2).*sin(x.2+y),-2,2,-1,1)例例triplequa

60、d(fun,xmin,xmax,ymin,ymax,zmin,zmaxtriplequad(fun,xmin,xmax,ymin,ymax,zmin,zmax, ,tol,methodtol,method) )例例 計(jì)算三重積分計(jì)算三重積分 triplequad(x,y,z)y.*sin(x)+z.*cos(x),0,pi,0,1,-1,1)ans =1.999999994362637 triplequad(inline(y.*sin(x)+z.*cos(x),0,pi,0,1,-1,1)ans =1.999999994362637 triplequad(y.*sin(x)+z.*cos(x)