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

}