3

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

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

但我不太明白我需要改变什么ODEfun以使其符合 Matlab 想要的'vectorized' ODEfun外观。

从文档中我看到您可以将示例范德波尔系统从其正常形式更改：

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应该代表什么，以及如何定义第二维的大小。我几乎可以只添加“ ,:”而不考虑它，但我遇到了问题，因为我已经在我的代码中进行了一些向量操作。

这是我当前功能的简化示例，还没有vectorized。它模拟 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

显然，在真正的函数中还有更多的事情发生，无论是 forT还是X. 如您所见，dXdt(1)是一个边界条件，需要自己计算。

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

我正在使用 Matlab R2006a。

4

1 回答 1

2

From 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) ...].

For the van der Pol example:

without vectorization: 

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

with vectorization: 

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´ ...]

where the y1, y2, y3, etc. refer to different vectors y at the same time t that ode15s will use to compute the next step.

For your example, you have to take into account that the y you get passed is no longer a vector, but a matrix in which every column represents a different vector you need to compute the derivative of: 

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
于 2013-09-06T12:28:17.907 回答