4

私はPythonとscipyに比較的慣れておらず、MATLABからの変換です。scipy.integrate の odeint 関数の簡単なテストを行っていたところ、この潜在的なバグに遭遇しました。次のスニペットを検討してください。

from scipy import *
from scipy.integrate import odeint
from scipy.interpolate import interp1d
from pylab import *

# ODE system with forcing function u(t)
def sis(x,t,u):
    return [x[1], u(t)]

# Solution time span
t = linspace(0, 10, 1e3)

# Profile for forcing function u(t)
acel = lambda t: 3*(t<2)-3*(t>8)

# Interpolator for acelerator
acel_interp = interp1d(t, acel(t), bounds_error=False, fill_value=0)    

# ODE integration with u(t) = acel, a lambda function
x_1 = odeint(sis, [0,0], t, args=(acel,) )            # Correct result
# ODE integration with u(t) = acel_interp, an interpolator
x_2 = odeint(sis, [0,0], t, args=(acel_interp,) )     # Incorrect result

両方の結果の違いを示すプロットを作成しました。ここをクリックしてください。

少なくとも私にとっては、不当な結果の違いについてどう思いますか? Python 2.6.6 の上で NumPy バージョン 1.5.0 と SciPy バージョン 0.8.0 を使用しています

4

1 に答える 1

3

これはバグではありません。問題は、bound_errorFalseに切り替えて、それらの値をゼロで埋めたことにあります。元のコードでTrueに設定した場合bound_error、補間の範囲を超えているため、積分にゼロを入れていることがわかります（したがって、関数の外側のポイントで関数を評価した場合とは異なる値が生成されます。のラムダで行うのと同じ範囲x_1。

次のことを試してみてください。正常に動作していることがわかります。基本的にt、補間を使用している範囲をカバーするのに十分な大きさの値の範囲をカバーするように拡張しました。

from scipy import *
from scipy.integrate import odeint
from scipy.interpolate import interp1d
from pylab import *

# ODE system with forcing function u(t)
def sis(x,t,u):
    return [x[1], u(t)]

# Solution time span
t = linspace(0, 10, 1e3)
t_interp = linspace(0,20,2e3)

# Profile for forcing function u(t)
acel = lambda t: 3*(t<2)-3*(t>8)

# Interpolator for acelerator
acel_interp = interp1d(t_interp, acel(t_interp))    

# ODE integration with u(t) = acel, a lambda function
x_1 = odeint(sis, [0,0], t, args=(acel,) )            
# ODE integration with u(t) = acel_interp, an interpolator
x_2 = odeint(sis, [0,0], t, args=(acel_interp,) )
于 2011-04-13T14:14:26.293 に答える