講義のページ/Pythonのお勉強/シミュレーション-２Leap Frog法と4次のRunge-Kutta法 - 講義のページ

[ トップ ] [ 編集 | 凍結 | 差分 | 添付 | リロード ] [ 新規 | | 検索 | 最終更新 | ヘルプ | ログイン ]

メニュー †

春学期

秋学期
番外

Arduino †

Arduino

電子工作 †

研究 †

研究のページ

Top/講義のページ/Pythonのお勉強/シミュレーション-２Leap Frog法と4次のRunge-Kutta法

数値計算法 †

お題 †

シミューレートする方程式
$\begin{eqnarray} \frac{dx}{dt} & = & v \\ \frac{dv}{dt} & = & -k x \end{eqnarray}$

Euler 法 †

概要
普通に差分化して計算する方法はオイラー法と呼ばれている。

すべて開く

すべて閉じる

def Euler(x) :
    xn = x[0]
    vn = x[1]
    dt = 0.01
    k  = 3
 
    xnn = xn + vn * dt
    vnn = vn - k * dt * xn
 
    return [xnn, vnn]

ところが、これを実行すると発散する。

すべて開く

すべて閉じる

x0 = [2,3]
 
for i in range(10000):
    x1 <- Euler(x0)
    x0 <- x1
print x0

Google Colaboratory でのプログラムと実行結果

Leap Frog 法 †

概要

$\begin{eqnarray} x_{n+1} &=& x_n + v_{n+1/2} \times {\Delta} t + O(\Delta t^3) \end{eqnarray}$
となり、 $\Delta t^2$ の項が消える。これは、
$\begin{eqnarray} x_{n+1} &=& x_{n+1/2} + x'_{n+1/2} \frac{\Delta t}{2} + \frac{1}{2} x''_{n+1/2} \left(\frac{\Delta t}{2}\right)^2 + \cdots \\ x_{n} &=& x_{n+1/2} - x'_{n+1/2} \frac{\Delta t}{2} + \frac{1}{2} x''_{n+1/2} \left(\frac{\Delta t}{2}\right)^2 + \cdots \end{eqnarray}$
の2式の間で引き算することで得られる。より高次の項は、 $\Delta t^3$ から始まる。

コード
この計算では発散しない。

すべて開く

すべて閉じる

def LeapFrog(x) :
    xn = x[0]
    vn = x[1]
    dt = 0.01
    k  = 3
 
    xnn = xn + vn * dt
    vnn = vn - k * dt * xnn
 
    return [xnn, vnn]
 
x0 = [2,3]
 
for i in range(100000):
    x1 = LeapFrog(x0)
    x0 = x1
print(x0)

Google Colaboratory でのプログラムと実行結果

4次の Runge-Kutta法 †

参考URL
- 倭算数理研究所
概要

$\begin{eqnarray} \frac{dx}{dt} &=& f(x) \\ x_{n+1} & = & x_{n} + \frac{\Delta t}{6} ( k_1+2k_2 + 2k_3 + k_4)\\ \ & & \\ k_1 & = & f(x_n ) \\ k_2 & = & f\left(x_n + \frac{k_1}{2} \Delta t \right) \\ k_3 & = & f\left(x_n + \frac{k_2}{2} \Delta t \right)\\ k_4 & = & f\left(x_n + {k_3} \Delta t \right) \end{eqnarray}$

コード

すべて開く

すべて閉じる

def RungeKutta(x):
    xn = x[0]
    vn = x[1]
    dt = 0.01
    k = 3
 
    k1x = vn
    k1v = -k * xn
 
    k2x = vn + 0.5 * dt * k1v
    k2v = -k * (xn + 0.5 * dt * k1x)
 
    k3x = vn + 0.5 * dt * k2v
    k3v = -k * (xn + 0.5 * dt * k2x)
 
    k4x = vn + dt * k3v
    k4v = -k * (xn + dt * k3x)
 
    xnn = xn + (dt / 6.0) * (k1x + 2 * k2x + 2 * k3x + k4x)
    vnn = vn + (dt / 6.0) * (k1v + 2 * k2v + 2 * k3v + k4v)
 
    return [xnn, vnn]
 
x0 = [2, 3]
 
for i in range(10000):
    x1 = RungeKutta(x0)
    x0 = x1
 
print(x0)

Google Colaboratory でのプログラムと実行結果

添付ファイル:

IMG_3128.jpg 24件 [詳細]

Last-modified: 2023-06-10 (土) 17:46:21