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

数値ラプラス逆変換


  数値ラプラス逆変換 島村1243 2008/07/22 07:33:55  (修正2回)
数値ラプラス逆変換  返事を書く  ノートメニュー
島村1243 <bjllmpcujp> 2008/07/22 07:33:55 ** この記事は2回修正されてます
数値ラプラス逆変換は一般にFFTを使って行うとの事ですが、(FFTのコーディングをしたことが無いので)
FFTを使わないやり方で、数値ラプラス逆変換が出来るプログラムを作りました。
(1)これをFFT法で処理すると、処理時間はどの位短くなるか。
(2)このプログラムで、処理時間を短縮する手段はあるか。
を教えて頂けないでしょうか。

作成したプログラムは下記のとおりです。

option arithmetic complex
SET WINDOW -1,11,-1.2,1.2
DRAW grid
LET j=SQR(-1)

!実行したいラプラス変換を設定する
LET w0=2*PI*0.1 ! 0.1Hz波の角周波数
DEF F(s)=s/(s^2+w0^2)*exp(-s) !開始時間遅れのあるcos(w0*(t-1))のラプラス変換

!数値計算解の関数曲線
SET LINE COLOR "black"
LET n=100 !0から無限(2π)までの分割数設定
LET dw=2*PI/n !微小分割角周波数
LET a=0.1
FOR t=0 TO 10 STEP 0.05
LET ft=0
FOR i=0 TO n*300
LET w=dw*i
LET s=a+j*w
LET aw=Re(F(s))
LET bw=-Im(F(s))
LET fti=aw*EXP(a*t)*COS(w*t)*dw/PI+bw*EXP(a*t)*SIN(w*t)*dw/PI
IF w=0 THEN LET fti=fti/2
LET ft=ft+fti
NEXT I
PLOT LINES:t,ft;
NEXT T

PLOT lines

!厳密解の関数曲線(赤打点)
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
  (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 
     └良く分かりました 島村1243 2008/07/25 21:02:36 

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