新しく発言する  EXIT  インデックスへ

数値ラプラス逆変換


  数値ラプラス逆変換 島村1243 2008/07/22 07:33:55  (修正2回)
  (2)について 山中和義 2008/07/23 21:22:50 
  │└約3分の1になりました 島村1243 2008/07/24 19:50:08 
  できるだけ演算回数を減らしてみました。 荒田浩二 2008/07/23 23:03:12  (修正1回)
  │└更に速い 島村1243 2008/07/24 19:54:40 
  (1)について 山中和義 2008/07/24 19:50:54 
   ├凄い速さで驚きです。 島村1243 2008/07/24 20:14:16 
   └Ganmaと誤差 島村1243 2008/07/25 09:14:17 
    └精度を上げるのは難しいようです。 山中和義 2008/07/25 14:53:36 
Re: Ganmaと誤差  返事を書く  ノートメニュー
山中和義 <drdlxujciw> 2008/07/25 14:53:36
精度を上げるのは難しいようです。

>最終時刻付近で発散します

γを大きくすれば精度は上がりそうですが、あとでExp(γt)をかけるので、
tが大きいところで発散します。

ハニング関数で、かなり抑えています。



>そこでLET r=Ganma/n を LET r=Ganma/T/n と変更して

γ=Ganma/T(γt=3〜7程度)と決めていて、
Exp(γt)をかける箇所は、Exp(γ*k*T/n)、k=0〜n-1となります。

r=Ganma/n=Ganma/T*T/n





このように記述した方がわかりやすいと思います。(メイン処理のみ)


(前略)


LET T=10 !時間区間 [0,T]

LET Ganma=5 !※3〜7
!γを大きくすれば精度は上がりそうですが、
!あとでExp(γt)をかけるので、tが大きいところで発散する。

LET r=Ganma/T !γを決める

FOR k=0 TO N/2 !変換データの作成
LET Fs=F( COMPLEX(r, 2*PI*k/T) ) !変換したい関数
LET Hw=(COS(2*PI*k/N)+1)/2 !ハニング関数 ※精度の向上
LET d(k)=N/T*Fs * Hw !0〜N/2

IF NOT(k=0 OR k=N/2) THEN LET d(N-k)=COMPLEX(Re(d(k)),-Im(d(k))) !N/2+1〜N-1は、共役複素数
NEXT k


CALL IFFT(d) !逆フーリエ変換


LET dt=T/N !時間刻み幅

FOR k=0 TO N-1 !結果の表示
LET tt=k*dt !経過時間
PLOT LINES: tt, Re(d(k))*EXP(r*tt); !Exp(γt)をかける ※Exp(γ*k*T/n)、k=0〜n-1
PRINT Re(d(k))*EXP(r*tt) !debug
NEXT k


(後略)
     └良く分かりました 島村1243 2008/07/25 21:02:36 

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