使用OpenGL进行Mandelbrot集的可视化

Mandelbrot集是哪一集??

Mandelbrot集不是哪一集!! 啊不对……

Mandelbrot集是哪一集!! 好像也不对……

Mandelbrot集是数集!! 所以……他不是一集而是数集??……

所以这个M...dem...集到底是什么啊??

Mandelbrot集是一个数集

Mandelbrot集\(\mathbb{M}\)(简称曼集)是一个由二元复数构成的集合,也就是一个复数集:

\[\mathbb{M}\subset\mathbb{C}
\]

也就是说,曼集的元素都是复数,也就是如下的形式:

\[\forall m\in \mathbb{M},m=a+b\mathrm{i},\left(a,b\in \mathbb{R}\right)
\]

为什么叫这么复杂的名字??

本华·曼德勃罗特(Benoit B.Mandelbrot)是分形学(Fractals) 的鼻祖,是的,Fractals这个词就是他提出的。

不过,这个集合本身并不是他本人提出的,而是一个叫做阿德里安·杜阿迪的法国数学家提出的,至于为什么以曼德勃罗特命名,根据wiki的说法是提出者为了向曼德勃罗特致敬而放弃用自己的名字命名。

这是个什么样子的集合??

一个复数\(z\)在或不在曼集中,取决于如下数列是否收敛

\[\zeta:z,f(z;c),f(f(z;c);c),\cdots,f\circ^n(z;c) \tag{1}
\]

\[z,c\in\mathbb{C}
\]

其中,\(f(z;c)=z^2+c\),\(c\)是一个复参数。

如果把\(f\circ^n(z;c)\)展开写的话就是:

\[\begin{split}
\zeta_n&=f\circ^n(z;c)=z_n^2+nc \\
z_n&=z_{n-1}^2+c
\end{split}
\tag{2}
\]

由于在其中迭代的前一个\(z\)都是以平方的形式出现,这也就意味着

\[\left|\Re\left(z_n\right)\right|\ge\left|\Re\left(z_{n-1}\right)\right|\ge 0 \\
\left|\Im\left(z_n\right)\right|\ge\left|\Im\left(z_{n-1}\right)\right|\ge 0
\]

这样一来,若\(z_n\)收敛,则\(z_m(m<n)\)必然收敛。

综上描述,有下结论:

\[z\in\mathbb{M} \Leftrightarrow \lim_{n\rightarrow\infty}\zeta_n\in\mathbb{C} \tag{3}
\]

式3即曼集判定的充要条件或定义。

为了让问题变得简单,这里可以令初始值\(z=0\),这样式2变成了仅与参数\(c\)相关的形式:

\[\begin{split}
z_0 &=0 \\
z_1 &=c \\
z_2 &= z_1^2+c = c^2+c \\
z_3 &= z_2^2+c = c^4+2c^3+c^2+c \\
&\cdots \\
\end{split}
\]

复数??迭代??还收敛??听起来好像很复杂……

这个……是的,这些东西听起来很不直观,以至于我也花了好久看了大量资料才搞清楚是什么意思。

但是,这件事本身的中心思想并不复杂,其无非做了这样的一件事情:

  1. 给定一个复数\(c\)
  2. 计算\(z_+=z^2+c\)(\(z_+\)会成为下一次的\(z\))
  3. 就这样,反复执行2,直到山穷水尽之时
    1. \(z_+\)仍然没有飞升为\(\infty\),那\(c\in\mathbb{M}\)
    2. 否则存在某一次\(z_+=\infty\)则\(c\notin\mathbb{M}\)

山穷水尽之时??那我恐怕等不到那一天了……

不必沮丧,我们甚至没打算去推导出\(z_n\)的完整展开形式,而且,计算机也无法处理需要运行无数次才能解决的问题(停机问题)。

但是,在有限的,可预见的次数(迭代数)内,进行这样的判定还是可行的并且能得出结果的。

但是这样的话肯定会有纰漏吧

说的对,因为只要迭代数\(n\)有限,\(z_n\)是总能求出确定值。但正因为如此,迭代数的增长才会产生实际的意义。这就好像,没有人能求出\(\pi\)的所有小数位数,但是每当有人多求出正确的一位的时候,其结果与实际值就总会更加接近。也就是说,只要范围有限,我们求得的曼集就总比实际的曼集大,但随着次数的增多,结果集合会向曼集逐渐地靠近。

所幸,在有限的次数内,我们尽管不能立即断言哪些数一定属于曼集,但我们能断言哪些数一定不属于曼集,那正是:

还没上高速公路就要起飞的,那必然是要玩完的

在有限次计算就体现出明显的发散趋势的时候,那基本也就抢救无效了

好吧好吧,那有哪些值是可以放弃治疗的?

根据迭代函数的定义,\(f(z;c)=z^2+c\),我们对函数的范数(也就是,或者称为复数的绝对值)进行考察,至少如果\(f\)收敛,其模与带入z的模至少是相等的(当然,原则上幅角也应当相等),:

\[|f(z;c)|=\left|z^2+c\right|=|z| \\
|z|\ge0
\]

根据绝对值三角不等式(这个对复数域同样成立,实际上就是向量绝对值三角不等式的移植),若\(|z|=0\)则\(|c|\)必然为\(0\),否则

\[|z|^2-|c|=|z^2|-|c|\le|f(z;c)|=|z|\le|z^2|+|c|=|z|^2+|c| \\
1-|c|\le|z|\le1+|c|
\]

由于迭代初始条件\(z_0=0\),因此经过n次迭代后的\(z\)一定是

\[z_n=\sum_{m=0}^{(n-2)^2}{k_mc^{m}} \\
c\le|z_n|\le\sum_{m=0}^{(n-2)^2}{|k_m||c^{m}|}
\]

上面的不等数右侧产生了一个很熟悉的形式——幂级数,如果要幂级数\(\sum_{m=0}^{(n-2)^2}{|k_m||c^{m}|}\)收敛,则\(|c|<1\)是必要条件(由于\(k_m\)的大小这里没有多加考虑,即使\(|c|<1\)满足,发散的可能同样存在,这个在后续的迭代过程中会进一步的被发掘出来),否则根据幂级数的性质级数必然发散,这样得到了一个必要阈值:\(|c_\Theta|=1\),将这个阈值带入,我们就能得到\(|z|\)的基础发散判定值:

\[|z|\ge2\Rightarrow \lim_{n\rightarrow\infty}{\zeta_n}=\infty
\]

在后面的程序中,我们将使用这个值作为发散判定,意思就是说,如果算出\(|z|\ge 2\)了,那就不用考虑了,你不可能收敛的,至于剩下的值,我只假设但不断言它一定在集合内,因此我暂且将它认为是集合内的,如果经过更多的计算发现也达到这个红线条件了,那么同样不论这个值挣扎了多久,OUT!

但是,复数的收敛并不限于模收敛,当复数被整理为辐角形式\(r\angle\theta\)时,若\(\theta\)没有稳定值即便\(r\)已经收敛仍然不能断言复数收敛(毕竟啊,那个复平面上的向量如果一直在原地打转,分量就总在以一个有界值变化,宛如论证的丧钟在转动),但是对于这个函数而言根本不会不存在这种情形

哇,断言的这么轻松??何出此言??

其实很容易,找到这样的向量,利用迭代函数平方相加后保持模不变即可,所幸的是,这样的值非常有限,因为首先一个条件就将范围急剧缩减:

\[|c|=1
\]

是的,只有\(|c|=1\)的时候,平方计算后才可能保证模不变(依据辐角计算法若\(z=r_z\angle\theta\),则\(z^2=r_z^2\angle{2\theta}\)),当然我是指仅一次,因为还另外需要一件事保证\(z^2+c=c^2+c\)后的模长也是\(|c|\)或者说是\(1\):

\[\cos{<c^2,c>}=-\frac{1}{2}
\]

也就是说\(\arg{z^2}-\arg{c}=\frac{2k\pi}{3}\)才可以(这个式子写的不是很严谨,但你当成是两个夹角为\(120^\circ\)就好),当然你会反驳:

\(c\)是随意取得,也就是说幅角就可以是任意值,当然,\(z^2\)的取值决定于\(c\),也就是\(z^2\)可以和\(c\)构造连续的函数映射,那么你根本保证不了\(z^2\)一定能够避开这些值啊你个代数白痴!!

啊是的,但是你很快就会意识到一个事实:

这个丧钟,转不久的!!

什么??

确实,当满足上述两个条件之后,得到的\(z_+\)的结果模仍然是\(1\),只是转了个角度,但是架不住这个函数是迭代的啊,这次得到的\(z_+\)很明显位于这\(z^2\)和\(c\)的角平分线处……

然后,夹角条件被打破了!!

就像雪崩一样,这个微妙的平衡很快就被打破了,并且马上开始急剧的发散(或收敛)。

如丸走坂

- Ramp Rollerball

这样一来,通过多次的迭代运算,不合格的点一点一点被筛出,如果次数足够多,这个集合最后形成的图形被称为Mandelbrot分形

大概长出来是这样的:

使用OpenGL进行Mandelbrot集的可视化

这个M..d什么b什么的分形有什么特点么??

作为一个分形,它最大的特点之一就是:

德罗斯特效应(自相似性,Droste Effect)

如果对这个图形局部放大,你会看到这个与最开始看到一致的子图性,而且:

放大一次有,一直放大一直有!放大放大再放大!每一个子部都看的清清楚楚!

德罗斯特效应本质是递归式的图形体现,因为迭代函数通过\(z_+=f(z;c)\)进行了递推,这样的话,也就意味着,分形具有无限细微的结构,而且跟整体是相似的。

而具有无限精细的结构也就意味着:它具有有限的面积,而周长却很可能无穷大

当然由于曼集的边缘性质非常复杂,这里就不再详细讨论他的周长是否收敛了(面积收敛是一看就能看出来的),毕竟这篇文章是为了说明如何通过OpenGL实现这个集合的可视化,上面其实已经花了大量的篇幅去推导一个阈值了。

Mandelbrot集的兄弟——Julia集简介

之前是以\(c\)作为参数进行迭代,得到了曼集,那么,如果把\(c\)和\(z\)交换一下,让\(z\)做参数,就得到了Mandelbrot集的兄弟——Julia集(以下简称朱集),以法国数学家加斯顿·朱利亚(Gaston Julia)命名。

— Oh,Julia~

— Oh,Mandelbrot~

— W...What??

同样,因为递归仍然存在,Julia集得到的也会是Julia分形

使用OpenGL进行Mandelbrot集的可视化

迭代函数没变,因此它的发散判定也是和曼集是一致的,这里仅仅介绍一下这个亲戚,并不打算实现这个Julia集的可视化。

OpenGL可视化实现

其实上面的函数都已经推导完毕的情况下,并且也知道集合是如何生成的之后,代码实现就容易多了,绘制曼集的代码如下(只考虑整数点):

基本绘制(二值化)

void DrawMandelbrot(int maxIter)
{
complex<double> z(0,0); // 位于标准库<complex>头文件中
complex<double> c;
complex<double> f(0,0);
bool diverge = false;
glBegin(GL_POINTS);
for (GLint x = -400; x < 400; x++)
{
diverge = false;
for (GLint y = -300; y < 300; y++)
{
c = { static_cast<double>(x)/150,static_cast<double>(y)/150 };
int iter;
for (iter = 0; iter < maxIter; iter++)
{
z = z*z + c;
if ( abs(z) >= 2 ) //发散判定
{
diverge = true;
break;
}
}
if (!diverge)
{
glColor3f(0.0f, 0.0f, 0.0f); //集合内黑色
}
else
{
glColor3f(1.0f, 1.0f, 1.0f); //集合外白色
diverge = false;
}
z = { 0,0 };
glVertex2i(x, y);
}
}
glEnd();
}

效果图:

使用OpenGL进行Mandelbrot集的可视化

有问题,如果发散判定值不取2会有什么影响??

问得好,实际上我正想说,其实判定值最好取2,如果实在无法取得的话,也至少请保证它大于2

如果小于2,则有些本来应该在曼集内的点被剔除,这样形成的曼集图形是不完整的。

而如果大于2,当迭代次数非常大的时候,实际上得到的结果与2的时候相差不大(细微的内部可能会有所差别),因为次数非常大的时候,该发散的值一般早都散得没影了,用再大的有限值根本拦不住。但是如果迭代次数非常小的时候,得到的曼集会包含大量冗余的点,因为阈值没能称职地起到约束作用。

迭代着色

当然,无论如何,我们可以让上面的图更加漂亮一些,比如我们以迭代多少次就发散了作为着色的依据,这里对上面的图形进行着色处理。

实际上只需要把集合外的着色代码里面的颜色分量改成与迭代数iter相关的形式:

//...
else
{
glColor3f(static_cast<float>(iter) / maxIter, (0.5f*iter) / maxIter, 0.15f*iter / maxIter);
diverge = false;
}

绘制结果:

使用OpenGL进行Mandelbrot集的可视化

值得注意的是,由于迭代次数是一个离散的值,因此整个图片的着色显得并不是那般的丝滑和连续,在其他文章中还有提到过对数化将使着色更加柔和漂亮的方法,这里不再阐述。

当指定参数maxIter的次数不同,形成的图形也不同,并且,随着迭代次数的增大, 形状逐渐趋于稳定,更接近于标准的曼集图样:

使用OpenGL进行Mandelbrot集的可视化

上一篇:20145225 实验四《Andoid开发基础》


下一篇:Day Three(Beta)