自平方分形

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

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

所以令二次判别式1-(4z2)/λ=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;

}

 

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

z0=c,c为常量复数，

zk=z2k-1+z0；

也即f(z)=z2+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;

}