<jupyter_start><jupyter_text>ルンゲ・クッタ法の誤差<jupyter_code>import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
%matplotlib inline
sns.set('notebook', 'whitegrid', 'dark', font_scale=2, rc={"lines.linewidth": 2, 'grid.linestyle': '--'})<jupyter_output><empty_output><jupyter_text>オイラー法$x'=f(t,x)$の初期値問題に対して漸化式$$    x_{n+1} = x_n + h f(t_n, x_n)$$によって近似解を計算する<jupyter_code>def Euler(t, x, f, h):
    return x + h * f(t, x)<jupyter_output><empty_output><jupyter_text>$x'=f(t,x)$の右辺の定義<jupyter_code>def func(t,x):
    return x<jupyter_output><empty_output><jupyter_text>$x'=x, x(0)=1$に対してオイラー法で$t=1$における近似解を計算し，真の解$x(1)=e$との誤差を調べる<jupyter_code>x0 = 1.0 # 初期値
a = 0.0 # 初期時刻
b = 1.0 # 最終時刻
Steps = [2**n for n in range(4,17)]
H = [(b - a) / N for N in Steps]
# いろんな刻み幅でx(1)の近似を計算する
Err = []
for n, N in enumerate(Steps):
    h = H[n]
    t = a
    x = x0
    for n in range(N):
        x = Euler(t, x, func, h)
        t = a + (n + 1) * h
    err = abs(x - np.exp(1))
    Err.append(err)
fig = plt.figure(figsize=(8,4))
ax = fig.add_subplot(111)
ax.set_xscale('log')
ax.set_yscale('log')
ax.set_xlabel("$h$")
ax.set_ylabel("error")
ax.plot(H, Err, '-o', label="Euler")
ax.plot(H, H, '-', label="$h$")
ax.legend()<jupyter_output><empty_output><jupyter_text>刻み幅$h$のときの誤差を$e_h$とおき，$$p_h := \log_2 \frac{e_{h}}{e_{2h}}$$を観察する．1次精度なら$p_h$は$1$に近いはず．<jupyter_code>P1 = [np.log2(Err[k]/Err[k+1]) for k in range(0, len(Err)-1)]
fig = plt.figure(figsize=(8,4))
ax = fig.add_subplot(111)
ax.set_xscale('log')
ax.set_yscale('linear')
ax.plot(H[0:-1], P1, '-ob')<jupyter_output><empty_output><jupyter_text>ホイン法\begin{align}    k_1 &= f(t_n, x_n),\\    k_2 &= f(t_n + h, x_n + h k_1),\\    x_{n+1} &= x_n + \frac{h}{2} (k_1 + k_2)\end{align}<jupyter_code>def Heun(t, x, f, h):
    k1 = f(t, x) 
    k2 = f(t + h, x + h * k1)
    return x + 0.5 * h * (k1 + k2)
Err2 = []
for n, N in enumerate(Steps):
    h = H[n]
    t = a
    x = x0
    for n in range(N):
        x = Heun(t, x, func, h)
        t = a + (n + 1) * h
    err = abs(x - np.exp(1))
    Err2.append(err)
fig = plt.figure(figsize=(8,4))
ax = fig.add_subplot(111)
ax.set_xscale('log')
ax.set_yscale('log')
ax.set_xlabel("$h$")
ax.set_ylabel("error")
ax.plot(H, Err, '-o', label="Euler")
ax.plot(H, H, '-', label="$h$")
ax.plot(H, Err2, '-o', label="Heun")
ax.plot(H, np.array(H)**2, '-', label="$h^2$")
ax.legend()
P2 = [np.log2(Err2[k]/Err2[k+1]) for k in range(0, len(Err)-1)]
fig = plt.figure(figsize=(8,4))
ax = fig.add_subplot(111)
ax.set_xscale('log')
ax.set_yscale('linear')
ax.plot(H[0:-1], P2, '-ob')<jupyter_output><empty_output><jupyter_text>古典的ルンゲ・クッタ法$$\begin{aligned}k_1 &= f(t_n, x_n),\\k_2 &= f(t_n + \tfrac{h}{2}, x_n + \tfrac{h}{2}k_1),\\k_3 &= f(t_n + \tfrac{h}{2}, x_n + \tfrac{h}{2}k_2),\\k_4 &= f(t_n + h, x_n + h k_3),\\x_{n+1} &= x_n + \frac{h}{6}(k_1 + 2 k_2 + 2 k_3 + k_4)\end{aligned}$$<jupyter_code>def RK4(t, x, f, h):
    k1 = f(t, x)
    k2 = f(t+0.5*h, x+0.5*h*k1)
    k3 = f(t+0.5*h, x+0.5*h*k2)
    k4 = f(t+h, x+h*k3)
    return x + h * (k1 + 2 * k2 + 2 * k3 + k4) / 6
Err3 = []
for n, N in enumerate(Steps):
    h = H[n]
    t = a
    x = x0
    for n in range(N):
        x = RK4(t, x, func, h)
        t = a + (n + 1) * h
    err = abs(x - np.exp(1))
    Err3.append(err)
fig = plt.figure(figsize=(8,6))
ax = fig.add_subplot(111)
ax.set_xscale('log')
ax.set_yscale('log')
ax.set_xlabel("$h$")
ax.set_ylabel("error")
ax.plot(H, Err, '-o', label="Euler")
ax.plot(H, H, '-', label="$h$")
ax.plot(H, Err2, '-x', label="Heun")
ax.plot(H, np.array(H)**2, '--', label="$h^2$")
ax.plot(H, Err3, '-v', label="RK4")
ax.plot(H, np.array(H)**4, '-.', label="$h^4$")
ax.legend()
# plt.savefig("rk_errors.pdf", bbox_inches='tight')
P3 = [np.log2(Err3[k]/Err3[k+1]) for k in range(0, len(Err3)-1)]
fig = plt.figure(figsize=(8,4))
ax = fig.add_subplot(111)
ax.set_xscale('log')
ax.set_yscale('linear')
ax.plot(H[0:-1], P3, '-ob')<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.24473216552424473, 'fraction_numerical': 0.04112718964204113, 'mean_word_length': 4.077319587628866, 'pattern_counts': {'":': 1, '<': 29, '<?xml version=': 0, '>': 29, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '7', 'dir': 'jupyter-structured-clean-dedup', 'id': '130033', 'license': 'BSD-3-Clause', 'n_tokens_mistral': 2123, 'n_tokens_neox': 2097, 'n_words': 354, 'path': 'notebooks/simulation/Runge-Kutta_error_estimates.ipynb', 'repo_name': 'tmiyaji/sgc164', 'size': '3939'}