Building a wrapper around script

He guys,

For research I'm running a simulation whereby the Tspan of the simulation varies with the driver frequency. The driver is a signal that is used as input, the simulation needs to simulate approx. 200 periods of the driver. This explains why the Tspan differs per driver frequency. Due to a large time resolution (Fs= 1/1000) i've tried to build a wrapper around the function to simulate this whole serie in chunks of 20 periods of the driver. Subsequently, i was planning on pasting the chunks of timeseries aswell as the outcome together so that I get a full signal simulated in chunks.

However, when I tried to save the outcome (per frequency) my matrix dimensions are conflicting with each other.

This is the wrapper which calls for fre_BG

tau=10; % in ms!!
Delta=1;
J=0;
g=3;
eta=1.5;% parameters used in Fig. 4a
a_mean0= 3.839807185522777e+02;        %a_mean0= 2.935035693817312e+02;
freq0= 37.938323789736410;%freq0= 30.275110904329726;


frequency=[0:1:(4.5*freq0)];         % driving frequency in Hz [note that 30 Hz is not exactly the osc. frequency for eta=1]


%variables wrapper
periods= 20;
total_i= 200/periods;

for i=1:total_i
    
    if i== 1;
        t_per= periods*1000;
        t_span= [0 t_per];
    else
        t_per=(periods*i)*1000;
        t_span= [((periods*(i-1))*1000) t_per];
    end
    
    int= num2str(i);
    
[t,v,r]= fre_BG(t_span,freq0,tau,Delta,J,g,eta,frequency);
save(['output_' int],'t','v','r')
end

And here is fre_BG posted

function [t,v,r]= fre_BG(t_span,freq0,tau,Delta,J,g,eta,frequency)
drive=@(eta,d_eta,f,t) eta+d_eta*sin(2*pi*f/1000*t); % recall that we use ms
d_eta=0.5;                              % driving amplitude
driven_fre=@(t,y,tau,Delta,eta,g,J)[ ...
    Delta/(tau*pi)+2*y(1,:).*y(2,:)-g*y(1,:); ...
    y(2,:).^2+eta(t)-(pi*tau*y(1,:)).^2+tau*J.*y(1,:) ...
    ]/tau;

r= zeros(200001,171);
v= zeros(200001,171);
    
for j=1:numel(frequency)

        % define the integration time dependent on the driving frequency
         if frequency(j)
           T= t_span/frequency(j); % roughly 20 periods of the driver
         else
            T=t_span/freq0;
         end
         
        [t,y]=ode45(@(t,y)driven_fre(t,y,tau,Delta,...
            @(t)drive(eta,d_eta,frequency(j),t),g,J),T(1):1/100:T(end),[0.1,1]);
        
          % convert into the 'real' r and v variables
        r(:,j)=y(:,1).*1000;      % this is to generate 'proper' Hz because time is in ms;
        v(:,j)=y(:,2).*1000;
        
      
end

The problem is the last bit, in order to make it a long signal (per frequency) I need to save the r and v for each frequency it is runned for. So far I did not manage to fix this. Does any of you might a solution for this problem?

Answer 1

[Star Strider]

This appears to run correctly, saving ‘t’, ‘v’, and ‘r’ as cell arrays. (I limited ‘total_i’ and ‘j’ each to 2 in my test.) I slightly changed the save call to make it more efficient, however I did not execute it because I have no need to save those results.

tau=10; % in ms!!
Delta=1;
J=0;
g=3;
eta=1.5;% parameters used in Fig. 4a
a_mean0= 3.839807185522777e+02;        %a_mean0= 2.935035693817312e+02;
freq0= 37.938323789736410;%freq0= 30.275110904329726;


frequency=[0:1:(4.5*freq0)];         % driving frequency in Hz [note that 30 Hz is not exactly the osc. frequency for eta=1]


%variables wrapper
periods= 20;
total_i= 200/periods

for i=1:total_i
    
    if i== 1;
        t_per= periods*1000;
        t_span= [0 t_per];
    else
        t_per=(periods*i)*1000;
        t_span= [((periods*(i-1))*1000) t_per];
    end
    
    int= num2str(i);
    
[t,v,r]= fre_BG(t_span,freq0,tau,Delta,J,g,eta,frequency)
save(sprintf('output_%02d.mat', int),'t','v','r')
end


function [t,v,r]= fre_BG(t_span,freq0,tau,Delta,J,g,eta,frequency)
drive=@(eta,d_eta,f,t) eta+d_eta*sin(2*pi*f/1000*t); % recall that we use ms
d_eta=0.5;                              % driving amplitude
driven_fre=@(t,y,tau,Delta,eta,g,J)[ ...
    Delta/(tau*pi)+2*y(1,:).*y(2,:)-g*y(1,:); ...
    y(2,:).^2+eta(t)-(pi*tau*y(1,:)).^2+tau*J.*y(1,:) ...
    ]/tau;

r= {zeros(200001,171)};
v= {zeros(200001,171)};
    
for j=1:numel(frequency)
    
        % define the integration time dependent on the driving frequency
         if frequency(j)
           T= t_span/frequency(j); % roughly 20 periods of the driver
         else
            T=t_span/freq0;
         end
         
        [t{j},y]=ode45(@(t,y)driven_fre(t,y,tau,Delta,...
            @(t)drive(eta,d_eta,frequency(j),t),g,J),T(1):1/100:T(end),[0.1,1]);
        
          % convert into the 'real' r and v variables
        r{j}=y(:,1).*1000;      % this is to generate 'proper' Hz because time is in ms;
        v{j}=y(:,2).*1000;
        
%         figure
%         plot(t{j}, y, '.-')
%         grid
      
end
end

The figure calls are optional. I put them in because it makes it easier to see the results of the integration.

.