2

我正在尝试使用牛顿迭代显示 Julia 集,但得到的结果如下所示。可能是什么问题呢?这是我的编辑:固定,工作代码:

#version 130

in vec3 vs_out_col;
in vec3 vs_out_pos;

out vec4 fs_out_col;

vec2 cplx_mul(vec2 z1, vec2 z2)
{
    return vec2(z1.x * z2.x - z1.y * z2.y, z1.y * z2.x + z1.x * z2.y); 
}

vec2 cplx_div(vec2 z1, vec2 z2)
{
    float denom = z2.x * z2.x + z2.y * z2.y;
    return vec2(
                (z1.x * z2.x + z1.y * z2.y) / denom,
                (z1.y * z2.x - z1.x * z2.y) / denom
            );
}

vec2 f(vec2 z)
{
    vec2 res = cplx_mul(cplx_mul(z, z), z);
    return vec2(res.x - 1.0f, res.y);
}

vec2 f_der(vec2 z)
{
    return 3 * cplx_mul(z, z);
}

void main()
{
    vec2 z = vs_out_pos.xy;

    for(int i = 0; i < 30; ++i){
        z = z - cplx_div(f(z), f_der(z));
    }

    vec2 root1 = vec2( 1.0f     ,  0.0f);
    vec2 root2 = vec2(-1.0f/2.0f,  1.0f/2.0f * sqrt(3.0f));
    vec2 root3 = vec2(-1.0f/2.0f, -1.0f/2.0f * sqrt(3.0f));

    if(abs(length(z - root1)) < 0.5f){
        fs_out_col = vec4 (1, 0, 0, 1);
    }
    else if(abs(length(z - root2)) < 0.01f){
        fs_out_col = vec4 (0, 1, 0, 1);
    }
    else if(abs(length(z - root3)) < 0.01f){
        fs_out_col = vec4 (0, 0, 1, 1);
    }
    else{
        fs_out_col = vec4 (0, 0, 0, 1);
    }
}

结果如下:在此处输入图像描述 固定结果在此处输入图像描述

4

1 回答 1

2

您应该正确计算 f(z)=z^3-1 的导数,导数为 f'(z)=3*z^2,缺少因子 3。

常数 1.0f 的减法真的那么简单吗?

于 2014-02-26T16:24:06.857 回答