File:Discrete dynamic in case of complex quadratic polynomial.gif
No higher resolution available.
Discrete_dynamic_in_case_of_complex_quadratic_polynomial.gif (400 × 300 pixels, file size: 307 KB, MIME type: image/gif, looped, 55 frames, 28 s)
 |
This is a file from the Wikimedia Commons |
Summary
| Description |
English: discrete dynamic in case of complex quadratic polynomial fc(z)=z^2 +c where c = -0.75 is a parabolic fixed point. Exterior of filled Julia set is white, known interior is gray, unknown pixels are red. Max Iteration number is in the upper left corner of the image. |
| Date | |
| Source | Own program[1] which uses code by Wolf Jung[2] and Claude Heiland-Allen[3] ( thx also for help with OMP code ) |
| Author | Adam majewski |
Long description
What can be seen ?
- in every image there is a number in the left upper corner of the image. It is a maximal number of the iteratrions ( maxiter variable in the c src code ), which was used to generate this image
- on first image ( maxiter=0 so no iterations was made ) one can see that :
- almost all pixels are red ( because : complex double zcenter = 0.0; double radius = 1.3; so all corners are inside cirle of center z=0 and radius= 2 )
- only a few pixels are gray. This is a small neighbourhood of alfa fixed point ( d2Min = pixel_spacing * 15; d2 = GiveDistance(z, z_alfa); )
- on images 1-400 escaping set is growing and attracting set is the same. The Filled Julia set ( now red ) is more and more detailed. On the rest of images ( >400) the escaping set also grows, but it can hardly be seen. It is because of Level Set method used to find escaping set
- on images > 400 attracting set is growing and fill ( almost) all the filled Julia set
- on images 500- 700 attracting petals (dark gray , two around parabolic fixed point alfa and its preimages around preimages of parabolic fixed point ) can be seen
- on images 718-1200 repelling petals can be seen ( red, two around parabolic fixed point alfa and its preimages around preimages of parabolic fixed point ) can be seen
- on the last image allmost all interior of filled Julia set is dark gray. Only few red points are seen ( unknown). It is because slow dynamics around parabolic fixed point.[4]
Algorithm
Making animated gif
- choose parabolic parameter c
- compute parabolic parameter c ( see GiveC function in c code )
- compute parabolic fixed point ( see GiveAlfaFixedPoint function in c code )
- for max iteration ( maxmaxiter ) from 0 to 1200 ( see main function in the c code ) compute and save ppm images ( render and image_save_ppm functions)
- convert ppm images to gif and add Max Iteration number in the upper left corner of the image ( see Bash and image Magic src code )
- convert series of gif static images into one animated gif ( see Bash and image Magic src code )
Color of pixel is computed with this pseudocede :
if (abs(z) > escape_radius) component= iExterior;
else { if (abs(z) < dMin) component= iInterior;
else component= iUnknown; }
Components
From math point of view there are :
Fatou set consist of components ( subsets) :
- attracting components ( here exterior of Julia set )
- parabolic compoent ) here interior of Julia set )
From programming point of view ( because of numerical methods and limited precision ) there are 3 components of image :
- known exterior = escaping set, all pixels that escape to infinity
- known interior = attracting set, all pixels which are attracted by alfa fixed point after limited iterations ( Max iter )
- unknown pixels, all pixels which do not escape or not attracted by alfa fixed point after limited iterations ( Max iter )
| Sets | color | description |
|---|---|---|
| escaping set | white ( light gray ) | |
| attracting set | dark gray | |
| unknown | red | |
| Julia set | black | not seen on this image |
C src code
/*
code from :
http://mathr.co.uk/blog/2014-11-02_practical_interior_distance_rendering.html
and by other people ;
see description of the functions
license GPL3+ <http://www.gnu.org/licenses/gpl.html>
---------- how to use ? -----------------------------------------------
to compile with gcc :
gcc -std=c99 -Wall -Wextra -pedantic -O3 -fopenmp d.c -lm
./a.out
gcc -std=c99 -Wall -Wextra -pedantic -O3 -fopenmp d.c -lm
time ./a.out
Text file saved.
real 19m4.124s
user 58m2.001s
sys 0m7.987s
------------- gnuplot --------------------
set terminal postscript portrait enhanced mono dashed lw 1 "Helvetica" 14
set output "data.ps"
set title "Points on dynamical plane "
set ylabel "ratio"
set xlabel "iteratio max"
set xrange [0:150]
plot "data.txt" using 1:2 title 'exterior' lc 1, "data.txt" using 1:3 title 'interior' lc 2, "data.txt" using 1:4 title 'unknown' lc 3
*/
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <complex.h>
#include <omp.h> // OpenMP; needs also -fopenmp
/* --------------------- global variables --------------- */
const double pi = 3.14159265358979323846264338327950288419716939937510;
const double escape_radius = 2.0 ;
double escape_radius_2 ; // = escape_radius * escape_radius
double d2Min; // minimal distance from z to z_alfa
/* components
white = exterior of Julia set
black = boundary of Julia set ( using DEM )
red = glitches
yellow = known interior of Julia set
blue = unknown
*/
int iExterior = 0;
int iInterior = 1;
int iUnknown = 2;
// Number of pixels in the image
int ExteriorPixels = 0;
int InteriorPixels = 0;
int UnknownPixels = 0;
double AllPixels;
/* ------------------ functions --------------------------*/
/* find c in component of Mandelbrot set
uses code by Wolf Jung from program Mandel
see function mndlbrot::bifurcate from mandelbrot.cpp
http://www.mndynamics.com/indexp.html
*/
double complex GiveC(unsigned int numerator, unsigned int denominator, double InternalRadius, unsigned int ParentPeriod)
{
double t = 2.0*pi*numerator/denominator; // Internal Angle from turns to radians
double R2 = InternalRadius * InternalRadius;
double Cx, Cy; /* C = Cx+Cy*i */
switch ( ParentPeriod ) // of component
{
case 1: // main cardioid
Cx = (cos(t)*InternalRadius)/2-(cos(2*t)*R2)/4;
Cy = (sin(t)*InternalRadius)/2-(sin(2*t)*R2)/4;
break;
case 2: // only one component
Cx = InternalRadius * 0.25*cos(t) - 1.0;
Cy = InternalRadius * 0.25*sin(t);
break;
// for each ParentPeriod there are 2^(ParentPeriod-1) roots.
default: // higher periods : not works, to do
Cx = 0.0;
Cy = 0.0;
break; }
return Cx + Cy*I;
}
/*
http://en.wikipedia.org/wiki/Periodic_points_of_complex_quadratic_mappings
z^2 + c = z
z^2 - z + c = 0
ax^2 +bx + c =0 // ge3neral for of quadratic equation
so :
a=1
b =-1
c = c
so :
The discriminant is the d=b^2- 4ac
d=1-4c = dx+dy*i
r(d)=sqrt(dx^2 + dy^2)
sqrt(d) = sqrt((r+dx)/2)+-sqrt((r-dx)/2)*i = sx +- sy*i
x1=(1+sqrt(d))/2 = beta = (1+sx+sy*i)/2
x2=(1-sqrt(d))/2 = alfa = (1-sx -sy*i)/2
alfa : attracting when c is in main cardioid of Mandelbrot set, then it is in interior of Filled-in Julia set,
it means belongs to Fatou set ( strictly to basin of attraction of finite fixed point )
*/
// uses global variables :
// ax, ay (output = alfa(c))
double complex GiveAlfaFixedPoint(double complex c)
{
double dx, dy; //The discriminant is the d=b^2- 4ac = dx+dy*i
double r; // r(d)=sqrt(dx^2 + dy^2)
double sx, sy; // s = sqrt(d) = sqrt((r+dx)/2)+-sqrt((r-dx)/2)*i = sx + sy*i
double ax, ay;
// d=1-4c = dx+dy*i
dx = 1 - 4*creal(c);
dy = -4 * cimag(c);
// r(d)=sqrt(dx^2 + dy^2)
r = sqrt(dx*dx + dy*dy);
//sqrt(d) = s =sx +sy*i
sx = sqrt((r+dx)/2);
sy = sqrt((r-dx)/2);
// alfa = ax +ay*i = (1-sqrt(d))/2 = (1-sx + sy*i)/2
ax = 0.5 - sx/2.0;
ay = sy/2.0;
return ax+ay*I;
}
// https://en.wikipedia.org/wiki/Euclidean_distance
double GiveDistance( complex double z1, complex double z2)
{
double dx,dy;
dx = creal(z1)-creal(z2);
dy = cimag(z1)-cimag(z2);
return (sqrt(dx*dx +dy*dy ));
}
/*
code from :
http://mathr.co.uk/blog/2014-11-02_practical_interior_distance_rendering.html
*/
static inline double cabs2(complex double z) {
return creal(z) * creal(z) + cimag(z) * cimag(z);
}
static inline unsigned char *image_new(int width, int height) {
return malloc(width * height * 3);
}
static inline void image_delete(unsigned char *image) {
free(image);
}
static inline void image_save_ppm(unsigned char *image, int width, int height, const char *filename) {
FILE *f = fopen(filename, "wb");
if (f) {
fprintf(f, "P6\n%d %d\n255\n", width, height);
fwrite(image, width * height * 3, 1, f);
fclose(f);
} else {
fprintf(stderr, "ERROR saving `%s'\n", filename);
}
}
static inline void image_poke(unsigned char *image, int width, int i, int j, int r, int g, int b) {
int k = (width * j + i) * 3;
image[k++] = r;
image[k++] = g;
image[k ] = b;
}
static inline void colour_pixel(unsigned char *image, int width, int i, int j, int component) {
int r, g, b;
// color depends on component
if (component == iExterior ) { r = 245; g = 245; b = 245; } // light gray
if (component == iInterior ) { r = 180; g = 180; b = 180; } // dark gray
if (component == iUnknown ) { r = 200; g = 0; b = 0; } // red
image_poke(image, width, i, j, r, g, b);
}
static inline void render(unsigned char *image, int maxiters, int width, int height,complex double c, complex double center, double radius, complex double z_alfa) {
double pixel_spacing = radius / (height / 2.0);
d2Min = pixel_spacing * 15;
#pragma omp parallel for ordered schedule(auto) shared(ExteriorPixels, InteriorPixels, UnknownPixels)
for (int j = 0; j < height; ++j) {
for (int i = 0; i < width; ++i) {
double x = i + 0.5 - width / 2.0;
double y = height / 2.0 - j - 0.5;
complex double z = center + pixel_spacing * (x + I * y);
int component ;
double r2=0.0; // r = cabs(z) = radius or abs value of complex number
double d2 = GiveDistance(z, z_alfa); // d = distance( z, z_alfa)
if (d2 > d2Min)
// iteration
for (int n = 1; n <= maxiters; ++n) {
z = z * z + c;
r2 = cabs2(z); // r = cabs(z)
d2 = GiveDistance(z, z_alfa); // d = distance( z, z_alfa)
if (r2 > escape_radius_2) break;
if (d2 < d2Min) break;
}
if (r2 > escape_radius_2) {component= iExterior;
#pragma omp atomic
ExteriorPixels++; }
else { if (d2 < d2Min) { component= iInterior;
#pragma omp atomic
InteriorPixels++;}
else {component= iUnknown;
#pragma omp atomic
UnknownPixels++;} }
colour_pixel(image, width, i, j, component);
}
}
}
int main() {
int maxmaxiter =1200; //maximal number of iterations : z(i+1) = fc( zi)
// parameter of the function fc(z)= z*z +c
complex double c = -0.75;
// image :
// integer sizes of the image
int width = 2000;
int height = 1500;
// description of the image : center and radius
complex double zcenter = 0.0;
double radius = 1.3;
unsigned int length = 30;
// text file name for data about images
// compare with http://mathr.co.uk/blog/2014-11-02_practical_interior_distance_rendering.html
AllPixels= width*height;
FILE * fp;
fp = fopen("data.txt","wb"); /*create new file,give it a name and open it in binary mode */
fprintf(fp,"# maxiter ExteriorPixels InteriorPixels UnknownPixesl\n"); /*writ`e header to the file*/
complex double z_alfa; // alfa fixed point alfa = f(alfa)
// setup
z_alfa = GiveAlfaFixedPoint(c);
escape_radius_2 = escape_radius * escape_radius ;
for (int maxiter = 0; maxiter <= maxmaxiter; ++maxiter) {
// filename from maxiter
char filename[length] ;
snprintf(filename, length, "%d.ppm", maxiter);
// setup
unsigned char *image = image_new(width, height);
// render image and save
render(image, maxiter, width, height, c, zcenter, radius, z_alfa);
image_save_ppm(image, width, height, filename);
// free image
image_delete(image);
// save info about image
fprintf(fp," %d %f %f %f\n", maxiter, ExteriorPixels/AllPixels, InteriorPixels/AllPixels, UnknownPixels/AllPixels);
// printf(" %d ; %d ; %d ; %d ; %d ; %.0f\n", maxiter,ExteriorPixels, InteriorPixels, UnknownPixels, ExteriorPixels+ InteriorPixels+ UnknownPixels, AllPixels);
ExteriorPixels = 0;
InteriorPixels = 0;
UnknownPixels = 0;
}
printf("Text file saved. \n");
fclose(fp);
return 0;
}
Bash, Image Magic anf gifsicle code
#!/bin/bash
# script file for BASH
# which bash
# save this file as g.sh
# chmod +x g.sh
# ./g.sh
# for all ppm files in this directory
for file in *.ppm ; do
# b is name of file without extension
b=$(basename $file .ppm)
# convert from pgm to gif and add text ( level ) using ImageMagic
convert $file -resize 600x400 -pointsize 100 -annotate +10+100 $b ${b}.gif
echo $file
done
# convert gif files to animated gif
convert -delay 50 -loop 0 %d.gif[0-30] a.gif
echo OK
I had :
- manually removed some images, when changes were invisible
- manualaly add each image to initial animated gif . Last image ( frame ) was added a few times to be better visualised ( it looks like animation was stoped on the last image )
- added the comment
convert -delay 50 f1200.gif 1200.gif g1200.gif convert -comment "Julia set for fc(z)=z*z -0.75, exterior = white, interior = gray, unknown=red; Adam Majewski " g1200.gif gc.gif
Better way :[7]
convert `ls *.gif | sort -n` -resize 800x600 -delay 50 -loop 0 a.gif
Removed unnecesary frmes and optimised :
gifsicle a0c.gif --delete "#11" "#12" "#13" "#14" "#16" "#17" "#18" "#19" "#21" "#22" "#23" "#24" "#26" "#27" "#28" "#29" "#31" "#33" "#34" "#35" "#36" "#38" "#39"> a0e.gif gifsicle a0e.gif -O3 > a0e3.gif
Fragmentarium src code
See how important is tuning parameters ( iMax and ar2 ) to get wanted result.
First set all values to minimum. You will see red circe ( center 0 and radius = er = 2 ) with blue esterior.
Then increase iMax until you will see good aproximation of filled Julia set ( red) with blue exterior. Good iMax = 100. It is Level set method.
Now you can increase ar2 : you will see green dots ( around fixed point alfa and it's preimages). Green points will fill all the interior
#include "2D.frag"
#group Julia set
// maximal number of iterations = quality of image
// but also ability to fall into circlae with radius ar
// around alfa fixed point
// if to big then all not escaping points are unknown ( green)
uniform int iMax; slider[1,1000,10000]
// escape radius = er; er2= er*er >= 4.0
uniform float er2; slider[4.0,100.0,1000.0]
// attrating radius (around fixed point alfa) = ar ; ar2 = ar*ar
uniform float ar2; slider[0.000001,0.0004,0.008]
vec2 c = vec2(-0.75,0.0); // initial value of c
vec2 za = vec2(-0.5,0.0); // alfa fixed point
vec3 GiveColor( int type)
{
switch (type)
{
case 0: return vec3(1.0, 0.0, 0.0); break; //unknown
case 1: return vec3(0.0, 1.0, 0.0); break; // interior
case 2: return vec3(0.0, 0.0, 1.0); break; // exterior
default: return vec3(1.0, 0.0,0.0); break;}
}
// compute color of pixel = main function here
vec3 color(vec2 z0) {
vec2 z=z0;
int type=0; // 0 =unknown; interior=1; exterior = 2;
int i=0; // number of iteration
// iteration
for ( i = 0; i < iMax; i++) {
// escape test
if (dot(z,z)> er2) { type = 2; break;}// exterior
// attraction test
if (dot(z-za,z-za)< ar2) { type = 1; break;}// interior
z = vec2(z.x*z.x-z.y*z.y,2*z.x*z.y) + c; // z= z^2+c
}
return GiveColor(type);
}
Licensing
I, the copyright holder of this work, hereby publish it under the following license:
This file is licensed under the Creative Commons Attribution-Share Alike 4.0 International license.
- You are free:
- to share – to copy, distribute and transmit the work
- to remix – to adapt the work
- Under the following conditions:
- attribution – You must give appropriate credit, provide a link to the license, and indicate if changes were made. You may do so in any reasonable manner, but not in any way that suggests the licensor endorses you or your use.
- share alike – If you remix, transform, or build upon the material, you must distribute your contributions under the same or compatible license as the original.
References
- ↑ | git repositiorium
- ↑ | Program mandel by Wolf Jung
- ↑ Practical interior distance rendering by Claude Heiland-Allen
- ↑ parabolic dynamics in wikiboks
- ↑ Julia set in wikipedia
- ↑ Classification of Fatou components in wikipedia
- ↑ Stackoverflow - how-to-make-animated-gif-from-non-consecutive-images
File history
Click on a date/time to view the file as it appeared at that time.
| Date/Time | Thumbnail | Dimensions | User | Comment | |
|---|---|---|---|---|---|
| current | 12:25, 7 February 2015 | | 400 × 300 (307 KB) | Soul windsurfer | removed unnecessary frames, optimised |
| 12:42, 1 February 2015 |  | 400 × 300 (541 KB) | Soul windsurfer | User created page with UploadWizard |
File usage
The following page uses this file:
