在 Octave / Matlab 中矢量化 ODE

Question

如果我有一首颂歌并以两种方式编写，比如这里：

function re=rabdab()
x=linspace(0,2000,2000)';
tic;
[T,Y] = ode45(@fun,[x],[0 1 1]);
[T,Y] = ode45(@fun,[x],[0 1 1]);
[T,Y] = ode45(@fun,[x],[0 1 1]);
toc;

tic;
[A,B] = ode45(@fun2,[x],[0 1 1]);
[A,B] = ode45(@fun2,[x],[0 1 1]);
[A,B] = ode45(@fun2,[x],[0 1 1]);
toc;



function dy = fun(t,y)
dy = zeros(3,1);    % a column vector
dy = [y(2) * y(3);...
-y(1) * y(3);...
-0.51 * y(1) * y(2);];

function dy = fun2(t,y)
dy = zeros(3,1);    % a column vector
dy(1) = y(2) * y(3);
dy(2) = -y(1) * y(3);
dy(3) = -0.51 * y(1) * y(2);

几乎没有时间上的差异。一个和另一个一样长。但我认为fun 是fun2 的矢量化版本。还是我在这里弄错了？ 目的是加快我的代码速度。该示例取自 matlab 网页。 我想我还没有真正理解“矢量化”是什么意思。 如果这已经是矢量化的，那么非矢量化的代码会是什么样子？

Answer 1

Vectorization 是一个与functional programming 密切相关的概念。在 MATLAB 中，这意味着在不隐式编写循环的情况下对数组（向量或矩阵）执行操作。

例如，如果您要为 1 到 100 之间的每个整数 x 计算函数 f(x) = 2x，您可以这样写：

for x = 1:100
    f(x) = 2 * x;
end

这不是矢量化代码。矢量化版本是：

x = 1:100; %// Declare a vector of integer values from 1 to 100
f = 2 * x; %// Vectorized operation "*"

甚至更短：

f = 2 * (1:100);

MATLAB 是一种解释型语言，因此很明显，解释器会将其翻译成某种“在幕后”的循环，但它已经过优化并且通常比实际解释循环要快得多（请阅读 this question 以供参考）。嗯，有点 - 直到最近发布的 MATLAB 版本之前都是这样，其中集成了 JIT acceleration（阅读 here）。

现在回到您的代码：这里有两个矢量化版本的代码：一个垂直连接三个值，另一个直接将这些值分配到列向量中。这基本上是一样的。如果您使用显式的 for 循环来完成它，这将不会被“矢量化”。关于“矢量化”循环（即，将 for 循环转换为矢量化操作）的实际性能增益，这取决于首先由于 JIT 加速导致 for 循环的实际速度。

似乎没有太多工作可以加快您的代码速度。您的函数非常基本，因此归结为 ode45 的内部实现，您无法修改。

如果您有兴趣进一步阅读有关矢量化和编写更快的一般 MATLAB 代码的信息，请阅读一篇有趣的文章：Loren on the Art of MATLAB: "Speeding Up MATLAB Applications"。

编码愉快！

Answer 2

虽然前面的答案在一般意义上是正确的，但 vectorized 在 ODE 的上下文中意味着更具体的东西。简而言之，如果 f(t,[y1 y2 ...]) 返回 [f(t,y1) f(t,y2) ...]，对于 y1,y2 列向量，函数 f(t,y) 是向量化的。根据文档 [1]，“这允许求解器减少计算雅可比矩阵的所有列所需的函数评估次数。”

下面的函数 fun3 和 fun4 在 ODE 意义上已正确矢量化：

function re=rabdab()
x=linspace(0,20000,20000)';
opts=odeset('Vectorized','on');

tic;
[T,Y] = ode45(@fun,[x],[0 1 1]);
[T,Y] = ode45(@fun,[x],[0 1 1]);
[T,Y] = ode45(@fun,[x],[0 1 1]);
toc;

tic;
[A,B] = ode45(@fun2,[x],[0 1 1]);
[A,B] = ode45(@fun2,[x],[0 1 1]);
[A,B] = ode45(@fun2,[x],[0 1 1]);
toc;

tic;
[A,B] = ode45(@fun3,[x],[0 1 1],opts);
[A,B] = ode45(@fun3,[x],[0 1 1],opts);
[A,B] = ode45(@fun3,[x],[0 1 1],opts);
toc;

tic;
[A,B] = ode45(@fun4,[x],[0 1 1],opts);
[A,B] = ode45(@fun4,[x],[0 1 1],opts);
[A,B] = ode45(@fun4,[x],[0 1 1],opts);
toc;


function dy = fun(t,y)
dy = zeros(3,1);    % a column vector
dy = [y(2) * y(3);...
-y(1) * y(3);...
-0.51 * y(1) * y(2);];

function dy = fun2(t,y)
dy = zeros(3,1);    % a column vector
dy(1) = y(2) * y(3);
dy(2) = -y(1) * y(3);
dy(3) = -0.51 * y(1) * y(2);

function dy = fun3(t,y)
dy = zeros(size(y)); % a matrix with arbitrarily many columns, rather than a column vector
dy = [y(2,:) .* y(3,:);...
-y(1,:) .* y(3,:);...
-0.51 .* y(1,:) .* y(2,:);];

function dy = fun4(t,y)
dy = [y(2,:) .* y(3,:);... % same as fun3()
-y(1,:) .* y(3,:);...
-0.51 .* y(1,:) .* y(2,:);];

（附带说明：用zeros 省略不必要的内存分配让fun4 运行速度比fun3 稍快。）

• 问：关于第一个参数的向量化怎么样？
• 答：对于 ODE 求解器，ODE 函数仅针对第二个参数进行矢量化。然而，边界值问题求解器bvp4c 确实需要对第一个和第二个参数进行矢量化。 [1]

官方文档 [1] 提供了有关 ODE 特定向量化的更多详细信息（请参阅“雅可比属性描述”部分）。

[1]https://www.mathworks.com/help/releases/R2015b/matlab/ref/odeset.html