根軌跡図の作成

　　根軌跡図の作成 島村1243 2007/06/23 10:05:45  (修正3回)
　　├プログラムから推測して、その処理を検討し... 山中和義 2007/06/23 20:33:56  (修正1回)
　　│└runしたら完璧な作図 島村1243 2007/06/23 21:52:14 
　　└!根軌跡図の意味を考える。 SECOND 2007/06/25 05:05:43  (修正3回)
　　　├実係数代数方程式の根（複素数）を複素平面... 山中和義 2007/06/25 12:02:22 
　　　│├ディジタル走査の、宿命的な問題ですね。 SECOND 2007/06/25 19:53:01 
　　　││└!貴兄の３Ｄ部分に習って、書き直してみまし... SECOND 2007/06/26 13:04:43  (修正2回)
　　　││　└変数引数について。 SECOND 2007/06/27 23:10:08 
　　　│└３Ｄで見る代数方程式の解（複素数） 山中和義 2007/06/27 12:05:59  (修正2回)
　　　│　└プログラムはこちら 山中和義 2007/06/27 12:10:25 
　　　│　　└もう１つ追加しておきます。 山中和義 2007/06/27 15:54:01  (修正1回)
　　　├補足、ＤＫＡの構造。（例は、５次の根を求... SECOND 2007/06/28 09:09:32  (修正1回)
　　　└高次方程式の複素数根を求めるのが本来の姿... 島村1243 2007/06/29 23:26:33  (修正1回)
　　　　└これからも、掲示を頂ければ、よろしく御願... SECOND 2007/07/01 01:04:50  (修正1回)

　　根軌跡図の作成 島村1243 2007/06/23 10:05:45  (修正3回) ツリーへ

根軌跡図の作成 返事を書くノートメニュー

島村1243 <bjllmpcujp> 2007/06/23 10:05:45 ** この記事は3回修正されてます

自動制御システムの特性評価方法の一つに、特性方程式の根軌跡を調べると言うEvansの方法が
有ります。
自動制御のテキストが教えるところにより、軌跡図が大体どの様になるかの推測は出来ていましたが、
手作業での作図は膨大な手間がかかるので実際に軌跡図がどの様になるかを作図したことは無く、頭
の中で、大体のイメージを想像する範囲を脱しませんでしたので正しい軌跡を把握した事は有りませ
んでした。
ところが複素数計算を扱え、かつ、その結果を簡単に作図出来る１０進BASICを知ったので、早速
応用してみました。（Vine4.0,FedoraCore5,Devian3.1）

システムの特性方程式が「1+K/s/(s+1)/(s+4)=0」で表される根軌跡図の求め方は、複素数平面上
に任意の点zを打って
点(0,0)からzに向かうベクトルの位相角をkaku0を求める。
点(-1,0)からzに向かうベクトルの位相角をkaku1を求める。
点(-4,0)からzに向かうベクトルの位相角をkaku2を求める。
kaku0+kaku1+kaku2の合計が180度または-180度になれば、そのzは根軌跡上に存在する。
と言うものです。

（１）x=-1からx=（約）-0.46迄は y=0の赤い直線が描かれる。
（２）x=（約）-0.46からx=0迄は y>0の赤い曲線が描かれ、x=0でy=2になる。

（３）x=0からx=（約）-0.46迄は y=0の青い直線（原点から-0.46に向かう直線）
（４）x=（約）-0.46からx=0迄は y<0の青い曲線が描かれ、x=0でy=-2になる。

が理論上正解です。（この他に、x=-4からx=-∞迄y=0の直線もあるが、作図の必要が無いのでプロ
グラムは作成省略してある）。

----プログラム-------------
OPTION ARITHMETIC COMPLEX
SET WINDOW -1.5,0.5,-2.5,2.5
DRAW grid
LET j=complex(0,1)
LET z0=0
LET z1=-1
LET z2=-4
SET LINE COLOR "red"
FOR y=0 TO 2.5 STEP 0.001
FOR x=-1.5 TO 0.5 STEP 0.001
LET z=x+j*y
LET kaku0=arg(z-z0)
LET kaku1=arg(z-z1)
LET kaku2=arg(z-z2)
LET kaku=kaku0+kaku1+kaku2
IF kaku<=PI THEN EXIT FOR
NEXT x
PLOT LINES:x,y;
NEXT y
PLOT LINES
SET LINE COLOR "blue"
FOR y=-2.5 TO 0 STEP 0.001
FOR x=-1 TO 0.5 STEP 0.001
LET z=x+j*y
LET kaku0=arg(z-z0)
LET kaku1=arg(z-z1)
LET kaku2=arg(z-z2)
LET kaku=kaku0+kaku1+kaku2
IF kaku>=-PI THEN EXIT FOR
NEXT x
PLOT LINES:x,y;
NEXT y
END
--------------------------------

プログラムをrunして得られる（４）の曲線は（２）の曲線に対して、実軸に対して線対称になって
いれば正しい結果ですのでこのプログラムで得られる図は（１）（２）（４）は正しい結果を表して
います。

＜欠点＞
（ア）複素数平面上に任意のzを設定する方法として、For-Nextで小刻みな分割計算をして180度または
-180度を越えたらForを抜けて線画する、としているので、正確さに欠けるという問題があるのと、
作図に時間がかかるのが欠点です。
x方向の計算に挟み撃ち法等を使えば、より正確で計算時間を短く出来ると思うのですが適用していま
せん。y方向については計算短縮の上手い方法が考え付きません。
（イ）（３）については同一のx値に対してyの値が２個存在する事になるので作図が出来ない状態
で、解答図としては不備になります。

以上について何か上手い方法が有りましたらご教示頂けると有り難いです。

　　├プログラムから推測して、その処理を検討し... 山中和義 2007/06/23 20:33:56  (修正1回) ツリーへ

Re: 根軌跡図の作成 返事を書くノートメニュー

山中和義 <drdlxujciw> 2007/06/23 20:33:56 ** この記事は1回修正されてます

プログラムから推測して、その処理を検討してみました。

>複素数平面上に任意のzを設定する方法

グラフィックス画面のピクセル（ドット）数を対象にすればよいかと思います。
実際は、実数x,yについて計算する訳ですので、精度をあげるため２倍の数で行っています。

OPTION ARITHMETIC complex

LET Xa=-5 !X=[-5,1]
LET Xb=1
LET Ya=-3 !Y=[-3,3]
LET Yb=3

SET WINDOW Xa,Xb,Ya,Yb !描画範囲
DRAW grid

ASK PIXEL SIZE (Xa,Ya; Xb,Yb) a,b !ドット換算
LET c_Div=2 !分解数
LET dx=(Xb-Xa)/a/c_Div !描画範囲内の描画対象となるドット数を得る
LET dy=(Yb-Ya)/b/c_Div

LET cEps=1E-3 !誤差　※精度が気になる!?

LET j=SQR(-1) !虚数単位

SET POINT STYLE 5

LET z0=0 !極 0
PLOT POINTS: re(z0),im(z0)
LET z1=-1 !-1
PLOT POINTS: re(z1),im(z1)
LET z2=-4 !-4
PLOT POINTS: re(z2),im(z2)

SET POINT STYLE 1

FOR k=0 TO c_Div*b-1 !Ｙ座標　※ドットをもとに
LET y=Ya+dy*k ![Ya,Yb]
FOR i=0 TO c_Div*a-1 !Ｘ座標
LET x=Xa+dx*i ![Xa,Xb]

LET z=x+y*j !空間上の点を計算する

LET a0=GetAngl(z,z0) !極とのなす角を得る
LET a1=GetAngl(z,z1)
LET a2=GetAngl(z,z2)
LET sumA=ABS(a0+a1+a2)
!!!PRINT x;ABS(sumA)-PI

IF ABS(sumA-PI) < cEps THEN PLOT POINTS: x,y !点画　※(2k+1)*PI　k=0,±1,±2,…
IF ABS(sumA-3*PI) < cEps THEN PLOT POINTS: x,y
NEXT i
NEXT k

FUNCTION GetAngl(z,z9) !なす角を得る
WHEN EXCEPTION IN
LET GetAngl=arg(z-z9)
USE
LET GetAngl=0 !同一点なら
END WHEN
END FUNCTION

END

　　│└runしたら完璧な作図 島村1243 2007/06/23 21:52:14  ツリーへ

Re: プログラムから推測して、その処理を検討し... 返事を書くノートメニュー

島村1243 <bjllmpcujp> 2007/06/23 21:52:14

runしたら完璧な作図

作図に省略部分の無い完璧な解析プログラム！ですね。この様な作図が出来るプログラムが考えられる
なんて凄いですね。（私の書いたプログラムは余りにも単純なのでお恥ずかしい限りです）
山中様の作成されたコード内容を直ぐには理解できないレベルですが、どんなテクニックを使用して
いるのか、久し振りに胸をわくわくしながら読ませて頂きます。ご教示有難う御座いました。

　　└!根軌跡図の意味を考える。 SECOND 2007/06/25 05:05:43  (修正3回) ツリーへ

Re: 根軌跡図の作成 返事を書くノートメニュー

SECOND <cszcthjjdj> 2007/06/25 05:05:43 ** この記事は3回修正されてます

!根軌跡図の意味を考える。

!　1+K/s/(s+1)/(s+4)=0　を変形して→　伝達関数F(s)=　K/s/(s+1)/(s+4)　= -1　とすると、

!信号入出力間の、伝達関数F(s)が、-1 に固定される時の、（反転して、閉ループ接続される時の）
!反転増幅率　K（実数０～∞）に対する解としての、s(複素数) の軌跡です。

!この軌跡上の s(複素数)は、閉ループ内に、信号電流ｉ(t)= exp^st が、自立出来る事を示します。

!何故そう言えるか、伝達関数F(s) の構成素子としてのインピーダンス Z(s) は、

!定常的な、交流インピーダンス　jωL　や、1/jωC　のような　Z(jω) が、
!電流ｉ(t)= exp^jωt　を想定して成立つ電圧降下率、インピーダンス　Z(jω)　であるように

!Z(s) も同様に、過渡的変化を含んでは、いますが、
!電流ｉ(t)= exp^(σ+jω)t　=exp^ｓt　を想定して成立つ電圧降下率、インピーダンス　です。

!そして、インピーダンス　Z(s)で組み立てられた伝達関数F(s) も又、
!信号電流ｉ(t)= exp^ｓt　を想定したときの入出力間、信号電圧比、だからです。

!閉ループ内に、信号電流ｉ(t)= exp^st が、自立出来るという事は、そこに対応している
!反転増幅率 K（実数０～∞）において、入力無しの制御系が、勝手に　exp^st　で動く
!意味を持ちますので、s=σ+jω の実数部が、
!正であれば、端へぶつかって戻らなくなったり、角速度ωの周波数で、振動増大したりします。
!負であっても、０に近ければ、減衰振動、鈍重さ、等を意味します。

!根軌跡図の　s=σ+jω に対応している K で、ｓの範囲を制限できます。

!s*(s+1)*(s+4)+K=0　　　s^3 +5*s^2 +4*s +K=0

!　ＤＫＡ法で直接的に、根の近似値を求める。
!
!　K=0 ～ K=40 まで、根ｓのグラフは、Ｘ印（ｓ= 0,-1,-4) から出発し、対応している。
!-------------
OPTION ARITHMETIC COMPLEX

LET N=3
DIM A(N),Xr(N)
LET A(1)=5
LET A(2)=4

SET WINDOW -6,+1,-3,+3
DRAW grid

SET POINT STYLE 5 !Ｘ点
FOR K=0 TO 40 STEP 0.02
LET A(3)=K
CALL DKA_00
SET POINT COLOR 2 !青
PLOT POINTS:re(Xr(1)),im(Xr(1))
SET POINT COLOR 1 !黒
PLOT POINTS:re(Xr(2)),im(Xr(2))
SET POINT COLOR 4 !赤
PLOT POINTS:re(Xr(3)),im(Xr(3))
SET POINT STYLE 1 !・点
NEXT K

!DKA (Durand Kerner Aberth) 法による高次方程式の根
!--------------------------------------------
! X^n + A(1)*X^(n-1) + ... + A(n) = 0

! Xr( 1 ~ n ) <== root
!--------------------------------------------

SUB DKA_00
LET r=1
FOR j=2 TO N
LET rn=ABS(A(j))^(1/j)
if r<rn then LET r=rn
NEXT j
FOR j=1 TO N
LET Xr(j)=-A(1)/N+r*EXP( complex(0,1)*2*PI/N *(j-3/4) )
NEXT j
FOR m=0 TO 100
LET mfx=0
LET maj=0
FOR j=1 TO N
LET Xk=1
LET fx=1
FOR w=1 TO N
LET fx=fx*Xr(j)+A(w)
IF w<>j THEN LET Xk=Xk*(Xr(j)-Xr(w))
NEXT w
LET Xr(j)=Xr(j)-fx/Xk
IF mfx<ABS(fx) THEN LET mfx=ABS(fx)
IF maj<ABS(fx/Xk) THEN LET maj=ABS(fx/Xk)
NEXT j
IF mfx<.0000001 AND maj<.0000001 THEN EXIT FOR
NEXT m
END SUB

END

　　　├実係数代数方程式の根（複素数）を複素平面... 山中和義 2007/06/25 12:02:22  ツリーへ

Re: !根軌跡図の意味を考える。 返事を書くノートメニュー

山中和義 <drdlxujciw> 2007/06/25 12:02:22

実係数代数方程式の根（複素数）を複素平面でグラフ化する

付属サンプル（SAMPLEフォルダ内）3DPLOT.BASをベースに３Ｄ表示してみました。
実係数代数方程式の根（複素数）は、突起部分になります。

f(s)=s*(s+1)*(s+4) + K　で、K=0～60でアニメーションさせてみると、、、

根の軌跡は、このグラフを上から（ＸＹ平面を）見ることになります。

突起の頂点が高低するのは、
FOR x=-6 TO 1 STEP 0.05 !複素平面をスキャンする
FOR y=-3 TO 3 STEP 0.05
で、うまく頂点を算出できないためです。
刻み幅を小さくすれば、該当する可能性が高くなりますが、アニメーションに負荷がかかります。

OPTION ARITHMETIC complex

LET i=SQR(-1) !虚数単位

SUB rotx(x,y,z,a) !回転
LET y0=y*cos(a)-z*sin(a)
LET z0=y*sin(a)+z*cos(a)
LET y=y0
LET z=z0
END SUB
SUB roty(x,y,z,a)
LET x0=x*cos(a)+z*sin(a)
LET z0=-x*sin(a)+z*cos(a)
LET x=x0
LET z=z0
END SUB
SUB rotz(x,y,z,a)
LET x0=x*cos(a)-y*sin(a)
LET y0=x*sin(a)+y*cos(a)
LET x=x0
LET y=y0
END SUB

SUB convert(x,y,z) !投影
CALL rotz(x,y,z,-30*PI/180)
CALL rotx(x,y,z,-60*PI/180)
LET x=x+2
LET y=y-1
END SUB

SUB plotTo(x,y,z)
LET x1=x
LET y1=y
LET z1=z
CALL convert(x1,y1,z1)
PLOT LINES: x1,y1;
END SUB

DEF f(s)=s*(s+1)*(s+4) + K !実係数代数方程式
!DEF f(s)=s^3+5*s^2+4*s + K

SET WINDOW -4,4,-4,4

!　z　y　で表示する
!　│／
!　・─x

FOR K=0 TO 40 STEP 0.25
SET DRAW MODE HIDDEN !描画開始
CLEAR

FOR x=-6 TO 1 STEP 0.05 !複素平面をスキャンする　x=[-6,1]
FOR y=-3 TO 3 STEP 0.05 !y=[-3,3]
LET z=ABS(1/f(x+y*i)) !関数の逆数の絶対値を描く
CALL PlotTo(x,y,z) !突起部分が解（複素数）になる
NEXT y
PLOT LINES
NEXT x
PLOT TEXT,AT -2.5,-2.5: "K="&STR$(K)

SET DRAW MODE EXPLICIT !描画終了
NEXT K

END

　　　│├ディジタル走査の、宿命的な問題ですね。 SECOND 2007/06/25 19:53:01  ツリーへ

Re: 実係数代数方程式の根（複素数）を複素平面... 返事を書くノートメニュー

SECOND <cszcthjjdj> 2007/06/25 19:53:01

ディジタル走査の、宿命的な問題ですね。
私も、刀尽き、矢折れました。グリッドが根に重なる有効な手立ては？ない！
その末、先輩たちに、泣きついてしまいました。

３Ｄにすると高さ方向に、Ｋも表現できて、すごい図になりそうですね。

　　　││└!貴兄の３Ｄ部分に習って、書き直してみまし... SECOND 2007/06/26 13:04:43  (修正2回) ツリーへ

Re: ディジタル走査の、宿命的な問題ですね。 返事を書くノートメニュー

SECOND <cszcthjjdj> 2007/06/26 13:04:43 ** この記事は2回修正されてます

!貴兄の３Ｄ部分に習って、書き直してみました。

OPTION ARITHMETIC complex

! s^3 +5*s^2 +4*s +K =0　の根のグラフ
LET N=3
DIM A(N),Xr(N)
LET A(1)=5
LET A(2)=4 !A(3)=K

!------------------
!　ｚ　ｙ　で表示する
!　│／
!　・─ｘ
!------------------
SUB rotx(x,y,z,a)
LET w=y*COS(a)-z*SIN(a)
LET z=y*SIN(a)+z*COS(a)
LET y=w
END SUB

SUB rotz(x,y,z,a)
LET w=x*COS(a)-y*SIN(a)
LET y=x*SIN(a)+y*COS(a)
LET x=w
END SUB

SUB plots(x,y,z)
CALL rotz(x,y,z,-PI/6)
CALL rotx(x,y,z,-PI/3)
PLOT LINES:x,y;
END SUB

SUB plott(x,y,z,w$)
CALL rotz(x,y,z,-PI/6)
CALL rotx(x,y,z,-PI/3)
PLOT TEXT,AT x,y:w$
END SUB

SET WINDOW -6,2,-2.5,5.5

!-----目盛り
FOR x=-6 TO 1
IF x=0 THEN SET LINE STYLE 1 ELSE SET LINE STYLE 3
CALL plots(x+0,-3, 0)
CALL plots(x+0,+3, 0) !　x+0 は、SUB からｘの書き戻しを止める。
PLOT LINES
NEXT x
FOR y=-3 TO 3
IF y=0 THEN SET LINE STYLE 1 ELSE SET LINE STYLE 3
CALL plots(-6,y+0, 0)
CALL plots( 1,y+0, 0) !　y+0 は、SUB からｙの書き戻しを止める。
PLOT LINES
NEXT y
SET LINE STYLE 1

!-----Ｋと根のグラフ
SET TEXT BACKGROUND "OPAQUE"
FOR K=0 TO 40 STEP 0.02
LET A(3)=K
CALL DKA_00
FOR s=1 TO 3
SET LINE COLOR s+1
CALL plots(re(Xr(s)),im(Xr(s)), 0)
CALL plots(re(Xr(s)),im(Xr(s)),K/10)
PLOT LINES
NEXT s
!----Ｋの印字
CALL plott(-6.2, 0, 0.65+INT(K/10), "K="&USING$("##",K) )
NEXT K

!DKA (Durand Kerner Aberth) 法による高次方程式の根
!--------------------------------------------
! X^n + A(1)*X^(n-1) + ... + A(n) = 0

! Xr( 1 ~ n ) <== root
!--------------------------------------------
SUB DKA_00
LET r=1
FOR j=2 TO N
LET rn=ABS(A(j))^(1/j)
if r<rn then LET r=rn
NEXT j
FOR j=1 TO N
LET Xr(j)=-A(1)/N+r*EXP( 2*PI*complex(0,1)/N *(j-3/4) )
NEXT j
FOR m=0 TO 100
LET mfx=0
LET maj=0
FOR j=1 TO N
LET Xk=1
LET fx=1
FOR w=1 TO N
LET fx=fx*Xr(j)+A(w)
IF w<>j THEN LET Xk=Xk*(Xr(j)-Xr(w))
NEXT w
LET Xr(j)=Xr(j)-fx/Xk
IF mfx<ABS(fx) THEN LET mfx=ABS(fx)
IF maj<ABS(fx/Xk) THEN LET maj=ABS(fx/Xk)
NEXT j
IF mfx<.0000001 AND maj<.0000001 THEN EXIT FOR
NEXT m
END SUB

END

　　　││　└変数引数について。 SECOND 2007/06/27 23:10:08  ツリーへ

Re: !貴兄の３Ｄ部分に習って、書き直してみまし... 返事を書くノートメニュー

SECOND <cszcthjjdj> 2007/06/27 23:10:08

変数引数について。

プログラムで使用している方法、これは、ＢＡＳＩＣに負担がかかって、良くないです。
　（
　　）
　CALL plots(x+0,+3, 0) !　x+0 は、SUB からｘの書き戻しを止める。
　（
　　）

ヘルプ(H)　目次(C)　→　手続き定義　→　変数引数
　引数が変更されないようにしたいときは，
　引数を括弧で括るなどして変数として解釈されないようにする。
　例 130 CALL exchange((a),(b))

この方が正しく、ＢＡＳＩＣ側に演算負荷がかかりません。
　（
　　）
　CALL plots((x),+3, 0) !　(x) は、SUB からｘの書き戻しを止める。
　（
　　）

　　　│└３Ｄで見る代数方程式の解（複素数） 山中和義 2007/06/27 12:05:59  (修正2回) ツリーへ

Re: 実係数代数方程式の根（複素数）を複素平面... 返事を書くノートメニュー

山中和義 <drdlxujciw> 2007/06/27 12:05:59 ** この記事は2回修正されてます

３Ｄで見る代数方程式の解（複素数）

通常は、ＸＹ座標のグラフ（本プログラムでは左下から右上方向の）を描いて、Ｘ軸との交点が解　云々・・・
「奥ゆき」をＸ虚部とした３Ｄグラフにすると、、、

X=-1、X=0からの２つ解は、重根になって、一組の複素数（共役）へと変化します。
まるで、衝突して（実数世界から）、飛び散る（複素数世界へ）みたいです。

　　　│　└プログラムはこちら 山中和義 2007/06/27 12:10:25  ツリーへ

Re: ３Ｄで見る代数方程式の解（複素数） 返事を書くノートメニュー

山中和義 <drdlxujciw> 2007/06/27 12:10:25

プログラムはこちら

OPTION ARITHMETIC complex

LET i=SQR(-1) !虚数単位
LET cEps=1e-6 !誤差

LET N=3 !代数方程式　s*(s+1)*(s+4) + K = 0
DIM A(N)
LET A(1)=5 !多項式 X^N+A(1)*X^(N-1)+A(2)*X^(N-2)+ … +A(N-1)*X+A(N)の係数
LET A(2)=4
!!LET A(3)=0

DIM Xr(N)
SUB DKA_00 !DKA法（Durand Kerner Aberth）による代数方程式の解
LET r=1 !初期値を仮定する
FOR j=2 TO N
LET rn=ABS(A(j))^(1/j)
if r<rn then LET r=rn
NEXT j
FOR j=1 TO N !半径rの円に等間隔に配置する
LET Xr(j)=-A(1)/N+r*EXP( 2*PI*i/N *(j-3/4) ) !アーバスの初期値
NEXT j

FOR m=0 TO 100 !反復
LET mfx=0
LET maj=0
FOR j=1 TO N
LET Xk=1
LET fx=1
FOR w=1 TO N
LET fx=fx*Xr(j)+A(w)
IF w<>j THEN LET Xk=Xk*(Xr(j)-Xr(w))
NEXT w
LET Xr(j)=Xr(j)-fx/Xk
IF mfx<ABS(fx) THEN LET mfx=ABS(fx)
IF maj<ABS(fx/Xk) THEN LET maj=ABS(fx/Xk)
NEXT j
IF mfx<cEps AND maj<cEps THEN EXIT FOR !収束したら
NEXT m
END SUB

SUB rotx(x,y,z,a) !回転
LET y0=y*cos(a)-z*sin(a)
LET z0=y*sin(a)+z*cos(a)
LET y=y0
LET z=z0
END SUB
SUB roty(x,y,z,a)
LET x0=x*cos(a)+z*sin(a)
LET z0=-x*sin(a)+z*cos(a)
LET x=x0
LET z=z0
END SUB
SUB rotz(x,y,z,a)
LET x0=x*cos(a)-y*sin(a)
LET y0=x*sin(a)+y*cos(a)
LET x=x0
LET y=y0
END SUB

SUB convert(x,y,z) !投影
CALL rotz(x,y,z,-30*PI/180)
CALL rotx(x,y,z,-60*PI/180)
LET x=x+2
LET y=y-2
END SUB

SUB plotTo(x,y,z)
LET x1=x
LET y1=y
LET z1=z
CALL convert(x1,y1,z1)
PLOT LINES: x1,y1;
END SUB
SUB PlotText(x,y,z,s$)
CALL convert(x,y,z)
PLOT TEXT ,AT x,y: s$
END SUB

SET WINDOW -4,4,-4,4

FOR K=0 TO 40 STEP 0.2
SET DRAW MODE HIDDEN !描画開始
CLEAR

!　y　x虚部　で表示する
!　│／
!　・─x実部

LET A(3)=K !定数項

FOR x=-6 TO 1 !目盛り
CALL PlotText(x+0,0,0,STR$(x))
NEXT x
CALL PlotTo(-6,0,0)
CALL PlotTo(1,0,0)
PLOT LINES

FOR y=-3 TO 3 !y=[-3,3]
CALL PlotText(0,y+0,0,STR$(y))
NEXT y
CALL PlotTo(0,-3,0)
CALL PlotTo(0,3,0)
PLOT LINES

FOR x=-6 TO 1 STEP 0.05 !関数を描く　x=[-6,1]
LET y=1 !Horner法（多項式）
FOR j=1 TO N
LET y=y*x+A(j)
NEXT j
CALL PlotTo(x,0,y/3) !※高さ調整
NEXT x
PLOT LINES

CALL DKA_00 !根を描く
FOR j=1 TO N
LET x1=re(Xr(j))
LET y1=im(Xr(j))
LET z1=0
CALL convert(x1,y1,z1)
DRAW disk WITH SCALE(0.1)*SHIFT(x1,y1)
NEXT j

PLOT TEXT,AT -2.5,-2.5: "K="&USING$("###.#",K)

SET DRAW MODE EXPLICIT !描画終了
WAIT DELAY 0.1
NEXT K

END

　　　│　　└もう１つ追加しておきます。 山中和義 2007/06/27 15:54:01  (修正1回) ツリーへ

Re: プログラムはこちら 返事を書くノートメニュー

山中和義 <drdlxujciw> 2007/06/27 15:54:01 ** この記事は1回修正されてます

もう１つ追加しておきます。

３つの解は、まるで
　左側が凹にめくれあがった、右側が凸に覆い被さったような曲面を
転げ落ちるみたいです。

OPTION ARITHMETIC complex

LET i=SQR(-1) !虚数単位
LET cEps=1e-6 !誤差

LET N=3 !代数方程式　s*(s+1)*(s+4) + K = 0
DIM A(N) !多項式 X^N+A(1)*X^(N-1)+A(2)*X^(N-2)+ … +A(N-1)*X+A(N)
LET A(1)=5 !係数
LET A(2)=4
!!LET A(3)=0

　　　：
　　　：
サブルーチン部分は同じなので、省略します。
　　　：
　　　：

SET WINDOW -4,4,-4,4

FOR K=0 TO 40 STEP 0.2
SET DRAW MODE HIDDEN !描画開始
CLEAR

LET A(3)=K !定数項

!　y　x虚部　で表示する
!　│／
!　・─x実部

FOR x=-6 TO 1 !目盛り
CALL PlotText(x+0,0,0,STR$(x))
NEXT x
CALL PlotTo(-6,0,0) !軸
CALL PlotTo(1,0,0)
PLOT LINES

FOR y=-3 TO 3 !y=[-3,3]
CALL PlotText(0,y+0,0,STR$(y))
NEXT y
CALL PlotTo(0,-3,0)
CALL PlotTo(0,3,0)
PLOT LINES

FOR x=-6 TO 1 STEP 0.2 !複素平面をスキャンする　x=[-6,1]
FOR y=-3 TO 3 STEP 0.2 !y=[-3,3]
LET Z=1 !Horner法（多項式）
FOR j=1 TO N
LET Z=Z*(x+y*i)+A(j)
NEXT j
CALL PlotTo(x,y,re(Z/10)) !※高さを調整する
NEXT y
PLOT LINES
NEXT x

CALL DKA_00 !解を描く
FOR j=1 TO N
LET x1=re(Xr(j))
LET y1=im(Xr(j))
LET z1=0
CALL convert(x1,y1,z1)
DRAW disk WITH SCALE(0.1)*SHIFT(x1,y1)
NEXT j

PLOT TEXT,AT -2.5,-2.5: "K="&USING$("###.#",K)

SET DRAW MODE EXPLICIT !描画終了
WAIT DELAY 0.1
NEXT K

END

　　　├補足、ＤＫＡの構造。（例は、５次の根を求... SECOND 2007/06/28 09:09:32  (修正1回) ツリーへ

Re: !根軌跡図の意味を考える。 返事を書くノートメニュー

SECOND <cszcthjjdj> 2007/06/28 09:09:32 ** この記事は1回修正されてます

補足、ＤＫＡの構造。（例は、５次の根を求める場合。）

　　　　　　　　真のｆ(x)＝ (x-真1) (x-真2) (x-真3) (x-真4) (x-真5)　　　真の根1~5

　　　　　　　　近似ｆ(x)＝ (x-近1) (x-近2) (x-近3) (x-近4) (x-近5)　　　近似根1~5

　とするとき、近似ｆ(x)　を、真のｆ(x)　に近づけるために

　プログラムは、以下の操作を繰返している。　ニュートン法に似ては、いるが、違う。

　　　　　　　　　　　　　　(近1-真1) (近1-真2) (近1-真3) (近1-真4) (近1-真5)
次の近1　＝　近1　－　────────────────────────────
　　　　　　　　　　　　　　　　　　　　(近1-近2) (近1-近3) (近1-近4) (近1-近5)

　　　　　　　　　　　　　　(近2-真1) (近2-真2) (近2-真3) (近2-真4) (近2-真5)
次の近2　＝　近2　－　────────────────────────────
　　　　　　　　　　　　　　　(近2-近1)　　　　　(近2-近3) (近2-近4) (近2-近5)

　　　　　　　　　　　　　　(近3-真1) (近3-真2) (近3-真3) (近3-真4) (近3-真5)
次の近3　＝　近3　－　────────────────────────────
　　　　　　　　　　　　　　　(近3-近1) (近3-近2)　　　　　(近3-近4) (近3-近5)

　　　　　　　　　　　　　　(近4-真1) (近4-真2) (近4-真3) (近4-真4) (近4-真5)
次の近4　＝　近4　－　────────────────────────────
　　　　　　　　　　　　　　　(近4-近1) (近4-近2) (近4-近3)　　　　　(近4-近5)

　　　　　　　　　　　　　　(近5-真1) (近5-真2) (近5-真3) (近5-真4) (近5-真5)
次の近5　＝　近5　－　────────────────────────────
　　　　　　　　　　　　　　　(近5-近1) (近5-近2) (近5-近3) (近5-近4)

この部分が、上記で、ｊ＝１～５について、繰返される。

! (N=5)
LET Xk=1
LET fx=1
FOR w=1 TO N
　LET fx=fx*Xr(j)+A(w)　　　　　・・・・・・・分子側、真のf(近ｊ)を、計算
　IF w<>j THEN LET Xk=Xk*(Xr(j)-Xr(w))　　・・分母側、近似f(近ｊ)から、(近ｊ-近ｊ)を除く。
NEXT w
LET Xr(j)=Xr(j)-fx/Xk　　　　　　・・・次の近ｊ

ｊ＝１～５について、真のf( 近似根ｊ)≒０になったら出来上がり。
プログラムでは、回毎の差分の検査まで、併用しているが、省略可。

以下は、最初にセットされる、近似根1~5　の初期値。
複素平面( -A(1)/5, 0 )を中心とし、半径r= max( A(1), A(2)^(1/2), …, A(5)^(1/5) )
の円周上に、均等角で、配置されています。かなり、ラフでも、構わない事になっていますが、
なぜ、この初期値が、好ましいかは、Durand Kerner Aberth　で検索、御調べ下さい。

! (N=5)
LET r=1
FOR j=2 TO N
　LET rn=ABS(A(j))^(1/j)
　if r<rn then LET r=rn
NEXT j
FOR j=1 TO N
　LET Xr(j)=-A(1)/N+r*EXP( complex(0,1)*2*PI/N *(j-3/4) )
NEXT j

　　　└高次方程式の複素数根を求めるのが本来の姿... 島村1243 2007/06/29 23:26:33  (修正1回) ツリーへ

Re: !根軌跡図の意味を考える。 返事を書くノートメニュー

島村1243 <bjllmpcujp> 2007/06/29 23:26:33 ** この記事は1回修正されてます

高次方程式の複素数根を求めるのが本来の姿だった。

山中さん、SECONDさん、ご教示有り難うございます。

根軌跡図の作成は、一般に高次特性方程式の根を直接求めるのは難しいので、根が
どの様に複素平面上を動くのかを幾何学計算（位相角の算術計）から知る唯一の簡
便法、と思っていました。
そのため、特性方程式の根を計算で求めてそれを複素数平面上にプロットするのが
根軌跡作成の本筋、と言うのをすっかり見失っていましたが、山中さんとSECONDさ
んの解説を見て目が覚めた思いです。

山中さんとSECONDさんの解説で、高次方程式の解法にDKA法と言うのがある事を初
めて知りました。そして１０進BASICを使うと複素数根も容易に求められると言う
事も知ることが出来ました。

Kと根軌跡の関係が一目で分かるSECONDさんの３D作画結果は、図を媒体として
根軌跡の意味が良く理解できると言う感想を持ちました。

山中さんが作成された最後のプログラムが描く図は、複素数の世界が持っている幻
想的な何かを垣間見せている、と言う印象を持ちました。

高度な知識を必要とするプログラムを見せて頂き、大変勉強になりました。有難う
御座いました。

　　　　└これからも、掲示を頂ければ、よろしく御願... SECOND 2007/07/01 01:04:50  (修正1回) ツリーへ

Re: 高次方程式の複素数根を求めるのが本来の姿... 返事を書くノートメニュー

SECOND <cszcthjjdj> 2007/07/01 01:04:50 ** この記事は1回修正されてます

これからも、掲示を頂ければ、よろしく御願いします。
人のあら探しばかりして、自分では、何一つ、自立する掲示のできない時間が、
ずっと続いています。恐縮いたします。

インデックスへ EXIT
新規発言を反映させるにはブラウザの更新ボタンを押してください。

根軌跡図の作成	返事を書くノートメニュー
島村1243 <bjllmpcujp> 2007/06/23 10:05:45 ** この記事は3回修正されてます
自動制御システムの特性評価方法の一つに、特性方程式の根軌跡を調べると言うEvansの方法が有ります。自動制御のテキストが教えるところにより、軌跡図が大体どの様になるかの推測は出来ていましたが、手作業での作図は膨大な手間がかかるので実際に軌跡図がどの様になるかを作図したことは無く、頭の中で、大体のイメージを想像する範囲を脱しませんでしたので正しい軌跡を把握した事は有りませんでした。ところが複素数計算を扱え、かつ、その結果を簡単に作図出来る１０進BASICを知ったので、早速応用してみました。（Vine4.0,FedoraCore5,Devian3.1）システムの特性方程式が「1+K/s/(s+1)/(s+4)=0」で表される根軌跡図の求め方は、複素数平面上に任意の点zを打って点(0,0)からzに向かうベクトルの位相角をkaku0を求める。点(-1,0)からzに向かうベクトルの位相角をkaku1を求める。点(-4,0)からzに向かうベクトルの位相角をkaku2を求める。 kaku0+kaku1+kaku2の合計が180度または-180度になれば、そのzは根軌跡上に存在する。と言うものです。（１）x=-1からx=（約）-0.46迄は y=0の赤い直線が描かれる。（２）x=（約）-0.46からx=0迄は y>0の赤い曲線が描かれ、x=0でy=2になる。（３）x=0からx=（約）-0.46迄は y=0の青い直線（原点から-0.46に向かう直線）（４）x=（約）-0.46からx=0迄は y<0の青い曲線が描かれ、x=0でy=-2になる。が理論上正解です。（この他に、x=-4からx=-∞迄y=0の直線もあるが、作図の必要が無いのでプログラムは作成省略してある）。 ----プログラム------------- OPTION ARITHMETIC COMPLEX SET WINDOW -1.5,0.5,-2.5,2.5 DRAW grid LET j=complex(0,1) LET z0=0 LET z1=-1 LET z2=-4 SET LINE COLOR "red" FOR y=0 TO 2.5 STEP 0.001 FOR x=-1.5 TO 0.5 STEP 0.001 LET z=x+jy LET kaku0=arg(z-z0) LET kaku1=arg(z-z1) LET kaku2=arg(z-z2) LET kaku=kaku0+kaku1+kaku2 IF kaku<=PI THEN EXIT FOR NEXT x PLOT LINES:x,y; NEXT y PLOT LINES SET LINE COLOR "blue" FOR y=-2.5 TO 0 STEP 0.001 FOR x=-1 TO 0.5 STEP 0.001 LET z=x+jy LET kaku0=arg(z-z0) LET kaku1=arg(z-z1) LET kaku2=arg(z-z2) LET kaku=kaku0+kaku1+kaku2 IF kaku>=-PI THEN EXIT FOR NEXT x PLOT LINES:x,y; NEXT y END -------------------------------- プログラムをrunして得られる（４）の曲線は（２）の曲線に対して、実軸に対して線対称になっていれば正しい結果ですのでこのプログラムで得られる図は（１）（２）（４）は正しい結果を表しています。＜欠点＞（ア）複素数平面上に任意のzを設定する方法として、For-Nextで小刻みな分割計算をして180度または -180度を越えたらForを抜けて線画する、としているので、正確さに欠けるという問題があるのと、作図に時間がかかるのが欠点です。 x方向の計算に挟み撃ち法等を使えば、より正確で計算時間を短く出来ると思うのですが適用していません。y方向については計算短縮の上手い方法が考え付きません。（イ）（３）については同一のx値に対してyの値が２個存在する事になるので作図が出来ない状態で、解答図としては不備になります。以上について何か上手い方法が有りましたらご教示頂けると有り難いです。

Re: 根軌跡図の作成	返事を書くノートメニュー
山中和義 <drdlxujciw> 2007/06/23 20:33:56 ** この記事は1回修正されてます
プログラムから推測して、その処理を検討してみました。 >複素数平面上に任意のzを設定する方法グラフィックス画面のピクセル（ドット）数を対象にすればよいかと思います。実際は、実数x,yについて計算する訳ですので、精度をあげるため２倍の数で行っています。 OPTION ARITHMETIC complex LET Xa=-5 !X=[-5,1] LET Xb=1 LET Ya=-3 !Y=[-3,3] LET Yb=3 SET WINDOW Xa,Xb,Ya,Yb !描画範囲 DRAW grid ASK PIXEL SIZE (Xa,Ya; Xb,Yb) a,b !ドット換算 LET c_Div=2 !分解数 LET dx=(Xb-Xa)/a/c_Div !描画範囲内の描画対象となるドット数を得る LET dy=(Yb-Ya)/b/c_Div LET cEps=1E-3 !誤差　※精度が気になる!? LET j=SQR(-1) !虚数単位 SET POINT STYLE 5 LET z0=0 !極 0 PLOT POINTS: re(z0),im(z0) LET z1=-1 !-1 PLOT POINTS: re(z1),im(z1) LET z2=-4 !-4 PLOT POINTS: re(z2),im(z2) SET POINT STYLE 1 FOR k=0 TO c_Divb-1 !Ｙ座標　※ドットをもとに LET y=Ya+dyk ![Ya,Yb] FOR i=0 TO c_Diva-1 !Ｘ座標 LET x=Xa+dxi ![Xa,Xb] LET z=x+yj !空間上の点を計算する LET a0=GetAngl(z,z0) !極とのなす角を得る LET a1=GetAngl(z,z1) LET a2=GetAngl(z,z2) LET sumA=ABS(a0+a1+a2) !!!PRINT x;ABS(sumA)-PI IF ABS(sumA-PI) < cEps THEN PLOT POINTS: x,y !点画　※(2k+1)PI　k=0,±1,±2,… IF ABS(sumA-3*PI) < cEps THEN PLOT POINTS: x,y NEXT i NEXT k FUNCTION GetAngl(z,z9) !なす角を得る WHEN EXCEPTION IN LET GetAngl=arg(z-z9) USE LET GetAngl=0 !同一点なら END WHEN END FUNCTION END

Re: プログラムから推測して、その処理を検討し...	返事を書くノートメニュー
島村1243 <bjllmpcujp> 2007/06/23 21:52:14
runしたら完璧な作図作図に省略部分の無い完璧な解析プログラム！ですね。この様な作図が出来るプログラムが考えられるなんて凄いですね。（私の書いたプログラムは余りにも単純なのでお恥ずかしい限りです）山中様の作成されたコード内容を直ぐには理解できないレベルですが、どんなテクニックを使用しているのか、久し振りに胸をわくわくしながら読ませて頂きます。ご教示有難う御座いました。

Re: 根軌跡図の作成	返事を書くノートメニュー
SECOND <cszcthjjdj> 2007/06/25 05:05:43 ** この記事は3回修正されてます
!根軌跡図の意味を考える。 !　1+K/s/(s+1)/(s+4)=0　を変形して→　伝達関数F(s)=　K/s/(s+1)/(s+4)　= -1　とすると、 !信号入出力間の、伝達関数F(s)が、-1 に固定される時の、（反転して、閉ループ接続される時の） !反転増幅率　K（実数０～∞）に対する解としての、s(複素数) の軌跡です。 !この軌跡上の s(複素数)は、閉ループ内に、信号電流ｉ(t)= exp^st が、自立出来る事を示します。 !何故そう言えるか、伝達関数F(s) の構成素子としてのインピーダンス Z(s) は、 !定常的な、交流インピーダンス　jωL　や、1/jωC　のような　Z(jω) が、 !電流ｉ(t)= exp^jωt　を想定して成立つ電圧降下率、インピーダンス　Z(jω)　であるように !Z(s) も同様に、過渡的変化を含んでは、いますが、 !電流ｉ(t)= exp^(σ+jω)t　=exp^ｓt　を想定して成立つ電圧降下率、インピーダンス　です。 !そして、インピーダンス　Z(s)で組み立てられた伝達関数F(s) も又、 !信号電流ｉ(t)= exp^ｓt　を想定したときの入出力間、信号電圧比、だからです。 !閉ループ内に、信号電流ｉ(t)= exp^st が、自立出来るという事は、そこに対応している !反転増幅率 K（実数０～∞）において、入力無しの制御系が、勝手に　exp^st　で動く !意味を持ちますので、s=σ+jω の実数部が、 !正であれば、端へぶつかって戻らなくなったり、角速度ωの周波数で、振動増大したりします。 !負であっても、０に近ければ、減衰振動、鈍重さ、等を意味します。 !根軌跡図の　s=σ+jω に対応している K で、ｓの範囲を制限できます。 !s(s+1)(s+4)+K=0　　　s^3 +5s^2 +4s +K=0 !　ＤＫＡ法で直接的に、根の近似値を求める。 ! !　K=0 ～ K=40 まで、根ｓのグラフは、Ｘ印（ｓ= 0,-1,-4) から出発し、対応している。 !------------- OPTION ARITHMETIC COMPLEX LET N=3 DIM A(N),Xr(N) LET A(1)=5 LET A(2)=4 SET WINDOW -6,+1,-3,+3 DRAW grid SET POINT STYLE 5 !Ｘ点 FOR K=0 TO 40 STEP 0.02 LET A(3)=K CALL DKA_00 SET POINT COLOR 2 !青 PLOT POINTS:re(Xr(1)),im(Xr(1)) SET POINT COLOR 1 !黒 PLOT POINTS:re(Xr(2)),im(Xr(2)) SET POINT COLOR 4 !赤 PLOT POINTS:re(Xr(3)),im(Xr(3)) SET POINT STYLE 1 !・点 NEXT K !DKA (Durand Kerner Aberth) 法による高次方程式の根 !-------------------------------------------- ! X^n + A(1)X^(n-1) + ... + A(n) = 0 ! Xr( 1 ~ n ) <== root !-------------------------------------------- SUB DKA_00 LET r=1 FOR j=2 TO N LET rn=ABS(A(j))^(1/j) if r<rn then LET r=rn NEXT j FOR j=1 TO N LET Xr(j)=-A(1)/N+rEXP( complex(0,1)2PI/N (j-3/4) ) NEXT j FOR m=0 TO 100 LET mfx=0 LET maj=0 FOR j=1 TO N LET Xk=1 LET fx=1 FOR w=1 TO N LET fx=fxXr(j)+A(w) IF w<>j THEN LET Xk=Xk*(Xr(j)-Xr(w)) NEXT w LET Xr(j)=Xr(j)-fx/Xk IF mfx<ABS(fx) THEN LET mfx=ABS(fx) IF maj<ABS(fx/Xk) THEN LET maj=ABS(fx/Xk) NEXT j IF mfx<.0000001 AND maj<.0000001 THEN EXIT FOR NEXT m END SUB END

Re: !根軌跡図の意味を考える。	返事を書くノートメニュー
山中和義 <drdlxujciw> 2007/06/25 12:02:22
実係数代数方程式の根（複素数）を複素平面でグラフ化する付属サンプル（SAMPLEフォルダ内）3DPLOT.BASをベースに３Ｄ表示してみました。実係数代数方程式の根（複素数）は、突起部分になります。 f(s)=s(s+1)(s+4) + K　で、K=0～60でアニメーションさせてみると、、、根の軌跡は、このグラフを上から（ＸＹ平面を）見ることになります。突起の頂点が高低するのは、 FOR x=-6 TO 1 STEP 0.05 !複素平面をスキャンする FOR y=-3 TO 3 STEP 0.05 で、うまく頂点を算出できないためです。刻み幅を小さくすれば、該当する可能性が高くなりますが、アニメーションに負荷がかかります。 OPTION ARITHMETIC complex LET i=SQR(-1) !虚数単位 SUB rotx(x,y,z,a) !回転 LET y0=ycos(a)-zsin(a) LET z0=ysin(a)+zcos(a) LET y=y0 LET z=z0 END SUB SUB roty(x,y,z,a) LET x0=xcos(a)+zsin(a) LET z0=-xsin(a)+zcos(a) LET x=x0 LET z=z0 END SUB SUB rotz(x,y,z,a) LET x0=xcos(a)-ysin(a) LET y0=xsin(a)+ycos(a) LET x=x0 LET y=y0 END SUB SUB convert(x,y,z) !投影 CALL rotz(x,y,z,-30PI/180) CALL rotx(x,y,z,-60PI/180) LET x=x+2 LET y=y-1 END SUB SUB plotTo(x,y,z) LET x1=x LET y1=y LET z1=z CALL convert(x1,y1,z1) PLOT LINES: x1,y1; END SUB DEF f(s)=s(s+1)(s+4) + K !実係数代数方程式 !DEF f(s)=s^3+5s^2+4s + K SET WINDOW -4,4,-4,4 !　z　y　で表示する !　│／ !　・─x FOR K=0 TO 40 STEP 0.25 SET DRAW MODE HIDDEN !描画開始 CLEAR FOR x=-6 TO 1 STEP 0.05 !複素平面をスキャンする　x=[-6,1] FOR y=-3 TO 3 STEP 0.05 !y=[-3,3] LET z=ABS(1/f(x+y*i)) !関数の逆数の絶対値を描く CALL PlotTo(x,y,z) !突起部分が解（複素数）になる NEXT y PLOT LINES NEXT x PLOT TEXT,AT -2.5,-2.5: "K="&STR$(K) SET DRAW MODE EXPLICIT !描画終了 NEXT K END

Re: 実係数代数方程式の根（複素数）を複素平面...	返事を書くノートメニュー
SECOND <cszcthjjdj> 2007/06/25 19:53:01
ディジタル走査の、宿命的な問題ですね。私も、刀尽き、矢折れました。グリッドが根に重なる有効な手立ては？ない！その末、先輩たちに、泣きついてしまいました。３Ｄにすると高さ方向に、Ｋも表現できて、すごい図になりそうですね。

Re: ディジタル走査の、宿命的な問題ですね。	返事を書くノートメニュー
SECOND <cszcthjjdj> 2007/06/26 13:04:43 ** この記事は2回修正されてます
!貴兄の３Ｄ部分に習って、書き直してみました。 OPTION ARITHMETIC complex ! s^3 +5s^2 +4s +K =0　の根のグラフ LET N=3 DIM A(N),Xr(N) LET A(1)=5 LET A(2)=4 !A(3)=K !------------------ !　ｚ　ｙ　で表示する !　│／ !　・─ｘ !------------------ SUB rotx(x,y,z,a) LET w=yCOS(a)-zSIN(a) LET z=ySIN(a)+zCOS(a) LET y=w END SUB SUB rotz(x,y,z,a) LET w=xCOS(a)-ySIN(a) LET y=xSIN(a)+yCOS(a) LET x=w END SUB SUB plots(x,y,z) CALL rotz(x,y,z,-PI/6) CALL rotx(x,y,z,-PI/3) PLOT LINES:x,y; END SUB SUB plott(x,y,z,w$) CALL rotz(x,y,z,-PI/6) CALL rotx(x,y,z,-PI/3) PLOT TEXT,AT x,y:w$ END SUB SET WINDOW -6,2,-2.5,5.5 !-----目盛り FOR x=-6 TO 1 IF x=0 THEN SET LINE STYLE 1 ELSE SET LINE STYLE 3 CALL plots(x+0,-3, 0) CALL plots(x+0,+3, 0) !　x+0 は、SUB からｘの書き戻しを止める。 PLOT LINES NEXT x FOR y=-3 TO 3 IF y=0 THEN SET LINE STYLE 1 ELSE SET LINE STYLE 3 CALL plots(-6,y+0, 0) CALL plots( 1,y+0, 0) !　y+0 は、SUB からｙの書き戻しを止める。 PLOT LINES NEXT y SET LINE STYLE 1 !-----Ｋと根のグラフ SET TEXT BACKGROUND "OPAQUE" FOR K=0 TO 40 STEP 0.02 LET A(3)=K CALL DKA_00 FOR s=1 TO 3 SET LINE COLOR s+1 CALL plots(re(Xr(s)),im(Xr(s)), 0) CALL plots(re(Xr(s)),im(Xr(s)),K/10) PLOT LINES NEXT s !----Ｋの印字 CALL plott(-6.2, 0, 0.65+INT(K/10), "K="&USING$("##",K) ) NEXT K !DKA (Durand Kerner Aberth) 法による高次方程式の根 !-------------------------------------------- ! X^n + A(1)X^(n-1) + ... + A(n) = 0 ! Xr( 1 ~ n ) <== root !-------------------------------------------- SUB DKA_00 LET r=1 FOR j=2 TO N LET rn=ABS(A(j))^(1/j) if r<rn then LET r=rn NEXT j FOR j=1 TO N LET Xr(j)=-A(1)/N+rEXP( 2PIcomplex(0,1)/N (j-3/4) ) NEXT j FOR m=0 TO 100 LET mfx=0 LET maj=0 FOR j=1 TO N LET Xk=1 LET fx=1 FOR w=1 TO N LET fx=fxXr(j)+A(w) IF w<>j THEN LET Xk=Xk*(Xr(j)-Xr(w)) NEXT w LET Xr(j)=Xr(j)-fx/Xk IF mfx<ABS(fx) THEN LET mfx=ABS(fx) IF maj<ABS(fx/Xk) THEN LET maj=ABS(fx/Xk) NEXT j IF mfx<.0000001 AND maj<.0000001 THEN EXIT FOR NEXT m END SUB END

Re: !貴兄の３Ｄ部分に習って、書き直してみまし...	返事を書くノートメニュー
SECOND <cszcthjjdj> 2007/06/27 23:10:08
変数引数について。プログラムで使用している方法、これは、ＢＡＳＩＣに負担がかかって、良くないです。　（　　）　CALL plots(x+0,+3, 0) !　x+0 は、SUB からｘの書き戻しを止める。　（　　）ヘルプ(H)　目次(C)　→　手続き定義　→　変数引数　引数が変更されないようにしたいときは，　引数を括弧で括るなどして変数として解釈されないようにする。　例 130 CALL exchange((a),(b)) この方が正しく、ＢＡＳＩＣ側に演算負荷がかかりません。　（　　）　CALL plots((x),+3, 0) !　(x) は、SUB からｘの書き戻しを止める。　（　　）

Re: 実係数代数方程式の根（複素数）を複素平面...	返事を書くノートメニュー
山中和義 <drdlxujciw> 2007/06/27 12:05:59 ** この記事は2回修正されてます
３Ｄで見る代数方程式の解（複素数）通常は、ＸＹ座標のグラフ（本プログラムでは左下から右上方向の）を描いて、Ｘ軸との交点が解　云々・・・「奥ゆき」をＸ虚部とした３Ｄグラフにすると、、、 X=-1、X=0からの２つ解は、重根になって、一組の複素数（共役）へと変化します。まるで、衝突して（実数世界から）、飛び散る（複素数世界へ）みたいです。

Re: ３Ｄで見る代数方程式の解（複素数）	返事を書くノートメニュー
山中和義 <drdlxujciw> 2007/06/27 12:10:25
プログラムはこちら OPTION ARITHMETIC complex LET i=SQR(-1) !虚数単位 LET cEps=1e-6 !誤差 LET N=3 !代数方程式　s(s+1)(s+4) + K = 0 DIM A(N) LET A(1)=5 !多項式 X^N+A(1)X^(N-1)+A(2)X^(N-2)+ … +A(N-1)X+A(N)の係数 LET A(2)=4 !!LET A(3)=0 DIM Xr(N) SUB DKA_00 !DKA法（Durand Kerner Aberth）による代数方程式の解 LET r=1 !初期値を仮定する FOR j=2 TO N LET rn=ABS(A(j))^(1/j) if r<rn then LET r=rn NEXT j FOR j=1 TO N !半径rの円に等間隔に配置する LET Xr(j)=-A(1)/N+rEXP( 2PIi/N (j-3/4) ) !アーバスの初期値 NEXT j FOR m=0 TO 100 !反復 LET mfx=0 LET maj=0 FOR j=1 TO N LET Xk=1 LET fx=1 FOR w=1 TO N LET fx=fxXr(j)+A(w) IF w<>j THEN LET Xk=Xk(Xr(j)-Xr(w)) NEXT w LET Xr(j)=Xr(j)-fx/Xk IF mfx<ABS(fx) THEN LET mfx=ABS(fx) IF maj<ABS(fx/Xk) THEN LET maj=ABS(fx/Xk) NEXT j IF mfx<cEps AND maj<cEps THEN EXIT FOR !収束したら NEXT m END SUB SUB rotx(x,y,z,a) !回転 LET y0=ycos(a)-zsin(a) LET z0=ysin(a)+zcos(a) LET y=y0 LET z=z0 END SUB SUB roty(x,y,z,a) LET x0=xcos(a)+zsin(a) LET z0=-xsin(a)+zcos(a) LET x=x0 LET z=z0 END SUB SUB rotz(x,y,z,a) LET x0=xcos(a)-ysin(a) LET y0=xsin(a)+ycos(a) LET x=x0 LET y=y0 END SUB SUB convert(x,y,z) !投影 CALL rotz(x,y,z,-30PI/180) CALL rotx(x,y,z,-60PI/180) LET x=x+2 LET y=y-2 END SUB SUB plotTo(x,y,z) LET x1=x LET y1=y LET z1=z CALL convert(x1,y1,z1) PLOT LINES: x1,y1; END SUB SUB PlotText(x,y,z,s$) CALL convert(x,y,z) PLOT TEXT ,AT x,y: s$ END SUB SET WINDOW -4,4,-4,4 FOR K=0 TO 40 STEP 0.2 SET DRAW MODE HIDDEN !描画開始 CLEAR !　y　x虚部　で表示する !　│／ !　・─x実部 LET A(3)=K !定数項 FOR x=-6 TO 1 !目盛り CALL PlotText(x+0,0,0,STR$(x)) NEXT x CALL PlotTo(-6,0,0) CALL PlotTo(1,0,0) PLOT LINES FOR y=-3 TO 3 !y=[-3,3] CALL PlotText(0,y+0,0,STR$(y)) NEXT y CALL PlotTo(0,-3,0) CALL PlotTo(0,3,0) PLOT LINES FOR x=-6 TO 1 STEP 0.05 !関数を描く　x=[-6,1] LET y=1 !Horner法（多項式） FOR j=1 TO N LET y=yx+A(j) NEXT j CALL PlotTo(x,0,y/3) !※高さ調整 NEXT x PLOT LINES CALL DKA_00 !根を描く FOR j=1 TO N LET x1=re(Xr(j)) LET y1=im(Xr(j)) LET z1=0 CALL convert(x1,y1,z1) DRAW disk WITH SCALE(0.1)*SHIFT(x1,y1) NEXT j PLOT TEXT,AT -2.5,-2.5: "K="&USING$("###.#",K) SET DRAW MODE EXPLICIT !描画終了 WAIT DELAY 0.1 NEXT K END

Re: プログラムはこちら	返事を書くノートメニュー
山中和義 <drdlxujciw> 2007/06/27 15:54:01 ** この記事は1回修正されてます
もう１つ追加しておきます。３つの解は、まるで　左側が凹にめくれあがった、右側が凸に覆い被さったような曲面を転げ落ちるみたいです。 OPTION ARITHMETIC complex LET i=SQR(-1) !虚数単位 LET cEps=1e-6 !誤差 LET N=3 !代数方程式　s(s+1)(s+4) + K = 0 DIM A(N) !多項式 X^N+A(1)X^(N-1)+A(2)X^(N-2)+ … +A(N-1)X+A(N) LET A(1)=5 !係数 LET A(2)=4 !!LET A(3)=0 　　　：　　　：サブルーチン部分は同じなので、省略します。　　　：　　　： SET WINDOW -4,4,-4,4 FOR K=0 TO 40 STEP 0.2 SET DRAW MODE HIDDEN !描画開始 CLEAR LET A(3)=K !定数項 !　y　x虚部　で表示する !　│／ !　・─x実部 FOR x=-6 TO 1 !目盛り CALL PlotText(x+0,0,0,STR$(x)) NEXT x CALL PlotTo(-6,0,0) !軸 CALL PlotTo(1,0,0) PLOT LINES FOR y=-3 TO 3 !y=[-3,3] CALL PlotText(0,y+0,0,STR$(y)) NEXT y CALL PlotTo(0,-3,0) CALL PlotTo(0,3,0) PLOT LINES FOR x=-6 TO 1 STEP 0.2 !複素平面をスキャンする　x=[-6,1] FOR y=-3 TO 3 STEP 0.2 !y=[-3,3] LET Z=1 !Horner法（多項式） FOR j=1 TO N LET Z=Z(x+yi)+A(j) NEXT j CALL PlotTo(x,y,re(Z/10)) !※高さを調整する NEXT y PLOT LINES NEXT x CALL DKA_00 !解を描く FOR j=1 TO N LET x1=re(Xr(j)) LET y1=im(Xr(j)) LET z1=0 CALL convert(x1,y1,z1) DRAW disk WITH SCALE(0.1)SHIFT(x1,y1) NEXT j PLOT TEXT,AT -2.5,-2.5: "K="&USING$("###.#",K) SET DRAW MODE EXPLICIT !描画終了 WAIT DELAY 0.1 NEXT K END

Re: !根軌跡図の意味を考える。	返事を書くノートメニュー
SECOND <cszcthjjdj> 2007/06/28 09:09:32 ** この記事は1回修正されてます
補足、ＤＫＡの構造。（例は、５次の根を求める場合。）　　　　　　　　真のｆ(x)＝ (x-真1) (x-真2) (x-真3) (x-真4) (x-真5)　　　真の根1~5 　　　　　　　　近似ｆ(x)＝ (x-近1) (x-近2) (x-近3) (x-近4) (x-近5)　　　近似根1~5 　とするとき、近似ｆ(x)　を、真のｆ(x)　に近づけるために　プログラムは、以下の操作を繰返している。　ニュートン法に似ては、いるが、違う。　　　　　　　　　　　　　　(近1-真1) (近1-真2) (近1-真3) (近1-真4) (近1-真5) 次の近1　＝　近1　－　──────────────────────────── 　　　　　　　　　　　　　　　　　　　　(近1-近2) (近1-近3) (近1-近4) (近1-近5) 　　　　　　　　　　　　　　(近2-真1) (近2-真2) (近2-真3) (近2-真4) (近2-真5) 次の近2　＝　近2　－　──────────────────────────── 　　　　　　　　　　　　　　　(近2-近1)　　　　　(近2-近3) (近2-近4) (近2-近5) 　　　　　　　　　　　　　　(近3-真1) (近3-真2) (近3-真3) (近3-真4) (近3-真5) 次の近3　＝　近3　－　──────────────────────────── 　　　　　　　　　　　　　　　(近3-近1) (近3-近2)　　　　　(近3-近4) (近3-近5) 　　　　　　　　　　　　　　(近4-真1) (近4-真2) (近4-真3) (近4-真4) (近4-真5) 次の近4　＝　近4　－　──────────────────────────── 　　　　　　　　　　　　　　　(近4-近1) (近4-近2) (近4-近3)　　　　　(近4-近5) 　　　　　　　　　　　　　　(近5-真1) (近5-真2) (近5-真3) (近5-真4) (近5-真5) 次の近5　＝　近5　－　──────────────────────────── 　　　　　　　　　　　　　　　(近5-近1) (近5-近2) (近5-近3) (近5-近4) この部分が、上記で、ｊ＝１～５について、繰返される。 ! (N=5) LET Xk=1 LET fx=1 FOR w=1 TO N 　LET fx=fxXr(j)+A(w)　　　　　・・・・・・・分子側、真のf(近ｊ)を、計算　IF w<>j THEN LET Xk=Xk(Xr(j)-Xr(w))　　・・分母側、近似f(近ｊ)から、(近ｊ-近ｊ)を除く。 NEXT w LET Xr(j)=Xr(j)-fx/Xk　　　　　　・・・次の近ｊｊ＝１～５について、真のf( 近似根ｊ)≒０になったら出来上がり。プログラムでは、回毎の差分の検査まで、併用しているが、省略可。以下は、最初にセットされる、近似根1~5　の初期値。複素平面( -A(1)/5, 0 )を中心とし、半径r= max( A(1), A(2)^(1/2), …, A(5)^(1/5) ) の円周上に、均等角で、配置されています。かなり、ラフでも、構わない事になっていますが、なぜ、この初期値が、好ましいかは、Durand Kerner Aberth　で検索、御調べ下さい。 ! (N=5) LET r=1 FOR j=2 TO N 　LET rn=ABS(A(j))^(1/j) 　if r<rn then LET r=rn NEXT j FOR j=1 TO N 　LET Xr(j)=-A(1)/N+rEXP( complex(0,1)2PI/N (j-3/4) ) NEXT j

Re: !根軌跡図の意味を考える。	返事を書くノートメニュー
島村1243 <bjllmpcujp> 2007/06/29 23:26:33 ** この記事は1回修正されてます
高次方程式の複素数根を求めるのが本来の姿だった。山中さん、SECONDさん、ご教示有り難うございます。根軌跡図の作成は、一般に高次特性方程式の根を直接求めるのは難しいので、根がどの様に複素平面上を動くのかを幾何学計算（位相角の算術計）から知る唯一の簡便法、と思っていました。そのため、特性方程式の根を計算で求めてそれを複素数平面上にプロットするのが根軌跡作成の本筋、と言うのをすっかり見失っていましたが、山中さんとSECONDさんの解説を見て目が覚めた思いです。山中さんとSECONDさんの解説で、高次方程式の解法にDKA法と言うのがある事を初めて知りました。そして１０進BASICを使うと複素数根も容易に求められると言う事も知ることが出来ました。 Kと根軌跡の関係が一目で分かるSECONDさんの３D作画結果は、図を媒体として根軌跡の意味が良く理解できると言う感想を持ちました。山中さんが作成された最後のプログラムが描く図は、複素数の世界が持っている幻想的な何かを垣間見せている、と言う印象を持ちました。高度な知識を必要とするプログラムを見せて頂き、大変勉強になりました。有難う御座いました。

Re: 高次方程式の複素数根を求めるのが本来の姿...	返事を書くノートメニュー
SECOND <cszcthjjdj> 2007/07/01 01:04:50 ** この記事は1回修正されてます
これからも、掲示を頂ければ、よろしく御願いします。人のあら探しばかりして、自分では、何一つ、自立する掲示のできない時間が、ずっと続いています。恐縮いたします。