Recent Changes

Friday, March 11

  1. page Algorithms for Critical Curves edited The following algorithm (implemented in C++) has been developed by Lorenzo Cerboni Baiardi of the …
    The following algorithm (implemented in C++) has been developed by Lorenzo Cerboni Baiardi of the University of Urbino Carlo Bo.
    Computation of Critique Lines of a 2D map
    #include <iostream>
    #include <fstream>
    #include <errno.h>
    #include <string.h>
    #include <cmath>
    #include <sys/time.h>
    #include <boost/numeric/ublas/matrix.hpp>
    #include <boost/numeric/ublas/matrix_proxy.hpp>
    #include <boost/numeric/ublas/vector.hpp>
    #include <boost/numeric/ublas/io.hpp>
    #include <boost/numeric/ublas/assignment.hpp> assignment.hpp had to be modified to work
    #include <boost/random/uniform_int.hpp>
    #include <boost/random/uniform_real.hpp>
    #include <boost/random/normal_distribution.hpp>
    #include <boost/random/variate_generator.hpp>
    #include <boost/random/mersenne_twister.hpp>
    #include <boost/numeric/ublas/vector.hpp>
    #include <boost/numeric/ublas/vector_proxy.hpp>
    #include <boost/numeric/ublas/matrix.hpp>
    #include <boost/numeric/ublas/triangular.hpp>
    #include <boost/numeric/ublas/lu.hpp>
    #include <boost/numeric/ublas/io.hpp>
    #include "boost/multi_array.hpp"
    #include <cassert>
    using namespace std;
    namespace alg = boost::numeric::ublas;
    namespace ublas = boost::numeric::ublas;
    void stop(string a)
    {
    int i=0;
    cout<<a<<endl;
    cin>>i;
    }
    PARAMETERS + MAP
    double a,v,b,c,g;
    double mapSim(double x, double y, double mis)
    {
    return x + v*x*(a-2.*b*x - b*y-c/(1.+(g+mis)*y));
    }
    JACOBIAN OF THE MAP (GIVEN)
    double J(double x, double y)
    {
    double J11 = 1. + v*a - 4.*b*v*x - b*v*y - c*v/(1.+g*y);
    double J22 = 1. + v*a - 4.*b*v*y - b*v*x - c*v/(1.+g*x);
    double J12 = -b*v*x + x*c*v*g/pow(1.+g*y,2.);
    double J21 = -b*v*y + y*c*v*g/pow(1.+g*x,2.);
    return J11*J22-J12*J21;
    }
    FUNCTIONS TO COMPUTE THE N-th DERIVATIVE OF A FUNCTION WITH RESPECT TO ONE VARIABLE
    #define N 200
    double hn[N];
    double M=10;
    double sigma=0.000001;
    int R=2;
    double dy=sigma/2.;
    double delta_t=1.;
    double dt=0.1;
    double estremo_sx=-10.;
    double estremo_dx=10.;
    double D=1.;
    void riempi_Hn(double p)
    {
    hn[0]=1.;
    hn[1]=2.*p;
    for(int i=2; i<N; i++)
    {
    hn[i]=2.*p*hn[i-1]-2.*(i-1)*hn[i-2];
    }
    }
    double fattoriale(double n)
    {
    double fatt=1.;
    for(double s=1.;s<=n; s++)
    {
    fatt=fatt*s;
    }
    return fatt;
    }
    double deltaM_l(double x, double y, int l)
    {
    double somma=0.;
    riempi_Hn((x-y)/(sqrt(2.)*sigma));
    int r=0;
    for(double t=0.; t<M/2.; t++ )
    {
    somma=somma+pow(-1.,l)*pow(-1./4.,t)*(1./(sqrt(2.*M_PI)*fattoriale(t)))*hn[2*r+l];
    r=r+1;
    }
    return exp(-(x-y)*(x-y)/(2.*sigma*sigma))*somma/(pow(2.,l/2.)*pow(sigma,l+1.));
    }
    THIS FUNCTION COMPUTE THE DERIVATIVE OF
    ORDER "l" OF A FUNCTION "f(x, ... other ...)"
    WITH RESPECT TO THE VARIABLE "x"
    double derivata1(derivation variable,
    order of derivation,
    other variables or parameters
    eventually present in the function to be derived)
    double derivata1(double i, int l, double x)
    {
    double der=0.;
    for(double ar=i-100.*sigma; ar<i+100.*sigma; ar=ar+dy)
    der=der+dy*deltaM_l(i,ar,l)*J(x,ar);
    J(x,ar) is the function to be derived. In particular this
    routine compute the derivative of order "l" of the function
    J(x,y) with respect to the "y" variable at the point
    "(x,y) = (x,i)".
    To compute the derivative of J with respect to x rather than y
    replace "J(x,ar)" with "J(ar,y)" giving the value "y" as the
    third input variable of the function "derivata1".
    return der;
    }
    this function compute the point "y" in which J(x,y) is zero,
    where "x" is given as an argument of the function, using the Newton
    recurrence starting at "y0" (the other argument of the "zeriJ" function
    double zeriJ(double x, double y0)
    {
    double precision = pow(10.,-12.);
    double y = y0;
    int i = 0;
    while(fabs(J(x,y))>precision || i<1000)
    {
    y = y - J(x,y)/derivata1(y,1,x);
    i = i + 1;
    }
    return y;
    }
    int main()
    {
    ofstream S2;
    S2.open("s2.dat");
    ofstream S3;
    S3.open("s3.dat");
    ofstream S31;
    S31.open("s31.dat");
    ofstream S32;
    S32.open("s32.dat");
    ofstream S33;
    S33.open("s33.dat");
    ofstream S34;
    S34.open("s34.dat");
    ofstream S35;
    S35.open("s35.dat");
    ofstream S36;
    S36.open("s36.dat");
    ofstream S37;
    S37.open("s37.dat");
    a = 10.;
    v = 0.32;
    b = 0.5;
    c = 2.;
    g = 0.15;
    double X = 0.;
    double Y = 0.;
    double Xp = 3.;
    double Yp = 4.;
    IF YOU WANT YOU CAN PRINT THE TRAJECTORY
    int Uno = 0;
    cout<<"If you want to have a look in the trajectory of the dyamical system press 1, or press any other number if not"<<endl;
    cin>>Uno;
    if(Uno == 1)
    {
    cout<<"If yuo want you can specify there both parameters values and initial conditions"<<endl;
    Uno = 0;
    cout<<"press 1 you you want specify initial conditions"<<endl;
    cin>>Uno;
    if(Uno == 1)
    {
    double p = 0.;
    cout<<"set a = "<<endl;
    cin>>p;
    a = p;
    cout<<"set v = "<<endl;
    cin>>p;
    v = p;
    cout<<"set b = "<<endl;
    cin>>p;
    b = p;
    cout<<"set c = "<<endl;
    cin>>p;
    c = p;
    cout<<"set g = "<<endl;
    cin>>p;
    g = p;
    }
    ofstream S1;
    S1.open("s1.dat");
    for(int t = 0; t<50000; t++)
    {
    X = mapSim(Xp,Yp,0.);
    Y = mapSim(Yp,Xp,0.02);
    Xp = X;
    Yp = Y;
    if you want print only the trajectory at times multiples of T0
    /*int T0 = 2000;
    if(t%T0==0)
    {
    Xp = 3.5 + rand()/((double)RAND_MAX)-0.5;
    Yp = 3.7 + rand()/((double)RAND_MAX)-0.5;
    }*/
    ignore the first transient point T1 of the trajectory
    int T1 = 0;
    if(t>T1)
    S1<<X<<" "<<Y<<endl;
    }
    cout<<"The trajectory is printed in the file s1.dat"<<endl;
    S1.close();
    }
    ROUTINE TO COMPUT CRITIQUE LINES
    stop("Press any number to start the computation of the critique lines");
    SPECIFY THE X-INTERVAL IN WHICH YOU WANT TO FIND THE ZEROS OF YOUR 2-D SURFACE
    SPECIFY ALSO THE POINT "y0" AT WHICH, FOR THE FIRST "x" IN THE CHOSEN INTERVAL, THE NEWTON
    METHOD WILL START (CONVERGENCE TO THE DESIRED ZEROS MAY DEPEND ON THIS STARTING POINT)
    double y0 = 1.;
    the stepsize dx may be changed for the specific map you are invertigating
    double dx = 0.05;
    for(double x = 1.75; x<5.48; x = x + dx)
    {
    COMPUTATION OF THE ZERO GIVEN THE ABSCISSA AND THE FIRST ATTEMPT "y0" OF THE ORDINATE
    double y = zeriJ(x,y0);
    ONCE YOU HAVE FOUND THE FIRST "(X,Y)" PAIR IN WHICH YOUR FUNCTION J(X,Y) = 0
    THAN YOU CAN CHOOSE THE STARTING POINTS y0 FOR THE SUCCESSIVE COMPUTATIONS OF
    THE ZEROS AS YOUR LAST COMPUTATION. IN THIS WAY CONFIDENT THAT YOU ARE STILL
    IN THE BASIN OF ATTRACTION OF THE NEXT ZERO AND THAT THE NEWTON METHOD WILL
    CONVERGE
    y0 = y;
    YOU CAN SPECIFY THE RENGE OF ORDINATE AXIS IN WHICH YOU ARE INTERESTED TO LOOK AT THE ZEROS
    if(y>0. && y<20.)
    {
    PRINT THE ZEROS IN THE FILE "s2.dat"
    S2<<x<<" "<<y<<endl;
    PRINT THE LC0 LINE IN THE FILE "s3.dat"
    S3<<mapSim(x,y,0.)<<" "<<mapSim(y,x,0.)<<endl;
    PRINT THE LC1 LINE IN THE FILE "s31.dat"
    double xLC2 = mapSim(mapSim(x,y,0.),mapSim(y,x,0.),0.);
    double yLC2 = mapSim(mapSim(y,x,0.),mapSim(x,y,0.),0.);
    S31<<xLC2<<" "<<yLC2<<endl;
    PRINT THE LC2 LINE IN THE FILE "s32.dat"
    S32<<mapSim(xLC2,yLC2,0.)<<" "<<mapSim(yLC2,xLC2,0.)<<endl;
    PRINT THE LC3 LINE IN THE FILE "s33.dat"
    double xLC22 = mapSim(mapSim(xLC2,yLC2,0.),mapSim(yLC2,xLC2,0.),0.);
    double yLC22 = mapSim(mapSim(yLC2,xLC2,0.),mapSim(xLC2,yLC2,0.),0.);
    S33<<xLC22<<" "<<yLC22<<endl;
    PRINT THE LC4 LINE IN THE FILE "s34.dat"
    S34<<mapSim(xLC22,yLC22,0.)<<" "<<mapSim(yLC22,xLC22,0.)<<endl;
    PRINT THE LC5 LINE IN THE FILE "s35.dat"
    double xLC23 = mapSim(mapSim(xLC22,yLC22,0.),mapSim(yLC22,xLC22,0.),0.);
    double yLC23 = mapSim(mapSim(yLC22,xLC22,0.),mapSim(xLC22,yLC22,0.),0.);
    S35<<xLC23<<" "<<yLC23<<endl;
    PRINT THE LC6 LINE IN THE FILE "s36.dat"
    S36<<mapSim(xLC23,yLC23,0.)<<" "<<mapSim(yLC23,xLC23,0.)<<endl;
    double xLC24 = mapSim(mapSim(xLC23,yLC23,0.),mapSim(yLC23,xLC23,0.),0.);
    double yLC24 = mapSim(mapSim(yLC23,xLC23,0.),mapSim(xLC23,yLC23,0.),0.);
    PRINT THE LC7 LINE IN THE FILE "s37.dat"
    S37<<xLC24<<" "<<yLC24<<endl;
    IF YOU WANT MORE CRITIQUE LINES AUGMENT THE PRINTINGS,
    OTHERWISE DELETE SOME ERASE OR COMMENT THOSE THAT ARE MORE
    }
    }
    S2.close();
    S3.close();
    S31.close();
    S32.close();
    S33.close();
    S34.close();
    S35.close();
    S36.close();
    S37.close();
    return 0;
    }

    (view changes)
    11:01 am
  2. 10:46 am

Sunday, January 19

  1. page Algorithms Implementation edited ... end end ... dynamical systems. Example Examples taken from ... 546–571 (2013). Downl…
    ...
    end
    end
    ...
    dynamical systems. ExampleExamples taken from
    ...
    546–571 (2013).
    Download all the following files in a single fold in Matlab.
    {simulater.m}

    % ===================================================================================================
    % hybriddde23 Simulater
    % ROUTINE ROUTINE FOR SOLVING
    ...
    MODELS
    %
    %
    % Università della Calabria (UNICAL)
    ...
    Davide Radi
    % email: radidavide85@gmail.com
    % Software developed under PRIN project "Local interactions of global dynamics in economics
    ...
    %
    % ===================================================================================================
    %
    % Matlab-routine for solving hybrid models involving delay-differential
    % equations with constant delays.

    %
    % Matlab-routine for solving hybrid models involving delay-differential
    % equations with constant delays.
    %
    % ARGUMENTS
    ARGUMENTS
    %
    % --------------------------------------------------------------------------------------------------|
    % modelct: | function | function handle that
    ...
    |
    % | for | for an example
    ...
    comand window; | |
    % --------------------------------------------------------------------------------------------------|
    % modeldt: | function | function handle that
    ...
    equations (map), | |
    % | for | for an example
    ...
    comand window; | |
    % --------------------------------------------------------------------------------------------------|
    % XH0c: | initial | initial conditions for
    ...
    history is | |
    % | taken | taken into account); | |
    % --------------------------------------------------------------------------------------------------|
    % XH0d: | initial | initial conditions for
    ...
    in modeldt; | |
    % --------------------------------------------------------------------------------------------------|
    % XH0bh: | initial | initial conditions for
    ...
    in modelbh; | |
    % --------------------------------------------------------------------------------------------------|
    % lags: | vector | vector of constant,
    ...
    for more | |
    % | information); | | information); |
    % --------------------------------------------------------------------------------------------------|
    % deltat: | 1/deltat | 1/deltat indicates the
    ...
    given as | |
    % | results | results by "hybriddde23",
    ...
    this see | |
    % | (options | (options for dde23
    ...
    matlab's help); | |
    % --------------------------------------------------------------------------------------------------|
    % tspan: | interval | interval of integration
    ...
    < tf; | |
    % --------------------------------------------------------------------------------------------------|
    % s: | interation | interation interval for
    ...
    map "modeldt"; | |
    % --------------------------------------------------------------------------------------------------|
    % options: | set | set options for
    ...
    more information; | |
    % --------------------------------------------------------------------------------------------------|
    % par: | parameters' | parameters' values for
    ...
    hybrid model; | |
    % --------------------------------------------------------------------------------------------------|
    % nameTH: | name | name for the
    ...
    is stored; | |
    % --------------------------------------------------------------------------------------------------|
    % nameXHc: | name | name for the
    ...
    is stored; | |
    % --------------------------------------------------------------------------------------------------|
    % nameXHd: | name | name for the
    ...
    is stored; | |
    % --------------------------------------------------------------------------------------------------|
    % nameTbh: | name | name for the
    ...
    is stored; | |
    % --------------------------------------------------------------------------------------------------|
    % nameXbh: | name | name for the
    ...
    is stored; | |
    % --------------------------------------------------------------------------------------------------|
    % flag1: | 1 | 1 if the
    ...
    (for hybridde23), | |
    % | 2 | 2 if the
    % --------------------------------------------------------------------------------------------------|
    % flag2: | 1 | 1 if the
    ...
    dde23-matlab routine, | |
    % | 2 | 2 if the
    ...
    ode45-matlab routine, | |
    % | 0 | 0 if there
    ...
    to solve; | |
    % --------------------------------------------------------------------------------------------------|
    %
    % OUTPUTS

    % OUTPUTS

    %
    % --------------------------------------------------------------------------------------------------|
    % TH: | time | time structure of
    ...
    in "nameTH.dat"; | |
    % --------------------------------------------------------------------------------------------------|
    % XHc: | integral | integral solution of
    ...
    in "nameXHc.dat"; | |
    % --------------------------------------------------------------------------------------------------|
    % XHd: | iteration | iteration results of
    ...
    in "nameXHd.dat; | |
    % --------------------------------------------------------------------------------------------------|
    % Tbh: | time | time structure of
    ...
    in "nameTbh.dat"; | |
    % --------------------------------------------------------------------------------------------------|
    % Xbh: | integral | integral solution of
    ...
    in "nameXHch.dat"; | |
    % --------------------------------------------------------------------------------------------------|
    % flag: | when | when flag=1 appears
    ...
    to work; | |
    % --------------------------------------------------------------------------------------------------|
    %
    % ===================================================================================================
    % ===================================================================================================
    % EXAMPLE EXAMPLE 1
    % RUN RUN hybriddde23 ROUTINE
    ...
    IN
    % Lamantia Lamantia and Radi
    %
    % * Lamantia, F. and Radi, D. 2013, “Exploitation of renewable resources with differential
    ...
    % Define (first) arguments of simulater as follows:
    % --------------------------------------------------------------------------------------------------|
    ...
    = @fishmodeldde23ct (right (right side of
    ...
    delay-differential
    % equations equations defined in
    ...
    = @fishmodeldde23dt (right (right side of
    ...
    difference-equations
    % (maps) (maps) defined in
    % modelbh = '' (becouse benchmark case is not considered);
    % XH0c = = [5 0
    ...

    % XH0d = = [.9];
    % X0bh = = [] (becouse
    ...

    % 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.15, 0.5, 1, 2, 0, 0, 0.4, 3, 3, 10, 1.8, 1, 1, .5, 1,
    ...
    .005];
    % for for the meaning
    ...

    % nameTH = = 'fishmodelTHsol';
    % nameXHc = 'fishmodelXHcsol';
    % nameXHd = 'fishmodelXHdsol';
    % nameTbh = '' (becouse benchmark case is not considered);
    % nameXbh = '' (becouse benchmark case is not considered);
    % flag1 = = 1;
    % flag2 = = [] (becouse
    % --------------------------------------------------------------------------------------------------|
    ...
    following routine-call:
    % flag=simulater(modelct,modeldt,modelbh,XH0c,XH0d,X0bh,lags,deltat,tspan,s,options,par,nameTH,nameXHc,nameXHd,nameTbh,nameXbh,flag1,flag2);
    %
    ...
    % ===================================================================================================
    % ===================================================================================================
    % EXAMPLE EXAMPLE 2
    % RUN RUN hybridode45 ROUTINE
    ...
    IN
    % Bischi, Bischi, Lamantia and
    %
    % * Bischi, G.I., Lamantia, F. and Radi, D. 2013, “Multispecies Exploitation
    ...
    % Define (first) arguments of hybridode45 as follows:
    % --------------------------------------------------------------------------------------------------|
    ...
    = @fishmodelode45ct (right (right side of
    ...
    delay-differential
    % equations equations defined in
    ...
    = @fishmodelode45dt (right (right side of
    ...
    difference-equations
    % (jump-equations, (jump-equations, maps) defined
    ...
    = @fishmodelode45bh (right (right side of
    ...
    difference-equations
    % defined defined in fishmodelode45bh.m,
    ...

    % XH0c = = [10, 10,
    ...

    % XH0d = = [.5];
    % X0bh = = [10, 10
    ...

    % lags = = [] (because
    ...

    % 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 OR = []
    ...

    % par = = [90, 140,
    ...
    s];
    % for for the meaning
    ...

    % nameTH = = 'fishmodelTHsol';
    % nameXHc = 'fishmodelXHcsol';
    % nameXHd = 'fishmodelXHdsol';
    % nameTbh = 'fishmodelTbhsol';
    % nameXbh = 'fishmodelXbhsol';
    % flag1 = = 2;
    % flag2 = = 2;
    % --------------------------------------------------------------------------------------------------|
    ...
    following routine-call:
    % flag=simulater(modelct,modeldt,modelbh,XH0c,XH0d,X0bh,lags,deltat,tspan,s,options,par,nameTH,nameXHc,nameXHd,nameTbh,nameXbh,flag1,flag2);
    %
    ...
    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.m}
    % ===================================================================================================
    % 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.m}
    % ===================================================================================================
    % 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
    {fishmodelode45dt.m}
    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
    {fishmodelode45ct.m}
    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
    {fishmodelode45bh.m}
    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
    {fishmodeldde23dt.m}
    % 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
    {fishmodeldde23ct.m}
    % 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

    (view changes)
    4:49 pm
  2. 4:46 pm
  3. 4:46 pm
  4. 4:46 pm
  5. 4:45 pm
  6. 4:45 pm
  7. 4:44 pm

More