|
|
!２階微分方程式の、ルンゲクッタ描画。
!ベッセルの微分方程式の場合。
!x^2*(d2y/dx2)+x*(dy/dx)+(x^2-a^2)*y=0 …ベッセル
!x^2*(d2y/dx2)+x*(dy/dx)-(x^2+a^2)*y=0 …変形ベッセル
!-------------------------------
OPTION ARITHMETIC NATIVE
OPTION BASE 0
SET TEXT background "opaque"
DIM col(4), y0(4)
MAT READ col
DATA 4,10,2,14, 12 ! ０次 ～４次の色分け
!-----------------------------------------------------------------------
!0白 1黒 2青 3緑 4赤 5水色 6黄色 7赤紫
! 8灰色 9濃い青 10濃い緑 11青緑 12えび茶 13オリーブ色 14濃い紫 15銀色
!-----------------------------------------------------------------------
SUB Dji( oddy,ody, iddy,idy,y)
LET oddy=-( x*idy +(sj*x^2-a^2)*y)/x^2 ! ddy=(d2y/dx2) :sj=1(ベッセル)
LET ody=-(x^2*iddy+(sj*x^2-a^2)*y)/x ! dy= (dy/dx) :sj=-1(変形ベッセル)
END SUB
SUB RungeKutta
CALL Dji( oddy1,ody1, iddy, idy, y)
CALL Dji( oddy2,ody2, iddy, idy+oddy1*dx/2, y+ody1*dx/2)
CALL Dji( oddy3,ody3, iddy, idy+oddy2*dx/2, y+ody2*dx/2)
CALL Dji( oddy4,ody4, iddy, idy+oddy3*dx , y+ody3*dx )
LET y=y+(ody1+2*ody2+2*ody3+ody4)*dx/6
LET idy=idy+(oddy1+2*oddy2+2*oddy3+oddy4)*dx/6
LET iddy=oddy4
END SUB
LET dx=.0005 !pitch
LET y0(0)=1 !ｙの初期値
LET y0(1)= 0.708248*dx
LET y0(2)=-1.2017 *dx^2
LET y0(3)=-0.015 *dx^3
LET y0(4)= 0.00615 *dx^4
!-----
LET sj=-1
LET m$="変形"
CALL window_
CALL Rk_bessel !変形ベッセルの微分方程式( ルンゲクッタ描画 )
CALL fn_bessel !変形ベッセル関数１種 In(x)
WAIT DELAY 2
!-----
LET sj=1
LET m$=""
CALL window_
CALL Rk_bessel ! ベッセルの微分方程式( ルンゲクッタ描画 )
CALL fn_bessel ! ベッセル関数１種 Jn(x)
SUB window_
CLEAR
IF sj=1 THEN !----ベッセル
LET h=1
LET l=-1
LET xr=20
SET WINDOW -3,xr, l, h
DRAW grid(xr/4,h/5)
ELSEIF sj=-1 THEN !----変形ベッセル
LET h=4
LET l=-1.5
LET xr=4
SET WINDOW -0.3,xr, l, h
DRAW grid(xr/8,h/8)
END IF
ASK PIXEL SIZE(0,0;xr,0) px,py
LET ss=Xr/px
END SUB
!-------------------- 微分方程式のルンゲクッタ描画
SUB Rk_bessel
SET TEXT COLOR 4
PLOT TEXT,AT (.325+.225*sj)*xr,.9*h-(.2-.04*sj)*(h-l) :"ルンゲクッタ描画"
FOR a=0 TO 4
LET y=y0(a)
LET idy=0
LET iddy=0
SET LINE COLOR "gray" !col(a+4)
SET TEXT COLOR "gray" !col(a+4)
PLOT TEXT,AT .1*xr,.9*h-.04*(h-l)*a :m$& "ベッセルの微分方程式 "& STR$(a)& "次"
!------------
FOR x=dx TO xr+dx STEP dx
PLOT LINES: x,y; ! PEN-on
IF FP(x)< dx THEN PRINT USING"##.## ###.######":x,y
CALL RungeKutta
NEXT x
!------------
PLOT LINES !PEN-off
PRINT
NEXT a
END SUB
!********************* 解のベッセル関数を重ねて照合する。
SUB fn_bessel
SET TEXT COLOR 1
PLOT TEXT,AT (.325+.225*sj)*xr,.9*h-(.2-.04*sj)*(h-l) :"関数で、重ね書き"
FOR n=0 TO 4
SET LINE COLOR col(n)
SET TEXT COLOR col(n)
PLOT TEXT,AT .1*xr,.9*h-.04*(h-l)*n :m$& "ベッセル関数１種 "& STR$(n)& "次 "
!------------画素間隔のStep
FOR x=dx TO xr+ss STEP ss
PLOT LINES: x,bessel(n,x); ! PEN-on
IF FP(x)< ss AND dx< x THEN PRINT USING"##.## ###.######":IP(x),bessel(n,IP(x))
NEXT x
!------------
PLOT LINES !PEN-off
PRINT
NEXT n
END SUB
!------- ベッセル関数 変形ベッセル関数
FUNCTION bessel(n,x)
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
IF sj=1 THEN !---- ベッセル関数 Jn(x)
FOR k=1 TO m/2
LET w=w+Tk(k*2,x)
NEXT k
LET bessel=Tk(n,x)/(Tk(0,x)+2*w)
ELSEIF sj=-1 THEN !---- 変形ベッセル関数 In(x)
FOR k=1 TO m
LET w=w+Tk(k,x)
NEXT k
LET bessel=EXP(x)*Tk(n,x)/(Tk(0,x)+2*w)
END IF
END FUNCTION
FUNCTION Tk(i,x)
LET t2=0
LET t1=1e-9
LET t0=2*(m+1)/x*t1-sj*t2
FOR kp1=m TO i+1 STEP -1
LET t2=t1
LET t1=t0
LET t0=2*kp1/x*t1-sj*t2
NEXT kp1
LET Tk=t0
END FUNCTION
END
|
|