数値ラプラス逆変換
数値ラプラス逆変換 島村1243 2008/07/22 07:33:55 (修正2回)
├（２）について 山中和義 2008/07/23 21:22:50
│└約3分の１になりました 島村1243 2008/07/24 19:50:08
├できるだけ演算回数を減らしてみました。 荒田浩二 2008/07/23 23:03:12 (修正1回)
│└更に速い 島村1243 2008/07/24 19:54:40
└（１）について 山中和義 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
数値ラプラス逆変換 島村1243 2008/07/22 07:33:55 (修正2回) 
ツリーへ
|
数値ラプラス逆変換
|
 返事を書く  ノートメニュー
|
|
島村1243 <bjllmpcujp> 2008/07/22 07:33:55 ** この記事は2回修正されてます
|
数値ラプラス逆変換は一般にFFTを使って行うとの事ですが、（FFTのコーディングをしたことが無いので）
FFTを使わないやり方で、数値ラプラス逆変換が出来るプログラムを作りました。
（１）これをFFT法で処理すると、処理時間はどの位短くなるか。
（２）このプログラムで、処理時間を短縮する手段はあるか。
を教えて頂けないでしょうか。
作成したプログラムは下記のとおりです。
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 !０から無限（２π）までの分割数設定
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
|
├（２）について 山中和義 2008/07/23 21:22:50 ツリーへ
|
Re: 数値ラプラス逆変換
|
返事を書く ノートメニュー
|
|
山中和義 <drdlxujciw> 2008/07/23 21:22:50
|
（２）について
ftの算出部分が（内側のループ）、各tに対して30000回のループですから、
ここの処理を軽減するのが得策でしょう。
そこで
・係数部分aw,bwは、tに依存しないためtのループの外に出す。
・同様に、ftiの共通項EXP(a*t)*dw/PIをiのループの外に出す。
の改修をしてみました。
ただし、メモリ性能は悪くなります。
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))のラプラス変換
LET t0=TIME
!数値計算解の関数曲線
SET LINE COLOR "black"
LET n=100 !０から無限（２π）までの分割数設定
LET dw=2*PI/n !微小分割角周波数
LET a=0.1
DIM aw(0 TO n*300),bw(0 TO n*300) !<-----
FOR i=0 TO n*300 !<-----
LET w=dw*i !<-----
LET Fs=F(a+j*w) !<-----
LET aw(i)=Re(Fs) !<-----
LET bw(i)=-Im(Fs) !<-----
NEXT i !<-----
FOR t=0 TO 10 STEP 0.05
LET C=EXP(a*t)*dw/PI !<-----
LET ft=0
FOR i=0 TO n*300
LET w=dw*i
LET fti=( aw(i)*COS(w*t) + bw(i)*SIN(w*t) )*C !<-----
IF w=0 THEN LET fti=fti/2
LET ft=ft+fti
NEXT i
PLOT LINES:t,ft;
NEXT t
PLOT LINES
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
|
│└約3分の１になりました 島村1243 2008/07/24 19:50:08 ツリーへ
|
Re: （２）について
|
返事を書く ノートメニュー
|
|
島村1243 <bjllmpcujp> 2008/07/24 19:50:08
|
約3分の１になりました
山中さん、いつも素早いご教示有り難うございます。
私のPCで約３２秒かかっていたのですが、約９秒になりました。
配列を利用すれば良い、と言う発想が出ませんでした。
|
├できるだけ演算回数を減らしてみました。 荒田浩二 2008/07/23 23:03:12 (修正1回) ツリーへ
|
Re: 数値ラプラス逆変換
|
返事を書く ノートメニュー
|
|
荒田浩二 <knrztrhoel> 2008/07/23 23:03:12 ** この記事は1回修正されてます
|
できるだけ演算回数を減らしてみました。
また、FORループの外に出せるものはできるだけ出してみました。
とくに IF w=0 THEN は、実質 i=0 のときだけですからループに入れてはもったいないと・・・。
私のPCでは、24.5秒から6.1秒に改善されました。
OPTION ARITHMETIC COMPLEX
LET time0=TIME
SET WINDOW -1,11,-1.2,1.2
DRAW grid
LET j=SQR(-1)
!実行したいラプラス変換を設定する
LET w0=2*PI*0.1 ! 0.1Hz波の角周波数
LET w0_2=(2*PI*0.1)^2
!DEF F(s)=s/(s^2+w0_2)*EXP(-s) !開始時間遅れのあるcos(w0*(t-1))のラプラス変換
!数値計算解の関数曲線
SET LINE COLOR "black"
LET n=100 !０から無限（２π）までの分割数設定
LET dw=2*PI/n !微小分割角周波数
LET a=0.1
DIM Re_F(n*300),Im_F(n*300)
FOR i=1 TO n*300
LET s=a+j*dw*i
LET Fs=s/(s^2+w0_2)*EXP(-s)
LET Re_F(i)=Re(Fs)
LET Im_F(i)=Im(Fs)
NEXT i
LET ft0=(Re(a/(a^2+w0_2)*EXP(-a))*COS(0)-Im(a/(a^2+w0_2)*EXP(-a)*SIN(0)))/2
FOR t=0 TO 10 STEP 0.05
LET edwp=EXP(a*t)*2/n ! EXP(a*t)*dw/PI
LET ft=ft0*edwp ! i=0のとき
FOR i=1 TO n*300
LET wt=dw*i*t
LET ft=ft+Re_F(i)*COS(wt)*edwp-Im_F(i)*SIN(wt)*edwp
NEXT I
PLOT LINES:t,ft;
NEXT T
PLOT lines
!厳密解の関数曲線（赤打点）
SET POINT COLOR "red"
FOR t=0 TO 1-0.1 STEP 0.1
PLOT POINTS:t,0
NEXT t
FOR t=1 TO 10 STEP 0.1
PLOT POINTS:t,cos(w0*(t-1))
NEXT T
PLOT TEXT ,AT 10,-1 : STR$(TIME-time0)
END
|
│└更に速い 島村1243 2008/07/24 19:54:40 ツリーへ
└（１）について 山中和義 2008/07/24 19:50:54 ツリーへ
|
Re: 数値ラプラス逆変換
|
返事を書く ノートメニュー
|
|
山中和義 <drdlxujciw> 2008/07/24 19:50:54
|
（１）について
!逆ラプラス変換
! 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
|
├凄い速さで驚きです。 島村1243 2008/07/24 20:14:16 ツリーへ
|
Re: （１）について
|
返事を書く ノートメニュー
|
|
島村1243 <bjllmpcujp> 2008/07/24 20:14:16
|
凄い速さで驚きです。
山中さん、手間のかかるコードを作成していただき有り難うございます。
私のPCで約0.９秒でした。
やはりFFTを使う（だけでは無いようですが）と速いのですね。
私のプログラムは、元々は、ラプラス変換の物理的イメージを捉える方法を考えている中で気づいたものです。
具体的には、電気的歪み波の複素フーリェ係数(An-jBn)から、AnとBnを抜き出して時間関数を復元出来る手法
をフーリェ変換に拡張利用し、更にラプラス変換に拡張利用出来れば、ラプラス変換のAw、Bwは、複素フーリェ
係数のAn、Bnと同じような物理的イメージとして捉えることが出来ると考え付いたものです。
十進BASICを知る前から、もしかしたら複素フーリェ係数(An-jBn)と同じ様な考え方でフーリェ変換F(jw)やラ
プラス変換F(s)も扱えるのかな？というぼんやりした物理イメージは有ったのですが、テストする手段が無くて
具体的な式を考えるには至りませんでした。しかし十進BASICを知って、複素数を扱う計算が簡単なので、今回
作成を試したものです。結果的には、数値ラプラス逆変換を行っていると言う事をウェブで知りました。
私のプログラムは物理的なイメージを分かり易く示していると思うのですが、処理精度と速度を真の数値逆変換
プログラムを作る事が出来ないので比べる事が出来ずにいました。
ご教示頂いた、真の数値逆変換プログラムを、これから勉強したいと思います。
お忙しい中、本当に有難うございました。
|
└Ganmaと誤差 島村1243 2008/07/25 09:14:17 ツリーへ
|
Re: （１）について
|
返事を書く ノートメニュー
|
|
島村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 ツリーへ
|
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 ツリーへ
|
Re: 精度を上げるのは難しいようです。
|
返事を書く ノートメニュー
|
|
島村1243 <bjllmpcujp> 2008/07/25 21:02:36
|
良く分かりました
山中さん、γの件、大変良く分かりました。どうも有難うございました。
ハニング関数の意味も、ウェブで調べ、HW=1と仮設定してプログラムを走らせた結果の波形図を見て良く分かり
ました。
私の示したやり方は、フーリェ級数による復元方法を使った非常に単純な方法で、時間がかかる事が大欠点です
が、ハニング関数等を使わなくても発散することが無い（証明した訳では無いので言いきれませんが）ので、つ
まり、無垢の元関数のみで復元できるので、捨てたものでもないのかなと少し思いました。
|
インデックスへ
EXIT
新規発言を反映させるにはブラウザの更新ボタンを押してください。
新しく発言する