MATLAB求函式零點—fzero函式
5.6 函式的零點
5.6.2 一元函式的零點
5.6.2.2 任意一元函式零點的精確解
【 * 例 5.6.2 .2-1 】通過求 的零點,綜合敘述相關指令的用法。
(1)構造一個行內函數物件
被解函式 以 為自變數, 和 為數。假如在
fzero 中直接採用字串表示被解函式,容易出錯。因此先構造行內函數如下:
y=inline('sin(t)^2*exp(-a*t)-b*abs(t)','t','a','b'); %<1>
(2)作圖法觀察函式零點分佈
a=0.1;b=0.5;t=-10:0.01:10; % 對自變數取樣,取樣步不宜太大。
y_char=vectorize(y); % 為避免迴圈,把 y 改寫成適合陣列運算形式。 <4>
Y=feval(y_char,t,a,b); % 在取樣點上計算函式值。
clf,plot(t,Y,'r');hold on,plot(t,zeros(size(t)),'k'); % 畫座標橫軸
xlabel('t');ylabel('y(t)'),hold off
圖 5.6.2 .2-1 函式零點分佈觀察圖
(3)利用 zoom 和 ginput 指令獲得零點的初始近似值(在 MATLAB 指令窗中進行)
zoom on % 在 MATLAB 指令窗中執行,獲區域性放大圖
[tt,yy]=ginput(5);zoom off % 在 MATLAB 指令窗中執行,用滑鼠獲 5 個零點猜測值。
圖 5.6.2 .2-2 區域性放大和利用滑鼠取值圖
(6)觀察 fzero 所採用的 options 預設設定,並更改控制計算精度的相對誤差設定。
op=optimset('fzero') % 提取 fzero 所採用的 options 預設設定
op =
......
Display: 'final'
......
TolX: 2.2204e-016
TypicalX: []
op=optimset('tolx',0.01); % 把終止計算的相對誤差閾值設定得較大
op.TolX % 觀察新設定值。注意 TolX 字母的大小寫。
ans =
0.0100
(7)利用新設定的選項引數重新求 tt(4)附近的零點,以便比較。
[t4n,y4n,exitflag]=fzero(y,tt(4),op,a,b) % 採用新的 op 設定引數。
Zero found in the interval: [0.57094, 0.60418].
t4n =
0 。 6042
y4n =
0 。 0017
exitflag =
1
〖說明〗
1、本例是採用行內函數形式求取函式零點的。
2、 若採用如下 M 函式檔案(該檔案必須放在搜尋路徑上)
[y_M.m]
function y=y_M(t,a,b)
y=sin(t).^2.*exp(-a*t)-b*abs(t);
則相應的求零點指令是 [t 4m ,y 4m ,exitflag]=fzero('y_M',tt(4),[],a,b)• 若直接用字串表達函式,則應把被解函式自變數 t 改成 x ,引數 a 、 b 改成 P1 、 P2 。相應的具體指令如下
P1=0.1;P2=0.5;
y_C='sin(x).^2.*exp(-P1*x)-P2*abs(x)';
[t 4c ,y 4c ,exitflag]=fzero(y_C,tt(4),[],P1,P2)
5.6.3 多元函式的零點
【例 5.6.3 -1 】求解二元函式方程組
(1)從三維座標初步觀察兩函式圖形相交情況
x=-2:0.05:2;y=x;[X,Y]=meshgrid(x,y); % 產生 x-y 平面上網點座標
F1=sin(X-Y);F2=cos(X+Y);
F0=zeros(size(X));
surf(X,Y,F1),
xlabel('x'),ylabel('y'),
view([-31,62]),hold on,
surf(X,Y,F2),surf(X,Y,F0),
shading interp,
hold off
圖 5.6.3 -0 兩函式的三維相交圖
(2)在某區域觀察兩函式 0 等位線的交點情況
clear;
x=-2:0.5:2;y=x;[X,Y]=meshgrid(x,y); % 產生 x-y 平面上網點座標
F1=sin(X-Y);F2=cos(X+Y);
v=[-0.2, 0, 0.2]; % 指定三個等位值,是為了更可靠地判斷 0 等位線的存在。
contour(X,Y,F1,v) % 畫 F1 的三條等位線。
hold on,contour(X,Y,F2,v),hold off % 畫 F2 的三條等位線。
圖 5.6.3 -1 兩個二元函式 0 等位線的交點圖
(3)從圖形獲取零點的初始近似值
在圖 5.6.3 -1 中,用 ginput 獲取兩個函式 0 等位線(即三線組中間那條線)交點的座標。
[x0,y0]=ginput(2); % 在圖上取兩個點的座標
disp([x0,y0])
-0.7926 -0.7843
0.7926 0.7843
(4)利用 fsolve 求精確解。以求( 0.7926,7843 )附近的解為例。
本例直接用字串表達被解函式。注意:在此,自變數必須寫成 x(1), x(2) 。假如寫成 xy(1), xy(2) ,指令執行將出錯。
fun='[sin(x(1)-x(2)),cos(x(1)+x(2))]'; %<12>
xy=fsolve(fun,[x0(2),y0(2)]) %<13>
xy =
0.7854 0.7854
(5)檢驗
fxy1=sin(xy(1)-xy(2));fxy2=cos(xy(1)+xy(2));disp([fxy1,fxy2])
1.0e-006 *
-0.0994 0.2019
〖說明〗指令 <12><13> 可用以下任何一組指令取代。
(A)行內函數形式指令
fun=inline('[sin(x(1)-x(2)), cos(x(1)+x(2))]', 'x'); % 項 'x' 必須有。
xy=fsolve(fun,[x0(2), y0(2)]);
(B) M 函式檔案形式及指令
先用如下 fun.m 表示被解函式(並在搜尋路徑上)
[fun.m]
function ff=fun(x)
ff(1)=sin(x(1)-x(2));
ff(2)=cos(x(1)+x(2));
然後執行指令 xy=fsolve('fun',[x0(2),y0(2)]) 。
第四步檢驗中的結果表明:所找零點處的函式值小於 ,是一個十分接近零的小數。該精度由 options.TolFun 控制。 options.TolFun 的預設值是
1.0000e-006 。它可以用下列指令看到
options=optimset('fsolve');
options.TolFun
ans =
1.0000e-006