第4章數(shù)值計算

1、第四章 數(shù)值計算MATLAB提供高性能的數(shù)學運算功能，特別是在數(shù)值分析方面，例如數(shù)值微分、數(shù)值積分、插值、求極值、方程求根、付立葉變換、常微分方程數(shù)直解、留數(shù)計算等，在求解各種工程問題時得到了廣泛應用。4.1 多項式運算4.1.1多項式的表示法多項式以行向量表示，其中元素為按降冪排列的多項式的系數(shù)?！纠?4-1】將多項式以向量表示。解：將多項式看成一元三次多項式，表示為：1 2 1。 如果將多項式看成一元四次多項式，應表示為：0 1 2 1。4.1.2多項式運算多項式的加減運算通過與多項式對應的行向量之間的加減運算實現(xiàn)，要求參與運算的行向量規(guī)模相同。多項式的乘除運算通過conv和deconv實

2、現(xiàn)。conv（u,v）函數(shù)：求兩向量卷積，此例中用以求得兩多項式之積。q,r=deconv（u,v）函數(shù)：求兩多項式相除的結果，u為被除多項式?！纠?-2】求的“商”及“余”。解：編寫matlab程序：fz=conv(1,0,2,2,4)+0,3,0,1;%計算分子多項式,注意缺項補零fm=1 0 2 1;q,r=deconv(fz,fm);quotient='商多項式為 ' yu='余多項式為 'disp(quotient,poly2str(q,'x'),disp(yu,poly2str(r,'x')運算結果：商多項式為 2余多

3、項式為 7 x2 + 74.1.3多項式求根roots函數(shù)求多項式為0的根。表達形式：r= roots(p) p是多項式系數(shù)形成的行向量，多項式由一個行向量表示，它的系數(shù)按降序排列。r為多項式為0的根，是一個列向量。【例4-3】，計算其根。解：在MATLAB命令窗輸入：p = 1 9 23 15r = roots (p)得結果：r =-1.0000-3.0000-5.0000poly函數(shù)：求已知根所對應的多項式。【例4-4】 已知根-1、-3、-5，求多項式系數(shù)。解：在MATLAB命令窗輸入：S = poly ( -1 -3 -5 )S =1.0000 9.0000 23.0000 15.00

4、00 % 返回行向量是多項式系數(shù)向量。FX=poly2str(S,'x') % 以習慣的方式顯示多項式4.1.4 求多項式的值polyval函數(shù)：求多項式值。表達形式：polyval(p,x)。p是多項式系數(shù)行向量，返回多項式在x的值?！纠?-5】 求在 x = 2 時的值。解：在MATLAB命令窗輸入：p = 1 -3 -4 15;polyval(p, 2)ans =3polyvalm函數(shù)：將矩陣作為變量來計算多項式的值，使用格式：Y=polyvalm(p,X) -p為多項式，X為已知的矩陣，返回也是一個矩陣?！纠?-6】 已知，求。解：編寫MATLAB程序：X=1 2;3

5、4;p=3 4 15;FX=polyvalm(p,X)運算結果：FX = 32 2233 654.1.5 多項式展開Residue留數(shù)函數(shù)，可用于多項式展開，在電路的復頻域分析中，可以用它進行拉氏反變化。假定 （1）可化為： （2） 是真分式，上式可化為：給定B(s)、A(s)，用residue函數(shù)可求出r1，r2, . rn , p1, p2, .pn,和k1, k2 , .kn 。表達形式：r,p,k=residue(num,den)num是B(s)多項式子系數(shù)按降序排列的行向量，den是A(s)多項式子系數(shù)按降序排列的行向量。 命令num,den=residue(r,p,k)將（2）式轉

6、換成（1）式?！纠?-7】 假定，求部分分式展開。解：在MATLAB輸入如下命令：num=4 3 6 10 20;den=1 2 5 2 8;r,p,k=residue(num,den)得結果：r = -1.6970 + 3.0171i -1.6970 - 3.0171i -0.8030 - 0.9906i -0.8030 + 0.9906ip = -1.2629 + 1.7284i -1.2629 - 1.7284i 0.2629 + 1.2949i 0.2629 - 1.2949ik =4即：4.1.6 多項式求導polyder(A)函數(shù)用于求多項式的導數(shù)。【例4-8】求多項式的導數(shù)。解：

7、編寫MATLAB程序：p=3 4 15;q=polyder(p)運算結果：q =6 44.2 數(shù)值分析4.2.1 多項式擬合在實際工程中，經(jīng)常要采集大量數(shù)據(jù)，并用多項式去逼近它，這是一項非常繁重的工作。MATLAB函數(shù)polyfit是一個多項式擬合函數(shù)，可以指定多項式的階數(shù)，使用格式是：p=polyfit(x,y,n)【例4-9】已知一個端口封閉網(wǎng)絡的近似電路模型如下圖，給電路施加電壓，測出端口電壓和通過端口的電流，得出以下表格數(shù)據(jù)，求參數(shù)R和。I12345U44.5688.5表4-1圖4-1解：從圖中看出存在以下關系：，是自變量的一元一次函數(shù)。編寫matlab程序：I=1 2 3 4 5;U

8、=4 4.5 6 8 8.5;z=polyfit(I,U,1)解得：z =1.2500 2.4500寫出關系式：R=1.25，=2.45V。4.2.2 一元插值插值有一元插值、二元插值、多元插值等多種方法。一元插值方法有：線性插值（linear）、三次多項式插值（cubic）、最近點插值（nearest）、三次樣條插值（spline）。一元插值函數(shù)及調用方法：yi=interp1(x,y,xi,method)x,y是已知數(shù)據(jù)，以向量表示，計算xi處的函數(shù)值yi，method為插值方法，有：linear：線性插值；cubic: 三次多項式插值；nearest：最近點插值；spline：三次樣條插

9、值?！纠?-10】將某電壓波形見圖4-2，經(jīng)A/D轉換后采集進計算機，得下表：圖4-2 某電壓波形編號1234567t(ms)00.51.11.52.12.53.1Uo(V)01.43832.52442.99252.72791.79540.42336編號8910111213t(ms)3.544.455.56.2Uo(V)-1.0523-2.2704-2.9326-2.8768-2.1166-0.83825表4-2 試估算當t=1.7ms時的Ux值。解：編寫MATLAB程序，分別用線性插值、三次多項式插值、最近點插值、三次樣條插值求Ux：t=00.51.11.52.12.53.13.544.45

10、5.56.2;Uo=01.43832.52442.99252.72791.79540.42336-1.0523-2.2704-2.9326-2.8768-2.1166-0.83825;tx=1.7;disp(使用線性插值);Ux_l=interp1(t,Uo,tx,linear)disp(使用三次多項式插值);Ux_c=interp1(t,Uo,tx,cubicr)disp(使用最近點插值插值);Ux_c=interp1(t,Uo,tx,nearest)disp(使用三次樣條插值=);Ux_c=interp1(t,Uo,tx,spline)運算結果：使用線性插值，Ux_l = 2.9043使用

11、三次多項式插值Ux_c = 2.9584使用最近點插值插值Ux_c = 2.9925使用三次樣條插值Ux_c =3.0682除了一元插值函數(shù)外，還有二元插值函數(shù)griddata和interp2，三元插值函數(shù)interp3，多元插值函數(shù)interpn，基于FFT方法的插值函數(shù)interpft，可以自行查詢幫助練習。4.2.3 數(shù)值積分利用trapz、quad和quad8函數(shù)，可以求解定積分問題：1. trapz函數(shù)使用梯形公式積分，其調用格式如下：Z=trapz(Y)：Y為被積函數(shù)值，點距為1；Z=trapz(X,Y)：被積函數(shù)由X、Y向量確定；2. quad函數(shù)，采用辛普生法。Z=quad(f

12、un,a,b)：計算函數(shù)fun在a,b區(qū)間的定積分，函數(shù)fun必須返回一個在x上的函數(shù)值向量；Z=quad(fun,a,b,tol)：與上相同，tol指定了容差；【例4-11】求積分解：先定義函數(shù)fx：%fx.mfunction y=fx(x)y=2.*x.*sin(2.*x+30/360*2*pi).*cos(3.*x)再編寫主程序：x=0:0.1:4;plot(x,fx(x);Z=quad(fx,0,4)運算結果：Z= -3.84713. quad8函數(shù)：采用牛頓-科特斯原理計算定積分，調用格式和quad方式相同。4.2.5 數(shù)值微分函數(shù)diff用以計算數(shù)組中的差分，調用格式為：D=dif

13、f(x)設函數(shù)y=f(x)，可近似表示為： (h>0)因此，可以用diff(y)./diff(x)近似計算函數(shù)的微分?！纠?-12】已知RLC串聯(lián)電路電流mA，L=20mH，求電感上電壓，并畫出電壓、電流波形。解：因，可編寫MATLAB程序：t=1e-9:1e-6:0.600001e-3;io=50e-3*exp(-8000*t).*(cos(6000*t)+2*sin(6000*t); %注意，數(shù)組元素相乘要使用.*。UL=20e-6*diff(io)./diff(t); %diff：相鄰元素的差，diff(io)./diff(t)：數(shù)值微分UL=UL,0; %數(shù)值微分后VL數(shù)組維數(shù)減

14、1,為了保證與t維數(shù)一致,須補充1位plot(t,io,t,UL),text(0.31,0.09,'電感電壓','sc'),text(0.51,0.09,'電感電壓','sc'),運行結果：圖4-3 電感電壓、電流波形4.2.6 常微分方程的數(shù)值解MATLAB可用于求解常微分方程（ordinary differential equations，簡稱ODE）的數(shù)值解，求解器Solver為ode45、ode23、ode113、ode15s、ode23s、ode23t、ode23tb。任何高階常微分方程都可以轉化為一階常微分方程，簡化求

15、解過程。求解具體ODE的基本過程：1. 工程問題以微分方程與初始條件描述為：2. 運用變量代換：，把高階（大于2階）的方程（組）寫成一階微分方程組：寫成矩陣形式：，3. 編寫計算導數(shù)的M-函數(shù)文件，將初始條件傳遞給求解器Solver中的一個，運行后就可得到ODE的在指定區(qū)間上的解列向量y（其中包含y及不同階的導數(shù)）。ODE函數(shù)調用格式 T,Y = solver(odefun,tspan,y0) %在區(qū)間tspan=t0,tf上，從t0到tf，用初始條件y0求解顯式微分方程y=f(t,y)。對于標量t與列向量y，函數(shù)f=odefun(t,y)必須返回一f(t,y)的列向量f。解矩陣Y中的每一行對

16、應于返回的時間列向量T中的一個時間點。要獲得問題在其他指定時間點t0,t1,t2,上的解，則令tspan=t0,t1,t2,tf（要求是單調的）。T,Y = solver(odefun,tspan,y0,options) %用參數(shù)options（用命令odeset生成）設置的屬性（代替了缺省的積分參數(shù)），再進行操作。常用的屬性包括相對誤差值RelTol（缺省值為1e-3）與絕對誤差向量AbsTol（缺省值為每一元素為1e-6）。T,Y =solver(odefun,tspan,y0,options,p1,p2) 將參數(shù)p1,p2,p3,.等傳遞給函數(shù)odefun，再進行計算。若沒有參數(shù)設置，則

17、令options=。因為沒有一種算法可以有效地解決所有的ODE問題，為此，MATLAB提供了多種求解器Solver，對于不同的ODE問題，采用不同的Solver。求解器Solver特點說明ode45一步算法；4，5階Runge-Kutta方程；累計截斷誤差達(x)3大部分場合的首選算法ode23一步算法；2，3階Runge-Kutta方程；累計截斷誤差達(x)3使用于精度較低的情形ode113多步法；Adams算法；高低精度均可到10-310-6計算時間比ode45短ode23t采用梯形算法適度剛性情形ode15s多步法；Gears反向數(shù)值微分；精度中等若ode45失效時，可嘗試使用ode23

18、s一步法；2階Rosebrock算法；低精度當精度較低時，計算時間比ode15s短ode23tb梯形算法；低精度當精度較低時，計算時間比ode15s短表4-3 不同求解器Solver的特點4在計算過程中，用戶可以對求解指令solver中的具體執(zhí)行參數(shù)進行設置（如絕對誤差、相對誤差、步長等）。屬性名取值含義AbsTol有效值：正實數(shù)或向量缺省值：1e-6絕對誤差對應于解向量中的所有元素；向量則分別對應于解向量中的每一分量RelTol有效值：正實數(shù)缺省值：1e-3相對誤差對應于解向量中的所有元素。在每步(第k步)計算過程中，誤差估計為：e(k)<=max(RelTol*abs(y(k),Ab

19、sTol(k)NormControl有效值：on、off缺省值：off為on時，控制解向量范數(shù)的相對誤差，使每步計算中，滿足：norm(e)<=max(RelTol*norm(y),AbsTol)Events有效值：on、off為on時，返回相應的事件記錄OutputFcn有效值：odeplot、odephas2、odephas3、odeprint缺省值：odeplot若無輸出參量，則solver將執(zhí)行下面操作之一：畫出解向量中各元素隨時間的變化；畫出解向量中前兩個分量構成的相平面圖；畫出解向量中前三個分量構成的三維相空間圖；隨計算過程，顯示解向量OutputSel有效值：正整數(shù)向量缺省

20、值：若不使用缺省設置，則OutputFcn所表現(xiàn)的是那些正整數(shù)指定的解向量中的分量的曲線或數(shù)據(jù)。若為缺省值時，則缺省地按上面情形進行操作Refine有效值：正整數(shù)k>1缺省值：k = 1若k>1，則增加每個積分步中的數(shù)據(jù)點記錄，使解曲線更加的光滑Jacobian有效值：on、off缺省值：off若為on時，返回相應的ode函數(shù)的Jacobi矩陣Jpattern有效值：on、off缺省值：off為on時，返回相應的ode函數(shù)的稀疏Jacobi矩陣Mass有效值：none、M、M(t)、M(t,y)缺省值：noneM：不隨時間變化的常數(shù)矩陣M(t)：隨時間變化的矩陣M(t,y)：隨時間

21、、地點變化的矩陣MaxStep有效值：正實數(shù)缺省值：tspans/10最大積分步長表4-4 Solver中options的屬性【例4-13】求解范德坡方程（Van Der Pol Equation）：解：令初值：，1 先寫出函數(shù)文件function dydt=vdp1(t,y)dydt=y(2);(1-y(1)2)*y(2)-y(1);2寫出matlab主程序：t,y = ode45('vdp1',0 20,2; 0); plot(t,y(:,1),'-',t,y(:,2),'-')title('Solution of van der P

22、ol Equation, mu = 1');xlabel('time t');ylabel('solution y');legend('y_1','y_2')圖4-4 4.2.7 求最小值求函數(shù)最小值可通過求導找導數(shù)為0的點求出，但是很多函數(shù)的導數(shù)很難求出，必須使用數(shù)值分析法來求解。MATLAB提供了求局部最小值函數(shù)：fminbnd、fminsearch函數(shù)（在早期版本中為fmin、fmins），可以求一維甚至n維函數(shù)的局部最小值。fminbnd：單變量非線性函數(shù)局部最小值函數(shù)，調用格式：x=fminbnd(FUN,x1,

23、x2) %返回值x是FUN函數(shù)在(x1,x2)上的局部最小值，F(xiàn)UN為單值非線性函數(shù)。x=fminbnd(FUN,x1,x2,options)x=fminbnd(FUN,x1,x2,options,p1,p2,)options用來控制算法的參數(shù)，有：Display：如果為非0，顯示結果的中間步數(shù)，缺省值為0；TolX：x的容許誤差，缺省值為1e-4(10-4)；MaxFunEval： 結果fun(x)的容許誤差，缺省值為1e-4(10-4)；MaxIter：疊代最大步數(shù)，缺省值為500?！纠?-14】求函數(shù)在區(qū)間(2.5,8)之間的最小值。解：先定義函數(shù)fminfun，以文件fminfun.m

24、保存：function ff=fminfunff='2*exp(-x)*(sin(x)2+0.2)'編寫程序求最小值：xmin=fminbnd(fminfun,2.5,8)x=xmin;ymin=eval(fminfun) %求在最小值點的函數(shù)值運行結果xmin = 6.3896ymin = 7.0945e-004fminsearch為多變量無束縛非線性局部最小值函數(shù)，調用格式：xmin=fminsearch(Fun,x0,options,p1,p2,)fminsearch采用Nelder-Mead單純形法求解多變量函數(shù)局部最小值，x0是最小值的初始猜測值，其它參數(shù)與fminbnd一致。fminbnd、fminsearch函數(shù)也可用于求解函數(shù)局部最大值。在求局部最大值時，定義一個函數(shù)為原函數(shù)的負函數(shù)，用fminbnd或fminsearch求出局部最小值后，即可求出原函數(shù)的最大值。4.2.8 求零點 MATL