自平方分形

在二维平面中，复数可表示为z=x+yi，x和y是实数，i²=-1。自平方分形的变换函数是z₂=λz(1-z)，λ是一个常量复数。

由z₂=λz(1-z)得z=f^-1(z₂)=1/2(1+sqrt(1-4z₂/λ)，

所以令二次判别式1-(4z₂)/λ=discr,设z=re(z)+im(z)，re(z),im(z)分别为实部和虚部。

则:

re(z)=1/2*(1+sqrt((abs(discr)+re(discr))/2)，

im(z)=1/2*sqrt((abs(discr)+re(discr))/2。

 

当λ=1.5+0.8i时的图形如下：

 

程序如下：

/* OpenGL源码，自平方分形z2=λz(1-z)，λ=1.5+0.8i */ #include<GL/glut.h> #include<GL/gl.h> #include<iostream> #include<stdlib.h> #include<math.h> #include<time.h> using namespace std; struct complexNum { GLfloat x,y; }; void init() { glClearColor(1.0,1.0,1.0,0.0); } void plotPoint(complexNum z) { glBegin(GL_POINTS); glVertex2f(z.x,z.y); glEnd(); } void solveQuadraticEq(complexNum lambda,complexNum* z) { GLfloat lambdaMagSq,discrMag; complexNum discr; static complexNum fourOverLambda={0.0,0.0}; static GLboolean firstPoint=true; if(firstPoint) { lambdaMagSq=lambda.x*lambda.x+lambda.y*lambda.y; fourOverLambda.x=4.0*lambda.x/lambdaMagSq; fourOverLambda.y=-4.0*lambda.y/lambdaMagSq; firstPoint=false; } discr.x=1.0-(z->x*fourOverLambda.x-z->y*fourOverLambda.y); discr.y=z->x*fourOverLambda.y+z->y*fourOverLambda.x; discrMag=sqrt(discr.x*discr.x+discr.y*discr.y); z->x=sqrt((discrMag+discr.x)/2.0); z->y=0.5*sqrt((discrMag-discr.x)/2.0); if(rand()<RAND_MAX/2) { z->x=-z->x; z->y=-z->y; } if(discr.y<0) z->x=-z->x; z->x=0.5*(1-z->x); } void selfSqTransf(complexNum lambda,complexNum z,GLint numPoints) { GLint k; for(k=0;k<10;k++) solveQuadraticEq(lambda,&z); for(k=0;k<numPoints;k++) { solveQuadraticEq(lambda,&z); plotPoint(z); } } void display() { clock_t start,end; start=clock(); GLint numPoints=2000000; complexNum lambda={1.5,0.8}; complexNum z0={1.5,0.4}; glClear(GL_COLOR_BUFFER_BIT); glColor3f(0.0,0.0,1.0); selfSqTransf(lambda,z0,numPoints); glFlush(); } void reshape(int w,int h) { glViewport(0,0,w,h); glMatrixMode(GL_PROJECTION); glLoadIdentity(); gluOrtho2D(-0.25,1.25,-0.75,0.75); glClear(GL_COLOR_BUFFER_BIT); } int main(int argc,char **argv) { glutInit(&argc,argv); glutInitDisplayMode(GLUT_SINGLE|GLUT_RGB); glutInitWindowSize(300,300); glutInitWindowPosition(100,100); glutCreateWindow(argv[0]); init(); glutDisplayFunc(display); glutReshapeFunc(reshape); glutMainLoop(); return 0; }

 

若修改一下自平方函数，令:

z₀=c,c为常量复数，

z_k=z²_k-1+z₀；

也即f(z)=z²+c。则可以得到Mandelbrot集。

如图：

#include<GL/glut.h> #include<GL/gl.h> #include<iostream> #include<stdlib.h> #include<math.h> #include<time.h> using namespace std; GLfloat xComplexMin=-2.0,xComplexMax=0.5; GLfloat yComplexMin=-1.2,yComplexMax=1.2; GLfloat complexWidth=xComplexMax-xComplexMin; GLfloat complexHeight=yComplexMax-yComplexMin; class complexNum{ public: GLfloat x,y; }; struct color{GLfloat r,g,b;}; void init() { glClearColor(1.0,1.0,1.0,0.0); } void plotPoint(complexNum z) { glBegin(GL_POINTS); glVertex2f(z.x,z.y); glEnd(); } complexNum complexSquare(complexNum z) { complexNum zSquare; zSquare.x=z.x*z.x-z.y*z.y; zSquare.y=2*z.x*z.y; return zSquare; } GLint mandelSqTransf(complexNum z0,GLint maxIter) { complexNum z=z0; GLint count=0; while((z.x*z.x+z.y*z.y<=4.0)&&(count<maxIter)) { z=complexSquare(z); z.x+=z0.x; z.y+=z0.y; count++; } return count; } void mandelbrot(GLint nx,GLint ny,GLint maxIter) { complexNum z,zIncr; color ptColor; GLint iterCount; zIncr.x=complexWidth/GLfloat(nx); zIncr.y=complexHeight/GLfloat(ny); for(z.x=xComplexMin;z.x<xComplexMax;z.x+=zIncr.x) for(z.y=yComplexMin;z.y<yComplexMax;z.y+=zIncr.y) { iterCount=mandelSqTransf(z,maxIter); if(iterCount>=maxIter) ptColor.r=ptColor.g=ptColor.b=0.0; else if (iterCount>(maxIter/8)) { ptColor.r=1.0; ptColor.g=0.5; ptColor.b=0.0; } else if(iterCount>(maxIter/10)) { ptColor.r=1.0; ptColor.g=ptColor.b=0.0; } else if(iterCount>(maxIter/20)) { ptColor.r=0.5; ptColor.g=ptColor.b=0.0; } else if(iterCount>(maxIter/40)) { ptColor.r=ptColor.g=1.0; ptColor.b=0.0; } else if(iterCount>(maxIter/100)) { ptColor.r=ptColor.b=0.0; ptColor.g=0.3; } else { ptColor.r=0.0; ptColor.g=ptColor.b=1.0; } glColor3f(ptColor.r,ptColor.g,ptColor.b); plotPoint(z); } } void display() { GLint nx=1000,ny=1000,maxIter=1000; glClear(GL_COLOR_BUFFER_BIT); mandelbrot(nx,ny,maxIter); glFlush(); } void reshape(int w,int h) { glViewport(0,0,w,h); glMatrixMode(GL_PROJECTION); glLoadIdentity(); gluOrtho2D(xComplexMin,xComplexMax,yComplexMin,yComplexMax); glClear(GL_COLOR_BUFFER_BIT); } int main(int argc,char **argv) { glutInit(&argc,argv); glutInitDisplayMode(GLUT_SINGLE|GLUT_RGB); glutInitWindowSize(600,400); glutInitWindowPosition(100,100); glutCreateWindow(argv[0]); init(); glutDisplayFunc(display); glutReshapeFunc(reshape); glutMainLoop(); return 0; }