3.6 Python言語によるプログラム例

(ここは工事中。 前項の Julia 言語によるプログラムと同じことをするプログラムを実現する。)

rk1ex1.py

# rk1ex1.py --- dx/dt=x (0<t<1), x(0)=1 をRunge-Kutta法で解く
# JuliaプログラムのPythonバージョン
# rk1ex1.jl(http://nalab.mind.meiji.ac.jp/~mk/labo/text/intro-ode-simulation/node25.html)
# Runge-Kutta法であるが、応用が効きにくい形。参考プログラム。
# シェルから起動したときコマンド・ライン引数を見るあたりは、参考になるかも。

import sys

def testrungekutta(Tmax=1.0, N=10):
    t0=0.0
    x0=1.0
    print('#N=%d Tmax=%f' % (N,Tmax))
    print('# t x')
    #Runge-Kutta法
    dt=(Tmax-t0)/N
    t=t0
    x=x0
    for j in range(N):
        k1=dt*f(t,x)
        k2=dt*f(t+dt/2,x+k1/2)
        k3=dt*f(t+dt/2,x+k2/2)
        k4=dt*f(t+dt,x+k3)
        x += (k1+2*k2+2*k3+k4) / 6.0
        t = t0 + (j+1) * dt
        print('%.4f %.7f' % (t, x))

def f(t,x):
    return x

if __name__ == "__main__":
    argc = len(sys.argv)
    if argc == 1:
        testrungekutta()
    elif argc == 2:
        testrungekutta(float(sys.argv[1]))
    elif argc == 3:
        testrungekutta(float(sys.argv[1]), int(sys.argv[2]))

実行の仕方

% python rk1ex1.py

#N=10 Tmax=1.000000
# t x
0.1000 1.1051708
0.2000 1.2214026
0.3000 1.3498585
0.4000 1.4918242
0.5000 1.6487206
0.6000 1.8221180
0.7000 2.0137516
0.8000 2.2255396
0.9000 2.4596014
1.0000 2.7182797
%

次は van der Pol 方程式の問題のJuliaによるプログラムである。 Tmax, N を実行時に入力できるような工夫を入れてみた。 また計算結果を rk2ex1.data というファイルに出力している。

rk2ex1.py

# rk2ex1.py --- Runge-Kutta method for van der Pol equation
# JuliaプログラムのPythonバージョン
# rk2ex1.jl(http://nalab.mind.meiji.ac.jp/~mk/labo/text/intro-ode-simulation/node25.html)
# お勧めという訳でもないが、それなりに応用が効く。

import sys
import numpy as np

def rungekutta(f,t,x,dt,args=()):
    k1=dt*f(t,x,args=args)
    k2=dt*f(t+dt/2, x+k1/2,args=args)
    k3=dt*f(t+dt/2, x+k2/2,args=args)
    k4=dt*f(t+dt, x+k3,args=args)
    return x + (k1 + 2 * k2 + 2 * k3 + k4) / 6

# van der Pol方程式
def f(t, x, args=(1.0,)):
    # return np.array([x[1],-x[0]+mu*(1.0-x[0]*x[0])*x[1]])
    mu=args[0]
    y=np.empty(np.size(x))
    y[0]=x[1]
    y[1]=-x[0]+mu*(1.0-x[0]*x[0])*x[1]
    return y

def testrungekutta(Tmax=50.0, N=1000):
    mu=1.0
    t0=0.0
    x0=np.array([0.1,0.1])
    print('#N=%d Tmax=%f' % (N,Tmax))
    print('# t x')
    #Runge-Kutta法
    dt=(Tmax-t0)/N
    t=t0
    x=x0
    with open('rk2ex1.data', mode='w') as fout:
        print('%f %f %f' % (t, x[0], x[1]))
        print('%f %f %f' % (t, x[0], x[1]), file=fout)
        for j in range(N):
            x=rungekutta(f,t,x,dt,args=(mu,))
            t = t0 + (j+1) * dt
            print('%f %f %f' % (t, x[0], x[1]))
            print('%f %f %f' % (t, x[0], x[1]), file=fout)

if __name__ == "__main__":
    argc = len(sys.argv)
    if argc == 1:
        testrungekutta()
    elif argc == 2:
        testrungekutta(float(sys.argv[1]))
    elif argc == 3:
        testrungekutta(float(sys.argv[1]), int(sys.argv[2]))

実行の仕方

% python rk2ex1.py
% gnuplot
gnuplot> plot "rk2ex1.data" using 2:3 with l
gnuplot> quit

% python rk2ex1.jl 100 2000
(Tmax を 100, N を 2000 として実行)

gnuplot を使わず、解軌道を直接描くのも簡単である。

rk2ex1plot.py

# rk2ex1.py --- Runge-Kutta method for van der Pol equation
# JuliaプログラムのPythonバージョン
# rk2ex1.jl(http://nalab.mind.meiji.ac.jp/~mk/labo/text/intro-ode-simulation/node25.html)
# お勧めという訳でもないが、それなりに応用が効く。

import sys
import numpy as np
import matplotlib.pyplot as plt

def rungekutta(f,t,x,dt,args=()):
    k1=dt*f(t,x,args=args)
    k2=dt*f(t+dt/2, x+k1/2,args=args)
    k3=dt*f(t+dt/2, x+k2/2,args=args)
    k4=dt*f(t+dt, x+k3,args=args)
    return x + (k1 + 2 * k2 + 2 * k3 + k4) / 6

# van der Pol方程式
def f(t, x, args=(1.0,)):
    # return np.array([x[1],-x[0]+mu*(1.0-x[0]*x[0])*x[1]])
    mu=args[0]
    y=np.empty(np.size(x))
    y[0]=x[1]
    y[1]=-x[0]+mu*(1.0-x[0]*x[0])*x[1]
    return y

def testrungekutta(Tmax=50.0, N=1000):
    mu=1.0
    t0=0.0
    x0=np.array([0.1,0.1])
    print('#N=%d Tmax=%f' % (N,Tmax))
    print('# t x')
    #Runge-Kutta法
    dt=(Tmax-t0)/N
    t=np.linspace(t0,Tmax,N+1)
    x=np.empty((N+1,2))
    x[0]=x0
    print('%f %f %f' % (t[0], x[0,0], x[0,1]))
    for j in range(N):
        x[j+1]=rungekutta(f,t,x[j],dt,args=(mu,))
        print('%f %f %f' % (t[j+1], x[j+1,0], x[j+1,1]))

    plt.plot(x[:,0],x[:,1])
    plt.show()

if __name__ == "__main__":
    argc = len(sys.argv)
    if argc == 1:
        testrungekutta()
    elif argc == 2:
        testrungekutta(float(sys.argv[1]))
    elif argc == 3:
        testrungekutta(float(sys.argv[1]), int(sys.argv[2]))