|
投稿者：yoshipyuta
投稿日：2018年 1月 3日(水)23時03分46秒
|
|
|
|
しばっちさん、立て続けののKidsの質問をお許し下さい。z^8-17*Z^4+16のニュートン・フラクタルの図示が面白くありません。この式は８つの解＋－１、＋－i、＋－２、＋－２iを有します。
下記の図でおおよその８rootsの位置とは思いますが、単一色での表示が残念です。カラー化できますか。また(z^4-1)*(Z^2-(1+i)^2)のように虚数iが式に入る場合はどうなりますか。
|
|
|
|
投稿者：しばっち
投稿日：2018年 1月 6日(土)10時35分22秒
|
|
|
|
> No.4444[元記事へ]
yoshipyutaさんへのお返事です。
> しばっちさん、立て続けののKidsの質問をお許し下さい。z^8-17*Z^4+16のニュートン・フラクタルの図示が面白くありません。この式は８つの解＋－１、＋－i、＋－２、＋－２iを有します。
> 下記の図でおおよその８rootsの位置とは思いますが、単一色での表示が残念です。カラー化できますか。
下記のようにしてみてください。
注意点として、計算精度(計算回数)を調整してください。特に拡大表示で再計算させた場合は計算回数不足により
収束しきれずに打ち切ってしまっている場合があるようです。(黒色は背景色のままで何も描画されていません)
また、使用している数値微分式や、Hの値によりニュートン法の収束精度が変わります。
計算区間を表示させていますので、計算精度(KS)と領域(LEFT,RIGHT,BOTTOM,TOP)及び、収束判定精度(EPS)を書き換えて
再計算してみてください。
また、新たな試みとしてNEWTON法以外でも試してみました。収束の仕方が変わり図柄が(若干?)変わります。
(※どうしても実行が遅い時はBasicAcc or BasicAcc2(旧 ParactBasic)をご使用ください)
OPTION ARITHMETIC COMPLEX
LET XSIZE=800 !'画像サイズ
LET YSIZE=800
CALL GINIT(XSIZE,YSIZE)
LET LEFT=-1.5 !'計算領域
LET RIGHT=1.5
LET BOTTOM=1.5
LET TOP=-1.5
LET KS=100 !'計算精度(計算回数) ※要調整
LET EPS=1E-5 !'収束判定精度 ※要調整
DIM Z(8)
LET Z(1)=COMPLEX(1,0) !'Z^8-17*Z^4+16の解
LET Z(2)=COMPLEX(-1,0)
LET Z(3)=COMPLEX(0,1)
LET Z(4)=COMPLEX(0,-1)
LET Z(5)=COMPLEX(2,0)
LET Z(6)=COMPLEX(-2,0)
LET Z(7)=COMPLEX(0,2)
LET Z(8)=COMPLEX(0,-2)
DO
CLEAR
SET WINDOW LEFT,RIGHT,BOTTOM,TOP
LET DX=(RIGHT-LEFT)/XSIZE
LET DY=(TOP-BOTTOM)/YSIZE
FOR ZR=LEFT TO RIGHT STEP DX
FOR ZI=BOTTOM TO TOP STEP DY
LET XX=COMPLEX(ZR,ZI)
FOR K=1 TO KS
LET X=XX
WHEN EXCEPTION IN
LET XX=X-FUNC(X)/DIFF(X,1) !'NEWTON法(2次収束式)
!'LET XX=X-2*FUNC(X)*DIFF(X,1)/(2*DIFF(X,1)^2-FUNC(X)*DIFF(X,2)) !'HALLEY法(3次収束式)
!'LET XX=X-FUNC(X)/(DIFF(X,1)-(FUNC(X)*(3*DIFF(X,2)*DIFF(X,1)-DIFF(X,3)*FUNC(X)))/(3*(2*DIFF(X,1)^2-DIFF(X,2)*FUNC(X)))) !'KISS法(4次収束式)
!'LET XX=X+(4*FUNC(X)^3*DIFF(X,3)-24*FUNC(X)^2*DIFF(X,1)*DIFF(X,2)+24*FUNC(X)*DIFF(X,1)^3)/(FUNC(X)^3*DIFF(X,4)-8*FUNC(X)^2*DIFF(X,1)*DIFF(X,3)-6*FUNC(X)^2*DIFF(X,2)^2+36*FUNC(X)*DIFF(X,1)^2*DIFF(X,2)-24*DIFF(X,1)^4) !'5次収束式
USE
CALL PSET(ZR,ZI,255) !'発散!?
EXIT FOR
END WHEN
LET FL=0
FOR I=1 TO 8
IF ABS(XX-Z(I))<EPS THEN !'収束判定
CALL PSET(ZR,ZI,I)
LET FL=1 !'フラグ
EXIT FOR
END IF
NEXT I
IF FL=1 THEN EXIT FOR
NEXT K
NEXT ZI
NEXT ZR
PAUSE "拡大する範囲を指定してください"
CALL GETSQUARE(LEFT,TOP,RIGHT,BOTTOM)
PRINT "LET LEFT=";LEFT !'計算区間を表示 (コピペして再計算)
PRINT "LET RIGHT=";RIGHT
PRINT "LET BOTTOM=";BOTTOM
PRINT "LET TOP=";TOP
PRINT "LET KS=";KS
PRINT "LET EPS=";EPS
PRINT
IF LEFT=RIGHT THEN EXIT DO
LOOP
END
EXTERNAL FUNCTION FUNC(Z)
OPTION ARITHMETIC COMPLEX
LET FUNC=Z^8-17*Z^4+16
END FUNCTION
EXTERNAL FUNCTION DIFF(X,N) !'再帰呼び出しによる高階数値微分
OPTION ARITHMETIC COMPLEX
LET H=1/256 !' H=1/2^n の形
!'LET H=1/1024
IF N=0 THEN
LET DIFF=FUNC(X)
EXIT FUNCTION
ELSE
LET DIFF=(DIFF(X+H,N-1)-DIFF(X-H,N-1))/(2*H) !'3点微分
!'LET DIFF=(-DIFF(X+2*H,N-1)+8*DIFF(X+H,N-1)-8*DIFF(X-H,N-1)+DIFF(X-2*H,N-1))/(12*H) !'5点微分
!'LET DIFF=(DIFF(X+3*H,N-1)-9*DIFF(X+2*H,N-1)+45*DIFF(X+H,N-1)-45*DIFF(X-H,N-1)+9*DIFF(X-2*H,N-1)-DIFF(X-3*H,N-1))/(60*H) !'7点微分
!'LET DIFF=(-3*DIFF(X+4*H,N-1)+32*DIFF(X+3*H,N-1)-168*DIFF(X+2*H,N-1)+672*DIFF(X+H,N-1)-672*DIFF(X-H,N-1)+168*DIFF(X-2*H,N-1)-32*DIFF(X-3*H,N-1)+3*DIFF(X-4*H,N-1))/(840*H) !'9点微分
END IF
END FUNCTION
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 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
SET COLOR MIX(255) .5,.5,.5
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
> (z^4-1)*(Z^2-(1+i)^2)のように虚数iが式に入る場合はどうなりますか。
下記のようにCOMPLEX文 又はSQR(-1)を使ってください。
(Z^4-1)*(Z^2-COMPLEX(1,1)^2)
> Z＊(Z^3-1)の２枚葉の発芽の図を下記に添付します。x: 0.5～0.7 y: -0.12～0.12領域
計算区間にずれがあります。
X:0.7-0.5=0.2
Y:0.12-(-0.12)=0.24
この時の表示画面が正方形なのに対し、X,Yの区間が違うため歪みが生じています。
区間幅を同じにするか、表示ウィンドウの上下左右の比率と合わせてください。
但し、SUB GETSQUAREで求める領域は正方形ウィンドウです(上下左右の比率が1)
Y:-0.1～0.1 0.1-(-0.1)=0.2
とするか、
0.24/0.2=1.2
LET XSIZE=800 !'画像サイズ
LET YSIZE=XSIZE*1.2
CALL GINIT(XSIZE,YSIZE)
又は
0.2/0.24=5/6
LET YSIZE=800
LEY XSIZE=INT(YSIZE*5/6)
CALL GINIT(XSIZE,YSIZE)
>さて、さらに進んでTan(z)やCos(z)に展開できるようにできますか？
複素数モードでもEXP,LOG,SQR関数などが使用できますが、SIN,COS,TANなどは使えません。
下記のように定義すれば複素数関数が使用できるようになります。
OPTION ARITHMETIC COMPLEX
END
EXTERNAL FUNCTION CSIN(Z)
OPTION ARITHMETIC COMPLEX
LET I=SQR(-1)
LET CSIN=(EXP(Z*I)-EXP(-Z*I))/2/I
END FUNCTION
EXTERNAL FUNCTION CCOS(Z)
OPTION ARITHMETIC COMPLEX
LET I=SQR(-1)
LET CCOS=(EXP(Z*I)+EXP(-Z*I))/2
END FUNCTION
EXTERNAL FUNCTION CTAN(Z)
OPTION ARITHMETIC COMPLEX
LET CTAN=CSIN(Z)/CCOS(Z)
END FUNCTION
EXTERNAL FUNCTION CACOS(Z)
OPTION ARITHMETIC COMPLEX
LET I=SQR(-1)
LET CACOS=-I*LOG(Z+I*SQR(1-Z*Z))
!' LET CACOS=I*LOG(Z-I*SQR(1-Z*Z))
END FUNCTION
EXTERNAL FUNCTION CASIN(Z)
OPTION ARITHMETIC COMPLEX
LET I=SQR(-1)
LET CASIN=-I*LOG(I*Z+SQR(1-Z*Z))
END FUNCTION
EXTERNAL FUNCTION CATAN(Z)
OPTION ARITHMETIC COMPLEX
LET I=SQR(-1)
LET CATAN=I/2*LOG((I+Z)/(I-Z))
!' LET CATAN=I/2*(LOG(1-I*Z)-LOG(1+I*Z))
END FUNCTION
EXTERNAL FUNCTION CACOSEC(Z)
OPTION ARITHMETIC COMPLEX
LET I=SQR(-1)
LET CACOSEC=-I*LOG(I/Z+SQR(1-1/Z/Z))
END FUNCTION
EXTERNAL FUNCTION CASEC(Z)
OPTION ARITHMETIC COMPLEX
LET I=SQR(-1)
LET CASEC=-LOG((1+SQR(1-Z*Z))/Z)/I
!'LET CASEC=-I*LOG(I*SQR(1-1/Z/Z)+1/Z)
END FUNCTION
EXTERNAL FUNCTION CACOTAN(Z)
OPTION ARITHMETIC COMPLEX
LET I=SQR(-1)
LET CACOTAN=I/2*LOG((Z-I)/(Z+I))
!'LET CACOTAN=I/2*(LLOG(1-I/Z)-LLOG(1+I/Z))
END FUNCTION
EXTERNAL FUNCTION CSEC(Z)
OPTION ARITHMETIC COMPLEX
LET CSEC=1/CCOS(Z)
END FUNCTION
EXTERNAL FUNCTION CCOSEC(Z)
OPTION ARITHMETIC COMPLEX
LET CCOSEC=1/CSIN(Z)
END FUNCTION
EXTERNAL FUNCTION CSINH(Z)
OPTION ARITHMETIC COMPLEX
LET I=SQR(-1)
LET CSINH=-I*CSIN(Z*I)
!'LET CSINH=(EXP(Z)-EXP(-Z))/2
END FUNCTION
EXTERNAL FUNCTION CCOSH(Z)
OPTION ARITHMETIC COMPLEX
LET I=SQR(-1)
LET CCOSH=CCOS(Z*I)
!'LET CCOSH=(EXP(Z)+EXP(-Z))/2
END FUNCTION
EXTERNAL FUNCTION CTANH(Z)
OPTION ARITHMETIC COMPLEX
LET I=SQR(-1)
LET CTANH=-I*CTAN(Z*I)
END FUNCTION
EXTERNAL FUNCTION CCOTAN(Z)
OPTION ARITHMETIC COMPLEX
LET CCOTAN=1/CTAN(Z)
END FUNCTION
EXTERNAL FUNCTION CSECH(Z)
OPTION ARITHMETIC COMPLEX
LET CSECH=1/CCOSH(Z)
END FUNCTION
EXTERNAL FUNCTION CCOSECH(Z)
OPTION ARITHMETIC COMPLEX
LET CCOSECH=1/CSINH(Z)
END FUNCTION
EXTERNAL FUNCTION CCOTANH(Z)
OPTION ARITHMETIC COMPLEX
LET CCOTANH=1/CTANH(Z)
END FUNCTION
EXTERNAL FUNCTION CASINH(Z)
OPTION ARITHMETIC COMPLEX
LET CASINH=LOG(Z+SQR(Z*Z+1))
END FUNCTION
EXTERNAL FUNCTION CACOSH(Z)
OPTION ARITHMETIC COMPLEX
LET I=SQR(-1)
LET CACOSH=LOG(Z+SQR(Z+1)*SQR(Z-1))
END FUNCTION
EXTERNAL FUNCTION CATANH(Z)
OPTION ARITHMETIC COMPLEX
LET CATANH=LOG((1+Z)/(1-Z))/2
END FUNCTION
EXTERNAL FUNCTION CACOTANH(Z)
OPTION ARITHMETIC COMPLEX
LET CACOTANH=LOG((Z+1)/(Z-1))/2
END FUNCTION
EXTERNAL FUNCTION CACOSECH(Z)
OPTION ARITHMETIC COMPLEX
LET CACOSECH=LOG(1/Z+SQR(1/Z/Z+1))
END FUNCTION
EXTERNAL FUNCTION CASECH(Z)
OPTION ARITHMETIC COMPLEX
LET CASECH=LOG(1/Z+SQR(1/Z+1)*SQR(1/Z-1))
END FUNCTION
|
|
|
|
投稿者：しばっち
投稿日：2018年 1月 6日(土)10時38分54秒
|
|
|
|
更に先に挙げた以外の方法で色々と試してみました。
主要な部分だけを抜き出していますので、付け足してから
実行してください。
(※関数式によっては全く収束しない場合があります)
OPTION ARITHMETIC COMPLEX
LET XSIZE=800
LET YSIZE=800
CALL GINIT(XSIZE,YSIZE)
LET LEFT=-1.5
LET RIGHT=1.5
LET BOTTOM=1.5
LET TOP=-1.5
LET EPS=1E-5
LET KS=100
CLEAR
SET WINDOW LEFT,RIGHT,BOTTOM,TOP
LET DX=(RIGHT-LEFT)/XSIZE
LET DY=(TOP-BOTTOM)/YSIZE
FOR ZR=LEFT TO RIGHT STEP DX
FOR ZI=BOTTOM TO TOP STEP DY
LET XX=COMPLEX(ZR,ZI)
LET X0=1
FOR K=1 TO KS
LET X=XX
WHEN EXCEPTION IN
LET XX=X-((X0-X)/(FUNC(X0)-FUNC(X)))*FUNC(X) !'SCANT法
LET X0=X
USE
CALL PSET(ZR,ZI,255)
EXIT FOR
END WHEN
IF ABS(X-XX)<EPS THEN
LET C=MOD(K,7)
CALL PSET(ZR,ZI,C+1)
EXIT FOR
END IF
NEXT K
NEXT ZI
NEXT ZR
END
EXTERNAL FUNCTION FUNC(X)
OPTION ARITHMETIC COMPLEX
LET FUNC=X^8-X^6+2*X^5+X^4+3*X^2-X+1
END FUNCTION
------------------------------------------------------------------------------
FOR ZR=LEFT TO RIGHT STEP DX
FOR ZI=BOTTOM TO TOP STEP DY
LET XX=COMPLEX(ZR,ZI)
FOR K=1 TO KS
LET X=XX
WHEN EXCEPTION IN
LET XX=X-FUNC(X)^2/(FUNC(FUNC(X)+X)-FUNC(X)) !'STEFFENSEN法
USE
CALL PSET(ZR,ZI,255)
EXIT FOR
END WHEN
IF ABS(X-XX)<EPS THEN
LET C=MOD(K,7)
CALL PSET(ZR,ZI,C+1)
EXIT FOR
END IF
NEXT K
NEXT ZI
NEXT ZR
END
EXTERNAL FUNCTION FUNC(X)
OPTION ARITHMETIC COMPLEX
LET FUNC=X^3-1
!'LET FUNC=X^3-3^X
!'LET FUNC=EXP(X*LOG(X))-1
END FUNCTION
------------------------------------------------------------------------------
FOR ZR=LEFT TO RIGHT STEP DX
FOR ZI=BOTTOM TO TOP STEP DY
LET XX=COMPLEX(ZR,ZI)
LET X0=XX
LET X2=0
LET X1=(X0+X2)/2
FOR K=1 TO KS
WHEN EXCEPTION IN
LET TT=(X2-X1)/(X2-X0)*(FUNC(X2)-FUNC(X0))/(FUNC(X2)-FUNC(X1))*FUNC(X1)/FUNC(X0) !'OSTROWSKI法
LET X3=(X1-X0*TT)/(1-TT)
LET X0=X1
LET X1=X2
LET X2=X3
USE
CALL PSET(ZR,ZI,255)
EXIT FOR
END WHEN
IF ABS(FUNC(X3))<EPS THEN
LET C=MOD(K,7)
CALL PSET(ZR,ZI,C+1)
EXIT FOR
END IF
NEXT K
NEXT ZI
NEXT ZR
END
EXTERNAL FUNCTION FUNC(X)
OPTION ARITHMETIC COMPLEX
LET FUNC=X^3-1
END FUNCTION
------------------------------------------------------------------------------
FOR ZR=LEFT TO RIGHT STEP DX
FOR ZI=BOTTOM TO TOP STEP DY
LET XX=COMPLEX(ZR,ZI)
FOR K=1 TO KS
LET X=XX
WHEN EXCEPTION IN
LET XX=X-FUNC(X)*DIFF(X,1)/(DIFF(X,1)^2-DIFF(X,2)*FUNC(X))
USE
CALL PSET(ZR,ZI,255)
EXIT FOR
END WHEN
IF ABS(X-XX)<EPS THEN
LET C=MOD(K,7)
CALL PSET(ZR,ZI,C+1)
EXIT FOR
END IF
NEXT K
NEXT ZI
NEXT ZR
END
EXTERNAL FUNCTION FUNC(X)
OPTION ARITHMETIC COMPLEX
LET FUNC=X^3-1
END FUNCTION
EXTERNAL FUNCTION DIFF(X,N) !'再帰呼び出しによる高階数値微分
OPTION ARITHMETIC COMPLEX
LET H=1/256 !' H=1/2^n の形
!'LET H=1/1024
IF N=0 THEN
LET DIFF=FUNC(X)
EXIT FUNCTION
ELSE
LET DIFF=(DIFF(X+H,N-1)-DIFF(X-H,N-1))/(2*H) !'3点微分
!'LET DIFF=(-DIFF(X+2*H,N-1)+8*DIFF(X+H,N-1)-8*DIFF(X-H,N-1)+DIFF(X-2*H,N-1))/(12*H) !'5点微分
!'LET DIFF=(DIFF(X+3*H,N-1)-9*DIFF(X+2*H,N-1)+45*DIFF(X+H,N-1)-45*DIFF(X-H,N-1)+9*DIFF(X-2*H,N-1)-DIFF(X-3*H,N-1))/(60*H) !'7点微分
!'LET DIFF=(-3*DIFF(X+4*H,N-1)+32*DIFF(X+3*H,N-1)-168*DIFF(X+2*H,N-1)+672*DIFF(X+H,N-1)-672*DIFF(X-H,N-1)+168*DIFF(X-2*H,N-1)-32*DIFF(X-3*H,N-1)+3*DIFF(X-4*H,N-1))/(840*H) !'9点微分
END IF
END FUNCTION
------------------------------------------------------------------------------
FOR ZR=LEFT TO RIGHT STEP DX
FOR ZI=BOTTOM TO TOP STEP DY
LET XX=COMPLEX(ZR,ZI)
FOR K=1 TO KS
WHEN EXCEPTION IN
LET X=XX
LET H=1
DO
LET H=H/2
LOOP UNTIL ABS(H*FUNC(X))<3
LET XX=X-H*FUNC(X) !'原始反復法
USE
CALL PSET(ZR,ZI,255)
EXIT FOR
END WHEN
IF ABS(X-XX)<EPS THEN
LET C=MOD(K,7)
CALL PSET(ZR,ZI,C+1)
EXIT FOR
END IF
NEXT K
NEXT ZI
NEXT ZR
END
EXTERNAL FUNCTION FUNC(X)
OPTION ARITHMETIC COMPLEX
LET FUNC=X^3-1
END FUNCTION
------------------------------------------------------------------------------
FOR ZR=LEFT TO RIGHT STEP DX
FOR ZI=BOTTOM TO TOP STEP DY
LET XX=COMPLEX(ZR,ZI)
FOR K=1 TO KS
LET X=XX
WHEN EXCEPTION IN
LET L=DIFF(X,1)*(FUNC(X)-FUNC(X-FUNC(X)/DIFF(X,1)))
IF L=0 THEN
LET XX=X
ELSE
LET XX=X-FUNC(X)^2/L !'改良ニュートン法
END IF
USE
CALL PSET(ZR,ZI,255)
EXIT FOR
END WHEN
IF ABS(X-XX)<EPS THEN
LET C=MOD(K,7)
CALL PSET(ZR,ZI,C+1)
EXIT FOR
END IF
NEXT K
NEXT ZI
NEXT ZR
END
EXTERNAL FUNCTION FUNC(X)
OPTION ARITHMETIC COMPLEX
LET FUNC=X^3-1
END FUNCTION
------------------------------------------------------------------------------
FOR ZR=LEFT TO RIGHT STEP DX
FOR ZI=BOTTOM TO TOP STEP DY
LET XX=COMPLEX(ZR,ZI)
FOR K=1 TO KS
LET X=XX
WHEN EXCEPTION IN
LET U=2
LET COUNT=0
DO
LET U=U/2
LET COUNT=COUNT+1
WHEN EXCEPTION IN
LET XX=X-U*FUNC(X)/DIFF(X,1) !'減速ニュートン法
USE
EXIT FOR
END WHEN
LOOP UNTIL ABS(FUNC(XX))<(1-U/2)*ABS(FUNC(X)) OR COUNT>=100
USE
CALL PSET(ZR,ZI,255)
EXIT FOR
END WHEN
IF ABS(X-XX)<EPS THEN
LET C=MOD(K,7)
CALL PSET(ZR,ZI,C+1)
EXIT FOR
END IF
NEXT K
NEXT ZI
NEXT ZR
END
EXTERNAL FUNCTION FUNC(X)
OPTION ARITHMETIC COMPLEX
LET FUNC=X^3-1
END FUNCTION
------------------------------------------------------------------------------
FOR ZR=LEFT TO RIGHT STEP DX
FOR ZI=BOTTOM TO TOP STEP DY
LET XX=COMPLEX(ZR,ZI)
FOR K=1 TO KS
LET X=XX
WHEN EXCEPTION IN
LET XX=X-FUNC(X)/(DIFF(X,1)-DIFF(X,2)/2*FUNC(X)/DIFF(X,1)-DIFF(X,3)/6*(FUNC(X)/DIFF(X,1))^2-DIFF(X,4)/24*(FUNC(X)/DIFF(X,1))^3) !'拡張ニュートン法?
USE
CALL PSET(ZR,ZI,255)
EXIT FOR
END WHEN
IF ABS(X-XX)<EPS THEN
LET C=MOD(K,7)
CALL PSET(ZR,ZI,C+1)
EXIT FOR
END IF
NEXT K
NEXT ZI
NEXT ZR
END
EXTERNAL FUNCTION FUNC(X)
OPTION ARITHMETIC COMPLEX
LET FUNC=X^10-1
END FUNCTION
EXTERNAL FUNCTION DIFF(X,N)
OPTION ARITHMETIC COMPLEX
SELECT CASE N
CASE 1
LET DIFF=10*X^9
CASE 2
LET DIFF=10*9*X^8
CASE 3
LET DIFF=10*9*8*X^7
CASE 4
LET DIFF=10*9*8*7*X^6
END SELECT
END FUNCTION
------------------------------------------------------------------------------
FOR ZR=LEFT TO RIGHT STEP DX
FOR ZI=BOTTOM TO TOP STEP DY
LET XX=COMPLEX(ZR,ZI)
LET X0=XX
FOR K=1 TO KS
LET X=XX
WHEN EXCEPTION IN
LET XX=X-FUNC(X)/DIFF(X0,1) !'修正ニュートン法?
IF ABS(X-XX)<EPS THEN
LET C=MOD(K,7)
CALL PSET(ZR,ZI,C+1)
EXIT FOR
END IF
USE
CALL PSET(ZR,ZI,255)
EXIT FOR
END WHEN
NEXT K
NEXT ZI
NEXT ZR
END
EXTERNAL FUNCTION FUNC(X)
OPTION ARITHMETIC COMPLEX
LET FUNC=X^3-1
END FUNCTION
------------------------------------------------------------------------------
FOR ZR=LEFT TO RIGHT STEP DX
FOR ZI=BOTTOM TO TOP STEP DY
LET X0=COMPLEX(ZR,ZI)
LET X1=FUNC(X0)
FOR K=1 TO KS
WHEN EXCEPTION IN
LET X2=X1+(X1-X0)/((X0-FUNC(X0))/(X1-FUNC(X1))-1) !'WEGSTEIN法
LET X0=X1
LET X1=X2
USE
CALL PSET(ZR,ZI,255)
EXIT FOR
END WHEN
IF ABS(FUNC(X2)-X2)<EPS THEN
LET C=MOD(K,7)
CALL PSET(ZR,ZI,C+1)
EXIT FOR
END IF
NEXT K
NEXT ZI
NEXT ZR
END
EXTERNAL FUNCTION FUNC(X)
OPTION ARITHMETIC COMPLEX
WHEN EXCEPTION IN
LET FUNC=X^3+X-1 !' f(x)+x
USE
END WHEN
END FUNCTION
------------------------------------------------------------------------------
FOR ZR=LEFT TO RIGHT STEP DX
FOR ZI=BOTTOM TO TOP STEP DY
LET X0=COMPLEX(ZR,ZI)
FOR K=1 TO KS
WHEN EXCEPTION IN
LET X1=FUNC(X0)
LET X2=FUNC(X1)
LET X3=X0-(X1-X0)*(X1-X0)/(X2-2*X1+X0) !' エイトケン法
LET X0=X3
USE
CALL PSET(ZR,ZI,255)
EXIT FOR
END WHEN
IF ABS(FUNC(X0)-X0)<EPS THEN
LET C=MOD(K,7)
CALL PSET(ZR,ZI,C+1)
EXIT FOR
END IF
NEXT K
NEXT ZI
NEXT ZR
END
EXTERNAL FUNCTION FUNC(X)
OPTION ARITHMETIC COMPLEX
LET FUNC=X^3+X-1 !' f(x)+x
END FUNCTION
------------------------------------------------------------------------------
FOR ZR=LEFT TO RIGHT STEP DX
FOR ZI=BOTTOM TO TOP STEP DY
LET XX=COMPLEX(ZR,ZI)
FOR K=1 TO KS
LET X=XX
WHEN EXCEPTION IN
LET XX=(FUNC(X)+X)/2 !'相加平均
USE
CALL PSET(ZR,ZI,255)
EXIT FOR
END WHEN
IF ABS(X-XX)<EPS THEN
LET C=MOD(K,7)
CALL PSET(ZR,ZI,C+1)
EXIT FOR
END IF
NEXT K
NEXT ZI
NEXT ZR
END
EXTERNAL FUNCTION FUNC(X)
OPTION ARITHMETIC COMPLEX
LET FUNC=1/X^2 !' x^3-1=0 → x^3=1 → x=1/x^2
END FUNCTION
------------------------------------------------------------------------------
FOR ZR=LEFT TO RIGHT STEP DX
FOR ZI=BOTTOM TO TOP STEP DY
LET XX=COMPLEX(ZR,ZI)
FOR K=1 TO KS
LET X=XX
WHEN EXCEPTION IN
LET XX=X-H(X)/DH(X) !'ハリー法
USE
CALL PSET(ZR,ZI,255)
EXIT FOR
END WHEN
IF ABS(X-XX)<EPS THEN
LET C=MOD(K,7)
CALL PSET(ZR,ZI,C+1)
EXIT FOR
END IF
NEXT K
NEXT ZI
NEXT ZR
END
EXTERNAL FUNCTION FUNC(X)
OPTION ARITHMETIC COMPLEX
LET FUNC=X^3-1
END FUNCTION
EXTERNAL FUNCTION H(X)
OPTION ARITHMETIC COMPLEX
LET H=FUNC(X)/SQR(DIFF(X,1))
END FUNCTION
EXTERNAL FUNCTION DH(X)
OPTION ARITHMETIC COMPLEX
LET DH=SQR(DIFF(X,1))-FUNC(X)*DIFF(X,2)/(2*DIFF(X,1)*SQR(DIFF(X,1)))
END FUNCTION
------------------------------------------------------------------------------
FOR ZR=LEFT TO RIGHT STEP DX
FOR ZI=BOTTOM TO TOP STEP DY
LET XX=COMPLEX(ZR,ZI)
FOR K=1 TO KS
LET X=XX
WHEN EXCEPTION IN
LET XX=FUNC(X) !'逐次代入法
USE
CALL PSET(ZR,ZI,255)
EXIT FOR
END WHEN
IF ABS((X+1)*(X+1)*(X+1)-3)<EPS THEN
LET C=MOD(K,7)
CALL PSET(ZR,ZI,C+1)
EXIT FOR
END IF
NEXT K
NEXT ZI
NEXT ZR
END
EXTERNAL FUNCTION FUNC(X)
OPTION ARITHMETIC COMPLEX
LET FUNC=2/(X*X+3*X+3) !'X=3^(1/3)-1 → (X+1)^3=3 → X*(X^2+3*X+3)=2 → X=2/(X^2+3*X+3)
END FUNCTION
------------------------------------------------------------------------------
FOR ZR=LEFT TO RIGHT STEP DX
FOR ZI=BOTTOM TO TOP STEP DY
LET XX=COMPLEX(ZR,ZI)
FOR K=1 TO KS
LET X=XX
WHEN EXCEPTION IN
LET XX=X-FUNC(X)/DIFF(X,1)*(1+FUNC(X)*DIFF(X,2)/2/DIFF(X,1)^2) !'HOUSEHOLDER法
USE
CALL PSET(ZR,ZI,255)
EXIT FOR
END WHEN
IF ABS(X-XX)<EPS THEN
LET C=MOD(K,7)
CALL PSET(ZR,ZI,C+1)
EXIT FOR
END IF
NEXT K
NEXT ZI
NEXT ZR
END
EXTERNAL FUNCTION FUNC(X)
OPTION ARITHMETIC COMPLEX
LET FUNC=X^3-1
END FUNCTION
EXTERNAL FUNCTION DIFF(X,N) !'高階数値微分
OPTION ARITHMETIC COMPLEX
LET H=1/128
IF MOD(N,2)=1 THEN LET SIGN=-1 ELSE LET SIGN=1
FOR I=0 TO N
LET S=S+(-1)^I*COMB(N,I)*FUNC(X-(N-2*I)*H)*SIGN
NEXT I
LET DIFF=S/((2*H)^N)
END FUNCTION
|
|
|
|
投稿者：しばっち
投稿日：2018年 1月 6日(土)10時39分28秒
|
|
|
|
> No.4452[元記事へ]
続き
FOR ZR=LEFT TO RIGHT STEP DX
FOR ZI=BOTTOM TO TOP STEP DY
LET XX=COMPLEX(ZR,ZI)
FOR K=1 TO KS
LET X=XX
WHEN EXCEPTION IN
LET XX=FUNC(X) !'定点法
IF ABS(XX-X)<EPS THEN
LET C=MOD(K,7)
CALL PSET(ZR,ZI,C+1)
EXIT FOR
END IF
USE
CALL PSET(ZR,ZI,255)
EXIT FOR
END WHEN
NEXT K
NEXT ZI
NEXT ZR
END
EXTERNAL FUNCTION FUNC(X) !' X^8-X^6+2*X^5+X^4+3*X^2-X+1=0 → 3*X^2=-X^8+X^6-2*X^5-X^4+X-1
OPTION ARITHMETIC COMPLEX
LET FUNC=SQR((-X^8+X^6-2*X^5-X^4+X-1)/3)
END FUNCTION
------------------------------------------------------------------------------
FOR ZR=LEFT TO RIGHT STEP DX
FOR ZI=BOTTOM TO TOP STEP DY
LET X1=COMPLEX(ZR,ZI)
LET X2=(1+X1)/2
LET X3=1
FOR K=1 TO KS
WHEN EXCEPTION IN
LET Q=(X3-X2)/(X2-X1)
LET A=Q*FUNC(X3)-Q*(1+Q)*FUNC(X2)+Q^2*FUNC(X1)
LET B=(2*Q+1)*FUNC(X3)-(1+Q)^2*FUNC(X2)+Q^2*FUNC(X1)
LET C=(1+Q)*FUNC(X3)
LET DD=B^2-4*A*C
IF ABS(B+SQR(DD))>ABS(B-SQR(DD)) THEN
LET X4=X3-(X3-X2)*2*C/(B+SQR(DD)) !'Muller法
ELSE
LET X4=X3-(X3-X2)*2*C/(B-SQR(DD))
END IF
LET X1=X2
LET X2=X3
LET X3=X4
USE
CALL PSET(ZR,ZI,255)
EXIT FOR
END WHEN
IF ABS(FUNC(X4))<EPS THEN
LET C=MOD(K,7)
CALL PSET(ZR,ZI,C+1)
EXIT FOR
END IF
NEXT K
NEXT ZI
NEXT ZR
END
EXTERNAL FUNCTION FUNC(X)
OPTION ARITHMETIC COMPLEX
LET FUNC=X^8-X^6+2*X^5+X^4+3*X^2-X+1
END FUNCTION
------------------------------------------------------------------------------
OPTION ARITHMETIC COMPLEX
PUBLIC NUMERIC N
LET XSIZE=800
LET YSIZE=800
CALL GINIT(XSIZE,YSIZE)
LET LEFT=-1.5
LET RIGHT=1.5
LET BOTTOM=1.5
LET TOP=-1.5
LET EPS=1E-5
LET KS=100
LET N=3
CLEAR
SET WINDOW LEFT,RIGHT,BOTTOM,TOP
LET DX=(RIGHT-LEFT)/XSIZE
LET DY=(TOP-BOTTOM)/YSIZE
FOR ZR=LEFT TO RIGHT STEP DX
FOR ZI=BOTTOM TO TOP STEP DY
LET XX=COMPLEX(ZR,ZI)
FOR K=1 TO KS
LET X=XX
WHEN EXCEPTION IN
LET L=(N-1)*(N*H(X)-G(X)^2)
!'IF ABS(G(X)+L)>ABS(G(X)-L) THEN
LET XX=X-N/(G(X)+L) !' Laguerre法
!'ELSE
!' LET XX=X-N/(G(X)-L)
!'END IF
USE
CALL PSET(ZR,ZI,255)
EXIT FOR
END WHEN
IF ABS(X-XX)<EPS THEN
LET C=MOD(K,7)
CALL PSET(ZR,ZI,C+1)
EXIT FOR
END IF
NEXT K
NEXT ZI
NEXT ZR
END
EXTERNAL FUNCTION FUNC(X)
OPTION ARITHMETIC COMPLEX
LET FUNC=X^N-1
END FUNCTION
EXTERNAL FUNCTION DIFF(X)
OPTION ARITHMETIC COMPLEX
LET DIFF=N*X^(N-1)
END FUNCTION
EXTERNAL FUNCTION DIFF2(X)
OPTION ARITHMETIC COMPLEX
LET DIFF2=N*(N-1)*X^(N-2)
END FUNCTION
EXTERNAL FUNCTION G(X)
OPTION ARITHMETIC COMPLEX
LET G=FUNC(X)/DIFF(X)
END FUNCTION
EXTERNAL FUNCTION H(X)
OPTION ARITHMETIC COMPLEX
LET H=G(X)^2-DIFF2(X)/FUNC(X)
END FUNCTION
------------------------------------------------------------------------------
OPTION ARITHMETIC COMPLEX
PUBLIC NUMERIC N,R
LET XSIZE=800
LET YSIZE=800
CALL GINIT(XSIZE,YSIZE)
LET LEFT=-1.5
LET RIGHT=1.5
LET BOTTOM=1.5
LET TOP=-1.5
LET EPS=1E-5
LET KS=100
LET N=3
LET R=1
CLEAR
SET WINDOW LEFT,RIGHT,BOTTOM,TOP
LET DX=(RIGHT-LEFT)/XSIZE
LET DY=(TOP-BOTTOM)/YSIZE
FOR ZR=LEFT TO RIGHT STEP DX
FOR ZI=BOTTOM TO TOP STEP DY
LET XX=COMPLEX(ZR,ZI)
FOR K=1 TO KS
LET X=XX
WHEN EXCEPTION IN
LET XX=H(X) !'Lambert's法
USE
CALL PSET(ZR,ZI,255)
EXIT FOR
END WHEN
IF ABS(X-XX)<EPS THEN
LET C=MOD(K,7)
CALL PSET(ZR,ZI,C+1)
EXIT FOR
END IF
NEXT K
NEXT ZI
NEXT ZR
END
EXTERNAL FUNCTION FUNC(X)
OPTION ARITHMETIC COMPLEX
LET FUNC=X^N-R
END FUNCTION
EXTERNAL FUNCTION H(X)
OPTION ARITHMETIC COMPLEX
LET H=((N-1)*X^N+(N+1)*R)/((N+1)*X^N+(N-1)*R)*X
END FUNCTION
------------------------------------------------------------------------------
FOR ZR=LEFT TO RIGHT STEP DX
FOR ZI=BOTTOM TO TOP STEP DY
LET X1=COMPLEX(ZR,ZI)
LET X2=(1+X1)/2
LET X3=0
FOR K=1 TO KS
WHEN EXCEPTION IN
LET XX=X2+P(X1,X2,X3)/Q(X1,X2,X3) !'Brent's法
LET X3=X2
LET X2=X1
LET X1=XX
USE
CALL PSET(ZR,ZI,255)
EXIT FOR
END WHEN
IF ABS(FUNC(XX))<EPS THEN
LET C=MOD(K,7)
CALL PSET(ZR,ZI,C+1)
EXIT FOR
END IF
NEXT K
NEXT ZI
NEXT ZR
END
EXTERNAL FUNCTION FUNC(X)
OPTION ARITHMETIC COMPLEX
LET FUNC=X^3-1
END FUNCTION
EXTERNAL FUNCTION P(X1,X2,X3)
OPTION ARITHMETIC COMPLEX
LET P=S(X1,X2)*(T(X1,X3)*(R(X2,X3)-T(X1,X3))*(X3-X2)-(1-R(X2,X3))*(X2-X1))
END FUNCTION
EXTERNAL FUNCTION Q(X1,X2,X3)
OPTION ARITHMETIC COMPLEX
LET Q=(T(X1,X3)-1)*(R(X2,X3)-1)*(S(X1,X2)-1)
END FUNCTION
EXTERNAL FUNCTION R(X2,X3)
OPTION ARITHMETIC COMPLEX
LET R=FUNC(X2)/FUNC(X3)
END FUNCTION
EXTERNAL FUNCTION S(X1,X2)
OPTION ARITHMETIC COMPLEX
LET S=FUNC(X2)/FUNC(X1)
END FUNCTION
EXTERNAL FUNCTION T(X1,X3)
OPTION ARITHMETIC COMPLEX
LET T=FUNC(X1)/FUNC(X3)
END FUNCTION
------------------------------------------------------------------------------
FOR ZR=LEFT TO RIGHT STEP DX
FOR ZI=BOTTOM TO TOP STEP DY
LET X1=COMPLEX(ZR,ZI)
LET X2=1
FOR K=1 TO KS
WHEN EXCEPTION IN
LET XX=X1-(X2-X1)/(FUNC(X2)-FUNC(X1))*FUNC(X1) !'Method of False Position法
LET X2=XX
USE
CALL PSET(ZR,ZI,255)
EXIT FOR
END WHEN
IF ABS(FUNC(XX))<EPS THEN
LET C=MOD(K,7)
CALL PSET(ZR,ZI,C+1)
EXIT FOR
END IF
NEXT K
NEXT ZI
NEXT ZR
END
EXTERNAL FUNCTION FUNC(X)
OPTION ARITHMETIC COMPLEX
LET FUNC=X^8-X^6+2*X^5+X^4+3*X^2-X+1
END FUNCTION
|
|
|
|
投稿者：しばっち
投稿日：2018年 1月 6日(土)17時22分35秒
|
|
|
|
> No.4451[元記事へ]
更に2次元(複素数 x+yi (i^2=-1))ではなく、4次元(Quaternion x+yi+zj+wk (i^2=j^2=k^2=-1))でも求めてみました。
但し、3D表示は敷居が高く、表示方法や仕様設定等が未定なので表示は2Dです。
4次元では平面が6つで、X-Y平面 , X-Z平面 , X-W平面 , Y-Z平面 , Y-W平面 , Z-W平面があります。
残った2軸にはスライドバーにより、(-1.5 <= a <= 1.5)(-1.5 <= b <=1.5) を代入して定数とした切断面を描画します。
計算にクォータニオン(4元数)を用い、ニュートン法にて方程式 x^3-1=0 を解きながら結果を描画しています。
(※X-Y平面で Z=0 , W=0とした時の切断面は複素数による一般的なものと同じになります)
OPTION ARITHMETIC NATIVE
OPTION BASE 0
DIM XX(3),X(3),S$(5),ST(6)
LET XSIZE=800 !'画像サイズ
LET YSIZE=800
CALL GINIT(XSIZE,YSIZE)
MAT READ S$
LOCATE CHOICE(S$) : N
DATA "X-Y平面 Z=A,W=B"
DATA "X-Z平面 Y=A,W=B"
DATA "X-W平面 Y=A,Z=B"
DATA "Y-Z平面 X=A,W=B"
DATA "Y-W平面 X=A,Z=B"
DATA "Z-W平面 X=A,Y=B"
LOCATE VALUE(1) ,RANGE -1.5 TO 1.5 ,AT 0: A
LOCATE VALUE(2) ,RANGE -1.5 TO 1.5 ,AT 0: B
SELECT CASE N
CASE 1
PRINT "X-Y平面 Z=";A;", W=";B
CASE 2
PRINT "X-Z平面 Y=";A;", W=";B
CASE 3
PRINT "X-W平面 Y=";A;", Z=";B
CASE 4
PRINT "Y-Z平面 X=";A;", W=";B
CASE 5
PRINT "Y-W平面 X=";A;", Z=";B
CASE 6
PRINT "Z-W平面 X=";A;", Y=";B
END SELECT
LET LEFT=-1.5
LET RIGHT=1.5
LET BOTTOM=1.5
LET TOP=-1.5
LET KS=100
LET EPS=1E-5
SET WINDOW LEFT,RIGHT,BOTTOM,TOP
LET DX=(RIGHT-LEFT)/XSIZE
LET DY=(TOP-BOTTOM)/YSIZE
MAT READ ST
DATA 1,7,11,13,17,19,23
FOR I=6 TO 0 STEP -1
FOR ZR=LEFT TO RIGHT STEP DX*ST(I)
FOR ZI=BOTTOM TO TOP STEP DY*ST(I)
SELECT CASE N
CASE 1
CALL QSET(XX,ZR,ZI,A,B) !' X-Y平面 Z=A,W=B 切断面
CASE 2
CALL QSET(XX,ZR,A,ZI,B) !' X-Z平面 Y=A,W=B 切断面
CASE 3
CALL QSET(XX,ZR,A,B,ZI) !' X-W平面 Y=A,Z=B 切断面
CASE 4
CALL QSET(XX,A,ZR,ZI,B) !' Y-Z平面 X=A,W=B 切断面
CASE 5
CALL QSET(XX,A,ZR,B,ZI) !' Y-W平面 X=A,Z=B 切断面
CASE 6
CALL QSET(XX,A,B,ZR,ZI) !' Z-W平面 X=A,Y=B 切断面
END SELECT
FOR K=1 TO KS
MAT X=XX
CALL NEWTON(X,XX) !' xx=x-f(x)/f'(x) ニュートン法
IF QABS(XX,X)<EPS THEN !'収束判定
CALL PSET(ZR,ZI,MOD(K,7)+1)
EXIT FOR
END IF
NEXT K
NEXT ZI
NEXT ZR
NEXT I
END
EXTERNAL FUNCTION QABS(X(),Y())
OPTION ARITHMETIC NATIVE
FOR I=0 TO 3
LET A=A+(X(I)-Y(I))^2
NEXT I
LET QABS=SQR(A)
END FUNCTION
EXTERNAL SUB NEWTON (X(),XX())
OPTION ARITHMETIC NATIVE
OPTION BASE 0
DIM A(3,3),Y1(3),Y2(3)
RESTORE
FOR I=3 TO 0 STEP -1
FOR J=0 TO 3
READ A(I,J)
NEXT J
NEXT I
DATA 1,0,0,0 !'(1+0i+0j+0k)*X^3
DATA 0,0,0,0 !'(0+0i+0j+0k)*X^2
DATA 0,0,0,0 !'(0+0i+0j+0k)*X
DATA -1,0,0,0 !'(-1+0i+0j+0k)
CALL HORNER(3,A,X,Y1) !' f(X)=X^3-1
CALL DERIVATIVE(3,A) !' 微分
CALL HORNER(2,A,X,Y2) !' f'(X)=3*X^2
CALL QDIV(Y1,Y2) !' (X^3-1)/(3*X^2)
MAT XX=X
CALL QSUB(XX,Y1) !' XX=X-(X^3-1)/(3*X^2)
END SUB
EXTERNAL SUB QSET(XX(),X,Y,Z,W)
OPTION ARITHMETIC NATIVE
LET XX(0)=X
LET XX(1)=Y
LET XX(2)=Z
LET XX(3)=W
END SUB
EXTERNAL SUB HORNER(N,A(,),X(),Y())
OPTION ARITHMETIC NATIVE
OPTION BASE 0
DIM AA(3)
CALL QSET(Y,A(N,0),A(N,1),A(N,2),A(N,3))
FOR I=N-1 TO 0 STEP -1
CALL QMUL(Y,X) !'Y=Y*X
CALL QSET(AA,A(I,0),A(I,1),A(I,2),A(I,3))
CALL QADD(Y,AA) !' Y=Y+A
NEXT I
END SUB
EXTERNAL SUB QMUL(A(),B())
OPTION ARITHMETIC NATIVE
OPTION BASE 0
DIM SS(3)
LET SS(0)=A(0)*B(0)-A(1)*B(1)-A(2)*B(2)-A(3)*B(3)
LET SS(1)=A(0)*B(1)+A(1)*B(0)+A(2)*B(3)-A(3)*B(2)
LET SS(2)=A(0)*B(2)-A(1)*B(3)+A(2)*B(0)+A(3)*B(1)
LET SS(3)=A(0)*B(3)+A(1)*B(2)-A(2)*B(1)+A(3)*B(0)
MAT A=SS
END SUB
EXTERNAL SUB QDIV(A(),B())
OPTION ARITHMETIC NATIVE
OPTION BASE 0
DIM BB(3)
CALL QINV(BB,B)
CALL QMUL(A,BB)
MAT A=(1/(B(0)^2+B(1)^2+B(2)^2+B(3)^2))*A
END SUB
EXTERNAL SUB QADD(A(),B())
OPTION ARITHMETIC NATIVE
MAT A=A+B
END SUB
EXTERNAL SUB QSUB(A(),B())
OPTION ARITHMETIC NATIVE
MAT A=A-B
END SUB
EXTERNAL SUB QINV(ZZ(),Z())
OPTION ARITHMETIC NATIVE
LET ZZ(0)=Z(0)
LET ZZ(1)=-Z(1)
LET ZZ(2)=-Z(2)
LET ZZ(3)=-Z(3)
END SUB
EXTERNAL SUB DERIVATIVE(N,A(,)) !'微分
OPTION ARITHMETIC NATIVE
OPTION BASE 0
DIM B(N,3)
FOR I=N TO 1 STEP-1
FOR J=0 TO 3
LET B(I-1,J)=I*A(I,J)
NEXT J
NEXT I
FOR I=N TO 0 STEP -1
FOR J=0 TO 3
LET A(I,J)=B(I,J)
NEXT J
NEXT I
END SUB
EXTERNAL SUB PSET(X,Y,C)
OPTION ARITHMETIC NATIVE
SET POINT COLOR C
PLOT POINTS:X,Y
END SUB
EXTERNAL SUB GINIT(XSIZE,YSIZE)
OPTION ARITHMETIC NATIVE
SET BITMAP SIZE XSIZE,YSIZE
!'SET WINDOW 0 , XSIZE-1 , YSIZE-1, 0
SET POINT STYLE 1
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
SET COLOR MIX(255) .5,.5,.5
CLEAR
END SUB
|
|
|
|
投稿者：yoshipyuta
投稿日：2018年 1月10日(水)15時17分34秒
|
|
|
|
> No.4455[元記事へ]
しばっちさんへのお返事です。
> 更に2次元(複素数 x+yi (i^2=-1))ではなく、4次元(Quaternion x+yi+zj+wk (i^2=j^2=k^2=-1))でも求めてみました。
> 但し、3D表示は敷居が高く、表示方法や仕様設定等が未定なので表示は2Dです。
しばっちさんのレベルにわがKidsはついていけませんね。現在、KidsはNewton・フラクタルプログラム①の理解に邁進中です。初め当惑していたDiamond Neckless模様が正しいことを理解しつつあるようです。国内外でこのことにコメントしているのは大変少ない。
Newton・フラクタルのしばっちプログラム①を用いて現在、複素解の研究を続行しています。プログラムは苦手なKidsですがまずは、色使いを変更（ソフトに）しました。
４次元複素数の話はまた勉強させてください。それと十進BASICでは陰関数や複素関数の３D表示プログラムを見たことがないのですが、難しいのでしょうか？KidsはWolframAlphaを用いてf(z)=sin(z)などのReal PartやImaginary Partを表示させていますが、色使いなどには満足はしていないようです。
しばっちさん、下記のような３Dは十進BASICでは無理ですか？是非Challengeしてください。
|
|
|
|
投稿者：yoshipyuta
投稿日：2018年 1月10日(水)16時48分32秒
|
|
|
|
> No.4451[元記事へ]
しばっちさんへのお返事です。Kidsへの詳細なご返事感謝いたします。
>> > (z^4-1)*(Z^2-(1+i)^2)のように虚数iが式に入る場合はどうなりますか。
> 下記のようにCOMPLEX文 又はSQR(-1)を使ってください。
> (Z^4-1)*(Z^2-COMPLEX(1,1)^2)
Newton・しばっちプログラム②を良く見てKidsもCOMPLEX(1,1)表示を理解しました。z^8-17*z^4+16のプログラムの解の部分を±1,±i,±(1+1)、f(z)=(z^4-1)*(Z^2-(1+i)^2)に変えても正しい図（下記の図）は得られません？
LET BOTTOM=2.5
LET TOP=-2.5
の定義のため？±(1+1)が-1+i、+1-1に位置している。LET TOP=-2.5にしているので当たり前でしょうか。
> > Z＊(Z^3-1)の２枚葉の発芽の図を下記に添付します。x: 0.5～0.7 y: -0.12～0.12領域
> 計算区間にずれがあります。
> X:0.7-0.5=0.2
> Y:0.12-(-0.12)=0.24
>
> この時の表示画面が正方形なのに対し、X,Yの区間が違うため歪みが生じています。
> 区間幅を同じにするか、表示ウィンドウの上下左右の比率と合わせてください。
> 但し、SUB GETSQUAREで求める領域は正方形ウィンドウです(上下左右の比率が1)
> Y:-0.1～0.1 0.1-(-0.1)=0.2
> とするか、
SUB GETSQUAREの時の範囲に注意するようにします。
> >さて、さらに進んでTan(z)やCos(z)に展開できるようにできますか？
> 複素数モードでもEXP,LOG,SQR関数などが使用できますが、SIN,COS,TANなどは使えません。
> 下記のように定義すれば複素数関数が使用できるようになります。
f(z)=sin(z)外部関数はNewton・しばっちプログラム①で用いるようにしました（既に入れてありました）。Newton・しばっちプログラム②ではf(z)の複素数解が判明している場合のみ使えるものと理解しています。
|
|
|
戻る
