|
投稿者：SECOND
投稿日：2009年11月 9日(月)01時52分53秒
|
|
|
|
!ベッセル関数ですが、このリストは、ｘの範囲が、実数までで、
!複素数が使用できません。拡張できる方、お願いします。
!
!※十進BASIC にも ベッセル関数、変形ベッセル関数１種を内包出来ないでしょうか。
!-------------------------------
OPTION ARITHMETIC NATIVE
OPTION BASE 0
SET TEXT background "opaque"
DIM col(5)
MAT READ col
DATA 4,10,2, 8,8,8 !Red darkGreen Blue Gray Gray Gray
!-----
CALL window_i
CALL G_besseli !変形ベッセル関数１種 In(x)
WAIT DELAY .5
CALL window_j
CALL G_besselj !ベッセル関数１種 Jn(x)
SUB window_i
CLEAR
LET h=4
LET l=-1.5
LET xr=4
SET WINDOW -.1*xr,xr, l,h
DRAW grid(xr/8,h/8)
ASK PIXEL SIZE (0,0; xr,0) j,i
LET dx=xr/j !pitch
END SUB
SUB window_j
CLEAR
LET h=1
LET l=-1
LET xr=20
SET WINDOW -.1*xr,xr, l,h
DRAW grid(xr/4,h/5)
ASK PIXEL SIZE (0,0; xr,0) j,i
LET dx=xr/j !pitch
END SUB
!-------
SUB G_besseli
FOR n=0 TO 2
SET LINE COLOR col(n)
SET TEXT COLOR col(n)
PLOT TEXT,AT .13*xr,h-.08*(h-l)-.04*(h-l)*n :"変形ベッセル関数１種 "& STR$(n)& "次"
FOR x=dx TO xr+dx STEP dx
LET y=besseli(n,x)
PLOT LINES: x,y; ! PEN-on
IF FP(x)< dx THEN PRINT USING"##.## ###.######":x,y
NEXT x
PLOT LINES !PEN-off
PRINT
NEXT n
END SUB
SUB G_besselj
FOR n=0 TO 2
SET LINE COLOR col(n)
SET TEXT COLOR col(n)
PLOT TEXT,AT .13*xr,h-.08*(h-l)-.04*(h-l)*n :"ベッセル関数１種 "& STR$(n)& "次"
FOR x=dx TO xr+dx STEP dx
LET y=besselj(n,x)
PLOT LINES: x,y; ! PEN-on
IF FP(x)< dx THEN PRINT USING"##.## ###.######":x,y
NEXT x
PLOT LINES !PEN-off
PRINT
NEXT n
END SUB
!-------
FUNCTION besseli(n,x)
LET m=2*INT( (6+MAX(n,1.5*x)+9*1.5*x/(1.5*x+2))/2)
LET w=0
FOR k=1 TO m
LET w=w+Tki(k)
NEXT k
LET besseli=EXP(x)*Tki(n)/(Tki(0)+2*w)
END FUNCTION
FUNCTION Tki(i)
LET t2=0
LET t1=1e-9
LET t0=2*(m+1)/x*t1+t2
FOR kp1=m TO i+1 STEP -1
LET t2=t1
LET t1=t0
LET t0=2*kp1/x*t1+t2
NEXT kp1
LET Tki=t0
END FUNCTION
!-------
FUNCTION besselj(n,x)
LET m=2*INT( (6+MAX(n,1.5*x)+9*1.5*x/(1.5*x+2))/2)
LET w=0
FOR k=1 TO m/2
LET w=w+Tk(k*2)
NEXT k
LET besselj=Tk(n)/(Tk(0)+2*w)
END FUNCTION
FUNCTION Tk(i)
LET t2=0
LET t1=1e-9
LET t0=2*(m+1)/x*t1-t2
FOR kp1=m TO i+1 STEP -1
LET t2=t1
LET t1=t0
LET t0=2*kp1/x*t1-t2
NEXT kp1
LET Tk=t0
END FUNCTION
END
|
|
|
|
投稿者：山中和義
投稿日：2009年11月11日(水)11時19分29秒
|
|
|
|
> No.714[元記事へ]
SECONDさんへのお返事です。
サイトを検索してみましたが、複素変数での実例（説明、サンプルなど）がほとんどありません。
本格的な計算はわかりません、、、
!Jn(z)=(1/π)*∫[0,π]{cos(n*θ-z*sinθ)}dθ 積分表示より、数値積分で求める
!（検算に使った）参考サイト http://keisan.casio.jp/
OPTION ARITHMETIC COMPLEX
LET j=SQR(-1) !虚数単位
DEF SIN(z)=(EXP(j*z)-EXP(-j*z))/(2*j) !三角関数
DEF COS(z)=(EXP(j*z)+EXP(-j*z))/2
FUNCTION complexbessel(n,z) !１種ベッセル関数 Jn(x+i*y)
LET div=1000 !分割数
LET u=0
LET h=PI/div
LET a=0
FOR i=1 TO div !数値積分の台形公式
LET u=u+( COS(n*a-z*SIN(a)) + COS(n*(a+h)-z*SIN(a+h)) )/2
LET a=h*i
NEXT i
LET complexbessel=h*u/PI
END FUNCTION
SET WINDOW -0.1*20,20, -1,1 !表示範囲
DRAW grid(20/4,1/5)
LET dx=0.2 !グラフの描画間隔
FOR n=0 TO 2
SET LINE COLOR n+2
SET TEXT COLOR n+2
PLOT TEXT,AT .13*20,1-.08*2-.04*2*n :"ベッセル関数１種 "& STR$(n)& "次"
FOR x=0 TO 20 STEP dx
LET y=Re( complexbessel(n,COMPLEX(x,0)) ) !実数
PLOT LINES: x,y;
IF FP(x)< dx THEN PRINT USING"##.## ###.######": x,y
NEXT x
PLOT LINES
PRINT
NEXT n
END
|
|
|
|
投稿者：SECOND
投稿日：2009年11月11日(水)11時56分51秒
|
|
|
|
> No.714[元記事へ]
!前回のベッセル関数について、変数 m を、
!絶対値処理するだけで、ｘの範囲を、複素数まで拡張できる事が分った。
!
! Jn( complex(x,0) )= complex(0,1)^(-n) *In( complex(0,x) )
! In( complex(x,0) )= complex(0,1)^(-n) *Jn( complex(0,x) ) !←文はこのテスト状態。
!-------------------------------
OPTION ARITHMETIC COMPLEX
OPTION BASE 0
SET TEXT background "opaque"
DIM col(5)
MAT READ col
DATA 4,10,2, 8,8,8 !Red darkGreen Blue Gray Gray Gray
!-----
CALL window_i
CALL G_besselj !ベッセル関数１種 Jn(x) …虚数入力
WAIT DELAY .5
CALL G_besseli !変形ベッセル関数１種 In(x)
SUB window_i
CLEAR
LET h=4
LET l=-1.5
LET xr=4
SET WINDOW -.1*xr,xr, l,h
DRAW grid(xr/8,h/8)
ASK PIXEL SIZE (0,0; xr,0) j,i
LET dx=xr/j !pitch
END SUB
SUB window_j
CLEAR
LET h=1
LET l=-1
LET xr=20
SET WINDOW -.1*xr,xr, l,h
DRAW grid(xr/4,h/5)
ASK PIXEL SIZE (0,0; xr,0) j,i
LET dx=xr/j !pitch
END SUB
!-------
SUB G_besseli
FOR n=0 TO 2
SET LINE COLOR col(n)
SET TEXT COLOR col(n)
PLOT TEXT,AT .13*xr,h-.08*(h-l)-.04*(h-l)*n :"変形ベッセル関数１種 "& STR$(n)& "次"
FOR t=dx TO xr+dx STEP dx
LET x=t
LET y=besseli(n,x)
! LET x=COMPLEX(0,t)
! LET y=COMPLEX(0,1)^(-n)*besseli(n,x)
PLOT LINES: ABS(x),y; ! PEN-on
IF FP(t)< dx THEN PRINT x;y
NEXT t
PLOT LINES !PEN-off
PRINT
NEXT n
END SUB
SUB G_besselj
FOR n=0 TO 2
SET LINE COLOR col(n+3)
SET TEXT COLOR col(n+3)
PLOT TEXT,AT .13*xr,h-.08*(h-l)-.04*(h-l)*n :"ベッセル関数１種 "& STR$(n)& "次 "
FOR t=dx TO xr+dx STEP dx
! LET x=t
! LET y=besselj(n,x)
LET x=COMPLEX(0,t)
LET y=COMPLEX(0,1)^(-n)*besselj(n,x)
PLOT LINES: ABS(x),y; ! PEN-on
IF FP(t)< dx THEN PRINT x;y
NEXT t
PLOT LINES !PEN-off
PRINT
NEXT n
END SUB
!-------
FUNCTION besseli(n,x)
IF x=0 THEN LET x=1e-32 !０の保護( 場合により外す)
LET m=2*INT( (6+MAX(n,1.5*ABS(x))+9*1.5*ABS(x)/(1.5*ABS(x)+2))/2 )
LET w=0
FOR k=1 TO m
LET w=w+Tki(k,x)
NEXT k
LET besseli=EXP(x)*Tki(n,x)/(Tki(0,x)+2*w)
END FUNCTION
FUNCTION Tki(i,x)
LET t2=0
LET t1=1e-9
LET t0=2*(m+1)/x*t1+t2
FOR kp1=m TO i+1 STEP -1
LET t2=t1
LET t1=t0
LET t0=2*kp1/x*t1+t2
NEXT kp1
LET Tki=t0
END FUNCTION
!-------
FUNCTION besselj(n,x)
IF x=0 THEN LET x=1e-32 !０の保護( 場合により外す)
LET m=2*INT( (6+MAX(n,1.5*ABS(x))+9*1.5*ABS(x)/(1.5*ABS(x)+2))/2 )
LET w=0
FOR k=1 TO m/2
LET w=w+Tk(k*2,x)
NEXT k
LET besselj=Tk(n,x)/(Tk(0,x)+2*w)
END FUNCTION
FUNCTION Tk(i,x)
LET t2=0
LET t1=1e-9
LET t0=2*(m+1)/x*t1-t2
FOR kp1=m TO i+1 STEP -1
LET t2=t1
LET t1=t0
LET t0=2*kp1/x*t1-t2
NEXT kp1
LET Tk=t0
END FUNCTION
END
|
|
|
|
投稿者：SECOND
投稿日：2009年11月11日(水)12時05分34秒
|
|
|
|
> No.717[元記事へ]
山中さんへ
いれちがいになったようで… ありがとうございました。問題は解決していました。ほっとしました
|
|
|
戻る