• August 02, 2021, 08:23:13 AM

Login with username, password and session length

Author Topic:  rational function  (Read 1619 times)

0 Members and 1 Guest are viewing this topic.

Offline Adam Majewski

  • Fractal Frankfurter
  • *
  • Posts: 541
rational function
« on: June 17, 2021, 05:16:22 PM »
I want to draw   Julia set for rational function (image bo36.png)  by  Michael Becker
Quote
f(z)=1/((0,15+0,15i)z5+z3+(-3+3i)z), dargestellt auf [-3;3]x[-3;3].

It look like connected Julia set. To speed up computaions I want to find critical points and periodic points

Code: [Select]

kill(all);
remvalue(all);
display2d:false;


f(z):= 1/(a*z^5+ z^3 + b*z);


/*
converts complex number z = x*y*%i
to the list in a draw format: 
[x,y]
*/
d(z):=[float(realpart(z)), float(imagpart(z))]$

ToPoints(myList):= points(map('d,myList))$


/* give Draw List from one point*/
ToPoint(z):=points([d(z)])$





/* f(z) is used as a global function
   I do not know how to put it as a argument */

GiveOrbit(z0):=
block(
 [z,Orbit, iLength ],
 iLength : 100,
 z:z0,
 Orbit:[z0],
 for i:1 thru iLength step 1 do
        ( z:f(z),
          z: float(z),
          z:rectform(z),
          Orbit:endcons(z,Orbit)),
 
 Orbit:ToPoints(Orbit),       
 return(Orbit)

)$




/* find fixed points  returns a list */
GiveFixedPoints():= block
(
  [s,z],
  s:solve(f(z)=z),
  /* remove "z="  from list s */
  s:map('rhs,s),
  s:map('rectform,s),
  s:map('float,s),
  return(s)
)$




/*compile(all);*/


a: (0.15+0.15*%i);
b: (-3+3*%i);


dz: diff(f(z) ,z,1)$

dz:ratsimp(float(ev(dz)));
 
 /*
 dz = -((150*%i+150)*z^4+600*z^2+600*%i-600) /(9*%i*z^10+(60*%i+60)*z^8-160*z^6+(1200*%i-1200)*z^4-3600*%i*z^2)
*/

s: solve(dz=0)$
s:map('rhs,s)$
s:map('float,s)$
s:map('rectform,s)$
s:map('float,s);
/*


 [0.389374737779177*%i-0.9400337727919567,
        0.9400337727919567-0.389374737779177*%i,
        (-1.816005797447402*%i)-0.7522142306508819,
        1.816005797447402*%i+0.7522142306508819]
 
*/


/* r1:solve(f(z)=z); */




iLength:100;
orbit1 :GiveOrbit(s[1])$
orbit2 :GiveOrbit(s[2])$
orbit3 :GiveOrbit(s[3])$
orbit4 :GiveOrbit(s[4])$



path:"~/Dokumenty/ijon/3/36/"$  /*  if empty then file is in a home dir */

load(draw); /* ( interface to gnuplot ) by Mario Rodriguez Riotorto http://www.telefonica.net/web2/biomates */

draw2d(
    title = concat("Critical orbit for f(z)=z^2 +", string(c)),
    terminal  = png,
    user_preamble = "set size square", /*    */
    file_name = concat(path ,"8"),
    dimensions = [1500,1500],    /* Since Maxima 5.23, pic_width and pic_height are deprecated. */
    xrange = [-3,3],
    yrange = [-3,3],
    xlabel     = "z.re ",
    ylabel     = "z.im",

   
   
    point_type    = 7,
    points_joined = false,
    point_size    = 0.8,
    key=" critical orbit ",
    color             =red,
    orbit1,
    orbit2,
    orbit3,
    orbit4,
   
    point_size    = 1.2,
    key="",
    ToPoint(s[1]),
    ToPoint(s[2]),
    ToPoint(s[3]),
    key= "critical point",
    color           = blue,
    ToPoint(s[4])
   
    /*
   

    key= "fixed point",
    color           = black,
    alfa
*/
 );


My image of critical orbits  shows that critical orbits tend to z=0, but on the image it seems to be in the exterior
How to find periodic points ?







Linkback: https://fractalforums.org/index.php?topic=4279.0
« Last Edit: June 17, 2021, 06:15:18 PM by Adam Majewski »

Offline claude

  • Uploader
  • *
  • Posts: 1938
    • mathr.co.uk
Re: rational function
« Reply #1 on: June 17, 2021, 05:27:02 PM »
critical orbits tend to z=0
Are you sure about that?  Possibly tending to fixed points near (but not equal to) 0?  I have not investigated, but it looks like there is a gap in the middle, in both images.

EDIT fixed points and period 2 points:
Code: [Select]
(%i10) ToPoints(solve(z = f(z), z));

rat: replaced 0.15 by 3/20 = 0.15

rat: replaced 0.15 by 3/20 = 0.15
(%o10) points([[z = - 0.1728968618954043, 0.0 = 0.4517976427751508],
[z = 0.1728968618954043, 0.0 = - 0.4517976427751508],
[z = - 1.542012716025067, 0.0 = 0.5826223205347855],
[z = 1.542012716025067, 0.0 = - 0.5826223205347855],
[z = - 1.052616733765589, 0.0 = - 2.511067897204381],
[z = 1.052616733765589, 0.0 = 2.511067897204381]])
(%i18) map('allroots, solve(f(f(z))=z));

rat: replaced 0.15 by 3/20 = 0.15

rat: replaced 0.15 by 3/20 = 0.15

rat: replaced 0.15 by 3/20 = 0.15

rat: replaced 0.15 by 3/20 = 0.15

rat: replaced 0.15 by 3/20 = 0.15

rat: replaced 0.15 by 3/20 = 0.15

rat: replaced 0.15 by 3/20 = 0.15

rat: replaced 0.15 by 3/20 = 0.15
(%o18) [[z = 0.0], [z = 0.1004507524894826 - 0.04160806403172667 %i,
z = 0.04160806403172667 %i - 0.1004507524894826,
z = (- 0.1972126334382917 %i) - 0.4417257204225227,
z = 0.1728968618954057 - 0.4517976427751519 %i,
z = 0.1972126334382916 %i + 0.4417257204225231,
z = 0.4517976427751519 %i - 0.172896861895406,
z = (- 2.008427151852299 %i) - 0.8319177653355924,
z = 1.486804283840914 - 0.6158544989534238 %i,
z = 0.5826223205446455 %i - 1.542012716044921,
z = 0.6783914544446792 %i - 1.502343841878277,
z = 1.545643025048946 - 0.6402263035701801 %i,
z = 1.542012716021183 - 0.5826223205577208 %i,
z = 0.6158544989143758 %i - 1.48680428382238,
z = 0.6402263036121009 %i - 1.545643025068914,
z = 1.502343841903455 - 0.678391454434491 %i,
z = 2.511067896520314 %i + 1.052616752123668,
z = 2.008427151854269 %i + 0.83191776533551,
z = (- 2.519905568763936 %i) - 1.031280709377238,
z = (- 2.511067897025448 %i) - 1.052616732010429,
z = 2.519905568548988 %i + 1.03128068923336,
z = 2.496853258644555 %i + 1.034230487296836,
z = (- 2.532331263153173 %i) - 1.048925954202523,
z = (- 2.496853270435978 %i) - 1.034230487368577,
z = 2.532331275662722 %i + 1.048925954304978]]
« Last Edit: June 17, 2021, 05:41:13 PM by claude, Reason: points »

Offline Adam Majewski

  • Fractal Frankfurter
  • *
  • Posts: 541
Re: rational function
« Reply #2 on: June 17, 2021, 05:37:10 PM »
Are you sure about that?  Possibly tending to fixed points near (but not equal to) 0?  I have not investigated, but it looks like there is a gap in the middle, in both images.

Ok. It looks that z=0 will be inside ( ?center) of very small component which is smth like limit of white components.

So period 1 point will be not helpfull for me.

I ask about it because usually testing for attracting periodic points makes the code faster.

Offline claude

  • Uploader
  • *
  • Posts: 1938
    • mathr.co.uk
Re: rational function
« Reply #3 on: June 17, 2021, 05:42:42 PM »
see edit, period 2 points include 0 and +/-(0.04160806403172667 %i - 0.1004507524894826)

Offline Adam Majewski

  • Fractal Frankfurter
  • *
  • Posts: 541
Re: rational function
« Reply #4 on: June 17, 2021, 07:29:01 PM »
Code: [Select]

/*

https://web.archive.org/web/20161024194536/http://www.ijon.de/mathe/julia/some_julia_sets_3.html

a: (0,15+0,15i)
b: (-3+3i)
f: 1/(a*z^5+ b*z^3 + z);
 dz = -(5*a*z^4+3*b*z^2+1)/(a*z^5+z^3+b*z)^2


  Adam Majewski
  adammaj1 aaattt o2 dot pl  // o like oxygen not 0 like zero
 
 
 
  Structure of a program or how to analyze the program
 
 
  ============== Image X ========================
 
  DrawImageOf -> DrawPointOf -> ComputeColorOf ( FunctionTypeT FunctionType , complex double z) -> ComputeColor
 
 
  check only last function  which computes color of one pixel for given Function Type
 
 

   
  ==========================================

 
  ---------------------------------
  indent d.c
  default is gnu style
  -------------------



  c console progam
 
  export  OMP_DISPLAY_ENV="TRUE"
  gcc d.c -lm -Wall -march=native -fopenmp
  time ./a.out > b.txt


  gcc d.c -lm -Wall -march=native -fopenmp


  time ./a.out

  time ./a.out >i.txt
  time ./a.out >e.txt
 
 
 
 
 
 
  convert -limit memory 1000mb -limit disk 1gb dd30010000_20_3_0.90.pgm -resize 2000x2000 10.png

 
 
 
*/

#include <stdio.h>
#include <stdlib.h> // malloc
#include <string.h> // strcat
#include <math.h> // M_PI; needs -lm also
#include <complex.h>
#include <omp.h> // OpenMP
#include <limits.h> // Maximum value for an unsigned long long int



// https://sourceforge.net/p/predef/wiki/Standards/

#if defined(__STDC__)
#define PREDEF_STANDARD_C_1989
#if defined(__STDC_VERSION__)
#if (__STDC_VERSION__ >= 199409L)
#define PREDEF_STANDARD_C_1994
#endif
#if (__STDC_VERSION__ >= 199901L)
#define PREDEF_STANDARD_C_1999
#endif
#endif
#endif




/* --------------------------------- global variables and consts ------------------------------------------------------------ */


//FunctionType
typedef enum  {Fatou = 0, IntLSM =1 , ExtLSM = 2, LSM = 3, DEM = 4, Unknown = 5 , BD = 6, MBD = 7 , SAC = 8, DLD = 9, ND = 10 , NP= 11, POT = 12 , Blend = 13

} FunctionTypeT;
// FunctionTypeT FunctionType;

// virtual 2D array and integer ( screen) coordinate
// Indexes of array starts from 0 not 1
//unsigned int ix, iy; // var
static unsigned int ixMin = 0; // Indexes of array starts from 0 not 1
static unsigned int ixMax; //
static unsigned int iWidth; // horizontal dimension of array

static unsigned int iyMin = 0; // Indexes of array starts from 0 not 1
static unsigned int iyMax; //

static unsigned int iHeight = 1000; // 
// The size of array has to be a positive constant integer
static unsigned long long int iSize; // = iWidth*iHeight;

// memmory 1D array
unsigned char *data;
unsigned char *edge;
//unsigned char *edge2;

// unsigned int i; // var = index of 1D array
//static unsigned int iMin = 0; // Indexes of array starts from 0 not 1
static unsigned int iMax; // = i2Dsize-1  =
// The size of array has to be a positive constant integer
// unsigned int i1Dsize ; // = i2Dsize  = (iMax -iMin + 1) =  ;  1D array with the same size as 2D array



// see SetPlane

double radius = 3.0;
complex double center = 0.0 ;
double  DisplayAspectRatio  = 1.0; // https://en.wikipedia.org/wiki/Aspect_ratio_(image)
// dx = dy compare setup : iWidth = iHeight;
double ZxMin; //= -1.3; //-0.05;
double ZxMax;// = 1.3; //0.75;
double ZyMin;// = -1.3; //-0.1;
double ZyMax;// = 1.3; //0.7;
double PixelWidth; // =(ZxMax-ZxMin)/ixMax;
double PixelHeight; // =(ZyMax-ZyMin)/iyMax;

// dem
double BoundaryWidth ; //= 1.0*iWidth/2000.0  ; //  measured in pixels ( when iWidth = 2000)
double distanceMax ; //= BoundaryWidth*PixelWidth;


double ratio;


/*
  ER = pow(10,ERe);
  AR = pow(10,-ARe);
*/
//int ARe ; // increase ARe until black ( unknown) points disapear
//int ERe ;
double ER;
double ER2; //= 1e60;
double AR; // bigger values do not works
double AR2;

double AR_max;
//double AR12;



int IterMax = 100000;
int IterMax_LSM = 1000;
int IterMax_DEM = 100000;

/* colors = shades of gray from 0 to 255

   unsigned char colorArray[2][2]={{255,231},    {123,99}};
   color = 245;  exterior
*/
unsigned char iColorOfExterior = 245;
unsigned char iColorOfInterior = 99;
//unsigned char iColorOfInterior2 = 183;
unsigned char iColorOfBoundary = 0;
unsigned char iColorOfUnknown = 5;

// pixel counters
unsigned long long int uUnknown = 0;
unsigned long long int uInterior = 0;
unsigned long long int uExterior = 0;



// critical points




// critical points
//complex double zcr = -0.35*I; // only one critical point
//complex double zc2 = -2.2351741790771484375e-08+9.4296410679817199707e-09*I;


/*

  0.389374737779177*%i-0.9400337727919567,
        0.9400337727919567-0.389374737779177*%i,
        (-1.816005797447402*%i)-0.7522142306508819,
        1.816005797447402*%i+0.7522142306508819]

*/

const complex double zp2 = 0.389374737779177*I-0.9400337727919567; //-0.04160806403172667*I - 0.1004507524894826;

/*
a: (0.15+0.15*%i);
b: (-3+3*%i);
f: 1/(a*z^5+ b*z^3 + z);
 dz = -(5*a*z^4+3*b*z^2+1)/(a*z^5+b*z^3+z)^2

*/
const complex double a =  0.15 + 0.15*I;
const complex double b = -3.0  + 3.0*I;
//const int period = 1;

// periodic points = attractors
//complex double z1 =  0.0 ; //fixed point (period  1) =  attracting cycle



/* ------------------------------------------ functions -------------------------------------------------------------*/


// complex function
complex double f(const complex double z0) {

double complex z = z0;
z = 1/(a*z*z*z*z*z + z*z*z + b*z);
return  z;
}


//
complex double dfz(const complex double z0) {


// dz= -(5*a*z^4+3*z^2+b) / (a*z^5+z^3+b*z)^2
double complex z = z0;
complex double z2= z*z;
complex double numerator = -5.0*a*z2*z2 - 3.0*z2 - b ;
complex double denom = a*z2*z2*z + z2*z + b*z;
denom = denom*denom;


return  numerator/denom;

}









//------------------complex numbers -----------------------------------------------------





// from screen to world coordinate ; linear mapping
// uses global cons
double
GiveZx (int ix)
{
  return (ZxMin + ix * PixelWidth);
}

// uses globaal cons
double
GiveZy (int iy)
{
  return (ZyMax - iy * PixelHeight);
} // reverse y axis


complex double
GiveZ (int ix, int iy)
{
  double Zx = GiveZx (ix);
  double Zy = GiveZy (iy);

  return Zx + Zy * I;




}



double cabs2(complex double z){

  return creal(z)*creal(z)+cimag(z)*cimag(z);


}



// ****************** DYNAMICS = trap tests ( target sets) ****************************


/* -----------  array functions = drawing -------------- */

/* gives position of 2D point (ix,iy) in 1D array  ; uses also global variable iWidth */
unsigned int Give_i (unsigned int ix, unsigned int iy)
{
  return ix + iy * iWidth;
}















//
unsigned char ComputeColorOfFatou (complex double z)
{





  double r2;


  int i; // number of iteration
  for (i = 0; i < IterMax; ++i)
    {




        r2 =cabs2(z);

      if (r2 > ER2) // esaping = exterior
{
  uExterior += 1;
  return iColorOfExterior;
}


if (cabs2(zp2-z) < AR2) //
{
 
  return iColorOfInterior;
}


     
z = f(z); // complex iteration f(z)=z^3 + c

    }

 
  return iColorOfUnknown;


}









// f(z)=1+z−3z2−3.75z3+1.5z4+2.25z5
unsigned char ComputeColorOfLSM (complex double z)
{





  double r2;


  int i; // number of iteration
  for (i = 0; i < IterMax_LSM; ++i)
    {




      // complex iteration f(z)=z^3 + c
            r2 =cabs2(z);

      if (r2 > ER2) // esaping = exterior
{
  uExterior += 1;
  return 255- ((i*15) % 255);
}

       if (cabs2(zp2-z) < AR2)  // interior
{
 
  return iColorOfInterior;
}

z = f(z);

    }

  return iColorOfUnknown;


}













// ***************************************************************************************************************************
// ************************** DEM/J*****************************************
// ****************************************************************************************************************************

unsigned char ComputeColorOfDEMJ(complex double z){
  // https://en.wikibooks.org/wiki/Fractals/Iterations_in_the_complex_plane/Julia_set#DEM.2FJ


 
  int nMax = IterMax_DEM;
  complex double dz = 1.0; //  is first derivative with respect to z.
  double distance;
  double cabsz;

  int n;

  for (n=0; n < nMax; n++){ //forward iteration
    cabsz = cabs(z);
    if (cabsz > 1e60 || cabs(dz)> 1e60) break; // big values
    if (cabs2(zp2-z) < AR2)  return iColorOfInterior; // falls into finite attractor = interior
 
    dz = dfz(z)*dz;
    z = f(z) ; /* forward iteration : complex cubic polynomial */
  }
 
 
  distance = 2.0 * cabsz* log(cabsz)/ cabs(dz);
  if (distance <distanceMax) return iColorOfBoundary; // distanceMax = BoundaryWidth*PixelWidth;
  // else
 
  return iColorOfExterior;

 
}






/* ==================================================================================================
   ============================= Draw functions ===============================================================
   =====================================================================================================
*/
unsigned char ComputeColor(FunctionTypeT FunctionType, complex double z){

  unsigned char iColor;



  switch(FunctionType){
 
  case Fatou :{iColor = ComputeColorOfFatou(z); break;}
 
 
  // case IntLSM :{iColor = ComputeColorOfIntLSM(z); break;}

  // case ExtLSM :{iColor = ComputeColorOfExtLSM(z); break;}
 
  case LSM :{iColor = ComputeColorOfLSM(z); break;}
 
 

  case DEM : {iColor = ComputeColorOfDEMJ(z); break;}

/*
  case Unknown : {iColor = ComputeColorOfUnknown(z); break;}

  case BD : {iColor = ComputeColorOfBD(z); break;}

  case MBD : {iColor = ComputeColorOfMBD(z); break;}

  case SAC : {iColor = ComputeColorOfSAC(z); break;}
 
  case DLD : {iColor = ComputeColorOfDLD(z); break;}

  case ND : {iColor = ComputeColorOfND(z); break;}

  case NP : {iColor = ComputeColorOfNP(z); break;}

  case POT : {iColor = ComputeColorOfPOT(z); break;}

  case Blend : {iColor = ComputeColorOfBlend(z); break;}
  */

  default: {}


  }

  return iColor;



}


// plots raster point (ix,iy)
int DrawPoint ( unsigned char A[], FunctionTypeT FunctionType, int ix, int iy)
{
  int i; /* index of 1D array */
  unsigned char iColor;
  complex double z;


  i = Give_i (ix, iy); /* compute index of 1D array from indices of 2D array */
 
  z = GiveZ(ix,iy);
 

  iColor = ComputeColor(FunctionType, z);
  A[i] = iColor ; //
 
  return 0;
}




int DrawImage ( unsigned char A[], FunctionTypeT FunctionType)
{
  unsigned int ix, iy; // pixel coordinate

  fprintf (stdout, "compute image \n");
  // for all pixels of image
#pragma omp parallel for schedule(dynamic) private(ix,iy) shared(A, ixMax , iyMax, uUnknown, uInterior, uExterior)
  for (iy = iyMin; iy <= iyMax; ++iy)
    {
      fprintf (stdout, " %d from %d \r", iy, iyMax); //info
      for (ix = ixMin; ix <= ixMax; ++ix)
DrawPoint(A, FunctionType, ix, iy); // 
    }
  fprintf (stdout, "\n"); //info
  return 0;
}







int PlotPoint(complex double z, unsigned char A[]){


  unsigned int ix = (creal(z)-ZxMin)/PixelWidth;
  unsigned int iy = (ZyMax - cimag(z))/PixelHeight;
  unsigned int i = Give_i(ix,iy); /* index of _data array */


  A[i]= 0; //255-A[i]; // Mark point with inveres color


  return 0;

}






// ***********************************************************************************************
// ********************** edge detection usung Sobel filter ***************************************
// ***************************************************************************************************

// from Source to Destination
int ComputeBoundaries(unsigned char S[], unsigned char D[])
{
 
  unsigned int iX,iY; /* indices of 2D virtual array (image) = integer coordinate */
  unsigned int i; /* index of 1D array  */
  /* sobel filter */
  unsigned char G, Gh, Gv;
  // boundaries are in D  array ( global var )
 
  // clear D array
  memset(D, iColorOfExterior, iSize*sizeof(*D)); // for heap-allocated arrays, where N is the number of elements = FillArrayWithColor(D , iColorOfExterior);
 
  // printf(" find boundaries in S array using  Sobel filter\n");   
#pragma omp parallel for schedule(dynamic) private(i,iY,iX,Gv,Gh,G) shared(iyMax,ixMax)
  for(iY=1;iY<iyMax-1;++iY){
    for(iX=1;iX<ixMax-1;++iX){
      Gv= S[Give_i(iX-1,iY+1)] + 2*S[Give_i(iX,iY+1)] + S[Give_i(iX-1,iY+1)] - S[Give_i(iX-1,iY-1)] - 2*S[Give_i(iX-1,iY)] - S[Give_i(iX+1,iY-1)];
      Gh= S[Give_i(iX+1,iY+1)] + 2*S[Give_i(iX+1,iY)] + S[Give_i(iX-1,iY-1)] - S[Give_i(iX+1,iY-1)] - 2*S[Give_i(iX-1,iY)] - S[Give_i(iX-1,iY-1)];
      G = sqrt(Gh*Gh + Gv*Gv);
      i= Give_i(iX,iY); /* compute index of 1D array from indices of 2D array */
      if (G==0) {D[i]=255;} /* background */
      else {D[i]=0;}  /* boundary */
    }
  }
 
   
 
  return 0;
}



// copy from Source to Destination
int CopyBoundaries(unsigned char S[],  unsigned char D[])
{
 
  unsigned int iX,iY; /* indices of 2D virtual array (image) = integer coordinate */
  unsigned int i; /* index of 1D array  */
 
 
  //printf("copy boundaries from S array to D array \n");
  for(iY=1;iY<iyMax-1;++iY)
    for(iX=1;iX<ixMax-1;++iX)
      {i= Give_i(iX,iY); if (S[i]==0) D[i]=0;}
 
 
 
  return 0;
}

















// *******************************************************************************************
// ********************************** save A array to pgm file ****************************
// *********************************************************************************************

int SaveArray2PGMFile (unsigned char A[],  char * n, char *comment)
{

  FILE *fp;
  const unsigned int MaxColorComponentValue = 255; /* color component is coded from 0 to 255 ;  it is 8 bit color file */
  char name[100]; /* name of file */
  snprintf (name, sizeof name, "%s", n ); /*  */
  char *filename = strcat (name, ".pgm");
  char long_comment[200];
  sprintf (long_comment, "fc(z)=z^2+z*c  %s", comment);





  // save image array to the pgm file
  fp = fopen (filename, "wb"); // create new file,give it a name and open it in binary mode
  fprintf (fp, "P5\n # %s\n %u %u\n %u\n", long_comment, iWidth, iHeight, MaxColorComponentValue); // write header to the file
  size_t rSize = fwrite (A, sizeof(A[0]), iSize, fp); // write whole array with image data bytes to the file in one step
  fclose (fp);

  // info
  if ( rSize == iSize)
  {
  printf ("File %s saved ", filename);
  if (long_comment == NULL || strlen (long_comment) == 0)
    printf ("\n");
  else { printf (". Comment = %s \n", long_comment); }
  }
  else {printf("wrote %zu elements out of %llu requested\n", rSize,  iSize);}

  return 0;
}




int PrintCInfo ()
{

  printf ("gcc version: %d.%d.%d\n", __GNUC__, __GNUC_MINOR__, __GNUC_PATCHLEVEL__); // https://stackoverflow.com/questions/20389193/how-do-i-check-my-gcc-c-compiler-version-for-my-eclipse
  // OpenMP version is displayed in the console : export  OMP_DISPLAY_ENV="TRUE"

  printf ("__STDC__ = %d\n", __STDC__);
  printf ("__STDC_VERSION__ = %ld\n", __STDC_VERSION__);
  printf ("c dialect = ");
  switch (__STDC_VERSION__)
    { // the format YYYYMM
    case 199409L:
      printf ("C94\n");
      break;
    case 199901L:
      printf ("C99\n");
      break;
    case 201112L:
      printf ("C11\n");
      break;
    case 201710L:
      printf ("C18\n");
      break;
      //default : /* Optional */

    }

  return 0;
}


int
PrintProgramInfo ()
{


  // display info messages
  printf ("Numerical approximation of Julia set for F(z) =   z*c) \n");
  //printf ("parameter C = ( %.16f ; %.16f ) \n", creal (C), cimag (C));
 

  printf ("Image Width = %f in world coordinate\n", ZxMax - ZxMin);
  printf ("PixelWidth = %.16f \n", PixelWidth);
  //printf ("AR = %.16f = %f *PixelWidth = %f %% of ImageWidth \n", AR, AR / PixelWidth, AR / ZxMax - ZxMin);



  // image corners in world coordinate
  // center and radius
  // center and zoom
  // GradientRepetition
  printf ("Maximal number of iterations = iterMax = %d \n", IterMax);
  printf ("ratio of image  = %f ; it should be 1.000 ...\n", ratio);
  //




  return 0;
}



int SetPlane(complex double center, double radius, double a_ratio){

  ZxMin = creal(center) - radius*a_ratio;
  ZxMax = creal(center) + radius*a_ratio; //0.75;
  ZyMin = cimag(center) - radius; // inv
  ZyMax = cimag(center) + radius; //0.7;
  return 0;

}



// Check Orientation of z-plane image : mark first quadrant of complex plane
// it should be in the upper right position
// uses global var :  ...
int CheckZPlaneOrientation(unsigned char A[] )
{
 
double Zx, Zy; //  Z= Zx+ZY*i;
unsigned i; /* index of 1D array */
unsigned int ix, iy; // pixel coordinate

fprintf(stderr, "compute image CheckOrientation\n");
  // for all pixels of image
#pragma omp parallel for schedule(dynamic) private(ix,iy, i, Zx, Zy) shared(A, ixMax , iyMax)
for (iy = iyMin; iy <= iyMax; ++iy){
    //fprintf (stderr, " %d from %d \r", iy, iyMax); //info
    for (ix = ixMin; ix <= ixMax; ++ix){
    // from screen to world coordinate
    Zy = GiveZy(iy);
    Zx = GiveZx(ix);
  i = Give_i(ix, iy); /* compute index of 1D array from indices of 2D array */
  if (Zx>0 && Zy>0) A[i]=255-A[i];   // check the orientation of Z-plane by marking first quadrant */
    }
    }
   
   
  return 0;
}







// *****************************************************************************
//;;;;;;;;;;;;;;;;;;;;;;  setup ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
// **************************************************************************************

int setup ()
{

  fprintf (stderr, "setup start\n");






  /* 2D array ranges */

  iWidth = iHeight* DisplayAspectRatio ;
  iSize = iWidth * iHeight; // size = number of points in array
  // iy
  iyMax = iHeight - 1; // Indexes of array starts from 0 not 1 so the highest elements of an array is = array_name[size-1].
  //ix

  ixMax = iWidth - 1;

  /* 1D array ranges */
  // i1Dsize = i2Dsize; // 1D array with the same size as 2D array
  iMax = iSize - 1; // Indexes of array starts from 0 not 1 so the highest elements of an array is = array_name[size-1].

 
  SetPlane( center, radius,  DisplayAspectRatio );
  /* Pixel sizes */
  PixelWidth = (ZxMax - ZxMin) / ixMax; //  ixMax = (iWidth-1)  step between pixels in world coordinate
  PixelHeight = (ZyMax - ZyMin) / iyMax;
  ratio = ((ZxMax - ZxMin) / (ZyMax - ZyMin)) / ((double) iWidth / (double) iHeight); // it should be 1.000 ...

  ER = 200.0; //
  ER2 = ER*ER;
  //AR_max = 5*PixelWidth*iWidth/2000.0 ; // adjust first number
  // GiveTunedAR(const int i_Max, const complex double zcr, const double c, const double zp){
  AR = 0.003; // GiveTunedAR(200, C, AR_max);
  AR2 = AR * AR;
  //AR12 = AR/2.0;
 
 
    BoundaryWidth = 1.0*iWidth/2000.0  ; //  measured in pixels ( when iWidth = 2000)
  distanceMax = BoundaryWidth*PixelWidth;



  /* create dynamic 1D arrays for colors ( shades of gray ) */
  data = malloc (iSize * sizeof (unsigned char));

  edge = malloc (iSize * sizeof (unsigned char));
  if (data == NULL || edge == NULL)
    {
      fprintf (stderr, " Could not allocate memory");
      return 1;
    }





 


  fprintf (stderr, " end of setup \n");

  return 0;

} // ;;;;;;;;;;;;;;;;;;;;;;;;; end of the setup ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;




int end ()
{


  fprintf (stderr, " allways free memory (deallocate )  to avoid memory leaks \n"); // https://en.wikipedia.org/wiki/C_dynamic_memory_allocation
  free (data);
  free(edge);


  PrintProgramInfo ();
  PrintCInfo ();
  return 0;

}

// ********************************************************************************************************************
/* -----------------------------------------  main   -------------------------------------------------------------*/
// ********************************************************************************************************************

int main ()
{
  setup ();

   DrawImage (data, Fatou); // first
  SaveArray2PGMFile (data,  "Fatou_1" , "Fatou ");
 
  DrawImage (data, LSM); // first
  SaveArray2PGMFile (data,  "LSM_1" , "LSM ");
 
 
 
  DrawImage (data, DEM); // first
  SaveArray2PGMFile (data,  "DEM_1" , "DEM ");
 

 
  end ();

  return 0;
}



« Last Edit: June 19, 2021, 11:38:38 AM by Adam Majewski, Reason: b »

Offline Adam Majewski

  • Fractal Frankfurter
  • *
  • Posts: 541
Re: rational function
« Reply #5 on: June 17, 2021, 08:20:56 PM »

Offline Adam Majewski

  • Fractal Frankfurter
  • *
  • Posts: 541
Re: rational function
« Reply #6 on: June 17, 2021, 09:44:20 PM »
Wolfram alfa

Wolfram Alpha 

    solve{ 1/((0.15*I+0.15)*z^5 +z^3 + (3*I-3)*z)==z}
 
  z≈ ± (1.0526 + 2.5111 i)
  z≈ ± (0.17290 - 0.45180 i)
  z≈ ± (1.5420 - 0.5826 i)

So there are 6 fixed points. Stability index = abs(f'(z))

[1.028219433424695,
1.028219433424695,
31.94911214885403,
 31.94911214885403,
 235.8426231893345,
235.8426231893345]




So none of fixed point is attracting.


Period 2 points ( numerically):

period 2 points :

 [
 0.0, // 1

+ 0.1004507524894826 - 0.04160806403172667 %i,  //2
- 0.1004507524894826 + 0.04160806403172667 %i , //2



//
- 0.1972126334382917 %i - 0.4417257204225227,
+ 0.1972126334382916 %i + 0.4417257204225231,
//
- 0.172896861895406  + 0.4517976427751519 %i , // period 1
+ 0.1728968618954057 - 0.4517976427751519 %i, // period 1
//
- 1.542012716044921 + 0.5826223205446455 %i , // period 1
+ 1.542012716021183 - 0.5826223205577208 %i, // period 1
//
- 1.48680428382238  + 0.6158544989143758 %i ,
+ 1.486804283840914 - 0.6158544989534238 %i,

//
- 1.545643025068914 + 0.6402263036121009 %i ,
+ 1.545643025048946 - 0.6402263035701801 %i,

//
- 1.502343841878277 + 0.6783914544446792 %i ,
+ 1.502343841903455 - 0.678391454434491 %i,
//
- 2.008427151852299 %i - 0.8319177653355924,
+ 2.008427151854269 %i + 0.83191776533551,
//
+ 2.511067896520314 %i + 1.052616752123668, // period 1
- 2.511067897025448 %i - 1.052616732010429, // period 1
//
- 2.519905568763936 %i - 1.031280709377238,
+ 2.519905568548988 %i + 1.03128068923336,
//
  2.496853258644555 %i + 1.034230487296836,
- 2.496853270435978 %i - 1.034230487368577,
//
- 2.532331263153173 %i - 1.048925954202523,
  2.532331275662722 %i + 1.048925954304978]







Rational Julia Sets by  Mark McClure

use ( comma separated lists of coefficients in the order of ascending degree from 0 to degree of numerator/denominator ) :
1
0,  -3 +3*i,0,1, 0,0.15 + 0.15*i


or xaos user fomule : 1/((0.15*I+0.15)*z^5+(3*I-3)*z+z^3)
It seem the image is upside down in Xaos ????




« Last Edit: June 19, 2021, 03:39:56 PM by Adam Majewski, Reason: not fails »

Offline pauldelbrot

  • 3f
  • ******
  • Posts: 2994
Re: rational function
« Reply #7 on: June 19, 2021, 03:20:49 PM »
The morphology of the Julia shapes makes it clear that the finite attractors are 2-cycles, so it's not surprising you found no attracting finite fixed points.

Offline xenodreambuie

  • Fractal Friar
  • *
  • Posts: 149
Re: rational function
« Reply #8 on: June 20, 2021, 12:40:17 AM »
I did this one in Jux some years ago, with a small change to the z^5 coefficient to 0.1591+0.1591i for more circular shapes. It appears to actually be period 4, with odd cycles and even cycles in separate clusters (not affected by the coefficient change).

Offline gerson

  • Fractal Freak
  • **
  • Posts: 676
Re: rational function
« Reply #9 on: June 20, 2021, 02:42:59 AM »
@ Adam Majewski very interesting equation
@ xenodreambuie liked a lot the image

Offline Adam Majewski

  • Fractal Frankfurter
  • *
  • Posts: 541
Re: rational function
« Reply #10 on: June 20, 2021, 08:59:33 AM »
The morphology of the Julia shapes makes it clear that the finite attractors are 2-cycles, so it's not surprising you found no attracting finite fixed points.
Probably {0, infinity} is superattractive period 2 cycle.
Are ther any other attractive cycles ?
Is this Julia set connected ?



Here is my naive numerical anlysis of periodic points ( c code below). Comments are welcome
Code: [Select]
k = 0

i = 0 z = +0.000000 +0.000000 
 z = +inf -nan cabs(zn) > 3.0 = escaping to infinity or period 2 cycle = {z0, inf} or repelling period 2

k = 1

i = 0 z = +0.100451 -0.041608 
i = 1 z = -0.831918 -2.008427 cabs(z) = 2.173906 d= 2.176623
i = 2 z = +0.100451 -0.041608 cabs(z) = 0.108727 d= 2.176623
period 2 cycle probably attracting :  fabs(dp-dn) < 0.000001 = distance is not changing or slowly changing

k = 2

i = 0 z = -0.100451 +0.041608 
i = 1 z = +0.831918 +2.008427 cabs(z) = 2.173906 d= 2.176623
i = 2 z = -0.100451 +0.041608 cabs(z) = 0.108727 d= 2.176623
period 2 cycle probably attracting :  fabs(dp-dn) < 0.000001 = distance is not changing or slowly changing

k = 3

i = 0 z = -0.441726 -0.197213 
i = 1 z = +0.441726 +0.197213 cabs(z) = 0.483750 d= 0.967501
i = 2 z = -0.441726 -0.197213 cabs(z) = 0.483750 d= 0.967501
period 2 cycle probably attracting :  fabs(dp-dn) < 0.000001 = distance is not changing or slowly changing

k = 4

i = 0 z = +0.441726 +0.197213 
i = 1 z = -0.441726 -0.197213 cabs(z) = 0.483750 d= 0.967501
i = 2 z = +0.441726 +0.197213 cabs(z) = 0.483750 d= 0.967501
period 2 cycle probably attracting :  fabs(dp-dn) < 0.000001 = distance is not changing or slowly changing

k = 5

i = 0 z = -0.172897 +0.451798 
i = 1 z = -0.172897 +0.451798 cabs(z) = 0.483750 d= 0.000000
i = 2 z = -0.172897 +0.451798 cabs(z) = 0.483750 d= 0.000000
 z = -0.172897 +0.451798   dp < 0.000001 = probably period 1

k = 6

i = 0 z = +0.172897 -0.451798 
i = 1 z = +0.172897 -0.451798 cabs(z) = 0.483750 d= 0.000000
i = 2 z = +0.172897 -0.451798 cabs(z) = 0.483750 d= 0.000000
 z = +0.172897 -0.451798   dp < 0.000001 = probably period 1

k = 7

i = 0 z = -1.542013 +0.582622 
i = 1 z = -1.542013 +0.582622 cabs(z) = 1.648409 d= 0.000000
i = 2 z = -1.542013 +0.582622 cabs(z) = 1.648409 d= 0.000000
 z = -1.542013 +0.582623   dp < 0.000001 = probably period 1

k = 8

i = 0 z = +1.542013 -0.582622 
i = 1 z = +1.542013 -0.582622 cabs(z) = 1.648409 d= 0.000000
i = 2 z = +1.542013 -0.582622 cabs(z) = 1.648409 d= 0.000000
 z = +1.542012 -0.582623   dp < 0.000001 = probably period 1

k = 9

i = 0 z = -1.486804 +0.615854 
i = 1 z = +1.034230 +2.496853 cabs(z) = 2.702575 d= 3.145437
i = 2 z = -1.486804 +0.615854 cabs(z) = 1.609305 d= 3.145437
i = 3 z = +1.034215 +2.496862 cabs(z) = 2.702577 d= 3.145430
i = 4 z = -1.487475 +0.614659 cabs(z) = 1.609468 d= 3.146682
i = 5 z = +0.938842 +2.545302 cabs(z) = 2.712930 d= 3.100709
i = 6 z = -0.171898 -0.241483 cabs(z) = 0.296417 d= 2.999986
i = 7 z = +0.771033 -0.122995 cabs(z) = 0.780781 d= 0.950346
i = 8 z = -0.177912 -0.300455 cabs(z) = 0.349179 d= 0.965395
i = 9 z = +0.637916 -0.158362 cabs(z) = 0.657279 d= 0.828110
unknown
 or if distance is constant probably attracting period 2 cycle
 or repelling period 2 cycle if distance is changing 

k = 10

i = 0 z = +1.486804 -0.615854 
i = 1 z = -1.034230 -2.496853 cabs(z) = 2.702575 d= 3.145437
i = 2 z = +1.486805 -0.615855 cabs(z) = 1.609306 d= 3.145437
i = 3 z = -1.034234 -2.496871 cabs(z) = 2.702592 d= 3.145450
i = 4 z = +1.488172 -0.616124 cabs(z) = 1.610673 d= 3.146385
i = 5 z = -1.055143 -2.607916 cabs(z) = 2.813281 d= 3.230432
i = 6 z = -0.239382 +0.188797 cabs(z) = 0.304873 d= 2.913257
i = 7 z = +0.101878 +0.781313 cabs(z) = 0.787927 d= 0.683765
i = 8 z = -0.203446 +0.173630 cabs(z) = 0.267466 d= 0.680075
i = 9 z = +0.079584 +0.889911 cabs(z) = 0.893462 d= 0.770171
unknown
 or if distance is constant probably attracting period 2 cycle
 or repelling period 2 cycle if distance is changing 

k = 11

i = 0 z = -1.545643 +0.640226 
i = 1 z = -1.048926 -2.532331 cabs(z) = 2.740976 d= 3.211207
i = 2 z = -1.545643 +0.640227 cabs(z) = 1.672992 d= 3.211207
i = 3 z = -1.048899 -2.532346 cabs(z) = 2.740979 d= 3.211226
i = 4 z = -1.544285 +0.642759 cabs(z) = 1.672709 d= 3.213518
i = 5 z = -0.799562 -2.636953 cabs(z) = 2.755507 d= 3.363201
i = 6 z = +0.064716 +0.091997 cabs(z) = 0.112479 d= 2.862541
i = 7 z = -2.059373 +0.355826 cabs(z) = 2.089888 d= 2.140411
i = 8 z = -0.085504 +0.016463 cabs(z) = 0.087074 d= 2.002830
i = 9 z = +1.518575 +2.246207 cabs(z) = 2.711368 d= 2.746785
unknown
 or if distance is constant probably attracting period 2 cycle
 or repelling period 2 cycle if distance is changing 

k = 12

i = 0 z = +1.545643 -0.640226 
i = 1 z = +1.048926 +2.532331 cabs(z) = 2.740976 d= 3.211207
i = 2 z = +1.545643 -0.640226 cabs(z) = 1.672992 d= 3.211207
i = 3 z = +1.048933 +2.532362 cabs(z) = 2.741007 d= 3.211236
i = 4 z = +1.542793 -0.639610 cabs(z) = 1.670123 d= 3.210188
 z = +1.107056 +2.828193 cabs(zn) > 3.0 = escaping to infinity or period 2 cycle = {z0, inf} or repelling period 2

k = 13

i = 0 z = -1.502344 +0.678391 
i = 1 z = +1.502344 -0.678391 cabs(z) = 1.648409 d= 3.296818
i = 2 z = -1.502344 +0.678391 cabs(z) = 1.648409 d= 3.296818
period 2 cycle probably attracting :  fabs(dp-dn) < 0.000001 = distance is not changing or slowly changing

k = 14

i = 0 z = +1.502344 -0.678391 
i = 1 z = -1.502344 +0.678391 cabs(z) = 1.648409 d= 3.296818
i = 2 z = +1.502344 -0.678391 cabs(z) = 1.648409 d= 3.296818
period 2 cycle probably attracting :  fabs(dp-dn) < 0.000001 = distance is not changing or slowly changing

k = 15

i = 0 z = -0.831918 -2.008427 
i = 1 z = +0.100451 -0.041608 cabs(z) = 0.108727 d= 2.176623
i = 2 z = -0.831918 -2.008427 cabs(z) = 2.173906 d= 2.176623
period 2 cycle probably attracting :  fabs(dp-dn) < 0.000001 = distance is not changing or slowly changing

k = 16

i = 0 z = +0.831918 +2.008427 
i = 1 z = -0.100451 +0.041608 cabs(z) = 0.108727 d= 2.176623
i = 2 z = +0.831918 +2.008427 cabs(z) = 2.173906 d= 2.176623
period 2 cycle probably attracting :  fabs(dp-dn) < 0.000001 = distance is not changing or slowly changing

k = 17

i = 0 z = +1.052617 +2.511068 
i = 1 z = +1.052616 +2.511064 cabs(z) = 2.722764 d= 0.000004
i = 2 z = +1.051601 +2.511174 cabs(z) = 2.722473 d= 0.001022
i = 3 z = +1.082008 +2.772741 cabs(z) = 2.976379 d= 0.263328
i = 4 z = +0.070973 -0.059033 cabs(z) = 0.092315 d= 3.006848
i = 5 z = -0.236804 -2.546474 cabs(z) = 2.557461 d= 2.506409
i = 6 z = +0.049184 -0.002415 cabs(z) = 0.049243 d= 2.560084
 z = -3.214658 -3.549163 cabs(zn) > 3.0 = escaping to infinity or period 2 cycle = {z0, inf} or repelling period 2

k = 18

i = 0 z = -1.052617 -2.511068 
i = 1 z = -1.052617 -2.511068 cabs(z) = 2.722768 d= 0.000000
i = 2 z = -1.052715 -2.511071 cabs(z) = 2.722809 d= 0.000098
i = 3 z = -1.052554 -2.488103 cabs(z) = 2.701579 d= 0.022969
i = 4 z = +0.462528 -0.999277 cabs(z) = 1.101130 d= 2.124165
i = 5 z = +0.014954 -0.195738 cabs(z) = 0.196309 d= 0.919781
i = 6 z = +0.771090 -0.911883 cabs(z) = 1.194198 d= 1.041444
i = 7 z = -0.070214 -0.203606 cabs(z) = 0.215373 d= 1.099749
i = 8 z = +0.972784 -0.475440 cabs(z) = 1.082752 d= 1.077839
i = 9 z = -0.130633 -0.302244 cabs(z) = 0.329266 d= 1.116926
unknown
 or if distance is constant probably attracting period 2 cycle
 or repelling period 2 cycle if distance is changing 

k = 19

i = 0 z = -1.031281 -2.519906 
i = 1 z = +1.031281 +2.519906 cabs(z) = 2.722768 d= 5.445536
i = 2 z = -1.031188 -2.519888 cabs(z) = 2.722717 d= 5.445485
i = 3 z = +1.035915 +2.498359 cabs(z) = 2.704610 d= 5.427312
i = 4 z = -1.610454 +0.770400 cabs(z) = 1.785239 d= 3.160555
i = 5 z = +0.134719 -0.407770 cabs(z) = 0.429448 d= 2.105639
i = 6 z = +0.223892 -0.495799 cabs(z) = 0.544008 d= 0.125304
i = 7 z = +0.124084 -0.414032 cabs(z) = 0.432226 d= 0.129026
i = 8 z = +0.235618 -0.484984 cabs(z) = 0.539190 d= 0.132190
i = 9 z = +0.114235 -0.422698 cabs(z) = 0.437862 d= 0.136431
unknown
 or if distance is constant probably attracting period 2 cycle
 or repelling period 2 cycle if distance is changing 

k = 20

i = 0 z = +1.031281 +2.519906 
i = 1 z = -1.031281 -2.519901 cabs(z) = 2.722764 d= 5.445531
i = 2 z = +1.032305 +2.519975 cabs(z) = 2.723221 d= 5.445984
i = 3 z = -0.993050 -2.779869 cabs(z) = 2.951919 d= 5.673660
i = 4 z = -0.040063 +0.082568 cabs(z) = 0.091774 d= 3.016908
i = 5 z = -0.836507 +2.428865 cabs(z) = 2.568877 d= 2.477788
i = 6 z = -0.004093 +0.026235 cabs(z) = 0.026553 d= 2.542743
 z = -5.232646 +7.169638 cabs(zn) > 3.0 = escaping to infinity or period 2 cycle = {z0, inf} or repelling period 2

k = 21

i = 0 z = +1.034230 +2.496853 
i = 1 z = -1.486803 +0.615854 cabs(z) = 1.609305 d= 3.145436
i = 2 z = +1.034230 +2.496787 cabs(z) = 2.702513 d= 3.145396
i = 3 z = -1.481624 +0.615777 cabs(z) = 1.604491 d= 3.141293
i = 4 z = +1.011372 +2.146284 cabs(z) = 2.372637 d= 2.925317
i = 5 z = -0.106695 +0.068182 cabs(z) = 0.126620 d= 2.359785
i = 6 z = +0.404193 +1.823882 cabs(z) = 1.868133 d= 1.828522
i = 7 z = -0.085208 +0.044710 cabs(z) = 0.096226 d= 1.845256
i = 8 z = +0.731458 +2.343232 cabs(z) = 2.454744 d= 2.439293
i = 9 z = -0.120305 -0.002622 cabs(z) = 0.120334 d= 2.495703
unknown
 or if distance is constant probably attracting period 2 cycle
 or repelling period 2 cycle if distance is changing 

k = 22

i = 0 z = -1.034230 -2.496853 
i = 1 z = +1.486804 -0.615854 cabs(z) = 1.609305 d= 3.145437
i = 2 z = -1.034230 -2.496860 cabs(z) = 2.702580 d= 3.145440
i = 3 z = +1.487300 -0.615818 cabs(z) = 1.609749 d= 3.145860
i = 4 z = -1.031051 -2.535870 cabs(z) = 2.737463 d= 3.166811
i = 5 z = -0.369923 +1.279341 cabs(z) = 1.331749 d= 3.872070
i = 6 z = -0.032201 +0.134449 cabs(z) = 0.138252 d= 1.193664
i = 7 z = -0.884125 +1.454884 cabs(z) = 1.702458 d= 1.571408
i = 8 z = +0.039601 +0.087176 cabs(z) = 0.095749 d= 1.650423
i = 9 z = -2.299897 +0.862691 cabs(z) = 2.456372 d= 2.464685
unknown
 or if distance is constant probably attracting period 2 cycle
 or repelling period 2 cycle if distance is changing 

k = 23

i = 0 z = -1.048926 -2.532331 
i = 1 z = -1.545643 +0.640226 cabs(z) = 1.672992 d= 3.211207
i = 2 z = -1.048927 -2.532322 cabs(z) = 2.740968 d= 3.211198
i = 3 z = -1.546480 +0.640104 cabs(z) = 1.673719 d= 3.211207
i = 4 z = -1.058996 -2.455627 cabs(z) = 2.674243 d= 3.133879
i = 5 z = +0.278573 -0.455770 cabs(z) = 0.534162 d= 2.405934
i = 6 z = +0.075496 -0.441610 cabs(z) = 0.448017 d= 0.203570
i = 7 z = +0.280187 -0.431527 cabs(z) = 0.514510 d= 0.204939
i = 8 z = +0.068681 -0.461364 cabs(z) = 0.466448 d= 0.213600
i = 9 z = +0.276134 -0.407608 cabs(z) = 0.492335 d= 0.214304
unknown
 or if distance is constant probably attracting period 2 cycle
 or repelling period 2 cycle if distance is changing 

k = 24

i = 0 z = +1.048926 +2.532331 
i = 1 z = +1.545642 -0.640226 cabs(z) = 1.672991 d= 3.211207
i = 2 z = +1.048928 +2.532431 cabs(z) = 2.741068 d= 3.211305
i = 3 z = +1.536481 -0.640006 cabs(z) = 1.664446 d= 3.209682
 z = +0.900158 +3.704194 cabs(zn) > 3.0 = escaping to infinity or period 2 cycle = {z0, inf} or repelling period 2


here is c src code
Code: [Select]
#include <stdio.h>
#include <stdlib.h> // malloc
#include <math.h> // M_PI; needs -lm also
#include <complex.h>
/*


a: (0.15+0.15**I);
b: (-3+3**I);

define(f(z), 1/(a*z^5+ z^3 + b*z));



period 2 points :
*/

double complex zp2[25] = {
 0.0, // period 2 = found by analysis of the function
//
+ 0.1004507524894826 - 0.04160806403172667 *I,  //2
-0.1004507524894826 + 0.04160806403172667 *I , //2
//
- 0.1972126334382917 *I - 0.4417257204225227, //2
+ 0.1972126334382916 *I + 0.4417257204225231,
//
- 0.172896861895406  + 0.4517976427751519 *I , // period 1
+ 0.1728968618954057 - 0.4517976427751519 *I, // period 1
//
- 1.542012716044921 + 0.5826223205446455 *I , // period 1
+ 1.542012716021183 - 0.5826223205577208 *I, // period 1
//
- 1.48680428382238  + 0.6158544989143758 *I , // probably repelling period 2
+ 1.486804283840914 - 0.6158544989534238 *I, // probably repelling period 2

//
- 1.545643025068914 + 0.6402263036121009 *I ,  // probably repelling period 2
+ 1.545643025048946 - 0.6402263035701801 *I,  // probably repelling period 2

//
- 1.502343841878277 + 0.6783914544446792 *I , // period 2
+ 1.502343841903455 - 0.678391454434491 *I,
//
- 2.008427151852299 *I - 0.8319177653355924, //period 2
+ 2.008427151854269 *I + 0.83191776533551,
//
+ 2.511067896520314 *I + 1.052616752123668, // period 1
- 2.511067897025448 *I - 1.052616732010429, // period 1
//
- 2.519905568763936 *I - 1.031280709377238,
+ 2.519905568548988 *I + 1.03128068923336,
//
  2.496853258644555 *I + 1.034230487296836, // ??
- 2.496853270435978 *I - 1.034230487368577,
//
- 2.532331263153173 *I - 1.048925954202523,
  2.532331275662722 *I + 1.048925954304978
};




complex double a = 0.15+0.15*I;
complex double b = -3.0+3.0*I;

complex double f(complex double z){

complex double z3 = z*z*z;
complex double z5 = z3*z*z;

return  1.0/(a*z5+ z3 + b*z);
}



void TestZ(complex double z0){



complex double zm = z0;
complex double zn;
double dp;
double dn;

int i=0;



printf("i = %d \t z = %+f %+f  \n",  i, creal(zm), cimag(zm) );


for (i = 1; i < 10; ++i)
  {
   
    zn = f(zm);
    if (cabs(zn) > 3.0){printf(" z = %+f %+f \t cabs(zn) > 3.0 = escaping to infinity or period 2 cycle = {z0, inf} or repelling period 2 \n\n", creal (zn), cimag(zn)); return;}
    dn = cabs(zm - zn);
   
    if (i>2) {
   
    if ( dn > 6.0 ){ printf(" dp>6.0\n"); break;}
    if ( dn < 0.000001 ) { printf(" z = %+f %+f \t  dp < 0.000001 = probably period 1 \n\n", creal (zn), cimag(zn)); return;}
    if ( fabs(dp-dn) < 0.000001 ) { printf("period 2 cycle probably attracting :  fabs(dp-dn) < 0.000001 = distance is not changing or slowly changing \n\n"); return;}
    }
   
    //
    printf("i = %d \t z = %+f %+f \t cabs(z) = %f \t d= %f\n", i, creal(zn), cimag(zn), cabs(zn), dn );
    zm= zn;
    dp = dn;
   
  }


printf("unknown \n or if distance is constant probably attracting period 2 cycle \n or repelling period 2 cycle if distance is changing  \n\n");




}



int main (void){

int k;
int kMax = sizeof(zp2)/sizeof(zp2[0]);


for (k = 0; k < kMax; ++k)
{printf("k = %d\n\n",k); TestZ(zp2[k]);}

return 0;
}
« Last Edit: June 20, 2021, 04:06:17 PM by Adam Majewski, Reason: code »

Offline Adam Majewski

  • Fractal Frankfurter
  • *
  • Posts: 541
Re: rational function
« Reply #11 on: June 20, 2021, 09:13:35 PM »
the best image

Offline xenodreambuie

  • Fractal Friar
  • *
  • Posts: 149
Re: rational function
« Reply #12 on: June 20, 2021, 11:35:17 PM »
Jux detects this as two basins each with period 4. Here they are coloured with warm colours in one and cool in the other. That seems to be at odds with the analysis of period 2 points by everyone else. Any ideas as to what might cause this disparity?

Offline Chris_M_Thomasson

  • Fractal Feline
  • **
  • Posts: 154
  • Programmer by trade; love fractals and math! :^)
    • Fractal Life 247
Re: rational function
« Reply #13 on: June 21, 2021, 01:10:30 AM »
Wolfram alfa

Wolfram Alpha 

[...]

or xaos user fomule : 1/((0.15*I+0.15)*z^5+(3*I-3)*z+z^3)
It seem the image is upside down in Xaos ????

Back when I used Xaos, it always erroneously inverted the y axis where +y pointed down, and -y pointed up. This is wrong. Btw, I found some interesting bugs in Xaos:

https://groups.google.com/g/xaos-devel/c/fuFN3MJteoc

https://element90.wordpress.com/2013/10/09/another-unexpected-difference/
It's a Fractal Life, 247... ;^)

Offline Adam Majewski

  • Fractal Frankfurter
  • *
  • Posts: 541
Re: rational function
« Reply #14 on: June 21, 2021, 02:50:50 PM »
Critical points are centers of biggest components ( gray dots on a_Fatou_cr.png)

Periodic points ( without zere and infinity) are on the boundary of componnets ( gray black  dots on a_Fatou-zp2.png are hard to seee, sorry) so this is parabolic case.

Zero and infinity is a period 2 cycle in the exterior of Julia set, so here there is no  gap in the middle. See that components also end on 2 periodic points near zero.


« Last Edit: June 22, 2021, 03:08:15 PM by Adam Majewski, Reason: spell »


xx
Definitions of degree of rational function

Started by xenodreambuie on Fractal Mathematics And New Theories

3 Replies
339 Views
Last post December 06, 2020, 07:30:43 AM
by Adam Majewski
question
How to convert a function to a recursive function

Started by Deliberate Dendrite on Noob's Corner

2 Replies
481 Views
Last post December 13, 2018, 05:32:30 PM
by Deliberate Dendrite
xx
Rounding function

Started by mclarekin on Color Snippets

0 Replies
241 Views
Last post March 26, 2019, 01:46:13 AM
by mclarekin
clip
Lambert W function

Started by superheal on Share a fractal

3 Replies
808 Views
Last post December 06, 2018, 05:38:25 AM
by FractalDave
clip
Programming a recursive function

Started by Fraktalist on Fractal Humor

0 Replies
1105 Views
Last post March 21, 2018, 11:48:01 AM
by Fraktalist