How to accelerate the ode45 solver?

Hi, everyone, I'm working with this code (see below) but it takes hours to finish and thus, to show me the results, I want to know if there is a way that I can speed up this, because maybe I'll need to add it a greater tspan, I tried to adjust the step size (as you can see in the code) but I'm not sure it is working, is there anything wrong in here? This is just one of the five codes that I'll need to solve this way, this is the simplest (the other ones have more than 2 ODE to be solved) so I need it to be solve quicker, how can I do it?

Thanks

function [t,x]= cobellilee
global out1

      tspan=[0:5:180];
      x0=[80; 0];
      options=odeset('InitialStep',5,'MaxStep',5)
      [t,x]=ode45(@f,tspan,x0,options);
      out1=[t,x]
      assignin('base','out1',out1);
      figure
      subplot(2,1,1)
      plot(t,x(:,1));
      subplot(2,1,2)
      plot(t,x(:,2));

function dxdt= f(t,x,uuno)
global tiempo
tiempo=t
assignin('base','tiempo',tiempo)

dxdt=[
    (-4.9e-2-x(2))*x(1) + 4.42;
    -9.1e-2*x(2) + 8.96*U1(t)
    ];

function [uuno]=U1(t)
    global uu1 tiempo s1
    uu1=importdata('datossinmodi.mat')
    assignin('base','uu1',uu1)

      s1=size(uu1,1)
      assignin('base','s1',s1)

      tiempo=t
      assignin('base','tiempo',tiempo)

      tiempo=round(tiempo)
      if (s1>=tiempo)
          if (tiempo==0)
              uuno=0;
          else
              uuno=uu1(tiempo)
          end
      else

          uuno=0;
      end

Note: 'datossinmodi.mat' is a 288x1 vector that I need the code to read depending on time, is what I did fine? Thanks again!

Answer 1

[Walter Roberson]

To speed up the routine:

1. do not use global
2. do not use assignin
3. read in datossinmodi.mat once at the beginning and pass the value in, rather than reading it every iteration

http://www.mathworks.com/help/matlab/math/parameterizing-functions.html

[Annie]

The thing is that I'll be doing a GUI with this code, this will go as part of a push button, that is why I need it to be way quicker, and that is why I'm using global and assignin, also I tried to run it this way but didn't worked:

function [t,x,out1,uu1,s1]= cobellilee_lapgde
% global out1
uu1=importdata('datossinmodi.mat')
s1=size(uu1,1)
    tspan=[0:5:10];
    x0=[80; 0];
    options=odeset('InitialStep',5,'MaxStep',5)
    [t,x]=ode45(@f,tspan,x0,options);
    out1=[t,x]
%     assignin('base','out1',out1);
    figure
    subplot(2,1,1)
    plot(t,x(:,1));
    subplot(2,1,2)
    plot(t,x(:,2));

function dxdt= f(t,x,uuno)
% global tiempo
tiempo=t
% assignin('base','tiempo',tiempo)

dxdt=[
    (-4.9e-2-x(2))*x(1) + 4.42;
    -9.1e-2*x(2) + 8.96*U1(t)
    ];

function [uuno]=U1(t,uu1,s1)
%     global uu1 tiempo s1
%     uu1=importdata('datossinmodi.mat')
%     assignin('base','uu1',uu1)

%     s1=size(uu1,1)
%     assignin('base','s1',s1)

      tiempo=t
  %     assignin('base','tiempo',tiempo)

      tiempo=round(tiempo)
      if (s1>=tiempo)
          if (tiempo==0)
              uuno=0;
          else
              uuno=uu1(tiempo)
          end
      else

          uuno=0;
      end

%

[Walter Roberson]

global and assignin are guaranteed to make the code slower, not faster. If you search the blogs for information on which kind of data references are faster or slower, global is the second worst.

      [t,x]=ode45(@(t,y) f(t,y,uu1), tspan, x0, options);

Answer 2

[Kelly Kearney]

If I understand your code correctly, the U1 function is simply interpolating the uu1 vector to the appropriate time, right? If so, you can eliminate all the loading and assigning, as Walter suggests:

uu1 = importdata('datossinmodi.mat');
ufun = @(t) interp1(1:length(uu1), uu1, t, 'nearest', 0);
f = @(t,x) [(-4.9e-2-x(2))*x(1) + 4.42; ...
            -9.1e-2*x(2) + 8.96*ufun(t)];

tspan=[0:5:10];
x0=[80; 0];
options=odeset('InitialStep',5,'MaxStep',5)
[t,x]=ode45(f,tspan,x0,options);

Because it's a variable step solver, ode45 may still be slow for your set of equations if it has to reduce its step size a lot while solving. The MaxStep option only controls the largest possible step it takes, but not the smallest. If your solver code is still pretty slow after you've eliminated the unnecessary loads, try experimenting with other solvers.

[Sean de Wolski]

http://blogs.mathworks.com/loren/2015/09/23/ode-solver-selection-in-matlab/

[Annie]

@Kelly Kearner: Yes, U1 function interpolates uu1 into the proper time, and I tried it and it is working, I checked more detailed the page Walter cited in here, and also checking Sean's page, I just have a doubt about the GUI's code, if in there I'm using assignin and loading multiple times some data, should I change all the code as it to be a nested function? Or can I still have the |assignin|to get some information that the GUI is going to be given by the user and just in the part that solves the ODE have this code? I'm sorry, I'm just starting with this kind of programming, thanks, all!

[Walter Roberson]

Nested would be more efficient. However if you are creating the GUI using GUIDE then it does not create nested functions, so you should set something in the handles structure, then call guidata() to save the change in the master copy of handles, then access the value from the handles structure in the place you need it.

assignin() does work, but it is not recommended for various reasons.

[Kelly Kearney]

I'd have to see a small snippet of your GUI code to see what might be the best way to pass things back and forth. I do very little GUI programming myself, never with GUIDE, and usually just stick variables in the figure's application data when it needs to be shared between multiple scopes. That may or may not work for you.

[Annie]

Hi, and what about if now I'm trying to do the same but I have problems because of some auxiliar functions:

uu1 = importdata('datossinmodi.mat');
ufun = @(t) interp1(1:length(uu1), uu1, t, 'nearest', 0);

          F01= 0.0097; 
    EGP0= 0.0161;
    k12=0.066;
    DG=0;
    AG=0.8;
    tmaxG=40;
    VG= 0.16;
    tmaxI=55;
    Ke=0.138;
    VI=0.12;
    Ka1=0.006;
    Ka2=0.06;
    Ka3=0.03;
    Kb1= 51.2e-4*Ka1; 
    Kb2= 8.2e-4*Ka2;  
    Kb3= 520e-4*Ka3;  
  % Auxiliar Functions:
   G=x(1)/VG;

if G>=4.5
    F01_C=F01;
else
    F01_C=F01*G/4.5;
end

if G>=9
    FR=0.003*(G-9)*VG;
else
    FR=0;
end

   UI=x(4)/tmaxI;

   UG=(DG*AG*t*exp(-t/tmaxG))/(tmaxG)^2;

f = @(t,x) [   
        -((F01_C/VG*G)*x(1))+ x(6)*x(1) + k12*x(2) - FR + UG + EGP0*(1-x(8));
        x(6)*x(1)-(k12+x(7))*x(2);     %(2)
        ufun(t)-(x(3)/tmaxI);                   %(3)
        (x(3)/tmaxI)-(x(4)/tmaxI);             %(4)
         UI/VI - Ke*x(5);                      %(5)
        -Ka1*x(6)+ Kb1*x(5);                   %(6)
        -Ka2*x(7)+ Kb2*x(5);                   %(7)
        -Ka3*x(8)+ Kb3*x(5);                   %(8)
        ];

tspan=[0;100];
x0 = [0.1249; 0.0941; 0; 0;0; 0.7665;0.9519; 0.8473];

global out1
[t,x]=ode45(@f,tspan,x0);
out1=[t,x]
 assignin('base','out1',out1);

   figure
      subplot(2,1,1)
      plot(t,x(:,1));
      subplot(2,1,2)
      plot(t,x(:,5));

And this is the error I get:

>> Hovorka2004_Prueba2
Error: File: Hovorka2004_Prueba2.m
Line: 59 Column: 14
"f" was previously used as a
variable,
conflicting with its use here as the
name of a function or command.
See MATLAB Programming, "How MATLAB
Recognizes Function Calls That Use
Command Syntax" for details.

I tried it with another name instead of 'f' but still I get errors, thanks

[Walter Roberson]

[t,x]=ode45(f,tspan,x0);