新しく発言する  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 
Re: (1)について  返事を書く  ノートメニュー
島村1243 <bjllmpcujp> 2008/07/25 09:14:17
Ganmaと誤差

ご教示頂いたプログラムを早速走らせてみました。描かれた波形が1.0を僅かですが越えているので、
試しにGanma=10に増やして見ました。すると1.0を越える事が無くなりました。

次にラプラス変換式を
F(s)=1/s*EXP(-s)
に変更して走らせると、t=10付近に至ると、濃い黒線が描かれ激しく発散しているようです。
Ganma=5に戻すと発散はしますが薄い黒線に変わりましたが、1.0を越えた関数波形になります。

以上の試行により、波形の精度はganmaを10程度にすると良いが、単位ステップ関数のラプラス変換
に適用すると最終時刻付近で発散します。

*********************
発散の原因を突き止めようと、最初の説明書きに
! y(t)=Exp(γt)∫Y(γ,ω)*Exp(i*2π*ω)dω [-∞,∞]
と有りましたので、プログラムを見ていくと

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) ) !ラプラス変換 <----ここ

と言う箇所で使用している Ganma/T{Y(γ,ω)のγと推測}と

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

と言うところで使われるrの意味合いが異なっていました。理由は、
! y(t)=Exp(γt)∫Y(γ,ω)*Exp(i*2π*ω)dω [-∞,∞]
に使われているExp(γt)とY(γ,ω)のγは同じ値のはず、と思ったからです。

そこでLET r=Ganma/n を LET r=Ganma/T/n と変更して F(s)=1/s の逆変換を試行したら
f(t)=EXP(-αt) らしき波形に復元されてしまいましたので、私の原因推測は誤っている事が分かりました。

確かに解析的には、ラプラス変換時に使用するγはその値を指定しない(都合の良い正の定数)ので、逆変換時
のγは都合の良い値に変えて構わないはずですが、数値計算の時にもそれが活きて来るのかが理解不能です。
    └精度を上げるのは難しいようです。 山中和義 2008/07/25 14:53:36 
     └良く分かりました 島村1243 2008/07/25 21:02:36 

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