Pythonでふつうは解けない微分方程式を解けたふうに見てみる

##オイラー法

$\dot{x} =f(x)$のとき、理想的には$x$を時間で積分して、$x=g(t)$となるような関数$g$を見つけたい。しかし非線形系においては解析的な積分はできない。そこで、時間を離散化することで数値的に解くことを考える。$f(x)$を時間で積分したのと近いような数値を、それぞれのtについてもとめることができればよい。そこで、長方形を組み合わせた区分求積法を応用したものであるオイラー法がつかえる。

わずかな時間幅$\Delta{t}$をとって,高さ$f(x_0)$、幅$\Delta{t}$の長方形の面積をもとめれば、その時間幅における積分値(＝求める解$g(t)$）に近い値が得られる。なので、これを使って初期値から$\Delta{t}$だけ進んだところに対応した値$x_1$を考えると、$$x(t_0+\Delta{t})\approx x_0+f(x_0)\Delta{t}=x_1$$とできる。一般化すると、$$x_{n+1}=x_n+f(x_n)\Delta{t}$$で、これを繰り返すことによって$g(t)$に近いものを求めることができる。

しかし、これは近似値によって次の近似値を求めていくという方法だから誤差がどんどん蓄積していく。近似の精度を高めるために、長方形の面積を求めるところでもう少し細かくみていくことにする。

##修正オイラー法

かんたんに考えつくのが、長方形の高さとして$f(x_n)$だけではなくて$f(x_{n+1})$との平均を使うことである。これは修正オイラー法と呼ばれる。オイラー法より誤差の収束がはやいが、誤差の計算理論はよくわからない。$\Delta{t}$の区間で単調な関数をイメージすれば、こうすることによって長方形からのはみ出し、削れがすくなくなることは簡単にわかる。

##Runge-Kutta法

オイラー法では（若い方の）一端、修正オイラー法では両端を長方形の高さとして採用したが、採用する数値をさらに巧妙に選んだものがRunge-Kutta法である。
計算量とのバランスをみて、ちょうどよいのが四次のRunge-Kutta法である。四次のRunge-Kutta法では、区間内の4つの点について長方形の値をもとめ、重み付きで平均することで区間内の積分値をもとめる。具体的には以下の通り

$$
x_{n+1}=x_n+\frac{1}{6}(k_1+2k_2+2k_3+k_4)
$$
$$k_1 = f(x_n)\Delta t
$$
$$
k_2 = f(x_n+\frac{1}{2}k_1)\Delta t
$$
$$
k_3 = f(x_n+\frac{1}{2}k_2\Delta t
$$
$$
k_4 = f(x_n+\frac{1}{2}k_3)\Delta t
$$

(なぜかふつうに改行できなかったので段落マシマシでお送りします)(最近は様々な改良版が知られていて、たとえばnumpyのstreamplotでは二次の適応刻み幅制御のRunge-Kutta法が使われている。)
よく使われるのは、いちばん興味のある初期値を与えてみたときの解軌道をRunge-Kuttaを何度もやって関数のグラフみたいに描いてみて、他の初期値については色々与えてみてRunge-Kuttaをそれぞれについて一度だけやってみてベクトル場や傾き場として描くことで全体の雰囲気を掴む、というもの。
以下に実装してみる。

import numpy as np
import matplotlib.pyplot as plt

#調べたい関数を定義
def dotx(x):
  return x*(1-x) 

#Runge-Kutta法でt-x平面にプロット
def rk4th(x0, n, dt):#x0:初期値、n：ｔの範囲、dt：ステップ幅
  t=0
  x=x0
  xvalues=[]
  tvalues=np.arange(0,n,dt)
  scope=len(tvalues)#要素数はこれで統一しておく

  #RungeKuttaをdotxに適用してx0に対応した軌道を描く
  for i in range(scope):
    xvalues.append(x)
    k1 = dotx(x)*dt
    k2 = dotx(x+k1/2)*dt
    k3 = dotx(x+k2/2)*dt
    k4 = dotx(x+k3)*dt
    x = x + (k1+ 2*k2 + 2*k3 + k4)/6
    t += dt  

  #他の初期値から求めてかたむき場を描く
  T, X = np.meshgrid(
    np.arange(0,n,1),
    np.arange(0,2,0.1)
  )
  k1 = dotx(X)*dt
  k2 = dotx(X+k1/2)*dt
  k3 = dotx(X+k2/2)*dt
  k4 = dotx(X+k3)*dt
  V = (k1+2*k2+2*k3+k4)/6
  plt.quiver(T,X,dt,V)
  plt.plot(tvalues, xvalues)
  plt.show()

rk4th(2, 10, 0.1)

以下に結果をしめす。これはほんとうは横軸がtだからベクトルではなくたんに傾きである