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

数値ラプラス逆変換

　　数値ラプラス逆変換 島村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回)

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 
　　　├凄い速さで驚きです。 島村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
新規発言を反映させるにはブラウザの更新ボタンを押してください。