|
投稿者：しばっち
投稿日：2018年 1月23日(火)21時52分51秒
|
|
|
|
マンデルブロート(ジュリア)を描画します。
スライドバーを19個使用し、パラメータを変化させます(Ver.7.8.1以上)
TABキーを押すと解像度を変更します。
ESCキーで終了します。
左ドラックすると平行移動します。(引っ張る感じになります)
右クリック後に左ドラッグで拡大範囲を指定します。
OPTION ARITHMETIC COMPLEX
LET LEFT=-2.5 !'計算範囲
LET RIGHT=2.5
LET BOTTOM=-2.5
LET TOP=2.5
LET XSIZE=400
LET YSIZE=400
CALL GINIT(XSIZE,YSIZE)
DIM S$(8)
MAT READ S$
DATA "200*200","300*300","400*400","500*500","600*600","700*700","800*800","1000*1000"
LOCATE VALUE NOWAIT(1),RANGE 1 TO 200,AT 30:KS !'計算精度
LOCATE VALUE NOWAIT(2),RANGE -1.5 TO 1.5,AT .5:R1
LOCATE VALUE NOWAIT(3),RANGE -1.5 TO 1.5,AT .5:I1
LOCATE VALUE NOWAIT(4),RANGE -1.5 TO 1.5,AT .5:R2
LOCATE VALUE NOWAIT(5),RANGE -1.5 TO 1.5,AT .5:I2
LOCATE VALUE NOWAIT(6),RANGE -1.5 TO 1.5,AT .5:R3
LOCATE VALUE NOWAIT(7),RANGE -1.5 TO 1.5,AT .5:I3
LOCATE VALUE NOWAIT(8),RANGE -1.5 TO 1.5,AT 0:R4
LOCATE VALUE NOWAIT(9),RANGE -1.5 TO 1.5,AT 0:I4
LOCATE VALUE NOWAIT(10),RANGE -1.5 TO 1.5,AT 0:R5
LOCATE VALUE NOWAIT(11),RANGE -1.5 TO 1.5,AT 0:I5
LOCATE VALUE NOWAIT(12),RANGE -1.5 TO 1.5,AT 0:R6
LOCATE VALUE NOWAIT(13),RANGE -1.5 TO 1.5,AT 0:I6
LOCATE VALUE NOWAIT(14),RANGE -1.5 TO 1.5,AT 0:R7
LOCATE VALUE NOWAIT(15),RANGE -1.5 TO 1.5,AT 0:I7
LOCATE VALUE NOWAIT(16),RANGE -1.5 TO 1.5,AT 0:R8
LOCATE VALUE NOWAIT(17),RANGE -1.5 TO 1.5,AT 0:I8
LOCATE VALUE NOWAIT(18),RANGE -1.5 TO 1.5,AT 0:R9
LOCATE VALUE NOWAIT(19),RANGE -1.5 TO 1.5,AT 0:I9
DO
IF FF<>0 THEN
LOCATE CHOICE(S$) : N
SELECT CASE N
CASE 1
LET XSIZE=200
LET YSIZE=200
CASE 2
LET XSIZE=300
LET YSIZE=300
CASE 3
LET XSIZE=400
LET YSIZE=400
CASE 4
LET XSIZE=500
LET YSIZE=500
CASE 5
LET XSIZE=600
LET YSIZE=600
CASE 6
LET XSIZE=700
LET YSIZE=700
CASE 7
LET XSIZE=800
LET YSIZE=800
CASE 8
LET XSIZE=1000
LET YSIZE=1000
END SELECT
CALL GINIT(XSIZE,YSIZE)
PRINT "画像サイズ";XSIZE;"*";YSIZE
LET FF=0
END IF
LET DX=(RIGHT-LEFT)/XSIZE
LET DY=(TOP-BOTTOM)/YSIZE
SET WINDOW LEFT,RIGHT,BOTTOM,TOP
PRINT "横:";LEFT;"～";RIGHT
PRINT "縦:";BOTTOM;"～";TOP
PRINT
CLEAR
DO
LOCATE VALUE NOWAIT(1):KS
LOCATE VALUE NOWAIT(2):R1
LOCATE VALUE NOWAIT(3):I1
LOCATE VALUE NOWAIT(4):R2
LOCATE VALUE NOWAIT(5):I2
LOCATE VALUE NOWAIT(6):R3
LOCATE VALUE NOWAIT(7):I3
LOCATE VALUE NOWAIT(8):R4
LOCATE VALUE NOWAIT(9):I4
LOCATE VALUE NOWAIT(10):R5
LOCATE VALUE NOWAIT(11):I5
LOCATE VALUE NOWAIT(12):R6
LOCATE VALUE NOWAIT(13):I6
LOCATE VALUE NOWAIT(14):R7
LOCATE VALUE NOWAIT(15):I7
LOCATE VALUE NOWAIT(16):R8
LOCATE VALUE NOWAIT(17):I8
LOCATE VALUE NOWAIT(18):R9
LOCATE VALUE NOWAIT(19):I9
FOR X=LEFT TO RIGHT STEP DX
FOR Y=BOTTOM TO TOP STEP DY
IF GetKeyState(9)<0 THEN !'TABキー
LET FF=1
EXIT DO
ELSEIF GetKeyState(27)<0 THEN !'ESCキー
STOP
END IF
CALL PSET(X,Y,0)
LET Z=COMPLEX(X,Y)
FOR K=1 TO KS
LET Z=COMPLEX(R9,I9)*Z^8+COMPLEX(R8,I8)*Z^7+COMPLEX(R7,I7)*Z^6+COMPLEX(R6,I6)*Z^5+COMPLEX(R5,I5)*Z^4+COMPLEX(R4,I4)*Z^3+COMPLEX(R3,I3)*Z^2+COMPLEX(R2,I2)*Z+COMPLEX(R1,I1)
IF ABS(Z)>2 THEN
CALL PSET(X,Y,MOD(K,7)+1)
EXIT FOR
END IF
NEXT K
MOUSE POLL X1,Y1,L,R
IF R<>0 THEN EXIT DO !'右クリック
IF L<>0 THEN
SET DRAW MODE NOTXOR
SET LINE STYLE 2
DO
PLOT LINES:X1,Y1;X2,Y2
PLOT LINES:X1,Y1;X2,Y2
MOUSE POLL X2,Y2,LL,RR
LOOP WHILE LL<>0 !'左ドラッグ
LET RIGHT=RIGHT-(X2-X1) !'平行移動
LET LEFT=LEFT-(X2-X1)
LET TOP=TOP-(Y2-Y1)
LET BOTTOM=BOTTOM-(Y2-Y1)
PRINT "移動量 X:";X2-X1;"Y:";Y2-Y1
SET DRAW MODE OVERWRITE
EXIT DO
END IF
NEXT Y
NEXT X
LOOP
IF R<>0 THEN
PAUSE "拡大する範囲を指定してください"
CALL GETSQUARE(LEFT,TOP,RIGHT,BOTTOM) !'左ドラッグで範囲指定
IF LEFT=RIGHT THEN EXIT DO
END IF
LOOP
PRINT "終了します"
END
EXTERNAL SUB PSET(X,Y,C)
OPTION ARITHMETIC COMPLEX
SET POINT COLOR C
PLOT POINTS:X,Y
END SUB
EXTERNAL SUB GINIT(XSIZE,YSIZE)
OPTION ARITHMETIC COMPLEX
SET BITMAP SIZE XSIZE,YSIZE
!'SET WINDOW 0,XSIZE-1,YSIZE-1,0
SET POINT STYLE 1
SET COLOR MODE "REGULAR"
SET COLOR MIX(0) 0,0,0
SET COLOR MIX(1) 0,0,1
SET COLOR MIX(2) 1,0,0
SET COLOR MIX(3) 1,0,1
SET COLOR MIX(4) 0,1,0
SET COLOR MIX(5) 0,1,1
SET COLOR MIX(6) 1,1,0
SET COLOR MIX(7) 1,1,1
CLEAR
END SUB
EXTERNAL SUB GETSQUARE(L,T,R,B)
OPTION ARITHMETIC COMPLEX
ASK LINE STYLE LSTYLE
SET DRAW MODE NOTXOR
SET LINE STYLE 2
DO
MOUSE POLL L,T,I,J
LOOP WHILE I=0
LET L0=L
LET T0=T
LET R0=L0
LET B0=T0
PLOT LINES:L0,T0;L0,B0;R0,B0;R0,T0;L0,T0
DO WHILE I=1
MOUSE POLL R,B,I,J
LET W=R-L
LET H=T-B
IF ABS(H)<ABS(W) THEN
LET B=T-SGN(H)*ABS(W)
ELSE
LET R=L+SGN(W)*ABS(H)
END IF
IF L0<>L OR R0<>R OR B0<>B OR T0<>T THEN
PLOT LINES:L0,T0;L0,B0;R0,B0;R0,T0;L0,T0
PLOT LINES:L,T;L,B;R,B;R,T;L,T
LET L0=L
LET T0=T
LET R0=R
LET B0=B
END IF
LOOP
WAIT DELAY 1
PLOT LINES:L,T;L,B;R,B;R,T;L,T
SET DRAW MODE OVERWRITE
SET LINE STYLE LSTYLE
IF L>R THEN SWAP L,R
IF B>T THEN SWAP B,T
END SUB
|
|
|
|
投稿者：yoshipyuta
投稿日：2018年 1月24日(水)16時20分53秒
|
|
|
|
> No.4478[元記事へ]
しばっちさんへのお返事です。
> マンデルブロート(ジュリア)を描画します。
> スライドバーを19個使用し、パラメータを変化させます(Ver.7.8.1以上)
> TABキーを押すと解像度を変更します。
> ESCキーで終了します。
> 左ドラックすると平行移動します。(引っ張る感じになります)
> 右クリック後に左ドラッグで拡大範囲を指定します。
>
しばっちさんのJuria②プログラムの内容もわからないKidsが早速、下記の２つのJuliaで類似部分があることを理解したようです。
LET Z=z^3+COMPLEX(1.503,-0.8046)*z^2 下図１
LET Z=z^2+COMPLEX(-0.122561,0.744862) 下図２
Kidsは現在Newtonフラクタルの研究で頭が一杯でJuliaには手が出ないようです。スライドバー１９個とは？４個のみ表示ですが。またRe(R1,R2)、Im(I1,I2)のバーが２つありますが、意味がわからないKidsです。c点を固定したらバーは動かない位しかわかりません。
このｃ点導入（可変・固定）をニュートン・フラクタルに導入したらKidsはパニックでしょう。LET XX=X-FUNC(X)/DIFF(X,1)+cなのか。
|
|
|
|
投稿者：yoshipyuta
投稿日：2018年 1月26日(金)17時47分28秒
|
|
|
|
> No.4478[元記事へ]
しばっちさんへのお返事です。
> マンデルブロート(ジュリア)を描画します。
> スライドバーを19個使用し、パラメータを変化させます(Ver.7.8.1以上)
> TABキーを押すと解像度を変更します。
> ESCキーで終了します。
> 左ドラックすると平行移動します。(引っ張る感じになります)
> 右クリック後に左ドラッグで拡大範囲を指定します。
>
このプログラムを用いてKidsは終日、２次Julia集合の探索をしています。Newton・フラクタルはしばしの休みです。赤いサンマルコ寺院と呼ばれるJulia内部に興味を持ったようです。
しばっち・Juriaプログラム②でサンマルコ寺院の中身J(z^2+complex(-0.745432,0.11321))を見ますと、下図１になります。中央部分の小さな渦の会合が見えにくい状況です。
拡大すると明瞭に渦の会合が観察されます。この描画可能でしょうか。
|
|
|
|
投稿者：しばっち
投稿日：2018年 1月29日(月)20時10分36秒
|
|
|
|
> No.4481[元記事へ]
yoshipyutaさんへのお返事です。
>
> 拡大すると明瞭に渦の会合が観察されます。この描画可能でしょうか。
上記サイト内にソースリストがありますので、とりあえず移植を試みてみました。
しかしながら、この言語についての知識はなく、不明な点もありパラメータも不明なので
サイト内の画像のように中心部分の拡大を2～3度繰り返すとサイト内の画像とはなんだか違う画像になりますが
サイトの内容を理解しているわけでもないのであしからず。
OPTION ARITHMETIC COMPLEX
OPTION BASE 0
LET XSIZE=1000 !'画像サイズ
LET YSIZE=1000
CALL GINIT(XSIZE,YSIZE)
DIM COUNT(XSIZE,YSIZE),FL(100)
LET KS=800
LET C=COMPLEX(-0.74543,0.11301)
LET R=(1+SQR(1+4*ABS(C)))/2
LET LEFT=-R
LET RIGHT=R
LET TOP=R
LET BOTTOM=-R
SET TEXT BACKGROUND "OPAQUE"
SET TEXT JUSTIFY "LEFT" , "TOP"
DO
SET COLOR COLORINDEX(1,1,1)
SET TEXT HEIGHT (TOP-BOTTOM)/5
MAT FL=ZER
LET KMAX=0
LET KMIN=0
SET WINDOW LEFT,RIGHT,BOTTOM,TOP
CLEAR
LET DX=(RIGHT-LEFT)/XSIZE
LET DY=(TOP-BOTTOM)/YSIZE
FOR ZR=LEFT TO RIGHT STEP DX
LET P=INT((ZR-LEFT)/(RIGHT-LEFT)*100)
IF FL(INT(P/5))=0 THEN
PLOT TEXT ,AT LEFT,TOP:STR$(P)&"%"
LET FL(INT(P/5))=1
END IF
FOR ZI=BOTTOM TO TOP STEP DY
LET Z=COMPLEX(ZR,ZI)
FOR K=1 TO KS
LET Z=Z*Z+C
IF ABS(Z)>R THEN EXIT FOR
NEXT K
LET X=PIXELX(ZR)
LET Y=PIXELY(ZI)
LET COUNT(X,Y)=K
LET KMAX=MAX(KMAX,K)
NEXT ZI
NEXT ZR
FOR ZR=LEFT TO RIGHT STEP DX
FOR ZI=BOTTOM TO TOP STEP DY
LET X=PIXELX(ZR)
LET Y=PIXELY(ZI)
LET K=COUNT(X,Y)
LET Z=COMPLEX(ZR,ZI)
CALL PSET(ZR,ZI,255*(K-KMIN)/(KMAX-KMIN),255*(1-(K-KMIN)/(KMAX-KMIN)),255*MIN(1,ABS(Z)/R))
NEXT ZI
NEXT ZR
PAUSE "拡大する範囲を指定してください"
CALL GETSQUARE(LEFT,TOP,RIGHT,BOTTOM)
IF LEFT=RIGHT THEN EXIT DO
LOOP
END
EXTERNAL SUB GINIT(XSIZE,YSIZE)
OPTION ARITHMETIC COMPLEX
SET BITMAP SIZE XSIZE,YSIZE
SET COLOR MIX(0) 0,0,0
SET COLOR MODE "NATIVE"
CLEAR
SET POINT STYLE 1
SET WINDOW 0,XSIZE-1,YSIZE-1,0
END SUB
EXTERNAL SUB PSET(X,Y,R,G,B)
OPTION ARITHMETIC COMPLEX
LET RR=MIN(255,MAX(0,INT(R)))
LET GG=MIN(255,MAX(0,INT(G)))
LET BB=MIN(255,MAX(0,INT(B)))
SET COLOR COLORINDEX(RR/255,GG/255,BB/255)
PLOT POINTS:X,Y
END SUB
EXTERNAL SUB GETSQUARE(L,T,R,B)
OPTION ARITHMETIC COMPLEX
ASK LINE STYLE LSTYLE
SET DRAW MODE NOTXOR
SET LINE STYLE 2
DO
MOUSE POLL L,T,I,J
LOOP WHILE I=0
LET L0=L
LET T0=T
LET R0=L0
LET B0=T0
PLOT LINES:L0,T0;L0,B0;R0,B0;R0,T0;L0,T0
DO WHILE I=1
MOUSE POLL R,B,I,J
LET W=R-L
LET H=T-B
IF ABS(H)<ABS(W) THEN
LET B=T-SGN(H)*ABS(W)
ELSE
LET R=L+SGN(W)*ABS(H)
END IF
IF L0<>L OR R0<>R OR B0<>B OR T0<>T THEN
PLOT LINES:L0,T0;L0,B0;R0,B0;R0,T0;L0,T0
PLOT LINES:L,T;L,B;R,B;R,T;L,T
LET L0=L
LET T0=T
LET R0=R
LET B0=B
END IF
LOOP
WAIT DELAY 1
PLOT LINES:L,T;L,B;R,B;R,T;L,T
SET DRAW MODE OVERWRITE
SET LINE STYLE LSTYLE
IF L>R THEN SWAP L,R
IF B>T THEN SWAP B,T
END SUB
|
|
|
|
投稿者：yoshipyuta
投稿日：2018年 1月31日(水)16時54分42秒
|
|
|
|
> No.4485[元記事へ]
しばっちさんへのお返事です。
>
> 上記サイト内にソースリストがありますので、とりあえず移植を試みてみました。
> しかしながら、この言語についての知識はなく、不明な点もありパラメータも不明なので
> サイト内の画像のように中心部分の拡大を2～3度繰り返すとサイト内の画像とはなんだか違う画像になりますが
> サイトの内容を理解しているわけでもないのであしからず。
しばっちさんのプログラムを理解しているわけではありませんが、描画は成功と思います。％描画時間待ちはKidsをワクワクさせる良きアイディアですね。今、Kidsはしばっち・Newton プログラム②を用いてf(z)=cos(z)のフラクタルを作成しています。アトラクターが実数軸上に周期πで並んでいます。
プログラムを手直しして、解を４つほど並べて解決と思ったらしい。
LET Z1=COMPLEX(PI/2,0)
LET Z2=COMPLEX(-PI/2,0)
LET z3=COMPLEX(3*PI/2,0)
LET z4=COMPLEX(-3*PI/2,0)
EXTERNAL FUNCTION FUNC(X)
OPTION ARITHMETIC COMPLEX
LET FUNC=ccos(x)
END FUNCTION
ところが原点に近づくほど遠くのアトラクターの収束領域になるので、下図のようにのような黒領域になって正しく表示されません。すべての解nπ-π/2（n:整数）をLet文で並べるわけにもいかないので、どのように解決できますか？
|
|
|
|
投稿者：しばっち
投稿日：2018年 2月 1日(木)21時54分19秒
|
|
|
|
> No.4486[元記事へ]
yoshipyutaさんへのお返事です。
> プログラムを手直しして、解を４つほど並べて解決と思ったらしい。
>
> LET Z1=COMPLEX(PI/2,0)
> LET Z2=COMPLEX(-PI/2,0)
> LET z3=COMPLEX(3*PI/2,0)
> LET z4=COMPLEX(-3*PI/2,0)
>
> EXTERNAL FUNCTION FUNC(X)
> OPTION ARITHMETIC COMPLEX
> LET FUNC=ccos(x)
> END FUNCTION
>
> ところが原点に近づくほど遠くのアトラクターの収束領域になるので、下図のようにのような黒領域になって正しく表示されません。すべての解nπ-π/2（n:整数）をLet文で並べるわけにもいかないので、どのように解決できますか？
まず、COS(x)=0 の一般解 x=π/2+nπなのでxからπ/2を引いたものはπの倍数、つまりπで割り切れるので
MOD(X-PI/2,PI)=0
という式が出てきます。
但し、この時のX は複素数なのでこのままでは使えません。RE(X)やIM(X)関数を使い、許容誤差を含めて
IF ABS(IM(X))<EPS AND ABS(MOD(RE(X)-PI/2,PI))<EPS THEN
とすれば解決するかと思います。
なお、色については、一般解のNを使うのがいいと思います。
N=(RE(X)-PI/2)/PI
このままではNが整数にならないかもしれないので
N=INT((RE(X)-PI/2)/PI+.5)
と四捨五入したほうがいいかもしれません。
このNの範囲については全くわかりませんが、例えば100色ならMOD(N,100)+1とすれば1～100までの色番号を使用します。
ただし、Nは負数の場合もあるはずなので絶対値をとって MOD(ABS(N),100)+1 とするか
補数を足して MOD(N+100,100)+1 のようにすればいいかと思います。
|
|
|
|
投稿者：yoshipyuta
投稿日：2018年 2月 2日(金)10時24分53秒
|
|
|
|
しばっちさんへのお返事です。
>
> まず、COS(x)=0 の一般解 x=π/2+nπなのでxからπ/2を引いたものはπの倍数、つまりπで割り切れるので
>
> MOD(X-PI/2,PI)=0
>
> という式が出てきます。
> 但し、この時のX は複素数なのでこのままでは使えません。RE(X)やIM(X)関数を使い、許容誤差を含めて
>
> IF ABS(IM(X))<EPS AND ABS(MOD(RE(X)-PI/2,PI))<EPS THEN
>
> とすれば解決するかと思います。
懇切なご指導ありがとうございます。わがKidsは最新のしばっち・Juliaプログラム②に夢中です。
ご返事中の色の取り扱いが分からなくて（プログラム上で）、ガムシャラKidsはどうも安易に解を左右に１６個並べて、この問題をスルーするらしいのです（下図、未完成）。
まあ、それはそれでプログラムを触ると様々なことが可能になるということが分かれば、ひとつの成果であります。
|
|
|
|
投稿者：yoshipyuta
投稿日：2018年 2月 2日(金)10時46分50秒
|
|
|
|
> No.4485[元記事へ]
しばっちさんへのお返事です。
> 上記サイト内にソースリストがありますので、とりあえず移植を試みてみました。
> しかしながら、この言語についての知識はなく、不明な点もありパラメータも不明なので
> サイト内の画像のように中心部分の拡大を2～3度繰り返すとサイト内の画像とはなんだか違う画像になりますが
> サイトの内容を理解しているわけでもないのであしからず。
しばっちさん、このプログラムはKidsにはかなり面白いようです。朝から晩までDR.Michel BeckerとともにJuliaの美しさに堪能のようですが探索はかなりむずかしそうです。
まずはC=complex(0.100,0.0700),LET Z=(Z^5-C*Z+0.09)/Z^3です。正に宇宙におけるRandom Walkです（上段の図）。但しプログラム出力図は見えにくいのでPhotoshopで色を変えています。
これをC=complex(0.０100,0.070２),LET Z=(Z^5-C*Z+0.09)/Z^3では何やら暗青色の海が見えてきます（中断の図）。
Tune-upすると超新星星雲のなかになんとDouadyのウサギが出現します。
|
|
|
|
投稿者：yoshipyuta
投稿日：2018年 2月 8日(木)00時10分52秒
|
|
|
|
> No.4485[元記事へ]
しばっちさんへの質問です。
>
> >
> > 拡大すると明瞭に渦の会合が観察されます。この描画可能でしょうか。
>
> 上記サイト内にソースリストがありますので、とりあえず移植を試みてみました。
> しかしながら、この言語についての知識はなく、不明な点もありパラメータも不明なので
> サイト内の画像のように中心部分の拡大を2～3度繰り返すとサイト内の画像とはなんだか違う画像になりますが
> サイトの内容を理解しているわけでもないのであしからず。
移植に成功したJuliaプログラムでKidsが困っています。大半の関数形では正しく表示されるのですが、中には手ごわいものがあります。どうも次数が上がるとパターンがはっきりとでないのです。たとえば次の関数のジュリア集合の場合です。
LET Z=z^7+z^3+0.01/z
で上段の図が得られます。かなり薄くパターンが出ており、これはPhotoshopでもお手上げです。
Dr.MihaelBeckerによると下段のような明確なパターンであると主張しています。精度KSの問題でもないようです。白黒ラインでラインを強く出す方法は考えられませんか？
|
|
|
|
投稿者：しばっち
投稿日：2018年 2月 8日(木)20時18分45秒
|
|
|
|
> No.4496[元記事へ]
yoshipyutaさんへのお返事です。
> Dr.MihaelBeckerによると下段のような明確なパターンであると主張しています。精度KSの問題でもないようです。白黒ラインでラインを強く出す方法は考えられませんか？
一般的な方法ではないのなら、特殊な方法で収束判定しているとしか思えません。
どこかにその方法が記されていない限り、理論も理解していない素人がその方法を思いつくことはありえません。
素人の私には残念ながら分かりません。yoshipyutaさんの方が、このjulia集合に関して詳しいのではないでしょうか。
ネット上にその回答があるのかは分かりませんが、ネットで調べるなりDr.MihaelBeckerの文献等を調べるかしない限り、
その方法は分からないと思います。
|
|
|
|
投稿者：yoshipyuta
投稿日：2018年 2月 9日(金)10時27分49秒
|
|
|
|
> No.4500[元記事へ]
しばっちさんへのお返事です。
> yoshipyutaさんへのお返事です。
>
> > Dr.MihaelBeckerによると下段のような明確なパターンであると主張しています。精度KSの問題でもないようです。白黒ラインでラインを強く出す方法は考えられませんか？
>
> 一般的な方法ではないのなら、特殊な方法で収束判定しているとしか思えません。
> どこかにその方法が記されていない限り、理論も理解していない素人がその方法を思いつくことはありえません。
> 素人の私には残念ながら分かりません。yoshipyutaさんの方が、このjulia集合に関して詳しいのではないでしょうか。
> ネット上にその回答があるのかは分かりませんが、ネットで調べるなりDr.MihaelBeckerの文献等を調べるかしない限り、その方法は分からないと思います。
ご迷惑をおかけしています。この問題は恐らく収束法判定の問題ではなくてTechnial Colouringの選択か？と思います。Dr.Mihael BeckerはソースCodeを公開していないので分かりませんが、点Plotではなく、J境界を陰関数表示している？まさかトレースではないと思いますが。
類似問題は実はNewton Fractalでもあります。しばっち・Newtonプログラム②でf(z)=z^3-2*z+2,R = 0.98を描くと上段の図が得られます。骨格ラインのみを強調すると（どうやって？）下段の図となります。
次のサイトのDistance Shadingという手法かな？と思ったりします。Codeも公開しています。
function colourDistance
parameters: h, s, v: real: components of the HSV colour space;
distance: real: the distance estimate calculated during iteration;
pixelWidth: real: the width represented by a pixel on the screen;
boundaryFraction: real: value > 0 used to control prominence of the boundary
let:
dScale = log2(distance / pixelWidth),
factor = 0;
if (dScale > 0) {
factor = 1;
}
else if (dScale > -boundaryFraction) {
factor = (boundaryFraction + dScale) / boundaryFraction;
}
// Darken v when factor < 1
return hsv2rgb(h, s, v * factor);
私やKidsにはさっぱり分かりませんので、是非チャレンジされて、お示しください。次のRoot Finding法も興味深い。
Finding Complex Roots of a Function
Basin colouring for Newtonian fractals requires us to know the roots of the function before we generate the fractal.
We can also anticipate a future requirement of identifying critical points of a function for generating Mandelbrot sets. Therefore, we need a generic root finding algorithm. This algorithm uses - what else? - Newton's method.
definition rootObject
properties: x: the real component of z;
y: the imaginary component of z;
r: the modulus of z;
theta: the angle of z
function findRoots
parameters: f: a pointer to a function that takes a complex number, z, and computes f(z);
fprime: a pointer to a function that takes a complex number, z, and computes f?(z)
constants: XMAX, XMIN, YMAX, YMIN: the limits of the complex plane to be tested
let:
eV = 1e-8,
// For checking suspiciously small real / imaginary components
delta = eV * 100,
// For comparing roots found
tolerance = delta * 100,
rSpan = XMAX - XMIN,
iSpan = YMAX - YMIN,
// If we don't find a root after n iterations, give up and move on
iterations = 100,
// We'll test this many points in the x- and y- axes
points = 30,
rInt = rSpan / points,
iInt = iSpan / points,
roots = [];
// Test a reasonably large sample of points in the plane; assuming that we only have to
// do this when we regenerate a fractal, we can use a conservatively large value to minimize
// the chances of failing to find distinct roots that are close to each other
for (i = 0..points) {
for(j = 0..points) {
let:
x = XMIN + i * rInt,
y = YMIN + j * iInt,
diff = 1,
iter = 0;
while (iter < iterations AND diff > eV) {
// Get the numerator and denominator
[nomX, nomY] = f(x, y);
[denomX, denomY] = fprime(x, y);
// Divide f(z) by f'(z) - multiply through by conjugate of denominator
tmp = nomX * denomX + nomY * denomY;
nomY = nomY * denomX - nomX * denomY;
nomX = tmp;
denomX = denomX * denomX + denomY * denomY;
if (!denomX) {
// Try another number
break;
}
// Now we can divide
_x = x - nomX / denomX;
_y = y - nomY / denomX;
// Get the differences from last round
diffX = x - _x;
diffY = y - _y;
diff = sqrt(diffX * diffX + diffY * diffY);
// Set up the next guess
x = _x;
y = _y;
++iter;
}
if (diff < eV) {
let:
newRoot = TRUE;
// Zero suspiciously small real / imaginary components so that the roots can
// be sorted deterministically
if (abs(x) < delta) {
x = 0;
}
if (abs(y) < delta) {
y = 0;
}
for (k = 0..LENGTH(roots)) {
if (complexCompare(x, y, roots[k].x, roots[k].y, tolerance)) {
newRoot = FALSE;
break;
}
}
if (newRoot) {
push(roots, rootObject{ x : x, y : y, r : sqrt(x * x + y * y), theta : atan2(y, x) });
}
}
}
}
SORT roots;
return roots;
|
|
|
|
投稿者：yoshipyuta
投稿日：2018年 2月 9日(金)11時30分47秒
|
|
|
|
> No.4500[元記事へ]
しばっちさんへのお返事です。
> yoshipyutaさんへのお返事です。
>
> > Dr.MihaelBeckerによると下段のような明確なパターンであると主張しています。精度KSの問題でもないようです。白黒ラインでラインを強く出す方法は考えられませんか？
>
> 一般的な方法ではないのなら、特殊な方法で収束判定しているとしか思えません。
> どこかにその方法が記されていない限り、理論も理解していない素人がその方法を思いつくことはありえません。
> 素人の私には残念ながら分かりません。yoshipyutaさんの方が、このjulia集合に関して詳しいのではないでしょうか。
> ネット上にその回答があるのかは分かりませんが、ネットで調べるなりDr.MihaelBeckerの文献等を調べるかしない限り、その方法は分からないと思います。
カリンさんが２０１０年に公開されている以下のCodeにも、参考になるかも知れません。なんとNewton Fractalの境界を光り輝くラインで描画しています。是非、しばっち流に改善してFractal研究に使わせて下さい。
!１０進ＢＡＳＩＣ専用
!任意の範囲をドラッグ＋正方形クリックで拡大操作
!
!カラーコード 193色(RGB)
data 000000
data FF0000,FF0800,FF1000,FF1800,FF2000,FF2800,FF3000,FF3800,FF4000,FF4800,FF5000,FF5800,FF6000,FF6800,FF7000,FF7800
data FF8000,FF8800,FF9000,FF9800,FFA000,FFA800,FFB000,FFB800,FFC000,FFC800,FFD000,FFD800,FFE000,FFE800,FFF000,FFF800
data FFFF00,F8FF00,F0FF00,E8FF00,E0FF00,D8FF00,D0FF00,C8FF00,C0FF00,B8FF00,B0FF00,A8FF00,A0FF00,98FF00,90FF00,88FF00
data 80FF00,78FF00,70FF00,68FF00,60FF00,58FF00,50FF00,48FF00,40FF00,38FF00,30FF00,28FF00,20FF00,18FF00,10FF00,08FF00
data 00FF00,00FF08,00FF10,00FF18,00FF20,00FF28,00FF30,00FF38,00FF40,00FF48,00FF50,00FF58,00FF60,00FF68,00FF70,00FF78
data 00FF80,00FF88,00FF90,00FF98,00FFA0,00FFA8,00FFB0,00FFB8,00FFC0,00FFC8,00FFD0,00FFD8,00FFE0,00FFE8,00FFF0,00FFF8
data 00FFFF,00F8FF,00F0FF,00E8FF,00E0FF,00D8FF,00D0FF,00C8FF,00C0FF,00B8FF,00B0FF,00A8FF,00A0FF,0098FF,0090FF,0088FF
data 0080FF,0078FF,0070FF,0068FF,0060FF,0058FF,0050FF,0048FF,0040FF,0038FF,0030FF,0028FF,0020FF,0018FF,0010FF,0008FF
data 0000FF,0800FF,1000FF,1800FF,2000FF,2800FF,3000FF,3800FF,4000FF,4800FF,5000FF,5800FF,6000FF,6800FF,7000FF,7800FF
data 8000FF,8800FF,9000FF,9800FF,A000FF,A800FF,B000FF,B800FF,C000FF,C800FF,D000FF,D800FF,E000FF,E800FF,F000FF,F800FF
data FF00FF,FF00F8,FF00F0,FF00E8,FF00E0,FF00D8,FF00D0,FF00C8,FF00C0,FF00B8,FF00B0,FF00A8,FF00A0,FF0098,FF0090,FF0088
data FF0080,FF0078,FF0070,FF0068,FF0060,FF0058,FF0050,FF0048,FF0040,FF0038,FF0030,FF0028,FF0020,FF0018,FF0010,FF0008
declare sub shift_sfn
declare sub fpr_sfn
declare sub select_area
declare sub updateparam
declare sub draw_area
declare sub tdraw_area
declare sub flood_area
declare sub draw_point
declare sub pointdraw
declare sub draw_side
declare sub draw_inside
declare sub draw_srline
declare sub printdata
set point style 1
set color mode "native"
set draw mode explicit
option base 0
option arithmetic complex
declare function fcos !複素数拡張 cos() 関数
declare function fsin !複素数拡張 sin() 関数
let cnum = 193 !色数(64×3＋1)
let kcn = cnum-1
let kcolor = 1/255 !RGB指数変換定数(24bits)
dim ind(kcn) !色指標用配列
let undolimit = 100
let ul = undolimit
let lv = 0
dim graphics(ul)
dim realpartmin(ul),imaginarypartmin(ul)
dim realpartmax(ul),imaginarypartmax(ul)
dim repeatlimit(ul),uupdate(ul)
!━━━━━━━━━━━━━━━━━━━━━━━━━━━カラーコードを色指標に変換━━━━━━━━━━━━━━━━━━━━
for g1 = 0 to kcn
read cl$
let r = bval(cl$(1:2),16)*kcolor
let g = bval(cl$(3:4),16)*kcolor
let b = bval(cl$(5:6),16)*kcolor
let ind(g1) = colorindex(r,g,b)
next g1
!━━━━━━━━━━━━━━━━━━━━━━━━━━━境界走査―境界描画と境界判定━━━━━━━━━━━━━━━━━━━
sub draw_side
for gx = mgx to mgxp
let gxf = gx1 + krg*gx
let gy = mgy
let gyf = gy1 + krg*gy
call draw_point!<--------------------------------draw_point -> pointdraw で判定
let gy = mgyp
let gyf = gy1 + krg*gy
call draw_point!<--------------------------------draw_point -> pointdraw で判定
next gx
for gy = mgy1 to mgym
let gyf = gy1 + krg*gy
let gx = mgx
let gxf = gx1 + krg*gx
call draw_point!<--------------------------------draw_point -> pointdraw で判定
let gx = mgxp
let gxf = gx1 + krg*gx
call draw_point!<--------------------------------draw_point -> pointdraw で判定
next gy
end sub
!━━━━━━━━━━━━━━━━━━━━━━━━━━━境界走査―内部の描画━━━━━━━━━━━━━━━━━━━━━━━
sub draw_inside(grx,gry)
let mgy = mgyt*grc
let mgyp = mgy+gry-1
let mgy1 = mgy+1
let mgym = mgyp-1
let sign = 1
call draw_side!<------------------------------------境界の計算、描画と判定
if sign = 1 then!-----------------------------------境界点がすべて非収束か非発散の場合内部も同様とみなす
call flood_area(mgx,mgy,mgx+grx-1,mgy+gry-1,ind(0))
else!<----------------------------------------------sign = 0 のときの内部計算と描画
for gx = mgx1 to mgxm
let gxf = gx1 + krg*gx
for gy = mgy1 to mgym
let gyf = gy1 + krg*gy
call draw_point
next gy
next gx
end if
end sub
!━━━━━━━━━━━━━━━━━━━━━━━━━━━境界走査―各列描画━━━━━━━━━━━━━━━━━━━━━━━━
sub draw_srline(grx1,gry1,gry2)
set draw mode hidden
let mgx = mgxt*grc
let mgxp = mgx+grx1-1
let mgx1 = mgx+1
let mgxm = mgxp-1
call flood_area(mgx,0,mgxp,wx,colorindex(1,1,1))!<---塗りつぶしを上手くするために描画する部分を消す
for mgyt = 0 to div
call draw_inside(grx1,gry1)
next mgyt
!---------------------------------------------------画像サイズが境界サイズで割り切れない場合の余り部分の描画
if gmod = 1 then
call draw_inside(grx1,gry2)
end if
set draw mode explicit
end sub
!━━━━━━━━━━━━━━━━━━━━━━━━━━━収束回数による色分けと点の描画━━━━━━━━━━━━━━━━━━
sub pointdraw
!------------------------------------------------------収束しないまたは振動の場合
if k = pn then
let nb = 0
else
!---------------------------------------------------収束または発散した場合 k が収束か発散までの反復回数
let nb = mod(k,kcn) + 1
let sign = 0!<-----------------------------------収束か発散した点が存在したとき sign = 0 として内部計算も行う
end if
set color ind(nb)
plot points :gx,gy
end sub
!━━━━━━━━━━━━━━━━━━━━━━━━━━━数列 z(n) の収束、振動及び発散の判定と描画━━━━━━━━━━━━
sub draw_point
when exception in
let c = complex(gxf,gyf)
call fpr_sfn
let zd = 1
let k = 1
!------------------------------------------------収束、発散の判定文
do while k =< num and abs(zd) > mindt and abs(zd) < maxdt
let gz = z
call shift_sfn
let zd = z-gz
let k = k + 1
loop
!------------------------------------------------
call pointdraw
use
call pointdraw
end when
end sub
!━━━━━━━━━━━━━━━━━━━━━━━━━━━拡大範囲選択の四角の描画━━━━━━━━━━━━━━━━━━━━━
sub tdraw_area
call draw_area(x1,y1,x1+trangex,y1+trangey,colorindex(0,0,1))
end sub
!━━━━━━━━━━━━━━━━━━━━━━━━━━━描画データの表示━━━━━━━━━━━━━━━━━━━━━━━━━
sub printdata
print "実部 "&str$(gx1)&"≦Re≦"&str$(gx2)
print "虚部 "&str$(gy1)&"≦Im≦"&str$(gy2)
print "幅 "&str$(xrange)
print "中心座標 ("&str$(lpx1)&","&str$(lpy1)&")"
print "色数 "&str$(cnum)&"色"
print "画像サイズ "&str$(wx)&"pixel×"&str$(wx)&"pixel"
print "境界走査幅 "&str$(grc)&"pixel"
print "計算回数 "&str$(num)&"回"
end sub
!━━━━━━━━━━━━━━━━━━━━━━━━━━━拡大部分の選択と画像保存、終了の操作をする━━━━━━━━━━━━
sub select_area
let left = 0
let right = 0
do while left = 0 and right = 0
mouse poll dx1,dy1,left,right
loop
if right = 1 then
if lu <> 1 then
let lu = lu - 1
set draw mode hidden
gload "~table"&str$(lu)&".bmp"
set draw mode explicit
set draw mode notxor
let wx = graphics(lu)
set window 0,wx-1,0,wx-1
end if
end if
if left = 1 then
set draw mode notxor
let x1 = dx1
let y1 = dy1
let x2 = x1
let y2 = y1
let update = 1
let trangex = 0
let trangey = 0
call tdraw_area
set draw mode hidden
do while left = 1
mouse poll x3,y3,left,right
if x3 <> x2 or y3 <> y2 then
set draw mode hidden
call tdraw_area
let trange = max(abs(x1-x3),abs(y1-y3))
let trangex = sgn(x3-x1)*trange
let trangey = sgn(y3-y1)*trange
set draw mode explicit
call tdraw_area
let x2 = x3
let y2 = y3
end if
loop
if x1 = x2 and y1 = y2 then
if lu <> lv then
let lu = lu + 1
set draw mode hidden
gload "~table"&str$(lu)&".bmp"
set draw mode explicit
set draw mode notxor
let wx = graphics(lu)
set window 0,wx-1,0,wx-1
end if
exit sub
end if
let left = 0
let right = 0
do while left = 0 and right = 0
mouse poll x4,y4,left,right
loop
let gx41 = sgn(x4-x1)
let tx41 = sgn(x4-x1-trangex)
let gy41 = sgn(y4-y1)
let ty41 = sgn(y4-y1-trangey)
if left = 1 then
if gx41 <> tx41 and gy41 <> ty41 then
let cont = 1
else
let cont = 0
call tdraw_area
end if
end if
if right = 1 then
if gx41 <> tx41 and gy41 <> ty41 then
call updateparam
let cont = 1
else
let cont = 0
call tdraw_area
do while right = 1
mouse poll x5,y5,left,right
loop
end if
end if
end if
end sub
!━━━━━━━━━━━━━━━━━━━━━━━━━━━描画範囲変更に伴う各値の更新━━━━━━━━━━━━━━━━━━━
sub updateparam
input prompt "繰り返しの回数、画像サイズ、境界走査サイズ":num,wx,grc
let cont = 1
if mod(wx,grc) <> 0 then
let gmod = 1
let grm = mod(wx,grc)
else
let gmod = 0
end if
let uupdate(lu) = 1
end sub
!━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━
!------------------------------------------------------初期値
let wx = 556 !グラフィックス画面の横画素数
let wy = wx !グラフィックス画面の縦画素数
let num = 256 !反復計算回数の上限
let lpx1_ = 0.0000001!グラフィックス画面の中心の実座標
let lpy1_ = 0.0000001!グラフィックス画面の中心の虚座標
let xrange_ = 8 !グラフィックス画面に割り当てる複素平面の範囲（初回描画時）
let mindt = 1e-10 !|z(n)-z(n+1)|=<mindtの時収束と判定する
let maxdt = 1e+10 !|z(n)-z(n+1)|>=maxdtの時発散と判定する
let grc = 10 !描画中再描画するまで計算するライン数
let grt = 1 !描画ラインカウント変数
let grm = 0
let lpx1 = lpx1_
let lpy1 = lpy1_
let xrange = xrange_
let zd = 1 !|z(n)-z(n+1)|を代入する変数
let gmod = 0
set bitmap size wx,wx
let gx1 = lpx1 - xrange/2
let gy1 = lpy1 - xrange/2
let gx2 = lpx1 + xrange/2
let gy2 = lpy1 + xrange/2
let count = 0
let wx_ = wx+1
let i = complex(0,1) !虚数単位
if mod(wx,grc) <> 0 then
let gmod = 1
let grm = mod(wx,grc)
end if
call printdata
!━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━
do while count <> 1
set draw mode overwrite
if wx <> wx_ then set bitmap size wx,wx
set window 0,wx-1,0,wx-1
let wx_ = wx
let krg = xrange/wx
let pn = num + 1
let start_time = time
let div = int(wx/grc)-1
let lv = lv + 1
let graphics(lv) = wx
let realpartmin(lv) = gx1
let imaginarypartmin(lv) = gy1
let realpartmax(lv) = gx2
let imaginarypartmax(lv) = gy2
let repeatlimit(lv) = num
!---------------------------------------------------描画総括
for mgxt = 0 to div
call draw_srline(grc,grc,grm)
next mgxt
if gmod = 1 then
call draw_srline(grm,grc,grm)
end if
!---------------------------------------------------拡大操作
set draw mode explicit
set draw mode notxor
print "計算時間 "&str$(round(time-start_time,3))&"秒"
print
gsave "~table"&str$(lv)&".bmp"
let lu = lv
let cont = 0
do while cont = 0
do while left = 1 or right = 1
mouse poll dxu,dyu,left,right
loop
call select_area
loop
if lv <> lu and uupdate(lu) = 0 then
let lv = lu
let wx = graphics(lu)
let gx1 = realpartmin(lu)
let gy1 = imaginarypartmin(lu)
let gx2 = realpartmax(lu)
let gy2 = imaginarypartmax(lu)
let num = repeatlimit(lu)
let krg = abs(gx1-gx2)/wx
end if
!---------------------------------------------------描画範囲変更に伴う各値の更新
if update = 1 then
let gx1_= gx1
let gy1_= gy1
let gx1 = gx1_+min(x1,x1+trangex)*krg
let gy1 = gy1_+min(y1,y1+trangey)*krg
let gx2 = gx1_+max(x1,x1+trangex)*krg
let gy2 = gy1_+max(y1,y1+trangey)*krg
let lpx1 = (gx1+gx2)/2
let lpy1 = (gy1+gy2)/2
let xrange = abs(trangex)*krg
let update = 0
end if
!---------------------------------------------------
call printdata
let uupdate(lu) = 0
loop
!━━━━━━━━━━━━━━━━━━━━━━━━━━━長方形を作る━━━━━━━━━━━━━━━━━━━━━━━━━━━
sub draw_area(x1,y1,x2,y2,cind)
set color cind
plot lines : x1,y1 ; x1,y2 ; x2,y2 ; x2,y1 ; x1,y1
end sub
!━━━━━━━━━━━━━━━━━━━━━━━━━━━長方形塗りつぶし━━━━━━━━━━━━━━━━━━━━━━━━━
sub flood_area(x1,y1,x2,y2,cind)
set area color cind
plot area : x1,y1 ; x1,y2 ; x2,y2 ; x2,y1 ; x1,y1
end sub
!━━━━━━━━━━━━━━━━━━━━━━━━━━━複素三角関数━━━━━━━━━━━━━━━━━━━━━━━━━━━
function fsin(z_)
let fsin = (exp(i*z_)-exp(-i*z_))/(2*i)
end function
function fcos(z_)
let fcos = (exp(i*z_)+exp(-i*z_))/2
end function
!━━━━━━━━━━━━━━━━━━━━━━━━━━━数列z(n)の反復式━━━━━━━━━━━━━━━━━━━━━━━━━
sub shift_sfn
let z2 = z*z
LET z = z-(z3*z2+z*z2*c+c)/(5*z2*z2+3*z2*c)
end sub
!━━━━━━━━━━━━━━━━━━━━━━━━━━━数列z(n)の初期値と複素数の変換式━━━━━━━━━━━━━━━━━
sub fpr_sfn
! let c = c
let z = -c/2 !数列z(n)の初期値
end sub
end
|
|
|
戻る
