前回の応用として、
を解き、
を満たすような
真の解は
より
以下がスクリプト。
[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
