2

scipy を使用して常微分方程式系を解いています。簡単にするために、私のコードを次のようにします。

import scipy as sp
import numpy as np
from scipy.integrate import odeint
from numpy import array

def deriv(y,t): # return derivatives of the array y
 a = -2.0
 b = -0.1
 return array([ y[1], a*y[0]+b*y[1] ])
time = np.linspace(0.0,10.0,100)
yinit = array([0.0005,0.2]) # initial values
y = odeint(deriv,yinit,time)

しかし今、定数「a」のいくつかの値についてこのシステムを解きたいと思います。したがって、たとえば a = -2.0 だけではなく、次のようにします。

a = array([[-2.0,-1.5,-1.,-0.5]])

a の各値についてシステムを解きます。配列の各要素をループせずにこれを行う方法はありますか? 一度に全部できますか？

4

1 に答える 1

1

の各値に対して 2 つの方程式を書くことで、連立方程式を再定式化できますa。それを行う1つの方法は

from scipy.integrate import odeint
from numpy import array,linspace,tile,empty_like
a = array([-2.0,-1.5,-1.,-0.5])
b = array([-.1,-.1,-.1,-.1])

yinit = tile(array([0.0005,0.2]),a.size)

def deriv(y,_):
    dy = empty_like(y)
    dy[0::2]=y[1::2]
    dy[1::2]=y[::2]*a+b*y[1::2]
    return dy
time = linspace(0.0,10.0,100)
y = odeint(deriv,yinit,time)

y[:,0]for のソリューション、yforのソリューションなどで、a=-2forの導関数であることがわかります。y[:,2]a=-1.5y[:,-1]ya=-.5

とにかく、これをベクトル化して速度を上げたい場合は、Python コードを c に変換してからコンパイルするライブラリに興味があるかもしれません。pydelay時間遅延もシミュレートする必要があるため、個人的に使用しており、お勧めします。C コードを扱う必要さえありません。翻訳は完全に自動化されています。

于 2013-09-25T21:55:32.973 に答える