コーシーの積分公式

　　コーシーの積分公式 クスクス 2006/12/08 20:50:27 
　　├!コーシーの積分公式:----複素正接函数tan(... クスクス 2006/12/08 20:55:01 
　　│└つづき クスクス 2006/12/08 20:56:24 
　　│　└注：複素正弦関数、余弦関数の定義について... クスクス 2006/12/09 00:40:55 
　　│　　├!難しいので、この状態から始めていいですか... SECOND 2006/12/09 03:19:17 
　　│　　│└関数の定義周りがモジュールを使って不必要... クスクス 2006/12/09 13:15:27 
　　│　　└!殆んど違いが有りませんでした。 SECOND 2006/12/10 02:52:50 
　　│　　　└実際の計算ではほとんど差がないようですね... クスクス 2006/12/10 16:50:28 
　　│　　　　└話が違って恐縮ですが、４次のルンゲ・クッ... SECOND 2006/12/10 17:24:51 
　　│　　　　　└ルンゲ・クッタ法で1階の常微分方程式を解く... クスクス 2006/12/13 01:49:26 
　　│　　　　　　└全くその通りになりました。h=0.1とかやって... SECOND 2006/12/13 04:59:38 
　　└tan(z)のコーシーの積分公式の検証２。先の... クスクス 2006/12/09 14:26:50 
　　　└プログラムコメント クスクス 2006/12/09 14:29:32 
　　　　└プログラム本体。 クスクス 2006/12/09 14:30:19

　　コーシーの積分公式 クスクス 2006/12/08 20:50:27  ツリーへ

コーシーの積分公式 返事を書く

クスクス 2006/12/08 20:50:27

これは「関数引数」で作ったサンプルのおまけです。十進BASICでは
複素数の加減剰余が自由にできます。せっかくの機能ですので、
単なる実数値の区分求積ではなく、複素積分を複素スチルチェス
積分の定義により求めるサンプルを作って見ました。複素正接
関数 tan(z) の原点中心、半径1.5の円周に沿うコーシー積分の
値を、-8 から 8 までの実軸上、実部と虚部に分けてプロット
します。

　　├!コーシーの積分公式:----複素正接函数tan(... クスクス 2006/12/08 20:55:01  ツリーへ

Re: コーシーの積分公式 返事を書く

クスクス 2006/12/08 20:55:01

! コーシーの積分公式: ---- 複素正接函数 tan(z) のコーシー積分のプロット ----
!
! ・tan(z) の半径1.5, 原点中心の円周 C 上におけるコーシー積分
! を I(z) とする。Re(I(z)) と Im(I(z))-2 のグラフを区間[-8,8]で描く。
! 複素函数論は、円周 C の内部において tan(z) がコーシー積分に
! 一致し、外で0になることを保障する。プログラムにより、コーシー
! 積分が[-1.5,1.5] 外の実軸上零になること、(-1.5, 1.5) において
! tan(x) に一致することが確認できる。
!
! ・複素線積分は複素スチルチェス積分として表される。STIELTJES_SUM(A, B, M)は
! 複素スチルチェス積分の定義により、コーシー積分を数値的に計算する。
! 本質的にこれは区分求積である。
!
! ・tan(z) は -π/2, π/2 を特異点に持つから、原点中心の積分路 C の半径を
! π/2 より大きく取ることはできない。実際プログラム中の LET DEF_FN.R = 1.5を
! LET DEF_FN.R = 1.6 等として動かしてみれば納得できる。

!-------------- 主プログラム --------------
OPTION ARITHMETIC COMPLEX

DECLARE EXTERNAL NUMERIC DEF_FN.A, DEF_FN.R
DECLARE EXTERNAL FUNCTION DEF_FN.CAUCHY, DEF_FN.C

LET X1 = -8
LET X2 = 8
LET Y1 = -8
LET Y2 = 8
LET DS = 0.02 ! 円弧サイズ
LET DEF_FN.R = 1.5 ! 積分路半径
LET M = CEIL(2*DEF_FN.R*PI/DS) ! 半径 DEF_FN.R の円周を長さ DS の円弧 M 個に分割する

SET WINDOW X1,X2, Y1,Y2
DRAW GRID

PLOT TEXT, AT X1+(X2-X1)*0.01,Y2-(Y2-Y1)*0.05: "線: 実正接函数 tan(x) のグラフ,"
PLOT TEXT, AT X1+(X2-X1)*0.01,Y2-(Y2-Y1)*0.1: "I(z): tan(z) の0中心, 半径" & STR$(DEF_FN.R) & "の円周に沿うコーシー積分,"
PLOT TEXT, AT X1+(X2-X1)*0.01,Y2-(Y2-Y1)*0.15: "赤○: Re(I(z)),"
PLOT TEXT, AT X1+(X2-X1)*0.01,Y2-(Y2-Y1)*0.2: "青○: Im(I(z))-2."
SET POINT STYLE 4
FOR T=-8 TO 8 STEP 0.1
LET DEF_FN.A = T
LET Z = STIELTJES_SUM(0, 2*PI, M)
SET POINT COLOR 4
PLOT POINTS:T,RE(Z)
SET POINT COLOR 2
PLOT POINTS:T,IM(Z)-2
NEXT T

SET LINE COLOR 4
FOR T=-PI/2 TO PI/2 STEP 0.1
PLOT LINES: t,tan(t);
NEXT T

END
!---------------- 主プログラム終わり ----------------

　　│└つづき クスクス 2006/12/08 20:56:24  ツリーへ

Re: !コーシーの積分公式:----複素正接函数tan(... 返事を書く

クスクス 2006/12/08 20:56:24

つづき

! 函数定義モジュール
MODULE DEF_FN
MODULE OPTION ARITHMETIC COMPLEX
SHARE NUMERIC I_
PUBLIC NUMERIC A,R ! A: コーシー積分の変数, R: 積分路半径
PUBLIC FUNCTION SINE,COSINE,F,CAUCHY,C

LET R = 1.5 ! 積分路 C の半径
LET I_ = COMPLEX(0,1) ! 虚数単位

! 複素正弦函数 sin(z)
EXTERNAL FUNCTION SINE(Z)
LET SINE = (EXP(I_*Z)-EXP(-I_*Z))/(2*I_)
END FUNCTION

! 複素余弦関数 cos(z)
EXTERNAL FUNCTION COSINE(Z)
LET COSINE = (EXP(I_*Z)+EXP(-I_*Z))/2
END FUNCTION

! 目標函数 (tan(z))
EXTERNAL FUNCTION F(Z)
LET F = SINE(Z)/COSINE(Z)
END FUNCTION

! F(Z)×コーシーの積分核/(2πi) = (1/2πi)sin(z)/(z-a)
EXTERNAL FUNCTION CAUCHY(Z)
WHEN EXCEPTION IN
LET CAUCHY = (1/(2*PI*I_))*DEF_FN.F(Z)/(Z-A)
USE
END WHEN
END FUNCTION

! 積分路 c(t) (半径 Rの原点中心の円周)
EXTERNAL FUNCTION C(T)
LET C = R*EXP(I_*T)
END FUNCTION
END MODULE

! DEF_FN.CAUCHY(DEF_FUN.C(T)) のスティルチェス和の計算
! ・T0: 始時刻, T1: 終時刻, M: 分点数
! ・DEF_FN.FUNC: 被積分函数(モジュール DEF_FN で定義する)
EXTERNAL FUNCTION STIELTJES_SUM(T0, T1, M)
OPTION ARITHMETIC COMPLEX

LET H = (T1-T0)/M
LET T = T0+H/2
LET SUM = 0
FOR I=0 TO M-1
LET SUM = SUM + DEF_FN.CAUCHY(DEF_FN.C(T))*(DEF_FN.C(T+H)-DEF_FN.C(T))
LET T = T + H
NEXT I
LET STIELTJES_SUM = SUM
END FUNCTION

　　│　└注：複素正弦関数、余弦関数の定義について... クスクス 2006/12/09 00:40:55  ツリーへ

Re: つづき 返事を書く

クスクス 2006/12/09 00:40:55

注：複素正弦関数、余弦関数の定義について。指数関数を使って
定義してしまいましたが、十進BASICでは双曲線関数が組み込みで
ありますから、本当は公式

cos(x+i*y)=cos(x)cosh(y)-i*sin(x)sinh(y)
sin(x+i*y)=sin(x)cosh(y)+i*cos(x)sinh(y)

を使って桁落ちを防ぐべきです。

　　│　　├!難しいので、この状態から始めていいですか... SECOND 2006/12/09 03:19:17  ツリーへ

Re: 注：複素正弦関数、余弦関数の定義について... 返事を書く

SECOND 2006/12/09 03:19:17

!難しいので、この状態から始めていいですか。

! コーシーの積分公式:　複素正接函数 tan(z) のコーシー積分のプロット
!
! ・tan(z) の半径1.5, 原点中心の円周 C 上におけるコーシー積分
! を I(z) とする。Re(I(z)) と Im(I(z))-2 のグラフを区間[-8,8]で描く。
! 複素函数論は、円周 C の内部において tan(z) がコーシー積分に
! 一致し、外で0になることを保障する。プログラムにより、コーシー
! 積分が[-1.5,1.5] 外の実軸上零になること、(-1.5, 1.5) において
! tan(x) に一致することが確認できる。
!
! ・複素線積分は複素スチルチェス積分として表される。STIELTJES_SUM(A, B, M)は
! 複素スチルチェス積分の定義により、コーシー積分を数値的に計算する。
! 本質的にこれは区分求積である。
!
! ・tan(z) は -π/2, π/2 を特異点に持つから、原点中心の積分路 C の半径を
! π/2 より大きく取ることはできない。実際プログラム中の LET DEF_FN.R = 1.5を
! LET DEF_FN.R = 1.6 等として動かしてみれば納得できる。

!----
OPTION ARITHMETIC COMPLEX
LET X1=-8
LET X2= 8
LET Y1=-8
LET Y2= 8
LET DS= 0.02 !円弧サイズ
LET R = 1.5 !積分路 C 半径
LET M = CEIL(2*R*PI/DS)!半径 DEF_FN.R の円周を長さ DS の円弧 M 個に分割する

SET WINDOW X1,X2, Y1,Y2
DRAW GRID
PLOT TEXT, AT X1+(X2-X1)*0.01,Y2-(Y2-Y1)*0.05: "線: 実正接函数 tan(x) のグラフ,"
PLOT TEXT, AT X1+(X2-X1)*0.01,Y2-(Y2-Y1)*0.1: "I(z): tan(z) の0中心, 半径" & STR$(R) & "の円周に沿うコーシー積分,"
PLOT TEXT, AT X1+(X2-X1)*0.01,Y2-(Y2-Y1)*0.15: "赤○: Re(I(z)),"
PLOT TEXT, AT X1+(X2-X1)*0.01,Y2-(Y2-Y1)*0.2: "青○: Im(I(z))-2."

SET POINT STYLE 4
FOR T= -8 TO 8 STEP 0.1
LET A= T
LET Z= STIELTJES_SUM(0, 2*PI, M, A, R)
SET POINT COLOR 4
PLOT POINTS:T,RE(Z)
SET POINT COLOR 2
PLOT POINTS:T,IM(Z)-2
NEXT T

SET LINE COLOR 4
FOR T= -PI/2 TO PI/2 STEP 0.1
PLOT LINES: t,tan(t);
NEXT T
END

!---- CAUCHY(C(T)) のスティルチェス和の計算
! ・T0: 始時刻, T1: 終時刻, M: 分点数
! A: コーシー積分の変数, R: 積分路C半径

EXTERNAL FUNCTION STIELTJES_SUM(T0, T1, M, A, R)
OPTION ARITHMETIC COMPLEX
LET I_= COMPLEX(0,1)
LET H = (T1-T0)/M
LET T = T0+H/2
LET SUM = 0
FOR I =0 TO M-1
LET SUM= SUM + CAUCHY(C(T))*(C(T+H)-C(T))
LET T= T + H
NEXT I
LET STIELTJES_SUM= SUM

!DEF SINE(Z)= (EXP(I_*Z)-EXP(-I_*Z))/(2*I_)!複素正弦函数 sin(z)
!DEF COSINE(Z)= (EXP(I_*Z)+EXP(-I_*Z))/2 !　複素余弦関数 cos(z)
DEF F(Z)= SINE(Z)/COSINE(Z) !　　　　　　　目標函数 (tan(z))
DEF C(T)= R*EXP(I_*T) !　　　　　　　　　　積分路 c(t) (半径 Rの原点中心の円周)
DEF CAUCHY(Z)= (1/(2*PI*I_))*F(Z)/(Z-A) !　F(Z)×コーシーの積分核/(2πi)= (1/2πi)sin(z)/(z-a)
!----
DEF COSINE(Z)= COS(re(z))*COSH(im(z))-I_*SIN(re(z))*SINH(im(z)) !新cos(z)
DEF SINE(Z)= SIN(re(z))*COSH(im(z))+I_*COS(re(z))*SINH(im(z)) !新sin(z)
END FUNCTION

　　│　　│└関数の定義周りがモジュールを使って不必要... クスクス 2006/12/09 13:15:27  ツリーへ

Re: !難しいので、この状態から始めていいですか... 返事を書く

クスクス 2006/12/09 13:15:27

関数の定義周りがモジュールを使って不必要に難解になっている
のは、関数引数問題を解決しようとして試作したプログラムを使
い回したからです。「被積分関数名は決め打ちにして、関数定義の
本体は別モジュールDEF_FNで行う。DEF_FN を差し替えるだけで、
同じプログラムを使い回せるようにする。」という事だったんですが・・・
本プログラムの趣旨とは何の関係もありませんので、こんな書き方
はすべきではありません。SECOND さんの書き方が正しいです。

　　│　　└!殆んど違いが有りませんでした。 SECOND 2006/12/10 02:52:50  ツリーへ

Re: 注：複素正弦関数、余弦関数の定義について... 返事を書く

SECOND 2006/12/10 02:52:50

!殆んど違いが有りませんでした。
!新sin,cos と、旧sin,cos との比較テストプログラムです。
!オーバーフローまでの耐性は、|Z|=710 と |Z|=711 、
!互いの差の、相対比は、いずれも 1~6( x10^-15) でした。
!sinh,cosh が内部で、exp() を使用していると想像されます。
!元へ戻される方が・・
!----
OPTION ARITHMETIC COMPLEX
LET I_= COMPLEX(0,1)

DEF SINE(Z)= (EXP(I_*Z)-EXP(-I_*Z))/(2*I_)!複素正弦函数 sin(z)
DEF COSINE(Z)= (EXP(I_*Z)+EXP(-I_*Z))/2 !　複素余弦関数 cos(z)

DEF 新SINE(Z)= SIN(re(z))*COSH(im(z))+I_*COS(re(z))*SINH(im(z))
DEF 新COSINE(Z)= COS(re(z))*COSH(im(z))-I_*SIN(re(z))*SINH(im(z))

SET WINDOW -800,1350,-11,9
SET COLOR MIX(15) 0.6,0.6,0.6
DRAW grid(200,1)

LET w1=-18
LET w2=18
LET w3=-8
LET w4=8

SUB prin(w,w$,w1$)
PLOT TEXT,AT SGN(w)*r,ABS(w)-10:w1$&"-overflow"
PLOT TEXT,AT SGN(w)*r,ABS(w)-11:"at |Z|="&STR$(r)
PLOT TEXT,AT SGN(w)*r,ABS(w)-12:w$&w1$
PLOT TEXT,AT SGN(w)*r,ABS(w)-13:"x10^-15"
LET w=0
END SUB

LET r=1
DO
FOR θ=0 TO 2*PI STEP PI/10
LET z=r*EXP(I_*θ)
WHEN EXCEPTION IN
LET a=sine(z)
USE
IF w1<>0 THEN CALL prin(w1,"(sinZ-新sinZ)/","sinZ")
END WHEN
WHEN EXCEPTION IN
LET b=新sine(z)
USE
IF w2<>0 THEN CALL prin(w2,"(sinZ-新sinZ)/","新sinZ")
END WHEN
IF w1<>0 AND w2<>0 THEN
PLOT POINTS:-r,10^15*ABS((a-b)/a)
PLOT POINTS: r,10^15*ABS((a-b)/b)
END IF
WHEN EXCEPTION IN
LET c=cosine(z)
USE
IF w3<>0 THEN CALL prin(w3,"(cosZ-新cosZ)/","cosZ")
END WHEN
WHEN EXCEPTION IN
LET d=新cosine(z)
USE
IF w4<>0 THEN CALL prin(w4,"(cosZ-新cosZ)/","新cosZ")
END WHEN
IF w3<>0 AND w4<>0 THEN
PLOT POINTS:-r,10^15*ABS((c-d)/c)-10
PLOT POINTS: r,10^15*ABS((c-d)/d)-10
END IF
NEXT θ
LET r=r+1
LOOP UNTIL w1=0 AND w2=0 AND w3=0 AND w4=0

END

　　│　　　└実際の計算ではほとんど差がないようですね... クスクス 2006/12/10 16:50:28  ツリーへ

Re: !殆んど違いが有りませんでした。 返事を書く

クスクス 2006/12/10 16:50:28

実際の計算ではほとんど差がないようですね。ただ実引数に対しての
三角関数、双曲線関数の有効桁数がはっきり分かっている場合、
新しい方の定義では、複素引数に対する三角関数の有効桁数が正確に
計算できるメリットがあると思います。

　　│　　　　└話が違って恐縮ですが、４次のルンゲ・クッ... SECOND 2006/12/10 17:24:51  ツリーへ

Re: 実際の計算ではほとんど差がないようですね... 返事を書く

SECOND 2006/12/10 17:24:51

話が違って恐縮ですが、４次のルンゲ・クッタ法について教えて下さい。
このページの「２重振り子のカオス」の中で、
微分方程式を、差分化し、ｔを用いないで使用しています。
これで、正しいでしょうか。

　　│　　　　　└ルンゲ・クッタ法で1階の常微分方程式を解く... クスクス 2006/12/13 01:49:26  ツリーへ

Re: 話が違って恐縮ですが、４次のルンゲ・クッ... 返事を書く

クスクス 2006/12/13 01:49:26

ルンゲ・クッタ法で1階の常微分方程式を解く場合、時間間隔の幅
h は解こうとする方程式に応じて決めねばなりません。1 だから
大きすぎるとか、小さく取れば取るほど近似の精度が上がるとか
ということはありません。重力加速度GをMKS単位系で9.8に取ると
桁あふれを起こすのは、刻みの幅が1に取られているのと関係が
あると思います。刻み幅を固定してしまった場合、収束を改善
するには、Gなどの解こうとする方程式側のパラーメータを調整
せざるおえなくなるでしょう。

RK法の最適刻み幅は試行錯誤で問題ごとに決めるしかありません。小さく取りすぎれば、逆に精度が落ちたりしますし、大きく
取り過ぎれば、計算値がバラけて全く収束しなかったり、解とは
程遠いものを計算したりすることになります。刻み幅を決めなが
ら計算できる算法もあるようですが、それはRK法を越える別の算法です。

　　│　　　　　　└全くその通りになりました。h=0.1とかやって... SECOND 2006/12/13 04:59:38  ツリーへ

Re: ルンゲ・クッタ法で1階の常微分方程式を解く... 返事を書く

SECOND 2006/12/13 04:59:38

全くその通りになりました。h=0.1 とかやって、g=9.8 に
合せ様とすると、BASIC 側の桁落ちの累積誤差が増える為と
想像しますが、破たんと、暴走、オーバーフローに至るまでの
時間が、遥かに短くなりました。h=1 の現状では、乗算が
不要な為か？、３０分は、持ちこたえています。　が・・
やはり、オーバーフローは、時間の問題でやってきます。
これを防ぐには、どうするのでしょうか。

数学的解は、カオスだとしても、現実の計算は、桁落ちのために
周期を持っていると、思うのですが？

　　└tan(z)のコーシーの積分公式の検証２。先の... クスクス 2006/12/09 14:26:50  ツリーへ

Re: コーシーの積分公式 返事を書く

クスクス 2006/12/09 14:26:50

tan(z) のコーシーの積分公式の検証２。先のプログラムは公式
を検証すると言っても、実軸じょうを調べただけでした。これでは
不十分です。そこで、複素平面上の正方形[-L,L]×[-L,L]内の点z
をランダムに抜き出して、z において公式が精度ε内で成立するか
否かをプロットするプログラムを作りました。パラメータは変数
にしてあるので変更は容易です。SECOND さんのアドバイスにより
関数の定義を変えたので、多少は読みやすくなったと思います。

画面上の○印は公式の成立を肯定、×印は否定する情報です。

　　　└プログラムコメント クスクス 2006/12/09 14:29:32  ツリーへ

　　　　└プログラム本体。 クスクス 2006/12/09 14:30:19  ツリーへ

Re: プログラムコメント 返事を書く

クスクス 2006/12/09 14:30:19

プログラム本体。

OPTION ARITHMETIC COMPLEX

LET I_ = COMPLEX(0,1) ! 虚数単位
DEF SINE(Z) = SIN(RE(Z))*COSH(IM(Z)) + I_*COS(RE(Z))*SINH(IM(Z)) ! 複素正弦関数 sin(z)
DEF COSINE(Z) = COS(RE(Z))*COSH(IM(Z)) - I_*SIN(RE(Z))*SINH(IM(Z)) ! 複素余弦関数 cos(z)
DEF TANG(Z) = SINE(Z)/COSINE(Z) ! 複素正接関数 tan(z)
DEF CAUCHY(Z, A) = (1/(2*PI*I_))*TANG(Z)/(Z-A) ! F(Z)×コーシーの積分核/(2πi) = (1/2πi)sin(z)/(z-a)
DEF C(T) = R*EXP(I_*T) ! 積分路 = 半径R, 原点中心の円周
DEF I(A,M) = STIELTJES_SUM(0, 2*PI, M) ! tan(z) の C に沿うコーシー積分を分割数 M で近似したもの

LET DS = 0.01 ! 円弧サイズ
LET R = 1.55 ! 積分路半径
LET M = CEIL(2*R*PI/DS) ! 半径 R の円周を長さ DS の円弧 M 個に分割する
LET L = 5 ! 複素平面上の正方形 [-L,L]×[-L,L]。この領域をプログラムが走査する。
LET EPSILON = 0.006 ! 誤差の上限
LET N = 10000 ! サンプル数

LET X1 = -L
LET X2 = L
LET Y1 = -L
LET Y2 = L
SET WINDOW X1,X2, Y1,Y2
DRAW GRID

SET LINE WIDTH 2 ! 積分路を黄色の円で描く
SET LINE COLOR "YELLOW"
DRAW CIRCLE WITH SCALE(R)
SET LINE COLOR "BLACK"
DRAW CIRCLE WITH SCALE((X2-X1)/300)*SHIFT(-PI/2,0) ! tan(z) の特異点を黒丸で描く
DRAW CIRCLE WITH SCALE((X2-X1)/300)*SHIFT(PI/2,0)
SET TEXT JUSTIFY "right", "bottom"
PLOT TEXT, AT -PI/2,0: "-π/2"
SET TEXT JUSTIFY "left", "bottom"
PLOT TEXT, AT PI/2,0: "π/2"

SET POINT COLOR 4
SET POINT STYLE 4
RANDOMIZE
DO WHILE N>0
LET X = X1+(X2-X1)*RND
LET Y = Y1+(Y2-Y1)*RND
LET A = complex(X,Y)
LET Z = STIELTJES_SUM(0, 2*PI, M)
IF ABS(A)>R THEN
IF ABS(Z)<EPSILON THEN
SET POINT COLOR 2
SET POINT STYLE 4
ELSE
SET POINT COLOR 1
SET POINT STYLE 5
END IF
PLOT POINTS: X,Y
ELSEIF ABS(A)<R THEN
IF ABS(TANG(A)-Z)<EPSILON THEN
SET POINT COLOR 4
SET POINT STYLE 4
ELSE
SET POINT COLOR 1
SET POINT STYLE 5
END IF
PLOT POINTS:X,Y
END IF
LET N=N-1
LOOP

! CAUCHY(C(T)) のスティルチェス和の計算
! ・T0: 始時刻, T1: 終時刻, M: 分点数
FUNCTION STIELTJES_SUM(T0, T1, M)

LET H = (T1-T0)/M
LET T = T0+H/2
LET SUM = 0
FOR J=0 TO M-1
LET SUM = SUM + CAUCHY(C(T),A)*(C(T+H)-C(T))
LET T = T + H
NEXT J
LET STIELTJES_SUM = SUM
END FUNCTION

END

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

コーシーの積分公式	返事を書く
クスクス 2006/12/08 20:50:27
これは「関数引数」で作ったサンプルのおまけです。十進BASICでは複素数の加減剰余が自由にできます。せっかくの機能ですので、単なる実数値の区分求積ではなく、複素積分を複素スチルチェス積分の定義により求めるサンプルを作って見ました。複素正接関数 tan(z) の原点中心、半径1.5の円周に沿うコーシー積分の値を、-8 から 8 までの実軸上、実部と虚部に分けてプロットします。

Re: コーシーの積分公式	返事を書く
クスクス 2006/12/08 20:55:01
! コーシーの積分公式: ---- 複素正接函数 tan(z) のコーシー積分のプロット ---- ! ! ・tan(z) の半径1.5, 原点中心の円周 C 上におけるコーシー積分 ! を I(z) とする。Re(I(z)) と Im(I(z))-2 のグラフを区間[-8,8]で描く。 ! 複素函数論は、円周 C の内部において tan(z) がコーシー積分に ! 一致し、外で0になることを保障する。プログラムにより、コーシー ! 積分が[-1.5,1.5] 外の実軸上零になること、(-1.5, 1.5) において ! tan(x) に一致することが確認できる。 ! ! ・複素線積分は複素スチルチェス積分として表される。STIELTJES_SUM(A, B, M)は ! 複素スチルチェス積分の定義により、コーシー積分を数値的に計算する。 ! 本質的にこれは区分求積である。 ! ! ・tan(z) は -π/2, π/2 を特異点に持つから、原点中心の積分路 C の半径を ! π/2 より大きく取ることはできない。実際プログラム中の LET DEF_FN.R = 1.5を ! LET DEF_FN.R = 1.6 等として動かしてみれば納得できる。 !-------------- 主プログラム -------------- OPTION ARITHMETIC COMPLEX DECLARE EXTERNAL NUMERIC DEF_FN.A, DEF_FN.R DECLARE EXTERNAL FUNCTION DEF_FN.CAUCHY, DEF_FN.C LET X1 = -8 LET X2 = 8 LET Y1 = -8 LET Y2 = 8 LET DS = 0.02 ! 円弧サイズ LET DEF_FN.R = 1.5 ! 積分路半径 LET M = CEIL(2DEF_FN.RPI/DS) ! 半径 DEF_FN.R の円周を長さ DS の円弧 M 個に分割する SET WINDOW X1,X2, Y1,Y2 DRAW GRID PLOT TEXT, AT X1+(X2-X1)0.01,Y2-(Y2-Y1)0.05: "線: 実正接函数 tan(x) のグラフ," PLOT TEXT, AT X1+(X2-X1)0.01,Y2-(Y2-Y1)0.1: "I(z): tan(z) の0中心, 半径" & STR$(DEF_FN.R) & "の円周に沿うコーシー積分," PLOT TEXT, AT X1+(X2-X1)0.01,Y2-(Y2-Y1)0.15: "赤○: Re(I(z))," PLOT TEXT, AT X1+(X2-X1)0.01,Y2-(Y2-Y1)0.2: "青○: Im(I(z))-2." SET POINT STYLE 4 FOR T=-8 TO 8 STEP 0.1 LET DEF_FN.A = T LET Z = STIELTJES_SUM(0, 2*PI, M) SET POINT COLOR 4 PLOT POINTS:T,RE(Z) SET POINT COLOR 2 PLOT POINTS:T,IM(Z)-2 NEXT T SET LINE COLOR 4 FOR T=-PI/2 TO PI/2 STEP 0.1 PLOT LINES: t,tan(t); NEXT T END !---------------- 主プログラム終わり ----------------

Re: !コーシーの積分公式:----複素正接函数tan(...	返事を書く
クスクス 2006/12/08 20:56:24
つづき ! 函数定義モジュール MODULE DEF_FN MODULE OPTION ARITHMETIC COMPLEX SHARE NUMERIC I_ PUBLIC NUMERIC A,R ! A: コーシー積分の変数, R: 積分路半径 PUBLIC FUNCTION SINE,COSINE,F,CAUCHY,C LET R = 1.5 ! 積分路 C の半径 LET I_ = COMPLEX(0,1) ! 虚数単位 ! 複素正弦函数 sin(z) EXTERNAL FUNCTION SINE(Z) LET SINE = (EXP(I_Z)-EXP(-I_Z))/(2I_) END FUNCTION ! 複素余弦関数 cos(z) EXTERNAL FUNCTION COSINE(Z) LET COSINE = (EXP(I_Z)+EXP(-I_Z))/2 END FUNCTION ! 目標函数 (tan(z)) EXTERNAL FUNCTION F(Z) LET F = SINE(Z)/COSINE(Z) END FUNCTION ! F(Z)×コーシーの積分核/(2πi) = (1/2πi)sin(z)/(z-a) EXTERNAL FUNCTION CAUCHY(Z) WHEN EXCEPTION IN LET CAUCHY = (1/(2PII_))DEF_FN.F(Z)/(Z-A) USE END WHEN END FUNCTION ! 積分路 c(t) (半径 Rの原点中心の円周) EXTERNAL FUNCTION C(T) LET C = REXP(I_T) END FUNCTION END MODULE ! DEF_FN.CAUCHY(DEF_FUN.C(T)) のスティルチェス和の計算 ! ・T0: 始時刻, T1: 終時刻, M: 分点数 ! ・DEF_FN.FUNC: 被積分函数(モジュール DEF_FN で定義する) EXTERNAL FUNCTION STIELTJES_SUM(T0, T1, M) OPTION ARITHMETIC COMPLEX LET H = (T1-T0)/M LET T = T0+H/2 LET SUM = 0 FOR I=0 TO M-1 LET SUM = SUM + DEF_FN.CAUCHY(DEF_FN.C(T))*(DEF_FN.C(T+H)-DEF_FN.C(T)) LET T = T + H NEXT I LET STIELTJES_SUM = SUM END FUNCTION

Re: つづき	返事を書く
クスクス 2006/12/09 00:40:55
注：複素正弦関数、余弦関数の定義について。指数関数を使って定義してしまいましたが、十進BASICでは双曲線関数が組み込みでありますから、本当は公式 cos(x+iy)=cos(x)cosh(y)-isin(x)sinh(y) sin(x+iy)=sin(x)cosh(y)+icos(x)sinh(y) を使って桁落ちを防ぐべきです。

Re: 注：複素正弦関数、余弦関数の定義について...	返事を書く
SECOND 2006/12/09 03:19:17
!難しいので、この状態から始めていいですか。 ! コーシーの積分公式:　複素正接函数 tan(z) のコーシー積分のプロット ! ! ・tan(z) の半径1.5, 原点中心の円周 C 上におけるコーシー積分 ! を I(z) とする。Re(I(z)) と Im(I(z))-2 のグラフを区間[-8,8]で描く。 ! 複素函数論は、円周 C の内部において tan(z) がコーシー積分に ! 一致し、外で0になることを保障する。プログラムにより、コーシー ! 積分が[-1.5,1.5] 外の実軸上零になること、(-1.5, 1.5) において ! tan(x) に一致することが確認できる。 ! ! ・複素線積分は複素スチルチェス積分として表される。STIELTJES_SUM(A, B, M)は ! 複素スチルチェス積分の定義により、コーシー積分を数値的に計算する。 ! 本質的にこれは区分求積である。 ! ! ・tan(z) は -π/2, π/2 を特異点に持つから、原点中心の積分路 C の半径を ! π/2 より大きく取ることはできない。実際プログラム中の LET DEF_FN.R = 1.5を ! LET DEF_FN.R = 1.6 等として動かしてみれば納得できる。 !---- OPTION ARITHMETIC COMPLEX LET X1=-8 LET X2= 8 LET Y1=-8 LET Y2= 8 LET DS= 0.02 !円弧サイズ LET R = 1.5 !積分路 C 半径 LET M = CEIL(2RPI/DS)!半径 DEF_FN.R の円周を長さ DS の円弧 M 個に分割する SET WINDOW X1,X2, Y1,Y2 DRAW GRID PLOT TEXT, AT X1+(X2-X1)0.01,Y2-(Y2-Y1)0.05: "線: 実正接函数 tan(x) のグラフ," PLOT TEXT, AT X1+(X2-X1)0.01,Y2-(Y2-Y1)0.1: "I(z): tan(z) の0中心, 半径" & STR$(R) & "の円周に沿うコーシー積分," PLOT TEXT, AT X1+(X2-X1)0.01,Y2-(Y2-Y1)0.15: "赤○: Re(I(z))," PLOT TEXT, AT X1+(X2-X1)0.01,Y2-(Y2-Y1)0.2: "青○: Im(I(z))-2." SET POINT STYLE 4 FOR T= -8 TO 8 STEP 0.1 LET A= T LET Z= STIELTJES_SUM(0, 2PI, M, A, R) SET POINT COLOR 4 PLOT POINTS:T,RE(Z) SET POINT COLOR 2 PLOT POINTS:T,IM(Z)-2 NEXT T SET LINE COLOR 4 FOR T= -PI/2 TO PI/2 STEP 0.1 PLOT LINES: t,tan(t); NEXT T END !---- CAUCHY(C(T)) のスティルチェス和の計算 ! ・T0: 始時刻, T1: 終時刻, M: 分点数 ! A: コーシー積分の変数, R: 積分路C半径 EXTERNAL FUNCTION STIELTJES_SUM(T0, T1, M, A, R) OPTION ARITHMETIC COMPLEX LET I_= COMPLEX(0,1) LET H = (T1-T0)/M LET T = T0+H/2 LET SUM = 0 FOR I =0 TO M-1 LET SUM= SUM + CAUCHY(C(T))(C(T+H)-C(T)) LET T= T + H NEXT I LET STIELTJES_SUM= SUM !DEF SINE(Z)= (EXP(I_Z)-EXP(-I_Z))/(2I_)!複素正弦函数 sin(z) !DEF COSINE(Z)= (EXP(I_Z)+EXP(-I_Z))/2 !　複素余弦関数 cos(z) DEF F(Z)= SINE(Z)/COSINE(Z) !　　　　　　　目標函数 (tan(z)) DEF C(T)= REXP(I_T) !　　　　　　　　　　積分路 c(t) (半径 Rの原点中心の円周) DEF CAUCHY(Z)= (1/(2PII_))F(Z)/(Z-A) !　F(Z)×コーシーの積分核/(2πi)= (1/2πi)sin(z)/(z-a) !---- DEF COSINE(Z)= COS(re(z))COSH(im(z))-I_SIN(re(z))SINH(im(z)) !新cos(z) DEF SINE(Z)= SIN(re(z))COSH(im(z))+I_COS(re(z))SINH(im(z)) !新sin(z) END FUNCTION

Re: !難しいので、この状態から始めていいですか...	返事を書く
クスクス 2006/12/09 13:15:27
関数の定義周りがモジュールを使って不必要に難解になっているのは、関数引数問題を解決しようとして試作したプログラムを使い回したからです。「被積分関数名は決め打ちにして、関数定義の本体は別モジュールDEF_FNで行う。DEF_FN を差し替えるだけで、同じプログラムを使い回せるようにする。」という事だったんですが・・・本プログラムの趣旨とは何の関係もありませんので、こんな書き方はすべきではありません。SECOND さんの書き方が正しいです。

Re: 注：複素正弦関数、余弦関数の定義について...	返事を書く
SECOND 2006/12/10 02:52:50
!殆んど違いが有りませんでした。 !新sin,cos と、旧sin,cos との比較テストプログラムです。 !オーバーフローまでの耐性は、\|Z\|=710 と \|Z\|=711 、 !互いの差の、相対比は、いずれも 1~6( x10^-15) でした。 !sinh,cosh が内部で、exp() を使用していると想像されます。 !元へ戻される方が・・ !---- OPTION ARITHMETIC COMPLEX LET I_= COMPLEX(0,1) DEF SINE(Z)= (EXP(I_Z)-EXP(-I_Z))/(2I_)!複素正弦函数 sin(z) DEF COSINE(Z)= (EXP(I_Z)+EXP(-I_Z))/2 !　複素余弦関数 cos(z) DEF 新SINE(Z)= SIN(re(z))COSH(im(z))+I_COS(re(z))SINH(im(z)) DEF 新COSINE(Z)= COS(re(z))COSH(im(z))-I_SIN(re(z))SINH(im(z)) SET WINDOW -800,1350,-11,9 SET COLOR MIX(15) 0.6,0.6,0.6 DRAW grid(200,1) LET w1=-18 LET w2=18 LET w3=-8 LET w4=8 SUB prin(w,w$,w1$) PLOT TEXT,AT SGN(w)r,ABS(w)-10:w1$&"-overflow" PLOT TEXT,AT SGN(w)r,ABS(w)-11:"at \|Z\|="&STR$(r) PLOT TEXT,AT SGN(w)r,ABS(w)-12:w$&w1$ PLOT TEXT,AT SGN(w)r,ABS(w)-13:"x10^-15" LET w=0 END SUB LET r=1 DO FOR θ=0 TO 2PI STEP PI/10 LET z=rEXP(I_θ) WHEN EXCEPTION IN LET a=sine(z) USE IF w1<>0 THEN CALL prin(w1,"(sinZ-新sinZ)/","sinZ") END WHEN WHEN EXCEPTION IN LET b=新sine(z) USE IF w2<>0 THEN CALL prin(w2,"(sinZ-新sinZ)/","新sinZ") END WHEN IF w1<>0 AND w2<>0 THEN PLOT POINTS:-r,10^15ABS((a-b)/a) PLOT POINTS: r,10^15ABS((a-b)/b) END IF WHEN EXCEPTION IN LET c=cosine(z) USE IF w3<>0 THEN CALL prin(w3,"(cosZ-新cosZ)/","cosZ") END WHEN WHEN EXCEPTION IN LET d=新cosine(z) USE IF w4<>0 THEN CALL prin(w4,"(cosZ-新cosZ)/","新cosZ") END WHEN IF w3<>0 AND w4<>0 THEN PLOT POINTS:-r,10^15ABS((c-d)/c)-10 PLOT POINTS: r,10^15ABS((c-d)/d)-10 END IF NEXT θ LET r=r+1 LOOP UNTIL w1=0 AND w2=0 AND w3=0 AND w4=0 END

Re: !殆んど違いが有りませんでした。	返事を書く
クスクス 2006/12/10 16:50:28
実際の計算ではほとんど差がないようですね。ただ実引数に対しての三角関数、双曲線関数の有効桁数がはっきり分かっている場合、新しい方の定義では、複素引数に対する三角関数の有効桁数が正確に計算できるメリットがあると思います。

Re: 実際の計算ではほとんど差がないようですね...	返事を書く
SECOND 2006/12/10 17:24:51
話が違って恐縮ですが、４次のルンゲ・クッタ法について教えて下さい。このページの「２重振り子のカオス」の中で、微分方程式を、差分化し、ｔを用いないで使用しています。これで、正しいでしょうか。

Re: 話が違って恐縮ですが、４次のルンゲ・クッ...	返事を書く
クスクス 2006/12/13 01:49:26
ルンゲ・クッタ法で1階の常微分方程式を解く場合、時間間隔の幅 h は解こうとする方程式に応じて決めねばなりません。1 だから大きすぎるとか、小さく取れば取るほど近似の精度が上がるとかということはありません。重力加速度GをMKS単位系で9.8に取ると桁あふれを起こすのは、刻みの幅が1に取られているのと関係があると思います。刻み幅を固定してしまった場合、収束を改善するには、Gなどの解こうとする方程式側のパラーメータを調整せざるおえなくなるでしょう。 RK法の最適刻み幅は試行錯誤で問題ごとに決めるしかありません。小さく取りすぎれば、逆に精度が落ちたりしますし、大きく取り過ぎれば、計算値がバラけて全く収束しなかったり、解とは程遠いものを計算したりすることになります。刻み幅を決めながら計算できる算法もあるようですが、それはRK法を越える別の算法です。

Re: tan(z)のコーシーの積分公式の検証２。先の...	返事を書く
クスクス 2006/12/09 14:29:32
プログラムコメント ! コーシーの積分公式２: ---- 複素正接函数 tan(z) のコーシーの積分公式の検証 ---- ! ! ・tan(z) の半径R, 原点中心の円周 C 上におけるコーシー積分 ! を I(a) とする。tan(z) のコーシーの積分公式の I(a) について ! 誤差 ε>0 内での成立を正方形 [-L,L]×[-L,L] 内の点をランダムに ! 抜き出して検査する。 ! ! ・黒丸は tan(z) の特異点 -π/2,π/2。黄色の円は積分路。 ! \|z\|>R のとき ! a) \|I(z)\| < ε なら z の位置に青○ ! b) \|I(z)\| >= ε ならの位置× ! \|z\|<R のとき ! a) \|tan(z)-I(z)\| < ε なら z の位置に赤○ ! b) \|tan(z)-I(z)\| >= ε なら　z の位置に× ! \|z\|=R のとき ! 円周上の点は除外集合。なにもプロットしない。 ! ! ・要するに積分公式の成立を否定するネガティブ情報が×印でプロットされる。 ! ! ・特異点-π/2, π/2 の周りに×が密集するのは、I(a) の近似の精度が ! そこで大きく下がるから。

Re: ルンゲ・クッタ法で1階の常微分方程式を解く...	返事を書く
SECOND 2006/12/13 04:59:38
全くその通りになりました。h=0.1 とかやって、g=9.8 に合せ様とすると、BASIC 側の桁落ちの累積誤差が増える為と想像しますが、破たんと、暴走、オーバーフローに至るまでの時間が、遥かに短くなりました。h=1 の現状では、乗算が不要な為か？、３０分は、持ちこたえています。　が・・やはり、オーバーフローは、時間の問題でやってきます。これを防ぐには、どうするのでしょうか。数学的解は、カオスだとしても、現実の計算は、桁落ちのために周期を持っていると、思うのですが？

Re: コーシーの積分公式	返事を書く
クスクス 2006/12/09 14:26:50
tan(z) のコーシーの積分公式の検証２。先のプログラムは公式を検証すると言っても、実軸じょうを調べただけでした。これでは不十分です。そこで、複素平面上の正方形[-L,L]×[-L,L]内の点z をランダムに抜き出して、z において公式が精度ε内で成立するか否かをプロットするプログラムを作りました。パラメータは変数にしてあるので変更は容易です。SECOND さんのアドバイスにより関数の定義を変えたので、多少は読みやすくなったと思います。画面上の○印は公式の成立を肯定、×印は否定する情報です。

Re: プログラムコメント	返事を書く
クスクス 2006/12/09 14:30:19
プログラム本体。 OPTION ARITHMETIC COMPLEX LET I_ = COMPLEX(0,1) ! 虚数単位 DEF SINE(Z) = SIN(RE(Z))COSH(IM(Z)) + I_COS(RE(Z))SINH(IM(Z)) ! 複素正弦関数 sin(z) DEF COSINE(Z) = COS(RE(Z))COSH(IM(Z)) - I_SIN(RE(Z))SINH(IM(Z)) ! 複素余弦関数 cos(z) DEF TANG(Z) = SINE(Z)/COSINE(Z) ! 複素正接関数 tan(z) DEF CAUCHY(Z, A) = (1/(2PII_))TANG(Z)/(Z-A) ! F(Z)×コーシーの積分核/(2πi) = (1/2πi)sin(z)/(z-a) DEF C(T) = REXP(I_T) ! 積分路 = 半径R, 原点中心の円周 DEF I(A,M) = STIELTJES_SUM(0, 2PI, M) ! tan(z) の C に沿うコーシー積分を分割数 M で近似したもの LET DS = 0.01 ! 円弧サイズ LET R = 1.55 ! 積分路半径 LET M = CEIL(2RPI/DS) ! 半径 R の円周を長さ DS の円弧 M 個に分割する LET L = 5 ! 複素平面上の正方形 [-L,L]×[-L,L]。この領域をプログラムが走査する。 LET EPSILON = 0.006 ! 誤差の上限 LET N = 10000 ! サンプル数 LET X1 = -L LET X2 = L LET Y1 = -L LET Y2 = L SET WINDOW X1,X2, Y1,Y2 DRAW GRID SET LINE WIDTH 2 ! 積分路を黄色の円で描く SET LINE COLOR "YELLOW" DRAW CIRCLE WITH SCALE(R) SET LINE COLOR "BLACK" DRAW CIRCLE WITH SCALE((X2-X1)/300)SHIFT(-PI/2,0) ! tan(z) の特異点を黒丸で描く DRAW CIRCLE WITH SCALE((X2-X1)/300)SHIFT(PI/2,0) SET TEXT JUSTIFY "right", "bottom" PLOT TEXT, AT -PI/2,0: "-π/2" SET TEXT JUSTIFY "left", "bottom" PLOT TEXT, AT PI/2,0: "π/2" SET POINT COLOR 4 SET POINT STYLE 4 RANDOMIZE DO WHILE N>0 LET X = X1+(X2-X1)RND LET Y = Y1+(Y2-Y1)RND LET A = complex(X,Y) LET Z = STIELTJES_SUM(0, 2PI, M) IF ABS(A)>R THEN IF ABS(Z)<EPSILON THEN SET POINT COLOR 2 SET POINT STYLE 4 ELSE SET POINT COLOR 1 SET POINT STYLE 5 END IF PLOT POINTS: X,Y ELSEIF ABS(A)<R THEN IF ABS(TANG(A)-Z)<EPSILON THEN SET POINT COLOR 4 SET POINT STYLE 4 ELSE SET POINT COLOR 1 SET POINT STYLE 5 END IF PLOT POINTS:X,Y END IF LET N=N-1 LOOP ! CAUCHY(C(T)) のスティルチェス和の計算 ! ・T0: 始時刻, T1: 終時刻, M: 分点数 FUNCTION STIELTJES_SUM(T0, T1, M) LET H = (T1-T0)/M LET T = T0+H/2 LET SUM = 0 FOR J=0 TO M-1 LET SUM = SUM + CAUCHY(C(T),A)(C(T+H)-C(T)) LET T = T + H NEXT J LET STIELTJES_SUM = SUM END FUNCTION END