”しばっち”さん、z^8-17z^4+16のニュートン解

 投稿者:yoshipyuta  投稿日:2018年 1月 3日(水)23時03分46秒
  しばっちさん、立て続けののKidsの質問をお許し下さい。z^8-17*Z^4+16のニュートン・フラクタルの図示が面白くありません。この式は8つの解+-1、+-i、+-2、+-2iを有します。

下記の図でおおよその8rootsの位置とは思いますが、単一色での表示が残念です。カラー化できますか。また(z^4-1)*(Z^2-(1+i)^2)のように虚数iが式に入る場合はどうなりますか。
 

Re: ”しばっち”さん、z^8-17z^4+16のニュートン解

 投稿者:しばっち  投稿日:2018年 1月 6日(土)10時35分22秒
  > No.4444[元記事へ]

yoshipyutaさんへのお返事です。

> しばっちさん、立て続けののKidsの質問をお許し下さい。z^8-17*Z^4+16のニュートン・フラクタルの図示が面白くありません。この式は8つの解+-1、+-i、+-2、+-2iを有します。
> 下記の図でおおよその8rootsの位置とは思いますが、単一色での表示が残念です。カラー化できますか。

下記のようにしてみてください。
注意点として、計算精度(計算回数)を調整してください。特に拡大表示で再計算させた場合は計算回数不足により
収束しきれずに打ち切ってしまっている場合があるようです。(黒色は背景色のままで何も描画されていません)
また、使用している数値微分式や、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)の2枚葉の発芽の図を下記に添付します。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
 

Re: ”しばっち”さん、z^8-17z^4+16のニュートン解

 投稿者:しばっち  投稿日: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
 

Re: ”しばっち”さん、z^8-17z^4+16のニュートン解

 投稿者:しばっち  投稿日: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
 

Re: ”しばっち”さん、z^8-17z^4+16のニュートン解

 投稿者:しばっち  投稿日: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
 

Re: ”しばっち”さん、z^8-17z^4+16のニュートン解

 投稿者: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ですがまずは、色使いを変更(ソフトに)しました。

4次元複素数の話はまた勉強させてください。それと十進BASICでは陰関数や複素関数の3D表示プログラムを見たことがないのですが、難しいのでしょうか?KidsはWolframAlphaを用いてf(z)=sin(z)などのReal PartやImaginary Partを表示させていますが、色使いなどには満足はしていないようです。

しばっちさん、下記のような3Dは十進BASICでは無理ですか?是非Challengeしてください。
 

Re: ”しばっち”さん、z^8-17z^4+16のニュートン解

 投稿者: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)の2枚葉の発芽の図を下記に添付します。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)の複素数解が判明している場合のみ使えるものと理解しています。
 

戻る