Параллельный пул на MATLAB для бифуркации

Question

Я новичок в этой концепции параллельного объединения в MATLAB (я использую версию 2019 a) и кодирования. Этот код, которым я собираюсь поделиться с вами, был доступен на net с некоторыми изменениями, которые я сделал для моих требований.

Постановка задачи: у меня нелинейная система (уравнение Росслера), и я должен построить ее диаграмму бифуркации, я попытался сделать это обычно, используя для l oop, но его время вычисления было слишком большим и мой компьютер несколько раз зависал, поэтому я получил совет распараллелить мой код, чтобы выйти из этой проблемы. Я попытался узнать, как распараллеливать пул с помощью MATLAB на net, но все же я не смог решить свои проблемы, так как все еще есть некоторые проблемы, так как в моем коде есть 2 цикла parfor, у меня проблемы с индексированием и в присвоение глобального параметра (Обратите внимание: этот код написан для обычного выполнения без использования параллельного пула).

Я прилагаю свой код ниже, пожалуйста, извините, если я упомянул много строк кода ,

clc;
a = 0.2; b = 0.2; global c;
crange = 1:0.05:90; % Range for parameter c
k = 0; tspan = 0:0.1:500; % Time interval for solving Rossler system
xmax = []; % A matrix for storing the sorted value of x1

for c = crange
f = @(t,x) [-x(2)-x(3); x(1)+a*x(2); b+x(3)*(x(1)-c)];
    x0 = [1 1 0]; % initial condition for Rossler system 
    k = k + 1;
    [t,x] = ode45(f,tspan,x0); % call ode() to solve Rossler system
    count = find(t>100); % find all the t values which is >10  
    x = x(count,:);
    j = 1; 
    n = length(x(:,1)); % find the length of vector x1(x in our problem)

    for i=2 : n-1
       % check for the min value in 1st column of sol matrix
       if (x(i-1,1)+eps) < x(i,1) && x(i,1) > (x(i+1,1)+eps)  
           xmax(k,j)=x(i,1); % Sorting the values of x1 in increasing order
           j=j+1;
       end 
    end

    % generating bifurcation map by plotting j-1 element of kth row each time
    if j>1
        plot(c,xmax(k,1:j-1),'k.','MarkerSize',1);
    end 

    hold on;
    index(k)=j-1;    
end
xlabel('Bifuracation parameter c');
ylabel('x max');
title('Bifurcation diagram for c'); 

Answer 1

Это можно сделать совместимым с parfor, сделав несколько относительно простых шагов. Во-первых, parfor работники не могут создавать экранную графику, поэтому нам нужно что-то изменить, чтобы получить результат. В вашем случае это не совсем тривиально, так как ваш первичный результат xmax присваивается не совсем равномерным образом - вы назначаете разное количество элементов на разных итерациях l oop. Мало того, кажется, что невозможно заранее предсказать, сколько столбцов нужно xmax.

Во-вторых, вам нужно внести небольшие изменения в итерацию l oop, чтобы она была совместима с parfor, для которого требуется итерация последовательного целого числа l oop.

Итак, основное изменение заключается в том, чтобы l oop записывал отдельные строки результатов в массив ячеек, который я назвал xmax_cell. Вне parfor l oop преобразование этого обратно в матричную форму тривиально.

Собирая все это вместе, мы в итоге получаем это, которое работает правильно в R2019b, насколько я могу судить:

clc;
a = 0.2; b = 0.2;
crange = 1:0.05:90; % Range for parameter c
tspan = 0:0.1:500; % Time interval for solving Rossler system

% PARFOR loop outputs: a cell array of result rows ...
xmax_cell = cell(numel(crange), 1);
% ... and a track of the largest result row
maxNumCols = 0;

parfor k = 1:numel(crange)
    c = crange(k);
    f = @(t,x) [-x(2)-x(3); x(1)+a*x(2); b+x(3)*(x(1)-c)];
    x0 = [1 1 0]; % initial condition for Rossler system
    [t,x] = ode45(f,tspan,x0); % call ode() to solve Rossler system
    count = find(t>100); % find all the t values which is >10
    x = x(count,:);
    j = 1;
    n = length(x(:,1)); % find the length of vector x1(x in our problem)
    this_xmax = [];
    for i=2 : n-1
        % check for the min value in 1st column of sol matrix
        if (x(i-1,1)+eps) < x(i,1) && x(i,1) > (x(i+1,1)+eps)
            this_xmax(j) = x(i,1);
            j=j+1;
        end
    end

    % Keep track of what's the maximum number of columns
    maxNumCols = max(maxNumCols, numel(this_xmax));
    % Store this row into the output cell array.
    xmax_cell{k} = this_xmax;
end

% Fix up xmax - push each row into the resulting matrix.
xmax = NaN(numel(crange), maxNumCols);
for idx = 1:numel(crange)
    this_max = xmax_cell{idx};
    xmax(idx, 1:numel(this_max)) = this_max;
end

% Plot
plot(crange, xmax', 'k.', 'MarkerSize', 1)
xlabel('Bifuracation parameter c');
ylabel('x max');
title('Bifurcation diagram for c');