Подходящие начальные условия для набора ODE с использованием parfor

Question

В настоящее время я пытаюсь использовать parfor для охвата ряда начальных условий для набора дифференциальных уравнений, решаемых ode45. Код работает отлично, используя два вложенных цикла, но я надеялся, что parfor может сделать процесс более эффективным. К сожалению, я столкнулся с проблемой, когда решатель может решить одну из комбинаций в матрице, представляющую начальные условия в диапазоне переменных, но остальные, похоже, имеют свои начальные значения, все заданные в 0, вместо указанных значений по начальным условиям. Возможно, это связано с тем, что мне нужно создать матрицу нулей («P»), в которую будут записаны результаты, возможно, переписывая начальные условия (?). Любая помощь будет принята с большой благодарностью.

Спасибо, Kyle

function help(C, R) 

A = 0.01; 
B = 0.00; 
C = [0.001,0.01]; 
D = 0.00; 
R = [1e-10,1e-9]; 

[CGrid,RGrid] = meshgrid(C,R); 

parfor ij = 1 : numel(CGrid) 
     c2 = [A; B; CGrid(ij); D; RGrid(ij)]; 
     [t,c] = ode45('ode_sys',[0:1:300],c2); 
     for k=1:length([0:1:300]) 
      for l=1:length(c2) 
       if c(k,l)<0 
        c(k,l)=0; 
       end 
      end 
     end 

    P = zeros(301,5,numel(R),numel(C)); 
    temp = zeros(301,5);  
    temp(:,1) = c(:,1); 
    temp(:,2) = c(:,2); 
    temp(:,3) = c(:,3); 
    temp(:,4) = c(:,4); 
    temp(:,5) = c(:,5); 
    P(:,:,ij)=temp; 

    parsave('data.mat', P); 
    end 
end 

Answer 1

У вас есть одна ошибка, и несколько возможностей для упрощения кода.

В цикле parfor у вас есть эта строка P = zeros(301,5,numel(R),numel(C));, которая перезаписывает P со всеми нулями на каждой итерации. Поместите это до петли parfor.

Первый двойной цикл for, который делает отрицательные элементы c, может быть равен нулю, используя max(c,0), что должно быть более эффективным. Вы также можете сделать P(:,:,ij)=c(:,1:5) напрямую.

Таким образом, вы можете заменить parfor петлю

P = zeros(301,5,numel(R),numel(C)); 
for ij = 1 : numel(CGrid) 
    c2 = [A; B; CGrid(ij); D; RGrid(ij)]; 
    [t,c] = ode45('ode_sys',0:300,c2); 
    c = max(c,0); 
    P(:,:,ij) = c(:,1:5); 
    parsave('data.mat',P); 
end