牛顿迭代分形

牛顿迭代是求出多项式解的一种方法,也可以画出漂亮的分形图。

牛顿迭代

过程就不写了,可查看注释,里面分析很清楚。1
公式:
$$ z_{n+1}=z_n+\frac{f(z_n)}{f’(z_n)} $$

分形图则是根据以上公式在复平面上绘制的。

  • 确定想要绘制数据范围,如 $z=a+bi,a\in[-1,1],b\in[-1,1]$
  • 将复平面数据映射在画布上,如 x 轴为实部(a),y 轴为虚部(b)
  • 画布上每个点进行牛顿迭代,当 $|z_{n+1}-z_n| < 0.001$,即迭代到多项式一个解,根据迭代次数和靠近的解设置该点的颜色2

实现

刚发现 golang 语言原生支持复数运算,省了不少事。
以下以 $f(z)=z^2+1$ 多项式为例。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
const (
_NEWTON_E = 0.001
_NEWTON_MAP_MAX = 3.0
)

func newton_f(z complex128) complex128 {
f := z*z + 1
return f
}

func newton_df(z complex128) complex128 {
df := 2 * z
return df
}

func newton_iter(z complex128) int {
for i := 0; i < 255; i += 9 {
z1 := z - newton_f(z)/newton_df(z)
if cmplx.Abs(z1-z) < _NEWTON_E {
return i
}
z = z1
}
return 255
}

func newton(img *image.RGBA, limit int) {
for i := 0; i < limit; i++ {
zx := _NEWTON_MAP_MAX*float64(i)/float64(limit) - _NEWTON_MAP_MAX/2
for j := 0; j < limit; j++ {
zy := _NEWTON_MAP_MAX*float64(j)/float64(limit) - _NEWTON_MAP_MAX/2
gray := uint8(newton_iter(complex(zx, zy)))
point(img, i, j, color.Gray{gray})
}
}
}

测试

1
2
3
4
5
6
7
8
func Test_newton(t *testing.T) {
const max_len = 500
img := image.NewRGBA(image.Rect(0, 0, max_len, max_len))
newton(img, max_len)
f, _ := os.OpenFile("newton.png", os.O_WRONLY|os.O_CREATE, 0600)
defer f.Close()
png.Encode(f, img)
}

结果:

分形欣赏

$f(z)=z^3-1$

$f(z)=z^8+15z^4-16$

$f(z)=z^7+iz-1$

无心插柳

开始我不知道 golang 自带复数运算,所以自己写了个错误的,结果生成了意想不到的图片,也挺好看的。
颇有水墨之风 (≖‿ゝ≖)✧


牛顿迭代分形
https://wishlily.github.io/article/code/2017/03/30/undefined/
作者
Wishlily
发布于
2017年3月30日
许可协议