Julia 集和 Mandelbrot 集的图片 / 计算机程序

我们将展示两种制作可以绘制 Mandelbrot 集的程序的方法，其中迭代趋于 ∞：一个程序从头开始，但只绘制单张图片，另一个程序可以做任何你想要的事情（缩放、改变颜色、渲染为文件），但是这些功能交给了另一个程序，即 *Ultra Fractal*。

一张图片

Mandelbrot 集的一段，用于 $z^{2}+0.009/z+c$

让我们制作一个尽可能简单的程序：它只绘制单张图片，不做其他事情（完成工作后它会关闭）。它最初是用汇编语言编写的，但这里我们用伪代码编写。它由两部分组成：绘制过程和颜色尺度。

我们为族 $f(z)=z^{2}+{\frac {p}{z}}+c\,$ 绘制 Mandelbrot 集的一段，其中 p = 0.009。我们让窗口大小为 800x600 像素，并且我们想象这段的左下角位于坐标为 (ax, ay) 的点，并且它具有宽度 h。我们有 $f'(z)=2z-{\frac {p}{z^{2}}}\,$，因此我们需要有限的临界点是方程 $z^{3}=p/2$ 的解。我们选择实数解 cri = $(p/2)^{1/3}$。

绘制过程

p = 0.009

cri = $(p/2)^{1/3}$

r = 1.0E200（跳出半径的平方）

u = log(log(r))

v = 1/log(2)（2 是有理函数的次数）

ax = -0.7028233689263（左下角的 x 坐标）

ay = 0.1142331238418（左下角的 y 坐标）

h = 0.0092906805859（这段的宽度）

g = h / 800

m = 850（最大迭代次数）

thick = h * 0.0005（边界的厚度）

dens = 24.9（颜色的密度）

disp = 432.0（颜色尺度的偏移）

for i = 0 to 799 do

begin

cx = ax + i * g（像素 (i, j) 的 x 坐标）

for j = 0 to 599 do

begin

cy = ay + j * g（像素 (i, j) 的 y 坐标）

x = cri

y = 0

xd = 0

yd = 0

f = 0

n = 0

while (n < m) and (f < r) do

begin

n = n + 1

x1 = x * x - y * y

y1 = 2 * x * y

f = x * x + y * y

x2 = 2 * x - p * x1 / (f * f)

y2 = 2 * y + p * y1 / (f * f)

temp = xd

xd = x2 * xd - y2 * yd + 1

yd = x2 * yd + y2 * temp

fd = xd * xd + yd * yd

x = x1 + p * x / f + cx

y = y1 - p * y / f + cy

f = x * x + y * y

end

if (n = m) or (log(f) * sqrt(f) < thick * sqrt(fd)) then

setpixel(i, 599 - j, 0)

else

begin

s = n - v * (log(log(f)) - u)

n = round(dens * s + disp) mod 720

col = paletteRGB(col1[n], col2[n], col3[n])

setpixel(i, 599 - j, col)

end

end

end

解释

我们设置了 $c=cx+icy$， $z=x+iy$， $z'=xd+iyd$ 以及 $x1+iy1=(x+iy)^{2}$，并且我们使用了 $1/z=z^{-}/|z|^{2}$。

导数z'的连续计算为$xd+iyd=(x2+iy2)*(xd+iyd)+1$，其中$x2+iy2=2*(x+iy)-p*(x1-iy1)/(f*f)$ 以及 $f=x^{2}+y^{2}$。迭代中的下一个点为$x+iy=(x1+iy1)+p*(x-iy)/f+(cx+icy)$。距离函数为

log(√f)*√f/√fd   (= log(f)*√f/(2*√fd))

对于最后一个z值，这里$f=x^{2}+y^{2}$ 以及 $fd=xd^{2}+yd^{2}$。从迭代次数中减去的[0, 1[中的数字（为了得到实际的迭代次数）为$log(log(|z|)/log(r))/log(2)=v*(log(log(f))-u)$，对于最后一个z，这里$v=1/log(2)$ 以及 $u=log(log(r))$.

paletteRGB(r, g, b) 是一个整数，用于索引具有 RGB 值 (r, g, b) 的颜色，0 表示黑色。col、col2 和 col3 是从 0 到 719 的数组，包含从 0 到 255 的整数，它们将在下一部分构建。

颜色比例

颜色在计算机中由其三种原色红、绿和蓝的组合立即确定，它们的比例用 0 到 255 之间的整数测量。这三个数字的组合是该颜色的 RGB 值。颜色对应于边长为 256 的立方体中的所有点，我们通过在该立方体中沿着椭圆进行规律移动来构建循环颜色比例。空间中的椭圆由以下确定：

坐标为(a, b, c)的中心

长轴rmaj和短轴rmin

两个角度，角度g和h，确定椭圆平面的方向

一个角度q，确定椭圆在此平面中的方向

u = cos(g * pi / 180) * cos(h * pi / 180) {(u,v,w) 是垂直于该平面的单位向量}

v = cos(g * pi / 180) * sin(h * pi / 180)

w = sin(g * pi / 180)

x = -a {(x,y,z) 是该平面中的向量}

y = -b

z = a * u / w + b * v / w

e = sqrt(sqr(x) + sqr(y) + sqr(z))

x1 = x / e {(x,y,z) 对应的单位向量}

y1 = y / e

z1 = z / e

x2 = v * z1 - w * y1 {该平面中垂直于 (x1,y1,z1) 的单位向量}

y2 = w * x1 - u * z1

z2 = u * y1 - v * x1

e = cos(q * pi / 180)

f = sin(q * pi / 180)

x1 = e * x1 + f * x2 {该平面中以角度 q 为方向的单位向量}

y1 = e * y1 + f * y2

z1 = e * z1 + f * z2

x2 = v * z1 - w * y1 {该平面中垂直于此方向的单位向量}

y2 = w * x1 - u * z1

z2 = u * y1 - v * x1

for i = 0 to 719 do

begin

e = rmaj * cos(i * pi / 360)

f = rmin * sin(i * pi / 360)

col1[i] = round(a + e * x1 + f * x2) {椭圆上点的三个坐标}

col2[i] = round(b + e * y1 + f * y2)

col3[i] = round(c + e * z1 + f * z2)

end

在图片中，我们有：(a, b, c) = (176, 176, 176), rmaj = rmin = 88, g = 146, h = -32, q = 0。

Ultra Fractal

Ultra Fractal 是一款用于绘制分形的程序，用户可以在很大程度上编写自己的公式程序，而主程序执行所有分形程序共有的操作，即从像素到像素的遍历、缩放、着色、以及将大图片保存为文件。您的工作量降至最低，除了无法通过像素到像素的方式常规绘制的图片，例如吸引子和景观，您可以随心所欲地绘制任何您想要的图片：非复杂函数、四元数等等。您可以直接操作复数。

首先，您应该下载程序的免费试用版并学习其工作原理。然后，您应该阅读文章 使用 Ultra Fractal 绘制曼德勃罗集和朱莉娅集的迷人图片，并从该网站下载（免费）程序，以便与 Ultra Fractal 一起运行。

为了进行绘制，需要将两个程序结合起来：公式程序和着色程序。公式程序中有一个名为“loop”的部分，当这个循环完成工作后，迭代次数和复数 *z* 的最后一个值将被传递到着色程序，着色程序根据这两个数字计算颜色。然而，Ultra Fractal 是为“所有”类型的分形设计的，不幸的是，朱莉娅集和曼德勃罗集并不完全符合这种形式：对于趋向于循环的迭代，必须继续迭代才能找到循环的阶数和吸引力，而我们并没有使用序列中 *z* 的最后一个值来确定颜色。因此，我们将这个循环替换为一个虚拟循环（只运行一次），并在“init”部分执行所有操作。这意味着我们不能使用 Ultra Fractal 的默认着色程序，我们需要编写一个新的着色程序。

您还需要注意的是，在 Ultra Fractal 中，复数 z 的范数 |z| 是其长度的 **平方**: |z| = $x^{2}+y^{2}$.

曼德勃罗集的一部分，对应于 $z^{2}/2+0.26*z^{4}+c$

我们将展示一个程序，该程序绘制曼德勃罗集，对应于函数族 $f(z)+c$，其中 $f(z)$ 是一个多项式，其前两项为零，因此从 $z^{2}$ 或 *z* 的更高次幂开始，因为我们可以将 0 作为有限临界点。

这两个程序可以分别复制并粘贴到一个空的公式和着色文档中。我们已经插入了 4 次的多项式 $z^{2}/2+p*z^{4}+c$，其中 p 是一个实参数。图片显示了 p = 0.26 的曼德勃罗集的一部分。

公式程序

Mandelbrot {

global

float p = 0.26 ;parameter

float deg = 4 ;degree of the polynomial

float v = 1 / log(deg)

float g = 10 * log(10)

float r = exp(g) ;square of the radius of the bail-out circle

float u = log(g)

float tb = sqr(@thick / (1000 * #magn)) ;for the thickness of the boundary

float h = 1 / (1500 * #magn * @width) ;a very small real number

init

complex z = 0 ;critical point

complex zd = 0 ;the sequence of the derivatives

complex z1 = 0

float w = 0

int n = 0

while n < #maxit && |z| < r

n = n + 1

z1 = z^2/2 + p*z^4

zd = ((z+h)^2/2 + p*(z+h)^4 - z1) * zd / h + 1

z = z1 + #pixel

endwhile

if n == #maxit || sqr(log(|z|)) * |z| < tb * |zd|

w = -1

else

w = n - v * (log(log(|z|)) - u)

endif

;begin fictive loop

z = 0

n = 0

loop

n = n + 1

z = z + #pixel

if n == 1

z = w

endif

bailout

n < 1

;end fictive loop

default

title = "Mandelbrot"

maxiter = 100

param thick

caption = "boundary"

default = 1.0

endparam

param width

caption = "width"

default = 640

endparam

}

解释

我们将逃逸圆的半径 *r* 的平方设置为 $10^{10}$，并将边界厚度设置为“thick/(1000 * #magn)” ，因此，当我们放大时，它会变得更薄。

我们通过 $(f(z+h)-f(z))/h$ 计算了函数 $f(z)$ 的导数，其中 h 是一个很小的数字，即“h = 1/(1500*#magn*width)” ，其中“width” 是图片的宽度。默认值为窗口的宽度（以像素为单位，这里设置为 640），但如果您绘制了大型图片，则应输入图片的宽度。

导数的连续计算$z'_{k+1}=f'(z_{k})z'_{k}+1$ 为 "zd = (f(z+h) - f(z)) * zd / h + 1"。下一轮迭代$z_{k+1}=f(z_{k})+c$ 为 "z = f(z) + #pixel"。距离估计的平方$log|z_{k}|*|z_{k}|/|z'_{k}|$ 为 "sqr(log(|z|)) * |z| / |zd|"。从迭代次数中减去的数字$log(log(|z_{k}|)/log(r))/log(deg)$ 为 "v * (log(log(|z|)) - u)"，其中 "v = 1/log(deg)" 和 "u = log(log(r))"。

最终的复数 z（将被传输到着色程序），实际上是实数，当该点属于曼德布罗特集合或边界时设置为 -1，否则设置为实数迭代次数。它被传输到我们称为 Gradient 的着色程序，作为 #z，并在此被设置为实数 s。当 s 为负时，颜色设置为 "#solid"，否则我们将 s 乘以一个决定密度的数字 "dens"，并向该数字添加一个决定比例偏移的数字 "disp"，由于我们希望此偏移为百分比，因此我们将结果除以 100： "u = (dens * s + disp) / 100"。在 Ultra Fractal 中，循环颜色比例的颜色由区间 [0, 1[ 中的实数索引，我们让这个数字 "#index" 成为数字 u 的非整数部分： "#index = u - trunc(u)"。

着色程序

Gradient {

final

float s = real(#z)

float u = 0

if s < 0

#solid = true

else

u = (@dens * s + @disp) / 100

#index = u - trunc(u)

endif

default

title = "Gradient"

param disp

caption = "displace"

default = 0

endparam

param dens

caption = "density"

default = 1.0

endparam

}