<jupyter_start><jupyter_text>ルンゲクッタ法による常微分方程式の数値解析 ここでは，オイラー法よりも計算精度が良く，広く用いられている4次のルンゲ・クッタ（Runge-Kutta）法を学ぶ．オイラー法のときと同様に，2階常微分方程式を1階化した$$\begin{cases} \frac{dx(t)}{dt} = v(t), \\ \frac{dv(t)}{dt} = f(x(t), v(t))\end{cases}\tag{2, 再掲}$$を考える．式（2）の数値解 $x_i$, $v_i$，$i=0, 1, 2, \ldots,$は，４次のルンゲ・クッタ法により次のように計算される．$$\begin{eqnarray}x_{i+1} &= x_i + \frac{h}{6} \left(k_{x1} + 2k_{x2} + 2k_{x3} + k_{x4} \right),\tag{9}\\%v_{i+1} &= v_i + \frac{h}{6} \left(k_{v1} + 2k_{v2} + 2k_{v3} + k_{v4} \right)\tag{10}\end{eqnarray}$$ここで，$$\begin{eqnarray}&k_{x1} = v_i, \quad&k_{v1} = f(x_i, v_i),\\%&k_{x2} = v_i + \frac{h k_{v1}}{2}, \quad&k_{v2} = f\left(x_i + \frac{h k_{x1}}{2}, v_i + \frac{h k_{v1}}{2} \right),\\%&k_{x3} = v_i + \frac{h k_{v2}}{2}, \quad&k_{v3} = f\left(x_i + \frac{h k_{x2}}{2}, v_i + \frac{h k_{v2}}{2} \right),\\%&k_{x4} = v_i + h k_{v3}, \quad&k_{v4} = f(x_i + h k_{x3}, v_i + h k_{v3}),\end{eqnarray}$$であり，$h$は積分刻みである． では，ルンゲクッタ法をPythonで実装してみよう．ここでは$$\begin{align}  f = -x(t) \tag{5}\\\end{align}$$とする．また，初期条件を$x(0)=1.0, y(0)=0$とする．このとき，ルンゲクッタ法を実装すると下記のようになる．（TODO） 下記のコードを理解し，セルを実行せよ<jupyter_code># 必要なライブラリのインポート
import numpy as np
from matplotlib import pyplot as plt

# 関数定義
def f(x, v):
  return -x # dv/dt = f(x, y) = -x の実装

# 計算条件を変数として定義
x = 1.0 # xの初期値
v = 0 # yの初期値
t_ini = 0 # 初期時刻
h = 0.2 # 積分刻み
t_max = 10 + h # 終了時刻

# 結果保存用の配列定義
ts = np.arange(start=t_ini, stop=t_max, step=h) # tの時系列
xs = np.zeros(len(ts)) # xの時系列
vs = np.zeros(len(ts)) # yの時系列

# 初期値の代入
xs[0] = x
vs[0] = v

# ルンゲクッタ法
for i in range(1, len(ts)): # i=1 から i=len(ts) までループ
  kx1 = v
  kv1 = f(x, v)
  kx2 = v + h * kv1 / 2.0
  kv2 = f(x + h * kx1 / 2.0, v + h * kv1 / 2.0)
  kx3 = v + h * kv2 / 2.0
  kv3 = f(x + h * kx2 / 2.0, v + h * kv2 / 2.0)
  kx4 = v + h * kv3
  kv4 = f(x + h * kx3, v + h * kv3)
  # v, xの更新
  v += (kv1 + 2.0 * kv2 + 2.0 * kv3 + kv4) * h / 6.0
  x += (kx1 + 2.0 * kx2 + 2.0 * kx3 + kx4) * h / 6.0
  # 結果の保存
  xs[i] = x
  vs[i] = v<jupyter_output><empty_output><jupyter_text>ルンゲ・クッタ法で計算された数値解が配列`xs`, `ys`に格納される．また，数値計算の各時刻は`ts`に格納される．（TODO） 下記のコードを実行し，`ts`を横軸，`xs`を縦軸にプロットせよ．<jupyter_code>plt.plot(ts, xs, '.')
plt.xlabel('t')
plt.ylabel('x')<jupyter_output><empty_output><jupyter_text>（TODO） 下記のコードを実行し，`ts`を横軸，`vs`を縦軸にプロットせよ．<jupyter_code>plt.plot(ts, vs, '.')
plt.xlabel('t')
plt.ylabel('v')<jupyter_output><empty_output><jupyter_text>オイラー法のときと同様に，得られた数値解と厳密解を比較してみよう．（TODO） 下記のコードを実行し，得られた厳密解とオイラー法の数値解を同時にplotせよ．<jupyter_code>xs_exact = np.cos(ts) # 厳密解

# 厳密解とオイラー法の数値解をplot
plt.plot(ts, xs_exact, '-r', label='Exact')
plt.plot(ts, xs, '.g', label='RK') # 既にルンゲ・クッタ法で計算済み
plt.xlabel('t')
plt.ylabel('x')
plt.legend() # labelを表示する<jupyter_output><empty_output><jupyter_text>（TODO） 次のコードを実行し，ルンゲ・クッタ法で得た数値解の配列`xs`と，厳密解の配列`xs_exact`のMAEを計算せよ．<jupyter_code>np.mean(np.abs(xs - xs_exact))<jupyter_output><empty_output>
StarcoderdataJupyterStructuredCleanDedup