第四章數值計算功能4.1多項式計算4.2數值插值和曲線擬合4.3函數極值4.4非線性方程問題求解4.5常微分方程的數值求解4.6數值微分與積分4.7概率統計4.1多項式運算

1．多項式表示法2．多項式運算函數

多項式求根多項式求值多項式乘法和除法多項式的微積分1．多項式表示法

(1)直接法:多項式多項式系數用行向量表示，按降冪排列；缺某冪次項，其冪次項系數為零。P=[a0,a1,...an-1,an]符串形式,x多項式變量y=poly2str(P,'x')(2)字符串表示法：y=poly2sym(P,x)(3)符號多項式表示法：完整形式p=[1,3,0,4,5]y=poly2str(p,'x')y=x^4+3x^3+4x+5symsxy=poly2sym(p,x)y=x^4+3*x^3+4*x+5>>p=[1,3,0,4,5]p=13045>>y=poly2str(p,'x')y=x^4+3x^3+4x+5>>rp=roots(p)rp=-3.23460.5594+1.1980i0.5594-1.1980i-0.8843>>P=[1,8,3,4,-10];>>pp=poly2str(P,'X')pp=X^4+8X^3+3X^2+4X-10例2.多項式的運算函數多項式的值；根和微分；擬合曲線；部分分式polyvalpolyvalpolyvalmcovolution

（1）多項式求根（roots）（1）多項式的根=0

n次多項式n個根，或實根，或若干對共軛復根。

r=roots(P)P：多項式系數向量；r：n個根b(1),b(2),…,b(n)的向量。>>r=roots(P)r=-7.6998-0.5572+1.1335i-0.5572-1.1335i0.8141例求多項式x4+8x3+3x2+4x-10的根P=[1,8,3,4,-10];r=roots(P)多項式根的r向量：[r(1),r(2),…,r(n)]（2）根行向量創建多項式（poly）P=poly(r)>>PP=poly(r)PP=1.00008.00003.00004.0000-10.0000>>formatrat避免浮點顯示>>p=poly(r)p=1834-10多項式的向量系數：(3)多項式求值（polyval）Y=polyval(P,x)若x數值，則求多項式在該點x的值；若X是向量或矩陣，每元素x(i,j)多項式求值。代數多項式P求值數組乘方運算:Y=a0*X.^n+a1*X.^(n-1)…an+1>>p=[1,3,0,2];>>poly2str(p,'x')ans=x^3+3x^2+2>>a=[1,2;3,4]a=1234>>y=polyval(p,a)y=62256114>>y=polyval(p,2)y=22例例

Y=polyvalm(P,X)矩陣多項式求值（polyvalm）Y=a0*X^n+a1*X^n-1…an+1X必為方陣，作自變量代入多項式求值；結果階數還是保持方陣階數。相當于矩陣X乘方運算：設A為方陣，P代表多項式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,0]a=2410>>pp=polyval(p,a)pp=-14-18-6-2>>pm=polyvalm(p,a)pm=-18-8-2-14矩陣乘方運算:數組乘方運算:例例14(4)多項式加減運算相同次數多項式，加減其系數向量,

不同次數,低次多項式中高次項用0補足。

p2=0x3-

0x2+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例例(5)多項式乘法-向量卷積w=conv(P1,P2)多項式x4+8x3-10與2x2-x+3的乘積。>>p1=[1,8,0,0,-10];>>p2=[2,-1,3];>>p12=conv(p1,p2)p12=Columns1through5215-524-20Columns6through710-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('Timeindexn');%標坐標軸>>ylabel('Amplitude');>>title('OutputObtainedbyConvolution')%圖形標題grid;卷積是兩個變量在某范圍內相乘后求和的結果。卷積的變量是序列x(n)和h(n)，y=conv(x,h)來實現卷級的，輸出的結果個數等于x的長度與h的長度之和減去1。卷積公式：z(n)=x(n)*y(n)=∫x(m)y(n-m)dm.例(6)多項式除法-向量解卷積[Q,r]=deconv(P1,P2)Q：多項式P1除以P2的商式，r：P1除以P2的余式，Q和r仍是多項式系數向量。P1=conv(P2,Q)+rdeconv是conv的逆函數求多項式x4+8x3-10除以2x2-x+3的結果。>>[q,r]=deconv(p12,p2)q=1800-10r=Columns1through500000Columns6through700>>[q,r]=deconv(p2,p12)q=0r=2-13例(7)多項式的一階導數（polyder）

k=polyder(P)k=polyder(P,Q)多項式乘積P·Q的導函數：多項式P的導函數：導函數的分子存入p，分母存入q。[p,q]=polyder(P,Q)多項式除法P/Q的導函數：例：求有理分式的導數。命令如下：P=[1];Q=[1,0,5];[p,q]=polyder(P,Q)例例21(8)多項式積分(polyint)I=polyint([2,-1,0,3],5);例：

p(x)

=2x3-

x2+3求常數項5k=polyint(P,c)k=polyint(P)不定積分，常數項取c不定積分，常數項取零例

4.2曲線擬合和數據插值4.2.1曲線擬合polyfit

4.2.2數據插值interp

擬合成多項式推算未知點值4.2.1.多項式曲線擬合（polyfit）

最小二乘法(誤差平方和最小)；x、y已知數據的橫、縱坐標；n為多項式次數;p=polyfit(x,y,n)1.N次多項式曲線擬合S結構體數組(struct)，估計預測誤差，含R，df和normr。

R：先輸入x構建范德蒙矩陣V，后QR分解，得上三角矩陣。

df：自由度，

df=length(y)-(n+1)。

df>0時，超定方程組求解，擬合點數比未知數(p(1)~p(n+1))多。

normr：標準偏差、殘差范數，normr=norm(y-V*p)，

此處的p為求解之后的數值。

[p,S]=polyfit(x,y,n)年份1900191019201930194019501960197019801990人口75.9991.97105.71123.20131.66150.69179.32203.21226.50249.63對不同年份人口數據分別進行1次、2次和9次多項式擬合（polyfit），用poly2str表示多項式完整形式法；1次、2次和9次多項式估計2000年人口(polyval)，結合美國2000年人口普查截至2000年4月1日美國人口2.81421906億數據繪制1900~2000年間的時間-人口數(polyval)曲線，要求用plot不同線型（LineSpec）繪制原始數據點、擬合的1次、2次和9次多項式曲線，標注圖例，說明是否階數越到高越好美國從1900年~1990年每隔10年進行人口普查的數據如下表所示：（百萬）例線型說明數據點標記符說明顏色說明-實線(默認)+加號符r紅色--雙劃線*星號g綠色-.點劃線.實心圓b藍色:虛線o空心圓c青綠色x叉號符m洋紅色

s正方形y黃色

d菱形k黑色

^上三角形w白色

v下三角形

>

右三角形

<

左三角形

p五角星

h六邊形

plot(x,y,‘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次多項式擬合nrfs1=poly2str(nrf1,'n')%1次多項式完整字符串表達式nrf2=polyfit(n,r,2)%2次多項式擬合nrf9=polyfit(n,r,9)%9次多項式擬合nrf10=polyfit(n,r,10)%10次多項式擬合nrf1_2000=polyval(nrf1,2000)%一次擬合得到的2000年人口數nrf2_2000=polyval(nrf2,2000)%2次擬合得到的2000年人口數nrf9_2000=polyval(nrf9,2000)%9次擬合得到的2000年人口數n20=1900:4:2000;%1900年到2000年的線性等分數組nrfv1=polyval(nrf1,n20);%從1900年到2000年間1次擬合人口數nrfv2=polyval(nrf2,n20);%從1900年到2000年間2次擬合人口數nrfv9=polyval(nrf9,n20);%從1900年到2000年間9次擬合人口數

plot(n,r,'or',n20,nrfv1,'-<b',n20,nrfv2,'-*k')holdonplot(n20,nrfv9,'-dg')legend('原始數據’,‘1次曲線','2次多項式曲線','9次多項式曲線');xlabel('年');ylabel('人口數')2009版多項式次數要適當，過低誤差大，過高波動明顯2014版Goodnessoffit

適合度

SSE擬合誤差

（誤差平方和）

RMSErootmeansquareerror均方根誤差

Rsquare

方程的確定系數，0~1之間，越接近1，表明方程的變量對y的解釋能力越強。

2.曲線擬合工具箱cftool

R-square=1，說明擬合的結果很好CustomEquations：用戶自定義的函數類型

Exponential：指數逼近，2種：a*exp(b*x)、a*exp(b*x)+c*exp(d*x)

Fourier：傅立葉逼近，7種，基礎型是a0+a1*cos(x*w)+b1*sin(x*w)

Gaussian：高斯逼近，8種，基礎型是a1*exp(-((x-b1)/c1)^2)

Interpolant：插值逼近，4種，linear、nearestneighbor、cubicspline、shape-PreservingPolynomial：多形式逼近，9種，linear~、quadratic~、cubic~、4-9thdegree~

工具箱提供的擬合類型Power：冪逼近，有2種類型，a*x^b、a*x^b+c

Rational：有理數逼近，分子、分母共有的類型是linear~、quadratic~、cubic~、4-5thdegree~；此外，分子還包括constant型

SmoothingSpline：平滑逼近（翻譯的不大恰當，不好意思）

SumofSinFunctions：正弦曲線逼近，有8種類型，基礎型是a1*sin(b1*x+c1)

Weibull：只有一種，a*b*x^(b-1)*exp(-a*x^b)擬合工具cftool

x12345678910y1311122832457080104x=1:10;y=[1,3,11,12,28,32,45,70,80,104];%1次多項式擬合p1=polyfit(x,y,1)Ps1=poly2str(p1,'x')xx=1:0.05:10;y1=polyval(p1,xx);%5次多項式擬合p2=polyfit(x,y,5)Ps2=poly2str(p2,'x')y2=polyval(p2,xx);%9次多項式擬合p3=polyfit(x,y,9)Ps3=poly2str(p3,'x')y3=polyval(p3,xx);>>p4=polyfit(x,y,11)y4=polyval(p4,xx);Warning:Polynomialisnotunique;degree>=numberofdatapoints.例holdon;plot(x,y,'*');plot(xx,y1,'r-');plot(xx,y2,'g-');plot(xx,y3,'b-')多項式次數要適當，過低誤差大，過高波動明顯已知：10個實驗數據，分別采用2次、5次、9次和11次擬合，選出最佳擬合次數3.函數擬合>>clear>>x=[0:5];>>y=[0.2097,0.3523,0.4339,0.5236,0.7590,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)例例插值構造的函數必須通過已知數據點；擬合則不要求,只要均方差最小。4.數據擬合和插值相同點都需根據已知數據構造函數；可使用得到函數計算未知點的函數值。不同點4.2.2數據插值構造函數近似表達式的方法。常用多項式作插值函數,稱多項式插值。多項式插值基本思想：高次多項式或分段的低次多項式為被插值函數f(x)的近似解析表達式。1.插值法根據已知輸入/輸出數據集和當前輸入估計輸出值2.插值函數多項式插值法：拉格朗日插值法

牛頓插值法

埃爾米特插值法分段插值法樣條插值法等

3.一維插值函數y=f(x)一維插值原理：調用格式：

yi=interp1(x,y,xi,'method')已知X,Y兩個等長向量，采樣點和樣本值；Xi是向量或標量，欲插值點；Yi是一個與Xi等長的插值結果。采用差值法估計美國的人口數量2000年人口及1900~1996年每隔4的人口數(interp1)；將上題原始數據點、2階擬合曲線和插值曲線繪制在同一圖，標注坐標軸為‘時間’和‘人口’（xlabel、ylabel）、圖形標題（title）為‘插值和擬合’和圖例legend年份1900191019201930194019501960197019801990人口75.9991.97105.71123.20131.66150.69179.32203.21226.50249.63例clearn=1900:10:1990;r=[7599,9197,10571,12320,13166,15069,17932,20321,22650,24963];n4=1900:4:1996;%1900~1996年每隔4的線性等分數組nrf2=polyfit(n,r,2)%2次多項式擬合nrfv2=polyval(nrf2,n4)nr4=interp1(n,r,n4,'spline')%插值1900~1996年每隔4的人口數nr2000=interp1(n,r,2000)%

線性插值2000年人口數nr2000=interp1(n,r,2000,'spline')%

插值2000年人口數plot(n,r,'or',n4,nrfv2,'-*k',n4,nr4,'-dg')%繪圖legend('原始數據','2階多項式擬合曲線','插值曲線')%圖例

例分段線性插值linear：插值點樣本值落在兩相鄰數據點的連線上.

(缺省).最鄰近插值法nearest：兩點間插值點對應值就是離兩點最近那點值。

3次多項式插值cubic：已知數據求3次多項式，用多項式求插值。

3次樣條插值spline：每分段內構造一3次多項式，插值函數除滿足插值條件和節點處光滑條件，再按照樣條函數插值。插值方法method

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(x,y,'-*',xx,y4,'dr')X1取值范圍不能超出X給定范圍，否則會給出“NaN”錯誤。例38.025 22138.075 26738.125 31838.175 48738.225 81538.275 150238.325 316238.375 620738.425 841438.475 815638.525 597538.575 347338.625 160138.675 71038.725 38638.775 32538.825 24738.875 23238.925 20738.975 18239.025 157插值解法XRD：nano-sic/Alcomposiste例plot(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.025038.075038.125038.175038.225038.275038.325038.375038.425038.475038.525038.575038.625038.675038.725038.775038.825038.875038.925038.975039.0250>>y=xy(:,2)y=22126731848781515023162620784148156597534731601710386325247232207182157yys=spline(xrd1(:,1),xrd1(:,2),xx)4.二維插值函數z=f(x,y)二維插值的原理：向量的采樣點X,Y，與采樣點對應函數值Z；向量或標量的欲插值點Xi,Yi，插值結果Zi。指定插值方法method：

‘linear’‘nearest’‘cubic’‘spline’

Xi,Yi取值不超X,Y給定范圍，否則“NaN”錯誤。

zi=interp2(x,y,z,xi,yi,method)函數interp2二維插值：例運行結果如下圖所示。向量的采樣點X,Y，與采樣點對應函數值Z；向量或標量的欲插值點Xi,Yi，插值結果Zi。指定插值方法method：linear(缺省)nearestcubicv4griddata算法

Xi,Yi取值不超X,Y給定范圍，否則“NaN”錯誤。

[xi,yi,zi]=griddata(x,y,z,xi,yi,method)函數griddata二維柵格插值：輸入參量（XI,YI）通常是規則的格點>>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,'cubic')mesh(Xi,Yi,Zi)[Xi,Yi,Zi]=griddata(X,Y,Z,Xi,Yi,'v4')mesh(Xi,Yi,Zi)

隨機生成30*3矩陣例4.3函數極值1.函數最小值和零點調用格式一樣[x,fval,exitflag,output]=

fminbnd(fun,x1,x2,options)最小值的解或零點根x、該點的函數值fval、程序退出標志exitflag輸出結果output、待求根函數fun、初始值x0、左右邊界x1,x2fval,exitflag,output和options可省略[x,fval,exitflag,output]=fminsearch(fun,x0,options)[x,fval,exitflag,output]=fzero(fun,x0,options)迭代解法輸出量x：最小值的x解或零點的根(自變量值)；fval=f(x)：最小值或零點時函數值f(x)；fun：待求最小值或根函數f(x)：函數名(少用)；字符串表達式，匿名對象、函數句柄、內聯函數；

exitflag：程序退出標志,

1收斂到解x；output：選定輸出結果；x0：搜索初始值；x1，x2：左右邊界

output=

迭代次數iterations:24函數評價次數funcCount:48算法algorithm:‘bisection,interpolation’對分插值退出信息message:'Zerofoundintheinterval[-28,190.51]'option：配置優化參數,

optimset函數定義參數;最優化工具箱的(20多個)選項設定;optimset命令:將option顯示出來;改變其中某個選項;display:選項函數調用時中間結果顯示方式‘off’不顯示;‘iter’每步都顯示;‘final’只顯示最終結果。optimset(‘display’,‘off’)options=optimset('param1',value1...)[x,fval,exitflag,output]=

fminbnd(fun,x1,x2,options)2.求一元函數最小值求區間[x1x2]內函數f(x)最小值時x：X：返回最小值的解;fval=F(x)，即最小值；fun：用于定義需求解函數;x1，x2：范圍；option：最優化工具箱的選項設定x=2.5148f=-0.4993e=1o=iterations:13funcCount:14algorithm:'goldensectionsearch,parabolicinterpolation'message:'優化已終止:當前的x滿足使用1.000000e-04的OPTIONS.Tol...'functionf=mmfun(x)f=exp(-0.1*x)*(sin(x)^2)-0.5*(x+0.1)*sin(x)end[x,f,e,o]=fminbnd('mmfun',-10,10,optimset('display','iter'))求在（-10，10）最小值例[x,f,e,o]=fminbnd('mmfun',6,10,optimset('display','iter'))x=8.0236f=-3.5680e=1o=iterations:9funcCount:10algorithm:'goldensectionsearch,parabolicinterpolation'message:'優化已終止:當前的x滿足使用1.000000e-04的OPTIONS.Tol...'黃金分割搜索，拋物線插值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)Fminbnd在指定區域找到真解最小最小求在（6，10）最小值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.514797840754235fval=-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.5680Fminbnd在指定區域找到真解例[x,fval]=fminbnd('-exp(-0.1*x)*(sin(x)^2)+0.5*(x+0.1)*sin(x)',-8,-2)x=-4.830203748934343fval=-3.947274022275747最大值求解：-f(x)在區間(a,b)上的最小值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)最大例3.求多元函數的最小值在初始x0附近尋找局部最小值；使用options選項來指定優化器的參數；fval,exitflag,output和options可以省略。[x,fval,exitflag,output]=fminsearch(fun,x0,options)著名Rosenbrock's"Banana"測試函數,理論極小值x=1,y=1.求極小值點

x=11fval=0exitflag=1iterations:24funcCount:49algorithm:'Nelder-Meadsimplexdirectsearch'message:'優化已終止:

當前的x滿足使用1.000000e-04的OPTIONS.TolX的終止條件，...例[x,fval,exitflag,output]=fminsearch('fxy',[1;1])functionf=fxy(x)f=100*(x(2)-x(1).^2).^2+(1-x(1)).^2end[x,fval,exitflag,output]=fminsearch('100*(x(2)-x(1).^2).^2+(1-x(1)).^2',[1;1])方程代數方程微分方程線性方程非線性方程常微分方程偏微分方程4.5常微分方程的求解2.5線性方程組求解4.4非線性方程求解4.4非線性方程數值求解4.4.1單變量非線性方程求解(fzero)4.4.2非線性方程組的求解(fsolve)迭代解法4.4.1單變量非線性方程求解(fzero)[x,fval,exitflag,output]=fzero(fun,x1,x2,options)[x,fval,exitflag,output]=fzero(fun,x0,

options)2.調用格式：函數是否穿越橫軸決定零點數值解，求方程f(x)=0根

1.解方程原理：fval,exitflag,output和options可省略3.求根函數fun【f(x)】的調用求f(x)=x-10x+2=0在x0=0.5附近的根。函數名：fzero（

‘funname’，x0）例最優化的結構體output，最優化取下列域：

algorithm使用算法

funcCount函數個數評估

interval

iterations發現區間的迭代次數

iterations發現零值點迭代次數

message退出信息

[x,f,e,o]=fzero('fun',0.5,optimset('display','iter'))x=0.3758；f=0；e=1o=intervaliterations:8iterations:5funcCount:21algorithm:‘bisection,interpolation‘二分法’message:'在區間[0.34,0.613137]中發現零'建立函數文件fun.m。functionfv=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.0012[x,fval]=fzero('fun',-1)x=-1.989761447718557fval=0(1)建立函數文件fun.m。functionfv=fun(x)fv=x-10.^x+2;(2)調用fzero函數求根。

[x,fval,e,output]=fzero('fun',0.5)x=0.3758fval=0>>[x,y]=fzero('fun',-3,0.9)x=-1.9898y=0例多個根只求x0最近的根[x,f,e,o]=fzero('fun',8,optimset('display','iter'))

x=NaNfval=NaNexitflag=-3正在退出fzero:終止搜索包含符號變化的區間因為在搜索期間遇到NaN或Inf函數值。請檢查函數或使用其他起始值重試。o=intervaliter:22iterations:0funcCount:44algorithm:'bisection,interpolation'message:'正在退出fzero:終止搜索包含符號變化的區間因為在搜索期間遇到NaN或Inf函數值。例無根為Nan,函數句柄（Functionhandle)>>[x,y]=fzero(@fun,0.9)x=0.3758y=2.2204e-016匿名函數：F=@(變量表)函數表達式>>[x,fval]=fzero(@(x)x-10^x+2,-1,6)x=0.3758fval=2.2204e-016class(f)function_handlehfun=@+函數名fun：待求最小值或根函數f(x)：函數名(少用)；字符串表達式，匿名對象、函數句柄、內聯函數；

字符串表達式：fzero(‘expression’,x0)>>[x,y,e,o]=fzero('x-10.^x+2',-1,6,optimset('Display','iter'))x=-1.9898y=0e=1o=intervaliterations:12iterations:6funcCount:31algorithm:'bisection,interpolation'message:'Zerofoundintheinterval[0.28,-2.28]'例例inline本質和function一樣，它直接內嵌在命令行里，不用單獨定義function。函數表達式內聯函數inlineinline(‘expression’)>>f=inline('x-10.^x+2');>>f(3)ans=-995f=inline('x-10.^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+2‘fi=inline('x-10^x+2')NameSizeBytesClassfa1x116function_handlefi1x1832inlinefs1x816charhd1x116function_handleclass(f)ans=inline例4.多項式求的根（roots）>>roots([1,0,-2,-5])ans=2.0946-1.0473+1.1359i-1.0473-1.1359i>>[x,fval]=fzero('x^3-2*x-5',3)x=2.0946fval=-8.8818e-016求解多項式例4.4.2非線性方程組求解x向量；F(x)

函數向量。x：返回的向量解;fval=F(x)：函數值向量；Jacobian：解x處的Jacobian陣；fun：用于定義需求解函數;x0：初始估計值，列向量；option：最優化工具箱的選項設定方程組的標準形式：F(x)=0[x,fval,exitflag,output,jacobian]=fsolve(fun,x0,options)functionF=myfun(x)F=[2*x(1)-x(2)-exp(-x(1));-x(1)+2*x(2)-exp(-x(2))];end[x,fval]=fsolve(@myfun,[-5;-5])求方程組解例FunctionF=myfun(x)。%x為自變量所構成的數組。F=[表達式1；表達式2；…表達式m]在x1=-5,x2=-5解

(1)建立函數文件myfun.m

functiony=myfun(x)y(1)=x(1)-0.6*sin(x(1))-0.3*cos(x(2));y(2)=x(2)-0.6*cos(x(1))+0.3*sin(x(2));[x,fval=fsolve('myfun',[0.5,0.5]',optimset('Display','iter'))在(0.5,0.5)附近解。例(2)初值x0=0.5,y0=0.5下，調用fsolve求方程的根。>>x=fsolve('myfun',[0.5,0.5]',optimset(‘display','iter'))x=0.63540.3734將求得的解代回原方程，檢驗結果是否正確，

可見得到了較高精度的結果。q=myfun([0.6354,0.3734])q=1.0e-004*-0.2744-0.6564例4.5常微分方程微分方程：y’=f(t,y),t0≤t≤tn初始條件：y(t0)=y01.一階常微分方程：方程的初值問題數值解：求解y(t)在節點t0,t1,t2,t3….tn,處近似值y0,y1,y2，y3….yn；采用等距節點tn=t0+nh，n=0,1,..n,h是步長微分方程描述：(ODE)t：時間列向量;y:對應t的數值解(列向量);odefun：顯式方程y’=f(t,y),或含混合矩陣方程M(t,y),tspan：時間t范圍，形式t0,tf,或[t0,tf],

y0：初始條件的值,options：用odeset設置可選參數龍格－庫塔法

[t,y]=ode23(odefun,tspan,y0,options)[t,y]=ode45(odefun,tspan,y0,options)一階常微分方程數值解

設有初值問題，y’=(y^2-t-2)/4/(t+1);0≤t≤1;y0(t=0)=2試求其數值解，并與精確解(y(t)=sqrt(t+1)+1)比較。(1)建立函數文件funt.m。functionyp=funt(t,y)yp=(y^2-t-2)/4/(t+1);(2)求解微分方程。t0=0;tf=1;y0=2;[t,y]=ode23('funt',[t0,tf],y0);%求數值解y1=sqrt(t+1)+1;%求精確解t'y'y1'y為數值解，y1為精確值，顯然兩者近似。t=linspace(0,1,11)例求解著名的Rossler微分方程組a=b=0.2,c=5.7,x(0)=y(0)=z(0)=0functionddx=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,[],a,b,c)>>plot(t,y)>>figure>>plot3(y(:,1),y(:,2),y(:,3))例例2.二階常微分方程：x’’=f(t,x,x’)

求解振蕩器的經典Verderpol的微分方程令y(1)=x,y(2)=dx/dty=[y(1),y(2)],y’=[y(1)’,y(2)’]一階微分方程組:初始條件：x(t=0)=1，y(1)(t=0)=1x’(t=0)=0,y(2)(t=0)=0Y0=[1;0]globalMUMU=1[t,y]=ode23('verderpol1',[0,40],[1;0]);plot(t,y)functionyy=verderpol1(t,y)globalMUyy=[y(2);

MU*(1-y(1).^2).*y(2)-y(1)];MU=32.高于2階常微分方程的求解基本過程

（1）規律、定律、公式的描述形式：初始條件：微分方程：（2）高階方程(組)轉換一階：一階微分方程組：

初始條件：

列向量（3）根據(1)與(2)，編寫導數M函數文件；（4）將M文件與初始條件傳遞給求解器Solver;（5）運行后得到ODE的、在指定時間區間解列向量y（其中包含y及不同階的導數）。3.顯式常微分方程組解法對比[t,y]=solver(odefun,tspan,y0)solversolver求解器solver奇異矩陣奇異矩陣4.求解器Solver與方程組的關系函數指令含

義函

數含

義求解器Solverode23普通2-3階法解ODEodefile包含ODE的文件ode23s低階法解剛性ODE選項odeset創建、更改Solver選項ode23t解適度剛性ODEodeget讀取Solver的設置值ode23tb低階法解剛性ODE輸出odeplotODE的時間序列圖ode45普通4-5階法解ODEodephas2ODE的二維相平面圖ode15s變階法解剛性ODEodephas3ODE的三維相平面圖ode113普通變階法解ODEodeprint在命令窗口輸出結果5.不同求解器Solver的特點求解器SolverODE類型特點說明ode45非剛性一步算法；4，5階Runge-Kutta方程；累計截斷誤差達(△x)3大部分場合的首選算法ode23非剛性一步算法；2，3階Runge-Kutta方程；累計截斷誤差達(△x)3使用于精度較低的情形ode113非剛性多步法；Adams算法；高低精度均可到10-3～10-6計算時間比ode45短ode23t適度剛性采用梯形算法適度剛性情形ode15s剛性多步法；Gear’s反向數值微分；精度中等若ode45失效時，可嘗試使用ode23s剛性一步法；2階Rosebrock算法；低精度當精度較低時，計算時間比ode15s短ode23tb剛性梯形算法；低精度當精度較低時，計算時間比ode15s短6.常微分方程組解法器參數4.6數值積分和微分4.6.1數值積分

4.6.2數值微分將區間[a,b]等分n個子區間[xi,xi+1]，i=1…n，x1=a，xn+1=b。求定積分就求和問題。數值積分方法：簡單的梯形法trapz

辛普生(Simpson)法

牛頓－柯特斯(Newton-Cotes)法4.6.1數值積分1.數值定積分基本原理向量自變量X和應變量Y(1)梯形法數值積分trapz用trapz函數計算定積分。

X=1:0.01:2.5;Y=exp(-X);%生成函數關系數據向量trapz(X,Y)ans=0.285796824163932.數值積分的實現方法z=trapz(Y)z=28.5797Z=trapz(Y)Z=trapz(X,Y)

y向量和>>[q,n]=quad('exp(-x)',1,2.5)q=0.2858n=13例

fun：被積函數；n：被積函數調用次數；a和b：定積分下限和上限；tol：控制積分精度，缺省1e-6；trace[q,n]=quad(fun,a,b,tol,trace)q=quad(fun,a,b,tol,trace)

(2).變步長辛普生法（自適應Simpleson）非0展現積分過程；0不展現，缺省時trace=0；返回參數I即定積分值是否展現積分過程P為壓力kpa，V為體積，m3；n為摩爾數kmol；R為理想氣體常數8.314kpam3/kmolK;T為溫度。氣缸內1mol氣體，溫度為300k，氣體在整個過程恒溫，體積由V1=1m3膨脹到V2=5m3求解氣缸活塞做功n=input('n=');T=input('T=');R=8.314;pp=@(v)n*R*T./v;W=quad(pp,1,5)例

(1)建立被積函數fesin.m。functionf=fesin(x)f=exp(-0.5*x).*sin(x+pi/6);(2)調用數值積分函數quad

[S,n]=quad('fesin',0,3*pi)S=0.9008n=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)formatshort>>disp(['quad積分',blanks(4),'sum求積分',blanks(4),'trapz求積分'])disp([s,ss,st])holdonplot(x,y,'linewidth',2)>>stairs(x,y,'-r*')>>stem(x,y,'-gh')quad積分sum求積分trapz求積分

1.31421.35781.3059例參數含義和quad相似；tol缺省值10-6；調用步數明顯小于quad，高效率求值；積分精度更高。[I,n]=quadl(fun,a,b,tol,trace)(3)牛頓－柯特斯法

求定積分(1)被積函數文件fx.m。functionf=fx(x)f=x.*sin(x)./(1+cos(x).*cos(x));(2)調用函數quad8求定積分。I=quad8('fx',0,pi)I=2.4674

例分別用quad和quadl求定積分的近似值，并在相同的積分精度下，比較函數的調用次數。調用函數quad求定積分：formatlong;fx=inline('exp(-x)');[I,n]=quad(fx,1,2.5,1e-10)I=0.28579444254766n=65調用函數quad8求定積分：formatlong;fx=inline('exp(-x)');[I,n]=quad8(fx,1,2.5,1e-10)I=0.28579444254754n=33例在[xmin,xmax]×[ymin,ymax]區域二重定積分；參數tol，trace的用法與函數quad相同；tol或method可以忽略。3.二重定積分的數值求解q=dblquad(fun,xmin,xmax,ymin,ymax,tol,method)

(1)建立函數文件fxy.m：functionf=fxy(x,y)globalki;ki=ki+1;%ki用于統計被積函數的調用次數f=exp(-x.^2/2).*sin(x.^2+y);(2)調用dblquad函數globalki;ki=0;I=dblquad('fxy',-2,2,-1,1)kiI=1.57449318974494ki=1038

計算二重積分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)例4.三重積分數值求解triplequad(fun,xmin,xmax,ymin,ymax,zmin,zmax,tol,method)例計算三重積分>>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)',0,pi,0,1,-1,1)tol或method可以忽略例4.7.2數值微分向前差分⊿f(x)=f(X+h)-f(x)向后差分⊿f(x)=f(X)-f(x-h)中心差分⊿f(x)=f(X+h/2)-f(x-h/2)h充分小差分1.數值差分與差商xX+hX-hX-h/2X+h/2f(x+h

）f(x）f(x-h）差商導數2.數值微分函數diff向量X向前差分：Y=diff(X)向量X的n階向前差分：Y=diff(X,n)矩陣A的n階差分：dX(i)=X(i+1)-X(i)，i=1,2,…,n-1。diff(X,2)=diff(diff(X))Y=diff(A,n,dim)dim=1時(缺省狀態)，按行計算；dim=2按列計算。

生成以向量V=[1,2,3,4,5,6]為基礎的范得蒙矩陣，按列進行差分運算。V=vander(1:6)DV=diff(V)%計算V的一階差分V=1111116842181279312566416416251252551dv=15731065195101753771036961910例3.梯度函數gradient中心差分FX=gradient(F)[FX,FY]=gradient(F)一元函數：二元函數：Fx：向量F的中心梯度H：F相鄰兩點間的間距。F是二維矩陣，FX相當于dF/dx=dF/HXFY相當于dF/dy=dF/HYHX，HY各方向相鄰點距離V=vander(1:5);[FX,FY]=gradient(V)FX=00000-8-6-3-3/2-1-54-36-12-4-2-192-120-30-15/2-3-500-300-60-12-4FY=1573104013410120286102724981036961910例clearde=5*eps;d=pi/100;x=0:d:2*pi;y=sin(x);yde=(sin(x+de)-sin(x))/de;yd=(sin(x+d)-sin(x))/d;ydd=diff(sin(x))/d;yg=gradient(sin(x))/d;holdonplot(x,y,x,yde,'-b*',x,yd,'yh',x(1:end-1),ydd,'r>',x,yg,'+k')已知y=sin(x),采用diff和gradient計算該函數在[0,2π]區間中的近似導函數。eps是一個極小的數:2.2204e-16，為了避免0帶來運算錯誤用eps代替，用機器零閾值eps替代理論0計算極限例f=inline('sqrt(x.^3+2*x.^2-x+12)+(x+5).^(1/6)+5*x+2');x=-3:0.01:3;p=polyfit(x,f(x),5);%用5次多項式p擬合f(x)dp=polyder(p);%對擬合多項式p求導數dpdpx=polyval(dp,x);%求dp在假設點的函數值f=inline('sqrt(x.^3+2*x.^2-x+12)+(x+5).^(1/6)+5*x+2');dx=diff(f([x,3.01]))/0.01;%直接對f(x)求數值導數g=inline('(3*x.^2+4*x-1)./sqrt(x.^3+2*x.^2-x+12)/2+1/6./(x+5).^(5/6)+5');%

f(x)導數解析式gx=g(x);%求函數f的導函數g在假設點的導數plot(x,dpx,x,dx,'.',x,gx,'-');%作圖

用不同方法求f(x)的導數，并在同一坐標中做f'(x)的圖像。例4.7概率統計4.7.1二項分布的概率計算

4.7.2正態分布的概率計算

4.7.3專用概率函數列表

4.7.4數字特征

4.7.5交互界面統計

隨機現象統計規律的數學方法1.二項分布

N次隨機試驗的隨機現象滿足條件：1)重復N次相互獨立隨機試驗；2)每次試驗兩結果成功概率p失敗的概率1-p。4.7.1二項分布bino的概率計算離散隨機變量：重復N次試驗取得p成功的次數為k，k可以在N+1個數即[0,1,…,n]范圍取值的離散隨機變量bino