返回列表 回复 发帖

(z-1)(z+1)(z-b+a)(z-b-a) 的牛迭M集征解(越玩越有趣)

这两天,我造作P(z) = (z + 1) * (z - 1) * (z - b + a) * (z - b - a)
;; F(z) = z - P(z)/P'(z)的M集,折腾来折腾去,都以失败告终。这是对a着色,b定位在原点,阈值设置为10^(-8).软件中的效果如下:
Fractal1.jpg
2010-6-7 20:07

不胜感激。
这个分形不好做,条件很多。没看出几何画板构造的思路来,等空下来后再试试另找它法。
1# 柳烟
z定位何处?
你把源代码发上来。到现在基本上没有作不出来的复变分形了!只要知道算法。
4# 分形几何


这是源代码,代为上传:
global:
  int n_attract = 4   ; number of attractors
  
  int n_orbits        ; number of orbits that have to be checked
  if @giveup=="all done" || @giveup=="any done"
    n_orbits = 2
  else
    n_orbits = 1
  endif   
  
  ; turning 'reorder' on will switch attractors 1 and 2 (of 0..3)
  int ixA
  int ixB
  if @reorder
    ixA = 1
    ixB = 2
  else
    ixA = 2
    ixB = 1
  endif
  
init:
  complex a = #pixel
   
  complex attract[4]  ; n_attract
  complex orbit[2]    ; n_orbits
  
  attract[0] = -1
  attract[ixA] = 1
  
; pre-calculate some factors outside the loop
  complex tmps
  complex fac0
  complex fac2
  
  if @pchoice=="halfdist"
    attract[ixB] = @b-a
    attract[3] = @b+a
    fac0  = @b^2 - a^2
    fac2  = fac0 - 1
    tmps  = @b^2 + 2*a^2 + 2
  else ; @pchoice=="square"
    complex tmp = sqrt(-a)
    attract[ixB] = @b + tmp
    attract[3] = @b - tmp
    fac0  = @b^2 + a
    fac2  = fac0 - 1
    tmps  = @b^2 - 2*a + 2
  endif

  if @giveup == "custom orbit start"
    orbit[0] = @z_init
  else   
    tmps = sqrt(tmps/3) / 2
$ifdef UNUSED
    if real(tmps)*imag(@b) > -imag(tmps)*real(@b)
      tmps = - tmps
    endif
$endif
    orbit[0] = @b/2 - tmps
    orbit[1] = @b/2 + tmps
   
    if @giveup == "other"
      orbit[0] = orbit[1]
      ; orbit[1] = @b/2 - tmps ; not used in this case
    endif
  endif
  
  z = orbit[0]-attract[0]
  complex zprev = z  ; give it some value to avoid compiler warnings
  complex dist = 0   ; ditto
  
  bool continue[2]  ; n_orbits
  continue[0] = true
  continue[1] = true
  
  int n_live    ; how many orbits can bail out until we're done?
  if @giveup=="all done"
    n_live = n_orbits
  else
    n_live = 1
  endif
  
  int step=0    ; which attractor are we checking ? (with @stripes)
  
loop:

  int ix = 0    ; count through the orbits
  repeat
    int found=-1
    if continue[ix]   ;  orbit still live ?
      z = orbit[ix]
      if (@stripes)
         ; stripes: check only one attractor
       dist = z-attract[step]
       if |dist| < @delta
         found=step
         endif
      else
        ; no stripes: check all attractors
        int j=0 ; loop through the attractors
        repeat
          dist = z-attract[j]
          if |dist| < @delta
            found = j
          endif
          j=j+1
        until j>=n_attract || found>=0
      endif ; stripes
      if (found>=0)
        continue[ix] = false
        n_live = n_live - 1
      elseif (!@stripes || step==n_attract-1)
        ; here we step z (this is the core of the formula)
        zprev = z
        z = (3*z^4 - 4*@b*z^3 +  fac2*z^2  +  fac0) / ( 4*z^3 - 6*@b*z^2 + 2*fac2*z + 2*@b)
      orbit[ix] = z
      endif ; found
    endif ; continue
    ix = ix + 1
  until ix >= n_orbits || n_live <= 0
  
  if @stripes
    step = step+1
    if step >= n_attract
      step = 0
    endif
  endif

  ; if we're bailing out, and zmode is != z, adjust the final z (for coloring algorithm)
  if @zmode != "z"
    if n_live <= 0
      if @zmode=="distance"
        z = dist
      elseif @zmode=="step"
        z = z-zprev
      endif ; zmode
    endif ; found
  endif ; zmode
  
bailout:
; -- this isn't the bailout but the continuing condition
;    false = bail; true = continue
  n_live > 0
  
default:
  title = "Newton-Mandel(a), degree-4 polynomial" ; v1.0
  periodicity = 0   ; turn off periodicity checking by default as it interferes with 'stripes'
  
  heading
    caption = "Fractal shape"
  endheading
  param pchoice
    caption="param interpretation"
    enum="halfdist" "square"
    default=0
  endparam
  complex param b
    caption="b"
    default=0
  endparam
  
  heading
    caption = "Bailout"
  endheading
  float param delta
    caption="Bailout delta"
    hint="How close can z get to a fixed point? \
      Actually, this is the square of the distance. \
      Smaller values mean more iterations."
    default=1e-10
    min=0
  endparam
  param giveup
    caption = "Stop when"
    enum = "all done" "any done" "one" "other" "custom orbit start"
    default = 0
  endparam
  complex param z_init
    caption = "orbit start"
    enabled = @giveup == "custom orbit start"
    default = 0   
  endparam
  
  heading
    caption = "Coloring"
  endheading
  bool param stripes
    caption="Staggered iteration"
    hint="Checking this box means you can use one of the striped coloring methods. Use 4 stripes."
    default=false
  endparam
  bool param reorder
    ; enabled = @stripes
    default = false
    caption = "Reorder attractors"
    hint = "Exchange the attractors 1 and 2 (of 0..3). \
      Matters if you use staggered iteration, or in the rare case when two attractors coincide."
  endparam
  param zmode
    caption = "final z"
    enum = "z" "distance" "step"
    default = 1 ; "distance"
  endparam
  
switch:
  type = "N_Poly4_J"
  pchoice = pchoice
  a = #pixel
  b = b
  delta = delta
  reorder = reorder
  stripes = stripes
  zmode = zmode
}
再增添榕坚老师转来代码之前一段外文:
N_Poly4_Ma {
; newton's method applied to P(z) = (z + 1) * (z - 1) * (z - b + a) * (z - b - a)
; Parameter plane (M-set) for all 'a' and a specific 'b'
; F(z) = z - P(z)/P'(z)
; F(z) = (3*z^4 - 4*b*z^3 +  (b^2-a^2-1)*z^2  +  (b^2-a^2)) / ( 4*z^3 - 6*b*z^2 + 2*(b^2-a^2-1)*z + 2*b)
; Attractors: the 4 zeros of P(z): 1,-1, b+a, b-a
; Critical points: the 4 zeros, as well as
;   b/2  +/-  sqrt(3)*sqrt(b^2 + 2*a^2 + 2)/6
; For b=0, this reverts to the 'Newt_fang_xxx' formulas from pwc_convert.ufm
我造M没有成功,但按Z着色,按软体中的a\b定位,造的J集与软体中一致,说明我的计算没错.
Fractal1.png
2010-6-8 13:14

未命名.JPG
2010-6-8 14:33

1068(z+1)(z-1)(z-b+a)(z-b-a)牛顿J集.gsp (25.71 KB)
运行软件,看一下这个属性面板内的值是多少?这个分形基本上作出来了,只是与贴图还差一点。原因是定位值你没有给出来!能把这些formula发给我吗?我的邮箱:422161240@qq.com
Snap4.jpg
2010-6-8 13:42

Snap5.jpg
2010-6-8 14:32
运行软件,看一下这个属性面板内的值是多少?这个分形基本上作出来了,只是与贴图还差一点。原因是定位值你没有给出来!能把这些formula发给我吗?我的邮箱:422161240@qq.com
5077
5078
分形几何 发表于 2010-6-8 13:42
Z值的定位一直没读懂代码,面板中看出一些点的定位。胡老师说的formula是指软件附带的范例吗?
z值定位不在代码里,而在属性面板中:starting(Re)与starting(im)就是迭代初始值的横纵坐标。能把那些范例发到我的邮箱里吗?
返回列表