Scipyで常微分方程式を解く (2) 最適化

前回の応用として、

を解き、


を満たすようなを求める。
真の解は

より


以下がスクリプト。

[1] a が与えられた時に、(x,y)=(0,1)から常微分方程式を積分し、x=x1でのyの値を返す関数y_at_x(a,x1)を定義。
[2] y_at_x(a=1,x1=3)を真の解と比較。
[3] F(a,x1,y1)=y_at_x(a,x1)-y1 という関数を定義。　
[4]　F(a,x1=3,y1=4)=0　を解いてaを求める。aはscipyのbrentq関数を使って0から1の範囲で探す。最後に真の解を使って答え合わせ。

import numpy as np
from scipy import integrate
from scipy import optimize 

#solve dy/dx=a y
# x=0, y=1  
# x=3, y=4
# what is a?

#define derivative
def dy_dx(y,x,a):
    return a*y

#[1]
def y_at_x(a,x1):
    y1, infodict = integrate.odeint(dy_dx, 1., [0,x1],args=(a,),full_output=True)
    #print infodict["message"]
    return y1


# [2] test y_at_x(a,x1)
y_test= y_at_x(a=1.,x1=3.)
print "test y_at_ax(1,3)",y_test[1],np.exp(1.*3.)

# [3] define y_at_x(a,x1)-y1
def F(a,x1,y1):
    return y_at_x(a,x1)[1][0]-y1


# [4] Solve  y_at_x(a,3)-4=0
asolution=optimize.brentq(F,0.,1.,args=(3.,4.))
print "solution of a=",asolution
print "true solution", np.log(4)/3.
print "test solution1",y_at_x(a=asolution,x1=3.)[1][0], np.exp(asolution*3.)

出力は

test y_at_ax(1,3) [ 20.08553873] 20.0855369232
solution of a= 0.462098103279
true solution 0.462098120373
test solution1 4.0 3.99999979487