Computation Time with respect to local variable assignment

I came across a funny computation time increase, which was entirely unexpected for me.

For context: I am programming a (single) shooting algorithm in an object oriented toolbox (for an ODE solution branch continuation).

The relevant part of the code is the function shooting_single_auto_fun.

Version 1:

function g = shoot_single_auto_fun(obj,x,mu,DYN)

    s = x(1:(end-1),:);         %state vector (last entry is the (unknown) base frequency of the solution)

% Call ode45 to compute a curve section. If the algorithm converges (periodic solution is found), the
% curve is closed.
    [~,temp] = ode45(@(t,y)FCNwrapper(t,y,@(tau,z)DYN.rhs(tau,z,mu)),[0,2*pi],x,obj.odeOpts); 
    f = DYN.rhs(0,DYN.x0(1:(end-1),:),mu);

    g = temp(end,1:(end-1)).'-s;    %The curve is closed, if the first and last point of the curve are identical
    g(end+1,:) = f.'*(s-DYN.x0(1:(end-1),:)); %Autonomous system: Frequency is also unknown. Additional equation needed. This is the Poincare condition

end

%For the autonomous case a function wrapper is needed, which adds the
%equation T' = 0, which needs to be added to the function and multiplies the right hand side f of
%the original ODE with T/2*pi
function dzdt = FCNwrapper(t,y,Fcn)
    z = y(1:(end-1),:);
    base_frqn = y(end,:);       
    dzdt = 1./base_frqn.*Fcn(t,z);
    dzdt(end+1,:) = 0;
end

This function "shoots"/computes one period of a solution curve. The intitial coditions IC (including the base frequency/period duration of the solution) gets adapted by a Newont-type solver (fsolve) until the computed curve is closed (g becomes approximately 0 - a periodic solution is found).

In this context, shoot_single_auto_fun is a method of a general Shooting class and DYN ist an object of class in which all relevant parameters are defined.

An alternative code version reads:

Version 2:

function g = shoot_single_auto_fun(obj,x,mu,DYN)

%This time: allocate the two values
    Fcn = DYN.rhs;  %function handle to the right side of the first order ODE
    x0 = DYN.IC;    %Initial condition - gets adapted by the Newton-type solver
    s = x(1:(end-1),:);

    [~,temp] = ode45(@(t,y)FCNwrapper(t,y,@(tau,z)Fcn(tau,z,mu)),[0,2*pi],x,obj.odeOpts);
    f = Fcn(0,x0(1:(end-1),:),mu);

    g = temp(end,1:(end-1)).'-s;
    g(end+1,:) = f.'*(s-x0(1:(end-1),:)); %This is the Poincare condition

end

%For the autonomous case a function wrapper is needed, which adds the
%equation T' = 0, which needs to be added to the function and multiplies the right hand side f of
%the original ODE with T/2*pi
function dzdt = FCNwrapper(t,y,Fcn)
    z = y(1:(end-1),:);
    base_frqn = y(end,:);       
    dzdt = 1./base_frqn.*Fcn(t,z);
    dzdt(end+1,:) = 0;
end

Now: version 2 is up to 50% faster, than version 1.

Since I know, which one is faster - the problem is "solved". But I wanna know why!

Does anyone have any idea why this is happening?

Answer 1

[Maneet Kaur Bagga]

Hi Simon,

The difference in performance between Version 1 and Version 2 of the code can be attributed to the way function handles are used in MATLAB as:

• In Version 1, the function handle "@(tau,z)DYN.rhs(tau,z,mu)" is created inside the "ode45" call. This means that for each time step, a new function handle is created, which incurs some overhead.
• In Version 2, the function handle "Fcn = DYN.rhs" is assigned outside the "ode45" call. This means that the function handle is created only once, and then reused for each time step. This eliminates the overhead of creating a new function handle repeatedly.
• By avoiding the creation of multiple function handles in each time step, Version 2 reduces the computational overhead and improves the overall performance of the code. This might be a possible reason for the better efficiency and time complexity of version 2.

Hope this helps!

Thank You

Maneet Bagga