In this page implementations of the same algorithms with different languages can be found.

C++

Directional coloring

Here is the code for implementing the "Visualization of unstable periodic points embedded
in chaotic attractors by the directional coloring" by T. Ueta and H. Kawakami.
  //  get ranges from the Graphical Library
  rxmin=graphicR-Xmin();
  rxmax=graphicR-Xmax();
  rymin=graphicR-Ymin();
  rymax=graphicR-Ymax();
   transform ranges in pixel
  graphicR-PixelCoord(rxmin,rymin,scrn);
  pxmin=scrn[0];
  pymax=scrn[1]; // it is not a mistake, y coordinates are upside down
  graphicR-PixelCoord(rxmax,rymax,scrn);
  pxmax=scrn[0];
  pymin=scrn[1]; // it is not a mistake, y coordinates are upside down
   draw points
  for (int pxcur=pxmin+1; pxcurpxmax; pxcur++){   //  for each horizontal pixel
    for (int pycur=pymin+1; pycurpymax; pycur++){ //  for each vertical pixel
       find plane coordinates
      graphicR-XYCoord(pxcur,pycur,coord);
      xoldval=coord[0];                         //     pass the initial coordinates
      yoldval=coord[1];
      for (short int i=0; iEN-Text.ToInt(); i++){ //  perform the prescribed number of iterations
                          //  (given by EN->Text)
           xcurval=mygamma1-Iteration(yoldval); //  compute a map iteration
           ycurval=mygamma2-Iteration(xoldval);
        // update values
        xoldval=xcurval;
        yoldval=ycurval;
 
      } // end iterations
      // compute angle
      if (xcurval-coord[0]!=0){             //  if tangent is defined
        angle= atan2l(ycurval-coord[1],xcurval-coord[0]);
      }
      else{                                // else, if tangent is infinite
        if ((ycurval-coord[1])0)
         angle=PI2.0;
        if ((ycurval-coord[1])0)
         angle=-PI2.0;
      }
      sector=(angle+PI)(PI)180.030.0;   // transform the angle in sector
      switch (sector) {                 //    select color depending on the sector
                case 0  IRight-Canvas-Pixels[pxcur][pycur]=0x002631FD; break;
                case 1  IRight-Canvas-Pixels[pxcur][pycur]=0x00257CFE; break;
                case 2  IRight-Canvas-Pixels[pxcur][pycur]=0x0025EEFE; break;
                case 3  IRight-Canvas-Pixels[pxcur][pycur]=0x0026FDB7; break;
                case 4  IRight-Canvas-Pixels[pxcur][pycur]=0x0025FE30; break;
                case 5  IRight-Canvas-Pixels[pxcur][pycur]=0x0082FD26; break;
                case 6  IRight-Canvas-Pixels[pxcur][pycur]=0x00DDFE25; break;
                case 7  IRight-Canvas-Pixels[pxcur][pycur]=0x00FEC225; break;
                case 8  IRight-Canvas-Pixels[pxcur][pycur]=0x00FF7124; break;
                case 9  IRight-Canvas-Pixels[pxcur][pycur]=0x00FF242F; break;
                case 10  IRight-Canvas-Pixels[pxcur][pycur]=0x00FE2582; break;
                case 11  IRight-Canvas-Pixels[pxcur][pycur]=0x00FE25D2; break;
      } //  end switch
    }  // end for vertical pixels
  }  // end for horizontal pixels



Visual Basic

Directional Coloring

Here you can find the VB version of the code for implementing the Kawakami-Ueta algorithm.
The form where the directional colors will appear is named "Colori" while the initial mask of the program (the one containing input-output textboxes, etc.) is named "Gestione".


Sub ColoriDirezionali()
 
'Definition of local variables
 
Dim Xp As Double 'Coordinates of the evaluating poing
Dim Yp As Double
Dim XA As Double 'Coordinates of the point where the evaluating point is mapped
Dim YA As Double
Dim iteraColor As Long  'Number of iterations
Dim I As Long
 
'Number of iterations for each point of the phase plane
 
iteraColor = Val(Gestione.iteraColor.Text) 'Take the number from the box called "iteraColor" located on the form "Gestione"
'iteraColor = 1                            'You can also give here the number of iteraction (alternative the to previous line)
 
'Iterations
 
For X0 = xmin To xmax Step (xmax - xmin) / 500
    For Y0 = ymin To ymax Step (ymax - ymin) / 500
 
    Xp = X0
    Yp = Y0
 
    X = X0
    Y = Y0
    For I = 1 To iteraColor
    itera
    Next I
 
    XA = X
    YA = Y
 
    'Attribution of the color corresponding to each direction
 
    If (XA - Xp) > 0 And (YA - Yp) > 0 And (XA - Xp) / (YA - Yp) > 2 Then
    Colori.PSet (X0, Y0), RGB(180, 247, 12)
    GoTo 5
    End If
 
    If (XA - Xp) > 0 And (YA - Yp) > 0 And (XA - Xp) / (YA - Yp) < 2 And (XA - Xp) / (YA - Yp) > 0.5 Then
    Colori.PSet (X0, Y0), RGB(247, 212, 12)
    GoTo 5
    End If
 
    If (XA - Xp) > 0 And (YA - Yp) > 0 And (XA - Xp) / (YA - Yp) < 0.5 Then
    Colori.PSet (X0, Y0), RGB(247, 102, 12)
    GoTo 5
    End If
 
    If (Xp - XA) > 0 And (YA - Yp) > 0 And (Xp - XA) / (YA - Yp) < 0.5 Then
    Colori.PSet (X0, Y0), RGB(247, 12, 67)
    GoTo 5
    End If
 
    If (Xp - XA) > 0 And (YA - Yp) > 0 And (Xp - XA) / (YA - Yp) < 2 And (Xp - XA) / (YA - Yp) > 0.5 Then
    Colori.PSet (X0, Y0), RGB(247, 12, 208)
    GoTo 5
    End If
 
    If (Xp - XA) > 0 And (YA - Yp) > 0 And (Xp - XA) / (YA - Yp) > 2 Then
    Colori.PSet (X0, Y0), RGB(176, 12, 247)
    GoTo 5
    End If
 
    If (Xp - XA) > 0 And (Yp - YA) > 0 And (Xp - XA) / (Yp - YA) > 2 Then
    Colori.PSet (X0, Y0), RGB(63, 12, 247)
    GoTo 5
    End If
 
    If (Xp - XA) > 0 And (Yp - YA) > 0 And (Xp - XA) / (Yp - YA) < 2 And (Xp - XA) / (Yp - YA) > 0.5 Then
    Colori.PSet (X0, Y0), RGB(12, 43, 247)
    GoTo 5
    End If
 
    If (Xp - XA) > 0 And (Yp - YA) > 0 And (Xp - XA) / (Yp - YA) < 0.5 Then
    Colori.PSet (X0, Y0), RGB(12, 165, 247)
    GoTo 5
    End If
 
    If (XA - Xp) > 0 And (Yp - YA) > 0 And (XA - Xp) / (Yp - YA) < 0.5 Then
    Colori.PSet (X0, Y0), RGB(12, 247, 204)
    GoTo 5
    End If
 
    If (XA - Xp) > 0 And (Yp - YA) > 0 And (XA - Xp) / (Yp - YA) < 2 And (XA - Xp) / (Yp - YA) > 0.5 Then
    Colori.PSet (X0, Y0), RGB(12, 247, 47)
    GoTo 5
    End If
 
    If (XA - Xp) > 0 And (Yp - YA) > 0 And (XA - Xp) / (Yp - YA) > 2 Then
    Colori.PSet (X0, Y0), RGB(94, 247, 12)
    GoTo 5
    End If
    DoEvents
5    Next Y0
 
DoEvents
Next X0
End Sub
 


Mathematica



Matlab

The following code shows a vector implementation of a two-dimensional map (example taken from Bischi, G.I. and M. Kopel "Equilibrium Selection in a Nonlinear Duopoly Game with Adaptive Expectations" Journal of Economic Behavior and Organization, vol. 46/1, pp. 73-100 (2001)). Note that this implentation provides the n-th iterate of the map. This function has to be stored in a files called "bk2.m" and can be direcly invoked by any .m matlab files, as shown below.

function w=bk2(z,a1,a2,m1,m2,n)
%% example from oligopoly game by Bischi-Kopel jebo 2001
%
 
for i=1:1:n
w(:,:,1)=(1-a1)*z(:,:,1)+a1*m1*z(:,:,2).*(1-z(:,:,2));
w(:,:,2)=(1-a2)*z(:,:,2)+a2*m2*z(:,:,1).*(1-z(:,:,1));
 
z=w;
end
 
end
 
For instance, if we want to analitically calculate periodic points (fixed points, cycle 2, etc.) we can try to invoke the matlab command "fsolve". For this reason we can embed in our program the following code to be later used by fsolve. It takes the previous vectorial definition of the map and translates it for the fsolve command.
% function definition for fsolve
function fz=myfun2(x)
        zf(1,1,1)=x(1);
        zf(1,1,2)=x(2);
        fa=bk2(zf,a1,a2,m1,m2,ciclo);
        fz=[fa(1,1,1)-x(1); fa(1,1,2)-x(2)];
end
At this point we can easily write a routine that calculates periodic points. The following routine is composed by three main parts:
1) Parameter definition: as fsolve needs an initial condition to solve the nonlinear system of equations we define a grid of points and other parameters;
2) Periodic point calculation: this procedure is invoked inside a for cycle, so we can start by looking for fixed points and then periodic points using the same routine on the n-th iterate of the map
3) Eliminate redundant solutions: if matlab finds a solution for every point of the grid we want to take only one of these solutions; moreover we want to say that two periodic points are the same if their distance is lower than the threshold spf.

The final outcome is a matrix of periodic points, called solipf. This matrix has always 3 columns: the first two columns contain the coordinate of the periodic point, whereas the third column indicates the periodicity of the point. It is very useful to store also this information: for example we can plot points with different periodicity with a different color, which is specified by its third coordinate.
%%%%%%%%%%%% Parameters %%%%%%%%%%%%%
xmin=0;xmax=1.2; %grid
ymin=0;ymax=1.2; %grid
nciclo=2; % see below
spf=.0001 ; % see previous discussion
npupf=8; % we consider npupf x npupf points in the grid
% as starting condition for fsolve
pf=[]; % initialize pf
spf=0.0001; % threshold to differentiate among periodic points
%%%%%%%%%%% Main calculation %%%%%%%%
 
for pp=0:nciclo
ciclo=2^pp;% here we calculate periodic points up to 2^nciclo
% grid definition
dx = (xmax-xmin)/npupf; dy = (ymax-ymin)/npupf;
x = [xmin:dx:xmax];  y = [ymin:dy:ymax];
% calculate periodic points starting from every point of the grid
    for i=1:npupf
            for j=1:npupf
        options = optimset('Display', 'off');
        pf2=fsolve(@myfun2,[x(i),y(j)],options);
        pf2=[pf2, ciclo];
        pf=[pf; pf2];
        end
    end
end
 
%%%%%%% the following routine eliminates redundant solutions
pfrid=sortrows(pf,[1,2,3]);
solipf=[pfrid(1,1),pfrid(1,2),pfrid(1,3)]; % initialize solipf
if length(pfrid)>1 % to perform only if we have more than one solution
    for i=1:length(pfrid)-1
        if abs(pfrid(i+1,1)-pfrid(i,1))+abs(pfrid(i+1,2)-pfrid(i,2))> spf
            solipf=[solipf; pfrid(i+1,1),pfrid(i+1,2),pfrid(i+1,3)];
        end
    end
end




The following code can be used to simulate a general class of hybrid nonlinear dynamical systems. Examples taken from Lamantia, F. and D. Radi "Exploitation of renewable resources with differentiated technologies: an evolutionary analysis" Mathematics and Computers in Simulation, in press (2014) and Bischi G. I. and Lamantia, F. and D. Radi "A prey-predator fishery model with endogenous switching of harvesting strategy" Applied Mathematics and Computation,vol. 219/20, pp. 10123–10142 (2013) and Bischi G. I. and Lamantia, F. and D. Radi "Multi-species exploitation with evolutionary switching of harvesting strategies" Natural Resource Modeling, vol. 26/4, pp. 546–571 (2013).
Download all the following files in a single fold in Matlab.





% ===================================================================================================
%                                           Simulater
%                               ROUTINE FOR SOLVING HYBRID MODELS
%                                                                   
%
% Università della Calabria (UNICAL)
% Davide Radi 
% email: radidavide85@gmail.com
% Software developed under PRIN project "Local interactions of global dynamics in economics
% and finance: models and tools" MIUR (Italian Ministry of Education, University and Research)
%
% ===================================================================================================
%                                           
%                   Matlab-routine for solving hybrid models involving delay-differential
%                                   equations with constant delays.
%
%                                           ARGUMENTS
%
% --------------------------------------------------------------------------------------------------|
% modelct:   |   function handle that evaluates the right side of the delay-differential equations, |
%            |   for an example type "edit fishmodeldde23ct" in the comand window;                  |
% --------------------------------------------------------------------------------------------------|
% modeldt:   |   function handle that evaluates the right side of the difference equations (map),   |
%            |   for an example type "edit fishmodeldde23dt" in the comand window;                  |    
% --------------------------------------------------------------------------------------------------|
% XH0c:      |   initial conditions for the state variables in modelct (only constant history is    |
%            |   taken into account);                                                               |
% --------------------------------------------------------------------------------------------------|
% XH0d:      |   initial conditions for the state variables in modeldt;                             |
% --------------------------------------------------------------------------------------------------|
% XH0bh:     |   initial conditions for the state variables in modelbh;                             |
% --------------------------------------------------------------------------------------------------|
% lags:      |   vector of constant, positive delays (see dde23 in matlab's help for more           |
%            |   information);                                                                      |
% --------------------------------------------------------------------------------------------------|
% deltat:    |   1/deltat indicates the number of integration results in a unit of time given as    |
%            |   results by "hybriddde23", it is not related to integration error, for this see     |
%            |   (options for dde23 in matlab's help);                                              |
% --------------------------------------------------------------------------------------------------|
% tspan:     |   interval of integration from t0=tspan(1) to tf=tspan(end) with t0 < tf;            |
% --------------------------------------------------------------------------------------------------|
% s:         |   interation interval for the map "modeldt";                                         |
% --------------------------------------------------------------------------------------------------|
% options:   |   set options for ode45-Matlab routine, see Matlab-help for more information;        |
% --------------------------------------------------------------------------------------------------|
% par:       |   parameters' values for the hybrid model;                                           |
% --------------------------------------------------------------------------------------------------|
% nameTH:    |   name for the file in which the output-varaible TH is stored;                       |
% --------------------------------------------------------------------------------------------------|
% nameXHc:   |   name for the file in which the output-varaible XHc is stored;                      |
% --------------------------------------------------------------------------------------------------|
% nameXHd:   |   name for the file in which the output-varaible XHd is stored;                      |
% --------------------------------------------------------------------------------------------------|
% nameTbh:   |   name for the file in which the output-varaible Tbc is stored;                      |
% --------------------------------------------------------------------------------------------------|
% nameXbh:   |   name for the file in which the output-varaible Xbh is stored;                      |
% --------------------------------------------------------------------------------------------------|
% flag1:     |   1 if the hybrid model does involve delay-differential equations (for hybridde23),  |
%            |   2 if the hybrid model does involve ordinary-differential equations(for hybriode45);|
% --------------------------------------------------------------------------------------------------|
% flag2:     |   1 if the benchmark model has to be solved with dde23-matlab routine,               |
%            |   2 if the benchmark model has to be solved with ode45-matlab routine,               |
%            |   0 if there is not benchmark model to solve;                                        |
% --------------------------------------------------------------------------------------------------|
% 
%                                           OUTPUTS
%
% --------------------------------------------------------------------------------------------------|
% TH:        |   time structure of integral solution of modelct (XHc), saved in "nameTH.dat";       |
% --------------------------------------------------------------------------------------------------|
% XHc:       |   integral solution of modelct,saved in "nameXHc.dat";                               |
% --------------------------------------------------------------------------------------------------|
% XHd:       |   iteration results of modeldt with time structure as TH, saved in "nameXHd.dat;     |
% --------------------------------------------------------------------------------------------------|
% Tbh:       |   time structure of integral solution of modelbh (Xbh), saved in "nameTbh.dat";      |
% --------------------------------------------------------------------------------------------------|
% Xbh:       |   integral solution of modelbh, saved in "nameXHch.dat";                             |
% --------------------------------------------------------------------------------------------------|
% flag:      |   when flag=1 appears in the command windows, the routine has finished to work;      |
% --------------------------------------------------------------------------------------------------|
%
% ===================================================================================================
% ===================================================================================================
%                                            EXAMPLE 1
%                            RUN hybriddde23 ROUTINE FOR THE MODEL IN
%                                    Lamantia and Radi 2013
%
% * Lamantia, F. and Radi, D. 2013, “Exploitation of renewable resources with differential
% technologies: an evolutionary analysis.
%
%
% Define (first) arguments of simulater as follows:
% --------------------------------------------------------------------------------------------------|
% modelct = @fishmodeldde23ct  (right side of the delay-differential
%               equations defined in fishmodeldde23ct.m, edit tofishmodeldde23ct);
% modeldt = @fishmodeldde23dt  (right side of the difference-equations
%               (maps) defined in fishmodeldde23dt.m, edit tofishmodeldde23dt);
% modelbh = '' (becouse benchmark case is not considered);
% XH0c    = [5 0 0];
% XH0d    = [.9];
% X0bh    = [] (becouse benchmark case is not considered);
% lags    = 2;
% deltat  = 1.e-4;
% tspan   = [0, 100];
% s       = 1;
% options = [] no ode45-matlab routine involved;
% par     = [19,      0.15,   0.5,  1,  2,  0,  0, 0.4,  3, 10, 1.8,  1, .5, 1, 1, .005];
%           for the meaning of the parameters see paper *
% nameTH  = 'fishmodelTHsol';
% nameXHc = 'fishmodelXHcsol';
% nameXHd = 'fishmodelXHdsol';
% nameTbh = '' (becouse benchmark case is not considered);
% nameXbh = '' (becouse benchmark case is not considered);
% flag1   = 1;
% flag2   = [] (becouse benchmark case is not considered);
% --------------------------------------------------------------------------------------------------|
% type (then) in the command windows the following routine-call: 
% flag=simulater(modelct,modeldt,modelbh,XH0c,XH0d,X0bh,lags,deltat,tspan,s,options,par,nameTH,nameXHc,nameXHd,nameTbh,nameXbh,flag1,flag2);
%
% Plot the results:
% TH=load('fishmodelTHsol'); XHc=load('fishmodelXHcsol'); XHd=load('fishmodelXHdsol');
% figure(1); plot(XHc(:,1),XHd(:,1)); figure(2); plot(TH,XHc(:,1));
% figure(3); plot(TH,XHd); figure(4); plot(TH,XHc(:,2:3));
%
% ===================================================================================================
% ===================================================================================================
%                                            EXAMPLE 2
%                            RUN hybridode45 ROUTINE FOR THE MODEL IN
%                                 Bischi, Lamantia and Radi 2013
%
% * Bischi, G.I., Lamantia, F. and Radi, D. 2013, “Multispecies Exploitation
% with Evolutionary Switching of Harvesting Strategies", Natural Resources Modelling.
%
%
% Define (first) arguments of hybridode45 as follows:
% --------------------------------------------------------------------------------------------------|
% modelct = @fishmodelode45ct  (right side of the delay-differential
%               equations defined in fishmodelode45ct.m, edit tofishmodelode45ct);
% modeldt = @fishmodelode45dt  (right side of the difference-equations
%               (jump-equations, maps) defined in fishmodelode45dt.m, edit tofishmodelode45dt);
% modelbh = @fishmodelode45bh  (right side of the difference-equations
%               defined in fishmodelode45bh.m, edit tofishmodelode45bh);
% XH0c    = [10, 10, 0, 0];
% XH0d    = [.5];
% X0bh    = [10, 10 , .5];
% lags    = [] (because no delay are involved);
% deltat  = 1.e-4;
% tspan   = [0, 100];
% s       = 3;
% options = odeset('RelTol',1e-10,'AbsTol',[1e-10 1.e-10 1.e-10 1.e-10]);
%               OR = [] if no specification is require;
% par     = [90, 140, 80, 80, 40, 50, 50, 9, 9, s];
%           for the meaning of the parameters see paper *
% nameTH  = 'fishmodelTHsol';
% nameXHc = 'fishmodelXHcsol';
% nameXHd = 'fishmodelXHdsol';
% nameTbh = 'fishmodelTbhsol';
% nameXbh = 'fishmodelXbhsol';
% flag1   = 2;
% flag2   = 2;
% --------------------------------------------------------------------------------------------------|
% type (then) in the command windows the following routine-call: 
% flag=simulater(modelct,modeldt,modelbh,XH0c,XH0d,X0bh,lags,deltat,tspan,s,options,par,nameTH,nameXHc,nameXHd,nameTbh,nameXbh,flag1,flag2);
%
% Plot the results:
% TH=load('fishmodelTHsol'); XHc=load('fishmodelXHcsol');
% XHd=load('fishmodelXHdsol'); Thb=load('fishmodelTbhsol'); Xhb=load('fishmodelXbhsol');
% figure(1); plot(XHc(:,1),XHc(:,2)); figure(2); plot(TH,XHc(:,1:2));
% figure(3); plot(TH,XHc(:,3:4)); figure(4); plot(TH,XHd); figure(5); plot3(XHc(:,1),XHc(:,2),XHd);
%
% ===================================================================================================
 
 
function flag=simulater(modelct,modeldt,modelbh,XH0c,XH0d,X0bh,lags,deltat,tspan,s,options,par,nameTH,nameXHc,nameXHd,nameTbh,nameXbh,flag1,flag2);
    if flag1==1;
        % use hybriddde23 routine to solve hybrid models made of delay-differential equations and difference equations::
        [TH,XHc,XHd]=hybriddde23(modelct,modeldt,XH0c,XH0d,lags,deltat,tspan,s,par);
        dlmwrite(nameTH,TH);
        dlmwrite(nameXHc,XHc);
        dlmwrite(nameXHd,XHd);
    elseif flag2==2
        % use hybridode45 routine to solve hybrid models made of ordinary-differential equations and difference equations:
        [TH,XHc,XHd]=hybridode45(modelct,modeldt,XH0c,XH0d,deltat,tspan,s,par,options);
        dlmwrite(nameTH,TH);
        dlmwrite(nameXHc,XHc);
        dlmwrite(nameXHd,XHd);
    else
        disp('flag1 must be either 1 or 2');
    end
    if flag2==1;
        % use dde23-Matlab routine to solve the benchmark model made of delay-differential-equations:
        sol=dde23(modelbh,lags,X0bh,tspan,[],par);
        dlmwrite(nameTbh,sol.x);
        dlmwrite(nameXbh,so.y);   
    elseif flag2==2;
        % use ode45-Matlab routine to solve the benchmark model made of ordinary-differential-equations:
        [Tbh,Xbh]=ode45(modelbh,tspan,X0bh,[],par);
        dlmwrite(nameTbh,Tbh);
        dlmwrite(nameXbh,Xbh); 
    else
        disp('No-benchmark model is considered');
    end
flag=1



% ===================================================================================================
%                                           hybridode45
%                               ROUTINE FOR SOLVING HYBRID MODELS
%                                                                   
%
% Università della Calabria (UNICAL)
% Davide Radi 
% email: radidavide85@gmail.com
% Software developed under PRIN project "Local interactions of global dynamics in economics
% and finance: models and tools" MIUR (Italian Ministry of Education, University and Research)
%
% ===================================================================================================
%                                           
%                   Matlab-routine for solving hybrid models involving ordinary-differential
%                                   equations with constant delays.
%
%                                           ARGUMENTS
%
% --------------------------------------------------------------------------------------------------|
% modelct:   |   function handle that evaluates the right side of the ordinary-differential         |
%            |   equations, for an example type "edit fishmodelode45ct" in the comand window;       |
% --------------------------------------------------------------------------------------------------|
% modeldt:   |   function handle that evaluates the right side of the difference equations (map),   |
%            |   for an example type "edit fishmodelode45dt" in the comand window;                  |    
% --------------------------------------------------------------------------------------------------|
% XH0c:      |   initial conditions for the state variables in modelct;                             |
% --------------------------------------------------------------------------------------------------|
% XH0d:      |   initial conditions for the state variables in modeldt;                             |
% --------------------------------------------------------------------------------------------------|
% deltat:    |   1/deltat indicates the number of integration results in a unit of time given as    |
%            |   results by "hybridode45", it is not related to integration error, for this see     |
%            |   (options for ode45 in matlab's help);                                              |
% --------------------------------------------------------------------------------------------------|
% tspan:     |   interval of integration from t0=tspan(1) to tf=tspan(end) with t0 < tf;            |
% --------------------------------------------------------------------------------------------------|
% s:         |   interval of application of the jump-equation "modeldt";                            |
% --------------------------------------------------------------------------------------------------|
% par:       |   row vector of parameters' values for the hybrid model;                             |
% --------------------------------------------------------------------------------------------------|
% options:   |   optional integration argument; a structure you create using the ddeset function.   |
%            |   (see ddeset for detailsinformation); If no specification is required set = []      |
% --------------------------------------------------------------------------------------------------|
% 
% ===================================================================================================
%                                            EXAMPLE
%                            RUN hybridode45 ROUTINE FOR THE MODEL IN
%                                 Bischi, Lamantia and Radi 2013
%
% * Bischi, G.I., Lamantia, F. and Radi, D. 2013, “Multispecies Exploitation
% with Evolutionary Switching of Harvesting Strategies", Natural Resources Modelling.
%
%
% Define (first) arguments of hybridode45 as follows:
% --------------------------------------------------------------------------------------------------|
% modelct = @fishmodelode45ct  (right side of the delay-differential
%               equations defined in fishmodelode45ct.m, edit tofishmodelode45ct);
% modeldt = @fishmodelode45dt  (right side of the difference-equations
%               (jump-equations, maps) defined in fishmodelode45dt.m, edit tofishmodelode45dt);
% XH0c    = [10, 10, 0, 0];
% XH0d    = [.5];
% deltat  = 1.e-4;
% tspan   = [0, 100];
% s       = 3;
% par     = [90, 140, 80, 80, 40, 50, 50, 9, 9, s];
%           for the meaning of the parameters see paper *
% options = odeset('RelTol',1e-10,'AbsTol',[1e-10 1.e-10 1.e-10 1.e-10]);
%               OR = [] if no specification is require;
% --------------------------------------------------------------------------------------------------|
% type (then) in the command windows the following routine-call: 
% [TH,XHc,XHd]=hybridode45(modelct,modeldt,XH0c,XH0d,deltat,tspan,s,par,options)
%
% Plot the results:
% figure(1); plot(XHc(:,1),XHc(:,2)); figure(2); plot(TH,XHc(:,1:2));
% figure(3); plot(TH,XHc(:,3:4)); figure(4); plot(TH,XHd); figure(5); plot3(XHc(:,1),XHc(:,2),XHd);
%
% ===================================================================================================
 
function [TH,XHc,XHd]=hybridode45(modelct,modeldt,XH0c,XH0d,deltat,tspan,s,par,options);
niter = tspan(2)/s;
TH=tspan(1);
XHc=XH0c;
XHd=XH0d;
parC=[par,XH0d]; % the discrete dynamical variable of the hybrid model is passed to the ode45-matlab routin as a parameter;
if tspan(2)<s
    [tmpTH,tmpXHc]=ode45(modelct,[tspan(1):deltat:tspan(2)],XH0c,options,parC);
    TH=[TH;tmpTH];XHc=[XHc;tmpXHc];XHd=[XHd;repmat(XH0d,length(tmpXHc(:,1)),1)];
    clear tmpTH tmpXHc
else
    for k = 1:1:floor(niter);
        t=(k-1)*s;
        [tmpTH,tmpXHc]=ode45(modelct,[t, t+s],XH0c,options,parC);
        TH=[TH;tmpTH];XHc=[XHc;tmpXHc];XHd=[XHd;repmat(XH0d,length(tmpXHc(:,1)),1)];
        %tindex=length(tmpTH);
        XH0c=tmpXHc(end,:);
        [tmpXHd]=modeldt(XH0c,XH0d,par);
        XH0d=tmpXHd;
        parC=[par,XH0d];
        XH0c=tmpXHc(end,:);
        clear tmpTH tmpXHc tmpXHd
    end
    if (t+s+deltat)<tspan(2);
        [tmpTH,tmpXHc]=ode45(modelct,[t+s, tspan(2)],XH0c,options,parC);
        TH=[TH;tmpTH];XHc=[XHc;tmpXHc];XHd=[XHd;repmat(XH0d,length(tmpXHc(:,1)),1)];
        clear tmpTH tmpXHc
    end
end




% ===================================================================================================
%                                           hybriddde23
%                               ROUTINE FOR SOLVING HYBRID MODELS
%                                                                   
%
% Università della Calabria (UNICAL)
% Davide Radi 
% email: radidavide85@gmail.com
% Software developed under PRIN project "Local interactions of global dynamics in economics
% and finance: models and tools" MIUR (Italian Ministry of Education, University and Research)
%
% ===================================================================================================
%                                           
%                   Matlab-routine for solving hybrid models involving delay-differential
%                                   equations with constant delays.
%
%                                           ARGUMENTS
%
% --------------------------------------------------------------------------------------------------|
% modelct:   |   function handle that evaluates the right side of the delay-differential equations, |
%            |   for an example type "edit fishmodeldde23ct" in the comand window;                  |
% --------------------------------------------------------------------------------------------------|
% modeldt:   |   function handle that evaluates the right side of the difference equations (map),   |
%            |   for an example type "edit fishmodeldde23dt" in the comand window;                  |    
% --------------------------------------------------------------------------------------------------|
% XH0c:      |   initial conditions for the state variables in modelct (only constant history is    |
%            |   taken into account);                                                               |
% --------------------------------------------------------------------------------------------------|
% XH0d:      |   initial conditions for the state variables in modeldt;                             |
% --------------------------------------------------------------------------------------------------|
% Lags:      |   vector of constant, positive delays (see dde23 in matlab's help for more           |
%            |   information);                                                                      |
% --------------------------------------------------------------------------------------------------|
% deltat:    |   1/deltat indicates the number of integration results in a unit of time given as    |
%            |   results by "hybriddde23", it is not related to integration error, for this see     |
%            |   (options for dde23 in matlab's help);                                              |
% --------------------------------------------------------------------------------------------------|
% tspan:     |   interval of integration from t0=tspan(1) to tf=tspan(end) with t0 < tf;            |
% --------------------------------------------------------------------------------------------------|
% s:         |   interation interval for the map "modeldt";                                         |
% --------------------------------------------------------------------------------------------------|
% par:       |   parameters' values for the hybrid model;                                           |
% --------------------------------------------------------------------------------------------------|
% 
% ===================================================================================================
%                                            EXAMPLE
%                            RUN hybriddde23 ROUTINE FOR THE MODEL IN
%                                    Lamantia and Radi 2013
%
% * Lamantia, F. and Radi, D. 2013, “Exploitation of renewable resources with differential
% technologies: an evolutionary analysis.
%
%
% Define (first) arguments of hybriddde23 as follows:
% --------------------------------------------------------------------------------------------------|
% modelct = @fishmodeldde23ct  (right side of the delay-differential
%               equations defined in fishmodeldde23ct.m, edit tofishmodeldde23ct);
% modeldt = @fishmodeldde23dt  (right side of the difference-equations
%               (maps) defined in fishmodeldde23dt.m, edit tofishmodeldde23dt);
% XH0c    = [5 0 0];
% XH0d    = [.9];
% lags    = 2;
% deltat  = 1.e-4;
% tspan   = [0, 100];
% s       = 1;
% par     = [19,      0.15,   0.5,  1,  2,  0,  0, 0.4,  3, 10, 1.8,  1, .5, 1, 1, .005];
%           for the meaning of the parameters see paper *
% --------------------------------------------------------------------------------------------------|
% type (then) in the command windows the following routine-call: 
% [TH,XHc,XHd]=hybriddde23(modelct,modeldt,XH0c,XH0d,lags,deltat,tspan,s,par)
%
% Plot the results:
% figure(1); plot(XHc(:,1),XHd(:,1)); figure(2); plot(TH,XHc(:,1));
% figure(3); plot(TH,XHd); figure(4); plot(TH,XHc(:,1));
%
% ===================================================================================================
 
function [TH,XHc,XHd]=hybriddde23(modelct,modeldt,XH0c,XH0d,lags,deltat,tspan,s,par);
    %format long
    t=tspan(1);
    niter=ceil(tspan(2)-tspan(1));
    TH=tspan(1);
    if s>0
        parM=[par,XH0d];
        tspand = [t+deltat:deltat:t+s];
        opts=[];
        sol=dde23(modelct,lags,XH0c',tspand,opts,parM) % returns a column vector of yi, as sol.y=[y1,y2,y3,...yn]
        XHd=[repmat(XH0d,length(sol.y(1,:)),1)];
        for i=1:s:niter
            t=t+s;
           [tmpXHd]=modeldt(sol.y',XH0d,par);
            if isrow(tmpXHd)==0;
                tmpXHd=tmpXHd';
            end
            XH0d=tmpXHd;
            clear tmpXHd
            parM=[par,XH0d];
            Test=sol.y;
            opts=ddeset('Jumps', sol.x(end));
            sol=dde23(modelct,lags,sol,[sol.x(end):deltat:t+s],opts,parM);
            XHd=[XHd;repmat(XH0d,length(sol.y(1,:))-length(Test(1,:)),1)];
        end
        XHc=sol.y';
        TH=sol.x';
    end
end




function dX=fishmodelode45dt(y,XH0d,par);
    rho1=par(1);
    rho2=par(2);
    K1=par(3);
    K2=par(4);
    N=par(5);
    a1=par(6);
    a2=par(7);
    gamma1=par(8);
    gamma2=par(9);
    dX=[XH0d*y(3)/(XH0d*y(3)+(1-XH0d)*y(4))];
end




function dydt=fishmodelode45ct(t,y,par)
    rho1=par(1);
    rho2=par(2);
    K1=par(3);
    K2=par(4);
    N=par(5);
    a1=par(6);
    a2=par(7);
    gamma1=par(8);
    gamma2=par(9);
    s=par(10);
    XH0d=par(11); %it is r
    dydt=[  y(1)*rho1*(1-y(1)/K1)-N*XH0d*a1*y(1)/(2*gamma1)
            y(2)*rho2*(1-y(2)/K2)-N*(1-XH0d)*a2*y(2)/(2*gamma2)
            a1*a1*y(1)/(s*4*gamma1)
            a2*a2*y(2)/(s*4*gamma2)];
end




function dydt=fishmodelode45bh(t,y,par)
    rho1=par(1);
    rho2=par(2);
    K1=par(3);
    K2=par(4);
    N=par(5);
    a1=par(6);
    a2=par(7);
    gamma1=par(8);
    gamma2=par(9);
    dydt=[  y(1)*rho1*(1-y(1)/K1)-N*y(3)*a1*y(1)/(2*gamma1)
            y(2)*rho2*(1-y(2)/K2)-N*(1-y(3))*a2*y(2)/(2*gamma2)
            y(3)*(1-y(3))*(a1*a1*y(1)/(4*gamma1)-a2*a2*y(2)/(4*gamma2))];
end




% Example for modeldt:
function dX=fishmodeldde23dt(X,XH0d,par)
    %format long
    theta=par(16);
    dX=XH0d*(exp(theta*X(end,2)))/(XH0d*exp(theta*X(end,2))+(1-XH0d)*exp(theta*X(end,3)));
end




% Example for modelct:
% y is state matrix;
% par is the parameters' vector;
% deltat is the size of the integration step;
function dydt=fishmodeldde23ct(t,X,Z,par)
    alpha=par(1);
    beta=par(2);
    sigma=par(3);
    a1=par(4);
    a2=par(5);
    b1=par(6);
    b2=par(7);
    c1=par(8);
    c2=par(9);
    N=par(10);
    q1=par(11);
    q2=par(12);
    g=par(13);
    k=par(14);
    s=par(15);
    XH0d=par(17);
    lag1=Z(2);
    lag2=Z(3);
    dydt =[ X(1)*(alpha-beta*X(1)-N*(XH0d*a1*q1+(1-XH0d)*a2*q2)/(2*g))
            k/(1-exp(-k*s))*(X(1)*a1*a1*q1/(4*g)-c1 -exp(-k*s)*(lag1*a1*a1*q1/(4*g)-c1) -(1-exp(-k*s))*X(2))
            k/(1-exp(-k*s))*(X(1)*a2*a2*q2/(4*g)-c2 -exp(-k*s)*(lag2*a2*a2*q2/(4*g)-c2) -(1-exp(-k*s))*X(3))];
end