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

Scipyで常微分方程式を解く

Scipyのodeint関数を使い、常微分方程式を解いてみる。

簡単な例として、

を、初期値

で解く。
真の解は

である。

以下がpythonスクリプト。xの0から3まで積分している。数値計算で得られた値を＋で、真の解より得られる値を実線で図示している。

import numpy as np
from scipy import integrate
import matplotlib.pyplot as plt 

#solve dy/dx=y

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

#integration
x1=np.linspace(0,3, 10)
x=np.linspace(0,3, 100)
y1, infodict = integrate.odeint(dy_dx, 1., x1, full_output=True)
print infodict

# plot
plt.plot(x1,y1,'+',markersize=12)  #solusion 
plt.plot(x,np.exp(x))            #true solution 
plt.savefig("ode_ex.png")
plt.show()

参照:
LoktaVolterraTutorial
http://www.scipy.org/LoktaVolterraTutorial 
A Coupled Spring-Mass System
http://www.scipy.org/Cookbook/CoupledSpringMassSystem

bloggerで数式を書くテスト

こちらの記事を参照した。
LaTeX for BloggerをSafariと新しい編集インターフェイスのために改造する。
http://satomacoto.blogspot.com/2009/11/latex-for-bloggersafariblogger.html 

以下のような数式が書ける。

PyNGLとPyNIOをインストール

PyNGL(1.31)とPyNIO(1.4.0)をUbuntu10.04にインストールした。

PyNGLは NCAR Command Language (NCL)をもとにしたpythonベースのグラフィックモジュールで、PyNIOはいろいろな形式のファイルを読み込むモジュールである。

以下がインストール手順。

(1)　プレコンパイル済みのファイルを入手
PyNGL-1.3.1.linux-debian-i686-gcc432-py255-numpy130.tar.gz　と　PyNIO-1.4.0.linux-debian-i686-gcc432-py255-numpy130.tar.gz　をEarth System Gridより入手。
Pythonのバージョンは2.6.5を使っているが、numpyは1.3.0を使っているので、プレコンパイル済みのバイナリの中からこれらを選択した。

(2)　以下を実行する。

sudo tar -C /usr/local -xzf PyNIO-1.4.0.linux-debian-i686-gcc432-py255-numpy130.tar.gz
sudo tar -C /usr/local -xzf PyNGL-1.3.1.linux-debian-i686-gcc432-py255-numpy130.tar.gz　

(3)　パスを通す。

export PYTHONPATH=$PYTHONPATH:/usr/local/lib/python2.5/site-packages:/usr/local/lib/python2.5/si
te-packages/PyNIO:/usr/local/lib/python2.5/site-packages/PyNGL
export PYNGL_NCARG=/usr/local/lib/python2.5/site-packages/PyNGL/ncarg

(4) pynglexの書き換え
/usr/local/bin/pynglexの一行目を

#!/usr/bin/env python 

のように編集する。