Search code examples
matlabperformancefor-loopdifferential-equations

Speeding up matlab for loop

I have a system of 5 ODEs with nonlinear terms involved. I am trying to vary 3 parameters over some ranges to see what parameters would produce the necessary behaviour that I am looking for.
The issue is I have written the code with 3 for loops and it takes a very long time to get the output.
I am also storing the parameter values within the loops when it meets a parameter set that satisfies an ODE event.

This is how I have implemented it in matlab. 

function [m,cVal,x,y]=parameters()

b=5000;
q=0;
r=10^4; 
s=0;
n=10^-8;

time=3000;

m=[];
cVal=[];
x=[];
y=[];

val1=0.1:0.01:5;
val2=0.1:0.2:8;
val3=10^-13:10^-14:10^-11;


for i=1:length(val1)
    for j=1:length(val2)
        for k=1:length(val3)
        options = odeset('AbsTol',1e-15,'RelTol',1e-13,'Events',@eventfunction);
        [t,y,te,ye]=ode45(@(t,y)systemFunc(t,y,[val1(i),val2(j),val3(k)]),0:time,[b,q,s,r,n],options);

         if length(te)==1
             m=[m;val1(i)];
             cVal=[cVal;val2(j)];
             x=[x;val3(k)];
             y=[y;ye(1)];

         end   
        end

    end
end

Is there any other way that I can use to speed up this process?

Profile viewer results enter image description here

I have written the system of ODEs simply with the a format like 

function s=systemFunc(t,y,p)
        s= zeros(2,1);
        s(1)=f*y(1)*(1-(y(1)/k))-p(1)*y(2)*y(1)/(p(2)*y(2)+y(1));
        s(2)=p(3)*y(1)-d*y(2);
end

f,d,k are constant parameters.
The equations are more complicated than what's here as its a system of 5 ODEs with lots of non linear terms interacting with each other.

Solution

  • Tommaso is right. Preallocating will save some time.

    But I would guess that there is fundamentally not a lot you can do since you are running ode45 in a loop. ode45 itself may be the bottleneck.

    I would suggest you profile your code to see where the bottleneck is:

    profile on 
parameters(... )
profile viewer

    I would guess that ode45 is the problem. Probably you will find that you should actually focus your time on optimizing the systemFunc code for performance. But you won't know that until you run the profiler.

    EDIT

    Based on the profiler output and additional code, I see some things that will help

    1. It seems like the vectorization of your values is hurting you. Instead of

      @(t,y)systemFunc(t,y,[val1(i),val2(j),val3(k)])

    try 

    @(t,y)systemFunc(t,y,val1(i),val2(j),val3(k))

    where your system function is defined as 

    function s=systemFunc(t,y,p1,p2,p3)
        s= zeros(2,1);
        s(1)=f*y(1)*(1-(y(1)/k))-p1*y(2)*y(1)/(p2*y(2)+y(1));
        s(2)=p3*y(1)-d*y(2);
end

    1. Next, note that you don't have to preallocate space in the systemFunc, just combine them in the output: 

       function s=systemFunc(t,y,p1,p2,p3)
    s = [ f*y(1)*(1-(y(1)/k))-p1*y(2)*y(1)/(p2*y(2)+y(1)), 
          p3*y(1)-d*y(2) ];
 end

    Finally, note that ode45 is internally taking about 1/3 of your runtime. There is not much you will be able to do about that. If you can live with it, I would suggest increasing your 'AbsTol' and 'RelTol' to more reasonable numbers. Those values are really small, and are making ode45 run for a really long time. If you can live with it, try increasing them to something like 1e-6 or 1e-8 and see how much the performance increases. Alternatively, depending on how smooth your function is, you might be able to do better with a different integrator (like ode23). But your mileage will vary based on how smooth your problem is.