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

数値ラプラス逆変換

　　数値ラプラス逆変換 島村1243 2008/07/22 07:33:55  (修正2回)
　　├（２）について 山中和義 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 
　　├できるだけ演算回数を減らしてみました。 荒田浩二 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

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