使用 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， 您也可下载压缩文件。