ルンゲ=クッタ法の気持ちを理解したい(1)

ルンゲ=クッタ法 - Wikipedia に関する課題が出たが、授業をロクに聞いていなかったため式だけを眺めていても気持ちが全く分からずに苦心した。そこでブログの投稿の練習も兼ねて記事にしてみる。

ルンゲ=クッタ法とは

Wikipediaによると

ルンゲ=クッタ法(英: Runge–Kutta method)とは、数値解析において微分方程式の初期値問題に対して近似解を与える一連の方法である。

上のページと同じく、{ \frac{\mathrm{d}y}{\mathrm{d}t} = f(t,y),\hspace{4mm} y(t_0)=y_0 } という初期値問題が与えられていると以下では考えることにする。

現実的な時間で十分な近似精度を得やすいので、数値解析の分野では重宝されるようだ。

オイラー

ルンゲ=クッタ法の一番簡単なケースとして、(前進)オイラー法がある。

以降に出てくる近似解を求めるアルゴリズムは、十分小さな刻み幅{ h }を設定し、時刻 { t_i = t_0 + i \cdot h } における { y(t_i) = y(t_0 + ih) } の近似値 { y_i } を順番に求めていくという部分は同じであるので覚えておいてほしい。

オイラー法は以下のような単純なアルゴリズムである。

{ y_{i+1} = y_i + h \cdot f(t_i,y_i) }

 { \frac{\mathrm{d}y}{\mathrm{d}t} = f(t,y) } であることを思い出せば、テイラー展開の1次までの項 { y(t_{i+1}) = y(t_i+h)=y(t_i) +   h \frac{\mathrm{d}y} {\mathrm{d}t}(t_i)  } を元にした近似式であることが分かる。

実はこの計算では十分な精度が得られない。

直観的には、直線近似をしているので { y_i }{ y(t_i) } の誤差がそのまま累積されていくことから、刻み幅として直線とみなせる十分に狭い区間をとる必要があると理解できる。

正確に言うと各ステップでの誤差が { {\rm O}(h^2) } であり、時刻{ t }までのステップ数は{ \frac{t}{h} } であるので、(時刻を固定したとき)誤差は { {\rm O}(h) } となる。

これは1桁精度をよくするために10倍の計算が必要になるということである。余りにも過大であるのが分かるだろう。

オイラー法のサンプルプログラム

次の初期値問題をオイラー法を用いて解いて、解析解と比較してみる。

{ \frac{\mathrm{d}y}{\mathrm{d}t} = f(t,y) = 2ty }

 {  y(0)=1}

解析的に解を求めると { y(t) = e^{t^2} } となる。 { y(1) = e }オイラー法で近似的に求めた値の誤差を両対数でプロットしてみる。

Pythonという言語が流行っているらしいのでPythonを使って書いたが、Pythonらしくない書き方をしていたらコメントで教えてほしい。

import numpy as np
import matplotlib.pyplot as plt

def f(t,y):
    return 2.0 * t * y

def ans(t):
    return np.exp(t**2.0)

m = 2
M = 7
N = np.logspace(m,M,M-m+1,10,dtype=int) # 10^m ... 10^M
Y = np.zeros(M-m+1)
for n in N:
    h = 1.0 / float(n) # h = 1 / n
    y = 1.0
    for i in range(0,n):
        t = float(i)/float(n)
        y = y + h * f(t,y)
    Y[np.where(N==n)] = y

plt.xscale("log")
plt.yscale("log")
plt.title("Euler method")
plt.xlabel("step width")
plt.ylabel("error")
plt.grid(which="both")
plt.plot(1.0/N,abs(Y-ans(1.0)),"o-",markersize=13) # ( h , |y - e| )
plt.show()

f:id:buko1062000:20161225233604p:plain

まさかオイラー法だけで(1)が終わってしまうとは思わなかった。 (2)へ (そのうち)