（１）について
!逆ラプラス変換
! y(t)=1/(2π*i)∫Y(s)*Exp(s*t)ds [γ-i*ω,γ+i*ω]、iは虚数単位、γは実数。
!ωは実数、s=γ+i*2πωとすると、ds=i*2πdω。
!これを上式に代入して、変形すると、
! y(t)=Exp(γt)∫Y(γ,ω)*Exp(i*2π*ω)dω [-∞,∞]
!右辺の積分部分は逆フーリエ変換。
OPTION ARITHMETIC COMPLEX
DECLARE EXTERNAL SUB IFFT !逆フーリエ変換
LET t0=TIME
SET WINDOW -1,11,-1.2,1.2
DRAW grid
!実行したいラプラス変換を設定する
LET w0=2*PI*0.1 ! 0.1Hz波の角周波数
DEF F(s)=s/(s^2+w0^2)*EXP(-s) !開始時間遅れのあるcos(w0*(t-1))のラプラス変換
LET n = 1024 !データ総数
DIM d(0 TO n-1) !入出力用配列
LET T=10 !時間幅 0〜T
LET Ganma=5
LET dt=T/n !時間刻み幅
FOR k=0 TO n/2
LET Fs=F( COMPLEX(Ganma/T,k*2*PI/T) ) !ラプラス変換
LET Hw=(COS(k*2*PI/n)+1)/2 !ハニング関数 ※精度の向上
LET c0=n/T * Hw
LET d(k)=COMPLEX(Re(Fs)*c0,Im(Fs)*c0)
IF k<>0 AND k<>n/2 THEN LET d(n-k)=COMPLEX(Re(d(k)),-Im(d(k))) !複素共役
NEXT k
CALL IFFT(d) !逆フーリエ変換
LET r=Ganma/n
FOR k=0 TO n-1 !結果の表示
PLOT LINES: k*dt, Re(d(k))*EXP(k*r);
NEXT k
PRINT TIME-t0
!厳密解の関数曲線（赤打点）
SET POINT COLOR "red"
FOR t=0 TO 10 STEP 0.1
IF t<1 THEN
LET ft=0
ELSE
LET ft=COS(w0*(t-1))
END IF
PLOT POINTS:t,ft
NEXT t
END
EXTERNAL SUB IFFT(x()) !逆フーリエ変換 x() : 入力/出力データ
OPTION ARITHMETIC COMPLEX
LET nx = SIZE(x)
LET theta = 2 * PI / nx ! W = exp(-j * i*pi/N * -1) = exp(j * theta)とする
!DFTの計算式は、要素番号 K=1..N に対して
! xk = 1/N * Σ [n=0..N-1] (xn * W^(-k*n) )
! W = exp(-j * 2*PI/N)
DIM w(0 TO nx-1), xtmp(0 TO nx-1)
MAT xtmp = x
FOR k=0 TO nx-1
LET tmp = theta * k
FOR n=0 TO nx-1
LET w(n)=EXP( COMPLEX(0, tmp * n) )
NEXT n
LET x(k) = DOT(w, xtmp)
NEXT k
MAT x = (1/nx) * x
END SUB
|
新しく発言する
EXIT
インデックスへ
返事を書く 
ノートメニュー