Octave / Matlab 中的 ODE 向量化

Question

如果我有一首Ode并以两种方式写出它，就像这样：

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网页。我想我还没有真正理解"向量化"是什么意思。如果这已经被向量化了，那么一段未被向量化的代码会是什么样子呢？

- lyvic

Answer 1

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

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

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加速（请在here中阅读）。

现在回到您的代码：您的代码中有两个矢量化版本：一个将三个值垂直连接在一起，另一个则直接将这些值分配到列向量中。基本上是同样的事情。如果使用显式的for循环执行，这就不会被“向量化”。关于从“向量化”循环（即将for循环转换为矢量化操作）中获得的实际性能提升，这取决于由于JIT加速而实际上for循环有多快。

看起来没有太多可以做来加速你的代码。你的函数相当基础，所以归结为 ode45 的内部实现，你无法修改。

如果你对向量化和编写更快的 MATLAB 代码感兴趣，这里有一篇有趣的文章：Loren on the Art of MATLAB: "Speeding Up MATLAB Applications"。

编写愉快！

- Eitan T

---

谢谢你帮我解决了这个困惑。所以，使用ode23可能会更快，因为它不太准确，但向量化在这里并不是问题。但是，由于你的解释，我可以找出代码中其他可以向量化的部分。 - lyvic

Answer 2

虽然前面的回答在一般情况下是正确的，但在ODE的上下文中，“向量化”意味着更具体的含义。简而言之，如果函数f(t,y)满足f(t,[y1 y2 ...])返回[f(t,y1) f(t,y2) ...]，其中y1,y2是列向量，则该函数是向量化的。根据文档[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特定向量化的详细信息（请参见“Jacobian属性描述”部分）。

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

- Vicky