ベッセル関数

 投稿者：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 		 

Re: ベッセル関数

 投稿者：山中和義  投稿日：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
 

Re: ベッセル関数

 投稿者：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 		 

Re: ベッセル関数

 投稿者：SECOND  投稿日：2009年11月11日(水)12時05分34秒
  > No.717[元記事へ]

山中さんへ
いれちがいになったようで… ありがとうございました。問題は解決していました。ほっとしました 		 

戻る