使用 ODE45 积分的 if 和 for 循环的 Matlab 矢量化

Question

我对优化我的代码速度很感兴趣，在使用 "Run and Time" 时，这是我代码中对速度有很大影响的函数，但我很难概念化如何正确矢量化此函数，因为我通常只是循环，在尝试中我运行进入错误，因为它也用于集成，我的原始函数如下并且不会导致错误

function [dotStates] = ODEFunc(t,states,params)
%ODE function
% Loading in and assigning the variables from parameters
K = params(1);                                                             
N = params(2);
nn = params(3);
% magnitude of the coupling based on the number of neighbours
kn = K/nn;
w = params(4:end);
dotStates=states;
    % For each oscillator
    for i=1:N
        % Use the oscillators natural frequency
        dotStates(i) = w(i);
        % For j number of neighbours
        for j=(i-nn):(i+nn)
            % neighbour number is positive and shorter than # of oscilators
            if (j > 0) && (j < length(dotStates))
                dotStates(i) =  dotStates(i) + (kn * sin( states(j)-states(i) ));
            end
        end
    end
end

我已尝试遵循 mathworks 矢量化指南：https://se.mathworks.com/help/matlab/matlab_prog/vectorization.html

到目前为止，我的尝试一直是遵循他们使用的一些输入，例如使用掩码并生成了以下代码

function [dotStates] = ODEFunc(t,states,params)
%ODE function
% Loading in and assigning the variables from parameters
K = params(1);                                                             
N = params(2);
nn = params(3);
% magnitude of the coupling based on the number of neighbours
kn = K/nn;
w = params(4:end);
dotStates=states;
% Use the oscillators natural frequency
dotStates = w';
% Mask of j states
j = (i-nn):(i+nn);
% neighbours cannot exceed boundaries
j = j(j>0 & j <=length(dotStates));
jstate = states(j);
jstate(numel(states)) = 0;
dotStates =  dotStates + (kn * sin( jstate'-states ));
end

我最终得到了一个比写入的向量更短的向量，我的解决方案是在 "jstate" 变量中添加一串零来弥补差异，但这并没有感觉像是一个适当的矢量化，当我运行代码时，我得到以下与集成步骤相关的错误

Warning: Colon operands must be real scalars.

In RK_ODE_2411>ODEFunc (line 99) In RK_ODE_2411>@(t,states)ODEFunc(t,states,params) In ode45 (line 324) In RK_ODE_2411 (line 58)

该函数依次用于以下部分，用于使用 ODE45 进行积分

%% Integration via ODE45
for K = 0:.1:Klen
    params(1) = K;
    K_count = K_count+1;
    nn_count = 0;
    for nn = nnlen:nnlen
        params(3) = nn;
        % index counter
        nn_count = nn_count+1;
        % 6th order runge kutta

        sol(K_count,nn_count) = ode45(@(t,states) ODEFunc(t,states,params),tSpan,init,options);
    end
end

第 58 行是

sol(K_count,nn_count) = ode45(@(t,states) ODEFunc(t,states,params),tSpan,init,options);

编辑：ODEFunc 中的第 99 行是

j = (i-nn):(i+nn);

Answer 1

试试这个片段

% For each oscillator
for i=1:N
    % For j number of neighbours
    j=(i-nn):(i+nn);

    % neighbour number is positive and shorter than # of oscilators
    lg = (j > 0) & (j < length(dotStates));
    dotStates(i) =  w(i) + sum(kn * sin( states(lg)-states(i) ));
end

最重要的是 dotStates 不会大于 stats，因为这会迫使 matlab 重新排列其内存，从而极大地减慢代码速度。

使用 ODE45 积分的 if 和 for 循环的 Matlab 矢量化

Matlab vectorization with if and for loops using ODE45 to integrate

matlab

if-statement

vectorization

nested-loops

ode45