Matlab:为 odeset 的 'Vectorized' 属性调整 myODEfun

Question

我有一个耦合 ODE 的刚性系统，我正在将其提供给 MATLAB 的 ode15s 求解器。它运行良好，但现在我正在尝试优化集成速度。我在 N 不同的空间站点上对 5 个不同的变量进行建模，给出 5N 个耦合方程。目前，N=20 积分时间约为 25 秒，但我想使用更大的 N 值。

我使用探查器发现绝大多数时间都花在了评估 myODEfun 上。我尽了最大努力优化代码，但这并没有改变函数中发生了相当多的事情并且它正在被评估约 50,000 次的事实。我读到为 ODEfunction 使用 'Vectorized' 属性可以减少所需的计算次数。

但我不太明白我的 ODEfun 到底需要更改什么以使其符合 Matlab 想要的 'vectorized' ODEfun 看起来像。

From the documentation我看到您可以将示例范德波尔系统从其正常形式更改为:

function dydt = vdp1000(t,y)
dydt = [y(2); 1000*(1-y(1)^2)*y(2)-y(1)];

到矢量化形式:

function dydt = vdp1000(t,y)
dydt = [y(2,:); 1000*(1-y(1,:).^2).*y(2,:)-y(1,:)];

我不完全理解这个新的 y 矩阵应该表示什么，以及如何定义第二个维度的大小。我几乎可以忍受只添加“,:”而不考虑它，但我遇到了问题，因为我已经在我的代码中进行了一些向量操作。

这是我当前函数的一个简化示例，尚未向量化。它对 2 个变量建模，生成 2*N 方程。请不要试图理解此处生成的 ODE:它们没有意义。我说的是正在发生的操作。

function dydt = exampleODEfun(t,y,N)

dydt = zeros(2*N,1);
dTdt = zeros(N,1);
dXdt = zeros(N,1);

T = y(1:N);
X = y(N+1:2*N);

a = [T(2:N).^2 T(2:N) ones(N-1,1)];
b = [3 5 -2];

dTdt(1:N) = 0;
dXdt(1) = 0;
dXdt(2:N) = a*b';

dydt(1:N) = dTdt;
dydt(N+1:2*N) = dXdt;

end

很明显，在实际函数中，T 和 X 都发生了很多事情。如您所见，dXdt(1) 是一个边界条件，需要自己计算。

盲目传递 odeset 'Vectorized','on' 并仅向所有索引添加“,:”是行不通的。例如，我现在需要将dTdt 和dXdt 初始化到什么大小？我该如何处理 ones(N-1,1)？我需要做什么才能使 (a*b') 仍然有意义？

我正在使用 Matlab R2006a。

Answer 1

来自 help odeset:

Vectorized - Vectorized ODE function [ on | {off} ]

Set this property 'on' if the ODE function F is coded so that 
F(t,[y1 y2 ...]) returns [F(t,y1) F(t,y2) ...].

以 van der Pol 为例:

没有矢量化:

function dydt = vdp1000(t,y)             %// 'y' is passed as [y1 
                                         %//                   y2]

    dydt = [y(2);                        %// 'dydt' is computed as [y1´
            1000*(1-y(1)^2)*y(2)-y(1)]   %//                        y2´]
                                         %// where the ´ indicates d/dt

使用矢量化:

function dydt = vdp1000(t,y)            %// 'y' is passed as [y11 y21 y31 ...
                                        %//                   y12 y22 y32 ...] 

    dydt = [y(2,:);                             %// 'dydt' is computed as 
            1000*(1-y(1,:).^2).*y(2,:)-y(1,:)]; %//   [y11´ y21´ y31´ ...
                                                %//    y12´ y22´ y32´ ...]

其中y1、y2、y3等同时引用不同的向量y t ode15s 将用于计算下一步。

对于您的示例，您必须考虑到您传递的 y 不再是一个向量，而是一个矩阵，其中每一列代表一个不同的向量，您需要计算其导数:

function dydt = exampleODEfun(t,y,N)

    %// Adjust sizes to meet size of y
    dydt = zeros(2*N, size(y,2));
    dTdt = zeros(N, size(y,2));
    dXdt = zeros(N, size(y,2));

    %// Adjust indices to extract proper vales of ALL vectors
    T = y(1:N,:);
    X = y(N+1:2*N,:);

    %// This sort of section is usually where all the "thought" goes into:
    %// you can't use a*b' anymore, so I sum over the third dimension of the 
    % 3D array I built from your original vector
    b  = [3 5 -2];
    ab = sum(cat(3, b(1)*T(2:N,:).^2, b(2)*T(2:N,:), b(3)*ones(N-1, size(y,2))), 3);

    %// and finish it off
    dTdt(1:N,:) = 0;
    dXdt(1,:) = 0;
    dXdt(2:N,:) = ab;

    dydt(1:N,:) = dTdt;
    dydt(N+1:2*N,:) = dXdt;

end