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;
}