給大家看點好玩的，用Python畫牛頓分形

今天給大家看點好玩兒的，如何用Python來畫牛頓分形。

前排提示：知乎圖片壓縮嚴重，非常影響觀感，建議點擊圖片顯示原圖觀看。

分形圖是世界上最美妙的圖形之一，它結構精細、反常，像是一頭怪物難以捉摸。它的每個細節與它的整體如出一轍，在我看來，分形是數學中最優美、吸引人的結構之一。

牛頓分形就是分形的一種，它與解方程的牛頓法(跟優化中的牛頓法是不同的方法)緊密相關，下面講述如何畫出牛頓分形。假如需要求解方程

其中， 的定義域是整個複平面。如何求解這個方程的解呢？我們可以用牛頓法。牛頓法是一種數值解法，我們首先會估算一個「比較好」的初始值 ，然後使用迭代公式

牛頓法可以確保，如果初始猜測值在根附近，那麼迭代必然收斂。而且牛頓法是個二階方法，收斂速度相當的快。下圖是迭代一步的示意圖

在 處沿著方向 下降，與 軸的交點即為 ，循環往複就能得到方程的根。

學過中學數學的我們都知道， 次方程在複數域上有 個根，那麼用牛頓法收斂的根就可能有 個目標。牛頓法收斂到哪個根取決於迭代的起始值。根據最後的收斂結果，我們把所有收斂到同一個根的起始點畫上同一種顏色，最終就形成了牛頓分形圖。下圖中展示的是方程 的情形

圖中的三種顏色代表了收斂的三個根，分別為 和 。左上角都是紅色的，代表了如果把左上角的點作為牛頓法迭代的初始值，最終會收斂到 ，左下角是綠色，代表這些初始值會收斂到 ，右邊是藍色，代表會收斂到 。神奇的是，中間的三個帶狀區域，是紅綠藍交錯的，而且無限重複自己的細節。

下面看代碼。首先是牛頓法解方程

def newton_method(f, df, c, eps, max_iter):n """f是一個函數,我們要求解f(x)=0 n df是f的導數n c是迭代的初始值 n 當兩步迭代結果的距離之差小於eps的時候即認為收斂到根n max_iter是最大迭代次數"""n for i in range(max_iter):n c2 = c - f(c) / df(c) # 迭代公式n if abs(f(c)) > 1e10: # 溢出n return None, Nonen if abs(c2 - c) < eps: # 達到收斂條件n return c2, i # 返回根和收斂所需的迭代次數n c = c2n

下面就可以用牛頓法畫圖了，首先我們需要用到下面兩個包

import numpy as npnfrom PIL import Imagen

然後，我們定義一個取色板函數，根據方程根的次序選取顏色

def color(ind, level):n """每種顏色用RGB值表示: (R, G, B),n level是灰度，收斂所用的迭代次數"""n colors = [(180, 0, 30), (0, 180, 30), (0, 30, 180),n (0, 190, 180), (180, 0, 175), (180, 255, 0),n (155, 170, 180), (70, 50, 0),n (150, 60, 0), (0, 150, 60), (0, 60, 150),n (60, 150, 0), (60, 0, 150), (150, 0, 60),n (130, 80, 0), (80, 130, 0), (130, 0, 80),n (80, 0, 130), (0, 130, 80), (0, 80, 130),n (110, 100, 0), (100, 110, 0), (0, 110, 100),n (0, 100, 100), (110, 0, 100), (100, 0, 110),n (255, 255, 255)]n if ind < len(colors):n c = colors[ind]n else:n c = (ind % 4 * 4, ind % 8 * 8, ind % 16 * 16)n if max(c) < 210:n c0 = c[0] + leveln c1 = c[1] + leveln c2 = c[2] + leveln return (c0, c1, c2)n else:n return cn

最後，就是我們的畫圖函數啦。為了少寫點參數，我把牛頓法的定義寫在`draw`這個函數中，與上面的版本稍有不同。

def draw(f, df, size, name, x_min=-2.0, x_max=2., y_min=-2.0, y_max=2.0, eps=1e-6,n max_iter=40):n Given function f, its derivative df,n generate its newton fractal with size, save to image with namen f是一個函數，我們需要求解方程 f(x) = 0n df是f的導數n size是圖片大小,單位是像素n name是保存圖像的名字n x_min, x_max, y_min, y_max定義了迭代初始值的取值範圍n eps是判斷迭代停止的條件n n def newton_method(c):n "c is a complex number"n for i in range(max_iter):n c2 = c - f(c) / df(c)n if abs(f(c)) > 1e10:n return None, Nonen if abs(c2 - c) < eps:n return c2, in c = c2nn return None, Nonenn roots = [] # 記錄所有根n img = Image.new("RGB", (size, size)) # 把繪畫結果保存為圖片n for x in range(size):n print("%d in %d" % (x, size))n for y in range(size): # 嵌套循環，遍歷定義域中每個點，求收斂的根n z_x = x * (x_max - x_min) / (size - 1) + x_minn z_y = y * (y_max - y_min) / (size - 1) + y_minnn root, n_converge = newton_method(complex(z_x, z_y))n if root:n cached_root = Falsen for r in roots:n if abs(r - root) < 1e-4: # 判斷是不是已遇到過此根n root = rn cached_root = Truen breakn if not cached_root:n roots.append(root)nn if root:n img.putpixel((x, y), color(roots.index(root), n_converge)) # 上色n print(roots) # 列印所有根n img.save(name, "PNG") # 保存圖片n

好了，有了這幾個函數，我們就可以開始實驗了。上面的 的圖可以用下面幾行代碼作出

def f(x):n return x ** 3 - 1ndef df(x):n return 3 * x * xndraw(f, df, 1000, "x^3-1.png")n

類似的，我們可以得到其他函數的牛頓分形圖。下面是看圖時間。

不知道為啥，知乎圖片壓縮的好厲害，細節丟失嚴重。。。

註：這次我們只是簡要介紹了作牛頓分形圖的原理，畫的圖還是比較原始和醜陋的。要想作出更漂亮的圖，還有兩方面需要改進，一是根據迭代收斂次數加入shading，二是顏色選取。這兩部分也是有一定技巧的，所以改天另開文章再寫吧。

最後，題圖是 ，來自Wikipedia