この微分方程式は、非常に簡単に解析的に解くことができます。
dy/dt = 0.01 * y * (1-y)
並べ替えて反対側の y 項と t 項を集める
100 dt = 1/(y * (1-y)) dy
lhs は自明に 100 * t に統合されますが、rhs はもう少し複雑です。2 つの商 * いくつかの定数の和として、常に 2 つの商の積を書くことができます。
1/(y * (1-y)) = A/y + B/(1-y)
A と B の値は、rh を同じ分母に置き、両側の定数項と 1 次の y 項を比較することで計算できます。この場合は単純で、A=B=1 です。したがって、統合する必要があります
1/y + 1/(1-y) 日
第 1 項は ln(y) に統合され、第 2 項は変数 u = 1-y から -ln(1-y) への変更で統合できます。そのための統合された方程式は次のようになります。
100 * t + C = ln(y) - ln(1-y)
積分定数を忘れないでください (ここで lhs に書くと便利です)。2 つの対数項を組み合わせることができます。
100 * t + C = ln( y / (1-y) )
y の正確な値について t を解くには、まず C の値を計算する必要があります。これは、初期条件を使用して行います。y が 1 から始まる場合、dy/dt = 0 であり、y の値が変化しないことは明らかです。したがって、最初に y と t の値を挿入します
100 * 0 + C = ln( y(0) / (1 - y(0) )
これにより、C の値が得られ (y が 0 または 1 でないと仮定)、y=0.8 を使用して t の値が取得されます。対数と係数 100 により、y の初期値が信じられないほど小さい場合を除き、ty は t 値の比較的短い範囲内で 0.8 に達することに注意してください。もちろん、上記の式を並べ替えて y を t で表現するのも簡単で、関数をプロットすることもできます。
編集:数値積分
解析的に解けないより複雑な ODE については、数値的に試す必要があります。最初は、ゼロ時間 y(0) での関数の値 (関数の軌跡を一意に定義するために少なくともそれを知る必要があります) と勾配の評価方法しか知りません。数値積分の考え方は、勾配に関する知識 (関数がどのように変化しているかを教えてくれる) を使用して、開始点付近で関数の値がどうなるかを計算できるということです。これを行う最も簡単な方法は、オイラー積分です。
y(dt) = y(0) + dy/dt * dt
オイラー積分は、勾配が t=0 と t=dt の間で一定であると仮定します。y(dt) がわかると、そこで勾配を計算し、y(2 * dt) などの計算に使用して、関数の完全な軌跡を徐々に構築することができます。特定の目標値を探している場合は、軌跡がその値を超えるまで待ってから、最後の 2 つの位置の間を補間して正確な t を取得します。
オイラー積分 (および他のすべての数値積分法) の問題は、仮定が有効な場合にのみ結果が正確になることです。勾配は時点のペア間で一定ではないため、積分ステップごとに一定量のエラーが発生し、時間が経つにつれて、答えが完全に不正確になるまで蓄積されます。積分の質を向上させるためには、勾配に対してより高度な近似を使用する必要があります。たとえば、Runge-Kutta 法を確認してください。これは、計算時間の増加を犠牲にして誤差項の累進次数を除去する積分器のファミリーです。関数が微分可能である場合、2 番目または 3 番目の微分を知ることで、積分誤差を減らすこともできます。
幸いなことに、もちろん、他の誰かがここで大変な作業を行っているので、数値安定性などの問題を解決することについてあまり心配する必要はなく、すべての詳細を深く理解する必要はありません (何が起こっているのかを大まかに理解することは非常に役立ちますが) )。すぐに使用できるインテグレーター クラスの例については、http://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.ode.html#scipy.integrate.odeを参照してください。例えば
from scipy.integrate import ode
def deriv(t, y):
return 0.01 * y * (1 - y)
my_integrator = ode(deriv)
my_integrator.set_initial_value(0.5)
t = 0.1 # start with a small value of time
while t < 3000:
y = my_integrator.integrate(t)
if y > 0.8:
print "y(%f) = %f" % (t, y)
break
t += 0.1
このコードは、y が 0.8 を超えたときに最初の t 値を出力します (0.8 に達しない場合は何も出力しません)。t のより正確な値が必要な場合は、前の t の y も保持し、それらの間を補間します。