Sweeping initial conditions for a set of ODEs using parfor

Question

I am currently trying to use parfor to sweep across a range of initial conditions for a set of differential equations solved by ode45. The code works fine using two nested for loops but I was hoping parfor could make the process more efficient. Unfortunately, I have run into an issue where the solver is able to solve one of the combinations in the matrix representing initial conditions across a range of variables, but the others seem to have their initial values all set at 0, instead of the values specified by the initial conditions. It may have something to do with the fact that I need to create a matrix of zeros ('P') that the results will be written into, perhaps overwriting the initial conditions(?) Any help would be greatly appreciated.

Thanks, 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

score 0 · Answer 1 · answered Mar 11 '16 at 01:40

0

You have one error, and a few opportunities to simplify the code.

In the parfor loop, you have this line P = zeros(301,5,numel(R),numel(C)); which overwrites P with all zeros at each iteration. Put this before the parfor loop.

The first double-for loop that makes negative elements of c zero can be done using max(c,0), which should be more efficient. You can also do P(:,:,ij)=c(:,1:5) directly.

So you can replace your parfor loop with

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

answered Mar 11 '16 at 01:40

David

8,449
1
22
32

Thanks for the Help! Unfortunately, I need to use parfor to be able to solve my system of equations on reasonable timescales. Additionally, it seems matlab requires placing the P = zeros(... code within the parfor loop, seemingly forcing the overwriting of the variable initial conditions I'm attempting to solve for. – ksm Mar 11 '16 at 18:57
In that case, try `P=zeros(301,5,numel(R)*numel(C));` with a `parfor`. – David Mar 12 '16 at 05:15

Sweeping initial conditions for a set of ODEs using parfor

1 Answers1