使用 NAG MATLAB® 工具箱

首頁 / 技術文章 / 使用 NAG MATLAB^® 工具箱

NAG MATLAB^®工具箱應用範例

使用 NAG MATLAB^® 工具箱

NAG 工具箱能夠讓使用者在 MATLAB 的環境中，使用完整的 NAG 程式庫函數。其中一個很大的特點是，透過 MATLAB 呼叫 NAG 函數，許多函數的參數都可以省略，或者以可選項方式呈現，這樣可以讓程式變的更容易閱讀與維護。

所有 NAG 程式庫的手冊 - 幾近 17 大冊，如今已經全部轉換為 MATLAB 的說明擋了。使用者可以很輕鬆的在 MATLAB 的環境中查找相關函數內容。在每一個函數中，也包含了 MATLAB 呼叫此函數的範例程式。

為了展示如何能夠輕鬆的使用 NAG 工具箱，我們將示範如何呼叫 NAG 函數，以及如何使用 MATLAB 的繪圖函數顯示計算結果。

注意： 本文中的範例皆是由我們的展示程式中截取的，若您直接以文中的 MATLAB 程式執行，也許無法完整的獲得如文中的結果。文中範例的完整的程式，請您自此處下載。

S 章節函數 - Bessel 函數

我們先從最簡單的 NAG 函數開始談起。在其他許多的特殊函數中，S 章節中包含了 Bessel 函數 Y0(x), Y1(x), J0(x) 以及 J1(x)。依據 NAG 函數的命名原則，他們分別對應至 s17ac, s17ad, s17ae 與 s17af 函數。

以下是呼叫此函數的命令，以及畫圖的指令：

% 求解 x 資料點的 Bessel 函數
x = [0.125, 0.25 : 0.25 : 40];
for i = 1 : length(x)
   [y0(i), ifail] = s17ac(x(i));
   [y1(i), ifail] = s17ad(x(i));
   [j0(i), ifail] = s17ae(x(i));
   [j1(i), ifail] = s17af(x(i));
end

% 繪製結果。我們忽略 Y1(x) 前面的兩個值，因為他們的值太大會影響圖形坐標軸的呈現
hold on;
plot(x, y0, 'Color', 'red');
plot(x(3:length(x)), y1(3:length(x)), 'Color', 'green');
plot(x, j0, 'Color', 'blue');
plot(x, j1, 'Color', [1.0 1.0 0.0]);

輸出結果如下圖：

Graph of Bessel Functions Using NAG Toolbox for MATLAB®
圖 1：Bessel 函數

本範例的 MATLAB 腳本位於 NAGToolboxDemos/Special_functions/s17a_demo.m，也可下載壓縮檔。

E02 函數 - 採用 bicubic splines 方法進行曲面擬合

在很多科學領域中，找出一個能夠估計逼近所有資料點的函數是極為普遍的問題。這些資料也許包含許多混亂隨機的資料，例如量測的誤差等，我們便需要將其平滑處理的。

e02dc 函數會在給定 x-y 平面的網格中，透過 bicubic spline 演算法進行資料估計。其中最簡單的參數就是控制擬合的緊密與平滑程度。

在 e02dc 函數中的平滑參數 s，若給定為最小值將會計算出最靠近資料點的擬合結果，但是其確切意義會隨著資料特性的不同而有所改變。理論上，s = 0 的值將會提供內插的樣條，雖然其有可能在原始資料附近會有震盪問題產生。

% 任意產生平面中的資料點，其中包含明顯的極大值
x = [0; 0.5; 1; 1.5; 2; 2.5; 3; 3.5; 4; 4.5; 5];
y = [0; 0.5; 1; 1.5; 2; 2.5; 3; 3.5; 4];
f = [1; 1; 1; 1; 1; 1; 1; 1; 1; 1; 1;
     1; 1; 1.1; 1.2; 1.1; 1; 1; 1; 1; 1; 1;
     1; 1; 1; 1; 1; 1; 1; 1; 1; 1; 1;
      1; 1; 1; 1; 1; 1; 1; 1; 1; 1; 1;
      1; 1; 1; 1; 1; 1.075; 1.15; 1.075; 1; 1; 1;
      1; 1; 1; 1; 1; 1; 1; 1; 1; 1; 1;
      1; 1; 1; 1; 1; 1; 1; 1; 1; 1; 1; 1.1;
      1.2; 1.1; 1; 1; 1; 1; 1; 1; 1; 1;
      1; 1; 1; 1; 1; 1; 1; 1; 1; 1; 1];

% 平滑參數，給定最為適合的數值
s = 0.0175593;

% 設定網格的座標參數
px = 0 : 0.1 : 5.0;
py = 0 : 0.1 : 4.0;

% 初始化 e02dc 函數所需的參數
start = 'C';
nx = int32(0);
lamda = zeros(15,1);
mu = zeros(13,1);
ny = int32(0);
wrk = zeros(592, 1);
iwrk = zeros(51, 1, 'int32');

% 使用 e02dc 函數計算 bicubic spline 的估計參數值
[nxOut, lamdaOut, nyOut, muOut, c, fp, wrkOut, iwrkOut, ifail] = ...
   e02dc(start, x, y, f, s, nx, lamda, ny, mu, wrk, iwrk);

% 使用 e02df 求解格點上的 spline 值
[ff, ifail] = e02df(px, py, lamdaOut(1:nxOut), muOut(1:nyOut), c);

圖 2 與 3 顯示原始資料，以及透過 bicubic spline 平滑化後的估計值。

Graph of data to be fitted by bicubic spline
圖 2：透過 bicubic spline 平滑後的圖型

the bicubic spline evaluated on a rectangular grid

圖 3：格點資料的 bicubic spline 估計值

本範例的 MATLAB 腳本位於 NAGToolboxDemos/Curve_and_surface_fitting/e02dc_demo.m，也可下載壓縮檔。

E04 函數 - 最佳化函數

NAG e04wd 函數能夠找出任意平滑方程的最小值，其中可以加入線性或非線性的限制式。它採用的是序列二次規劃法 (Sequential Quadratic Pogramming, SQP) 演算法。可以參考 e04wd 說明文件了解此方法的詳細說明。

我們用 e04wd 函數嘗試找出 Rosenbrock 方程的最小值：

  f(x,y) = 100*(y-x*x)^2 + (1-x)^2

這個求解的 Rosenbrock 是一個 無限制性的問題 - 表示求解的 x 與 y 數值並沒有任何的限制條件，可在空間中任意搜尋。

由於 Rosenbrock 方程的等高線圖有如同香蕉外形的特質，如下圖所示，這個方程通常被用來測試最佳化函數。最簡單的演算法例如很基本的最速下降法 (Steepest Descent)，在追蹤最小值時，是很難隨著彎曲的低窪處進行追蹤的，但是 NAG 提供的函數不會有這個困難。

e04wd 函數的一個特點是，它不需要使用者提供目標方程的一階導數。如果沒有提供一階導數，或者很難求得其導數，e04wd 可以透過 difference approximations 估計出。然而若使用者能夠提供導數，則演算法的收斂將會更為快速且穩定。

提供給 e04wd 求解最佳化的函數是透過 MATLAB 的 .m 檔案 - 在這裡我們命名為 objfun.m。如果需要的話，這個檔案內也同時計算目標函數的一階導數。然而，若沒有提供導數呢？那就不需要計算它；e04wd 能夠很聰明的我們是否有提供，若未提供，它會自行計算求得。

以下是 objfun.m 的內容：

function [mode, objf, grad, user] = ...
          objfun(mode, n, x, grad, nstate, user)
objf = 100*(x(2)-x(1)^2)^2 + (1 - x(1))^2;
if (mode == 1 || mode == 2)
   % 導數：grad(i) 分別表示 x(i) 的導數，i = 1, 2.
   grad(1) = 2*(1-200*x(2))*x(1) + 400*x(1)*x(1)*x(1) - 2;
   grad(2) = 200*(x(2) - x(1)*x(1));
end

利用 e04wd 求解：

% 沒有線性或非線性限制式
a = [];
istate = zeros(2, 1, 'int32');
ccon = [];
cjac = [];
clamda = zeros(2,1);
hess = zeros(2,2);

% 任意設定變數的邊界
bl = [-10; -10];
bu = [10; 10];

% 初值
x = [-2.75; 2];

% 使用 e04wc 函數初始化 e04wd 
[iw, rw] = e04wc();

% user 變量包含我們想傳遞給 objfun 函數的數值
user = [0,0,0,0,0];

% 使用 e04wd 函數求解最佳化
[majits, istateOut, cconOut, cjacOut, clamdaOut, objf, grad, ...
       hessOut, xOut, iwOut, rwOut, user] = ...
   e04wd(a, bl, bu, 'confun', 'objfun', istate, ...
         ccon, cjac, clamda, hess, x, iw, rw, 'user', user);

請您注意在 objfun.m 函數中，e04wd 求解時也需要一個 confun.m 函數，提供非線性的限制式。在此例中，我們並未提供限制式，所以我們就給他一個空值。

為了簡化起見，所有 MATLAB 繪圖的命令我們都省略。

下圖顯示當提供導數時，求解 Rosenbrock 函數的最小值的路徑。當每一次目標函數評值時，我們就將得出的求解點繪出。初始點我們用藍框表示，最小值 (1,1) 用綠框顯示。您可以發現，函數如何追蹤並回朔及改正最小值求解。

Path to the minimum of Rosenbrock's function when derivatives are known
圖 4：當導數已知時，Rosenbrock 函數求解路徑

為了要顯示 e04wd 如何在沒有提供導數的情況下，還能正確求解，我們修改了 objfun.m 中的程式：

function [mode, objf, grad, user] = ...
          objfun(mode, n, x, grad, nstate, user)
objf = 100*(x(2)-x(1)^2)^2 + (1 - x(1))^2;
% 未提供一階導數資訊
end

Path to the minimum of Rosenbrock's function when derivatives are not supplied
圖 5：未提供導數，Rosenbrock 函數求解路徑

這次我們並未提供 e04wd 函數任何的導數。所以函數必須要使用額外計算求解導數。整體來說，我們需要 231 次求解，相反的，若我們提供導數，則只需要 54 次的求解。此外，獲得的最小值並不如有導數的，雖然其還是能取得正確解 0.0。

這說明了甚麼呢？如果你可以提供導數的話，請您提供；然而若沒有提供導數，則 NAG 函數會幫您計算！

本範例的 MATLAB 腳本位於 NAGToolboxDemos/Minimization/e04wd_demo.m， 您也可下載壓縮檔。