スポンサードリンク
今回は、pythonの勉強の一環でやってきた簡単な数値計算のソースコードを紹介していきます。
一次元減衰強制振動
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 |
# -*- coding: utf-8 -*- import numpy as np import matplotlib.pyplot as plt from mpl_toolkits.mplot3d import Axes3D import math #####################設定パラメーター################## k = 1.0 m = 1.0 w = 0.1 u = 5 ################################################### tr=[] xr=[] te=[] xe=[] def g(t, x): return np.array([x[1], (-k/m*x[0])-2*w*x[1]+u]) x = np.array([0.01, 0.0]) # 初期条件 x(0) = 1, v(0) = 0 dt = 0.1 ######################ルンゲクッタ######################### for i in range(800): t = i*dt k1 = g(t, x) k2 = g(t + dt/2, x + k1*dt/2) k3 = g(t + dt/2, x + k2*dt/2) k4 = g(t + dt, x + k3*dt) x += (k1 + 2*k2 + 2*k3 + k4)*dt/6 xr.append(x[0]) tr.append(t) #####################オイラー法############################# x1 = np.array([0.01, 0.0]) # 初期条件 x(0) = 1, v(0) = 0 for i in range(800): t = i*dt k1 = g(t, x1) k2 = g(t + dt, x1 + k1*dt) x1 += (k1 + k2)*dt/2 xe.append(x1[0]) te.append(t) ##################グラフセッティング###################### fig = plt.figure() ax = plt.gca() ax.plot(tr, xr, c='b', label="Runge-Kutta") ax.plot(te, xe, c='r', label="Eular") ax.set_xlabel("time") ax.set_ylabel("x(t)") ax.set_title("x-t graph") plt.legend() plt.show() |
x-t(振幅ー時間)グラフは以下のようになります。
両解法（オイラー法とルンゲクッタ法）において、以前にも記事に書いた通り、一次元で方程式が単純ということもあり、そんなに解に差は見られません。
python初心者の方は、設定パラメータを動かして、どんなグラフになるか、遊んでみてください。
以下に、一次元減衰強制振動のモデルを示します。
ばね定数：k、減衰係数：c、質量：m、強制外力：m/u
スポンサードリンク
ロトカヴォルテラ方程式
自然界における捕食者と被食者の個体数変動の数理モデルです。
以下の数式で定義されています。
\[
\left\{
\begin{array}{ll}
\frac{dx}{dt}= ax-bxy\\
\frac{dy}{dt}= cxy-dy
\end{array}
\right.
\]
xが被食者の個体数、yは捕食者の個体数、tは時間、4つの係数a、b、c、dは正の実数パラメータです。
以下は、a=8, b=3, c=4 ,d=18の時の計算結果。
周期関数になることがわかりますが、直感的に何が起きているか不明なので、定性的に述べていきます。ライオンとシマウマの例で考えると、
- 捕食者であるライオンが増加
- それにより被食者であるシマウマが減少
- エサ（シマウマ）が減少するので、ライオンは減少
- 天敵（ライオン）が減るので、シマウマは増加
- １以降を繰り返す
といったことを繰り返しています。よくできたモデルですが、簡易化されているため、実際はもっと複雑です。以下は上で述べたことを図にしたものです。赤が捕食者（ライオン）で、青が被食者（シマウマ）です。
ルンゲクッタ法とオイラー法の比較
オイラー法とルンゲクッタ法の誤差の確認として、被食者数（上で言えばシマウマ）をグラフにしてみました。
ぱっと見では、解に差が内容に思われますが、縮尺を変更すると
時刻t=122では,0.5くらいの差があります。
割合でいうと、0.5/12=4%くらいの差があります。この程度の計算回数で誤差４パーは論外ですね笑
まあルンゲクッタが正解かというとそうでもないですが。
ついでに、グラフが直線的になっているのは、計算する際の微分間隔が雑なためです。
以下にコードを載せますが、そこの「微分間隔」と書いてあるところを小さくしていけば、理論的な微分に近づき、滑らかな解（グラフ）が得られます。
計算に用いたプログラム
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 |
# -*- coding: utf-8 -*- import numpy as np import matplotlib.pyplot as plt from mpl_toolkits.mplot3d import Axes3D # rorenz equation def lotka(x, y, a=8, b=3, c=4 ,d=18): x_dot = a*x -b*x*y y_dot = c*x*y -d*y return np.array([x_dot, y_dot]) t=0 ##### 初期時間###### dt = 0.001 ##### 微分間隔###### stepCnt = 300 ######時間ステップ## # Need one more for the initial values xs = np.empty((stepCnt + 1,)) ys = np.empty((stepCnt + 1,)) # Setting initial values xs[0], ys[0] = (10, 7) k0=[0,0,0] k1=[0,0,0] k2=[0,0,0] k3=[0,0,0] listData=[0] listData1=[0] ########ルンゲクッタ#################### # Stepping through "time". for i in range(stepCnt): x,y=xs[i],ys[i] k0 = dt * lotka(x,y) k1 = dt * lotka(x+k0[0]/2., y+k0[1]/2.) k2 = dt * lotka(x+k1[0]/2., y+k1[1]/2.) k3 = dt * lotka(x+k2[0], y+k2[1]) dx = (k0[0]+2.0*k1[0]+2.0*k2[0]+k3[0])/6.0 dy = (k0[1]+2.0*k1[1]+2.0*k2[1]+k3[1])/6.0 xs[i+1] = xs[i] + dx ys[i+1] = ys[i] + dy listData.append(i) # Need one more for the initial values xs2 = np.empty((stepCnt + 1,)) ys2 = np.empty((stepCnt + 1,)) # Setting initial values xs2[0], ys2[0] = (10, 7) k02=[0,0,0] k12=[0,0,0] listData1=[0] ########オイラー法#################### # Stepping through "time". for i in range(stepCnt): x,y=xs2[i],ys2[i] k02 = dt * lotka(x,y) k12 = dt * lotka(x+k02[0], y+k12[1]) dx = (k02[0]+k12[0])/2.0 dy = (k02[1]+k12[1])/2.0 xs2[i+1] = xs2[i] + dx ys2[i+1] = ys2[i] + dy listData1.append(i) ###########all-time-datasiries(graph-setting)################## fig = plt.figure() ax = plt.gca() ax.plot(listData, xs, c='b', label="Runge-Kutta") ax.plot(listData1, xs2, c='r', label="Eular") ax.set_xlabel("time") ax.set_ylabel("Populazione") ax.set_title("Lotka-Volterra-model") plt.legend() plt.show() |
以上です。基礎的な微分方程式や運動方程式の数値解法の勉強する際には役立てみてください。
スポンサードリンク
