קובץ:Discrete dynamic in case of complex quadratic polynomial.gif
אין גרסה ברזולוציה גבוהה יותר.
Discrete_dynamic_in_case_of_complex_quadratic_polynomial.gif (400 × 300 פיקסלים, גודל הקובץ: 307 ק"ב, סוג MIME: image/gif, בלולאה, 55 תמונות, 28 שניות)
|
זהו קובץ שמקורו במיזם ויקישיתוף. תיאורו בדף תיאור הקובץ המקורי (בעברית) מוצג למטה. |
תקציר
| תיאור |
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. |
| תאריך יצירה | |
| מקור | Own program[1] which uses code by Wolf Jung[2] and Claude Heiland-Allen[3] ( thx also for help with OMP code ) |
| יוצר | 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);
}
רישיון
אני, בעל זכויות היוצרים על עבודה זו, מפרסם בזאת את העבודה תחת הרישיון הבא:
הקובץ הזה מתפרסם לפי תנאי רישיון קריאייטיב קומונז ייחוס-שיתוף זהה 4.0 בין־לאומי.
- הנכם רשאים:
- לשתף – להעתיק, להפיץ ולהעביר את העבודה
- לערבב בין עבודות – להתאים את העבודה
- תחת התנאים הבאים:
- ייחוס – יש לתת ייחוס הולם, לתת קישור לרישיון, ולציין אם נעשו שינויים. אפשר לעשות את זה בכל צורה סבירה, אבל לא בשום צורה שמשתמע ממנה שמעניק הרישיון תומך בך או בשימוש שלך.
- שיתוף זהה – אם תיצרו רמיקס, תשנו, או תבנו על החומר, חובה עליכם להפיץ את התרומות שלך לפי תנאי רישיון זהה או תואם למקור.
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
היסטוריית הקובץ
ניתן ללחוץ על תאריך/שעה כדי לראות את הקובץ כפי שנראה באותו זמן.
| תאריך/שעה | תמונה ממוזערת | ממדים | משתמש | הערה | |
|---|---|---|---|---|---|
| נוכחית | 15:25, 7 בפברואר 2015 | | 300 × 400 (307 ק"ב) | Soul windsurfer | removed unnecessary frames, optimised |
| 15:42, 1 בפברואר 2015 |  | 300 × 400 (541 ק"ב) | Soul windsurfer | User created page with UploadWizard |
שימוש בקובץ
אין בוויקיפדיה דפים המשתמשים בקובץ זה.
שימוש גלובלי בקובץ
אתרי הוויקי השונים הבאים משתמשים בקובץ זה:
- שימוש באתר en.wikibooks.org
