牛顿迭代分形

Wishlily 收录于开发

2017-03-30 约 478 字预计阅读 1 分钟

$/article/2017/03/30/newton-iter-fractal/newton_err_3.png$

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

牛顿迭代

过程就不写了，可查看注释，里面分析很清楚。¹
公式：
$$ 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$，即迭代到多项式一个解，根据迭代次数和靠近的解设置该点的颜色²

实现

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

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})
        }
    }
}

测试

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)
}