====== 【数値情報処理a】第6回 常微分方程式の数値解法 ======
===== 前回の課題 =====

[[ans5]]

===== 常微分方程式の初期値問題 =====

変数 $t$ と未知の関数 $x(t)$ が含まれる以下のような方程式を考える。
\[
\frac{dx(t)}{dt}=f(x(t),t)\qquad(a\le t\le b)
\]
ここで $f(x(t),t)$ は $x(t)$ と $t$ の任意の関数である。
このような一つの独立変数 $t$ の関数 $x(t)$ とその導関数 $dx(t)/dt$ が入った式を__**常微分方程式**__という。
常微分方程式を初期値 $x(a)=x_{0}$ の元で解けば、$x(t)$ がどのような関数か分かるので、
これを__**常微分方程式の初期値問題**__という。

----
==== 例 ====

ボールを時刻 $t=0 \mathrm{[s]}$ に鉛直上方向に初速度 $v_{0}$ で投げ上げる。
$x(t)$ をボールの地上からの高度とすると、ボールの速度 $\frac{dx(t)}{dt}$ が $t$ とともにどのように変化するかは次の式で与えられる。
\[
\frac{dx(t)}{dt}=v_{0}-gt
\]
これは $f(x(t),t)=v_{0}-gt$ と置いた常微分方程式である。
ここで $g=9.8\mathrm{[m/s^{2}]}$ は重力加速度である。
初期値は $x(0)=0 \mathrm{[m]}$ として $x(t)$ を求めることを考える。

この常微分方程式は次のように積分で解析的に解くことができる。
\[
\int_{x(0)}^{x(t)}dx=\int_{0}^{t}(v_{0}-gt)dt
\]
\[
[x]_{x(0)}^{x(t)}=\left[v_{0}t-\frac{1}{2}gt^{2}\right]_{0}^{t}
\]
\[
x(t)=v_{0}t-\frac{1}{2}gt^{2}
\]
得られた $x(t)$ を $t$ で微分すると元の常微分方程式が得られるので、これが解（解析解）ということが確認できる。

この解が $0$ になるところを探すと
\[
v_{0}t-\frac{1}{2}gt^{2}=t\left(v_{0}-\frac{1}{2}gt\right)=0
\]
なので、$t=0, \frac{2v_{0}}{g}$ のときに $x(t)=0$ になることがわかる。
つまり、時刻 $0$ に高度 $0$ で投げ上げたボールは頂点に達した後、時刻 $\frac{2v_{0}}{g}$ 後に再び高度 $0$ になる。

例えば、$v_{0}=4.9 \mathrm{[m/s]}$ とすると、元の高さに戻るのは $\frac{2v_{0}}{g}=\frac{2\times 4.9}{9.8}=1$ 秒後になる。
これは次のコマンドで $x(t)$ のグラフを描くことで確かめられる。これはいわゆる放物線である。
<code rsplus>
> curve(4.9 * t - 0.5 * 9.8 * t ^ 2, xname = "t", ylab = "x")
</code>

{{ parabola.png?nolink |放物線}}

物理学では常微分方程式の初期値問題に帰結する問題が良く出てくるが、解析的に解けない場合が多い。
そのような場合は数値計算を行う。

===== オイラー法 =====
常微分方程式の初期値問題は、数値微分のように微分の定義を利用して数値的に解くことができる。
$t$ を微小な $h$ ずつ進めて $x(t)$ を計算していく方法で、最も単純な方法を__**オイラー法**__という。

$t$ の区間 $[a,b]$ において等間隔に $n+1$ 個の分点をとり、$n$ 個の小区間に分けるとする。
$t$ の小さい方から $i$ 番目（$i=0,1,2,\cdots,n$）の分点は
\[
t_{i}=a+ih
\]
である。
ここで、$h=(b-a)/n$ は小区間の幅（刻み幅）で、$t_{0}=a$ と $t_{n}=b$ は端点である。

数値微分の前進差分を思い出してみると、関数 $f(x)$ の導関数は
\[
\frac{f(x)}{dx}=f^{\prime}(x)\simeq f_{1}(x, h)=\frac{f(x+h)-f(x)}{h}
\]
で近似していた。
これと同様に $x(t)$ の導関数を
\[
\frac{d x(t)}{dt}\simeq\frac{x(t+h)-x(t)}{h}
\]
で近似すれば、変形して
\begin{eqnarray}
x(t+h)
&\simeq&x(t)+\frac{d x(t)}{dt}h\\
&\simeq&x(t)+f(x(t),t)h
\end{eqnarray}
となる。
この式を使って $x(t)$ から $x(t+h)$ を近似的に求められる。

$t=t_{0}=a$ の初期値 $x(a)=x_{0}$ から始めて $n$ 回繰り返すと
\begin{eqnarray}
x_{1}&=&x_{0}+f(x_{0},t_{0})h\\
x_{2}&=&x_{1}+f(x_{1},t_{1})h\\
x_{3}&=&x_{2}+f(x_{2},t_{2})h\\
&\vdots&\\
x_{n-1}&=&x_{n-2}+f(x_{n-2},t_{n-2})h\\
x_{n}&=&x_{n-1}+f(x_{n-1},t_{n-1})h
\end{eqnarray}
と $t$ が $h$ ずつ後の $x_{i}$ $(i=1,2,3,\cdots,n)$ が逐次求められる。
最後の $x_{n}$ は $t=t_{n}=b$ に対応する $x(b)$ であるが、近似式を繰り返して求めた値なので $x(b)$ と同じ値になるとは限らない。

{{ :numericalipa:euler.png?nolink |オイラー法}}

オイラー法では以下の手順で数値的に計算する。
  - $i=1$ とする。
  - $x_{i}$ を $x_{i}=x_{i-1}+f(x_{i-1},t_{i-1})h$ で計算する。
  - $i$ を 1 増やす。
  - $i<n$ なら 2 に戻る。

----
==== 例 ====

常微分方程式
\[
\frac{dx(t)}{dt}=v_{0}-gt
\]
を初期値 $x(0)=0 \mathrm{[m]}$ としてオイラー法を使って $t=1$ の $x(t)$ を求めてみよう。
$v_{0}=4.9 \mathrm{[m/s]}$ とする。

次のプログラム euler.R は $t=a$ から $t=b$ までオイラー法で数値計算を行う。
比較のために赤い線で解析解を表示する。

<file rsplus euler.R[enable_line_numbers="true",highlight_lines_extra="30"]>
# オイラー法

# 関数
f <- function(x, t) {
  return(4.9 - 9.8 * t)
}

# tの範囲
a <- 0
b <- 1

# 初期値
x0 <- 0
t0 <- a

# 区間数
n <- 10

# 刻み幅
h <- (b - a) / n

# x,tのベクトル
x <- c(x0)
t <- c(t0)

# オイラー法による繰り返し
xi <- x0
ti <- t0
for (i in 1 : n) {
  xi <- xi + f(xi, ti) * h
  ti <- a + i * h
  x <- c(x, xi)
  t <- c(t, ti)
}

# 結果の出力
print(sprintf("h = %f", h), quote = FALSE)
print(sprintf("t = %f, x(t) = %f", t[n + 1], x[n + 1]), quote = FALSE)

# 解との比較
plot(t, x, type = "o", pch = 18)
curve(4.9 * t - 0.5 * 9.8 * t ^ 2, xname = "t", add = TRUE, col = "red")
</file>

区間数 $n$ (刻み幅 $h$) を変えて実行すると以下のようになる。

=== $n=10$ ($h=0.1$) の場合 ===

<code rsplus>
[1] h = 0.100000
[1] t = 1.000000, x(t) = 0.490000
</code>

{{ euler_n10.png?nolink |オイラー法（n=10）}}

=== $n=20$ ($h=0.05$) の場合 ===

<code rsplus>
[1] h = 0.050000
[1] t = 1.000000, x(t) = 0.245000
</code>

{{ euler_n20.png?nolink |オイラー法（n=20）}}

=== $n=40$ ($h=0.025$) の場合 ===

<code rsplus>
[1] h = 0.025000
[1] t = 1.000000, x(t) = 0.122500
</code>

{{ euler_n40.png?nolink |オイラー法（n=40）}}

$n$ が大きくなるにつれて（$h$ が小さくなるにつれて）解析解に近づいていく。
解析解では $x(1)=0$ なので、$n$ を2倍にすると誤差が $1/2$ になっていることが分かる。 

数値的に解を探すときに、$n$ を小さな数からはじめて徐々に $n$ を大きくしていくとよい。
ただし、オイラー法で精度を上げるにはかなり大きな $n$ にしないといけない。
そうすると丸め誤差の問題も出てくる。

===== 小テスト① =====

[[KMS>|Moodle Server（非公式）]]で第6回の小テスト①を受験しなさい。

===== ルンゲ・クッタ法 =====
オイラー法を改良したのが__**ルンゲ・クッタ法**__で、一般的によく使われる方法である。
一般的にオイラー法よりもかなり精度が高い。

ルンゲ・クッタ法では以下の手順で計算する。
  - $i=1$ とする。
  - $x_{i}$ を $x_{i}=x_{i-1}+\frac{1}{6}(k_{1}+2k_{2}+2k_{3}+k_{4})h$ で計算する。
  - $i$ を 1 増やす。
  - $i<n$ なら 2 に戻る。

ここで
\begin{align*}
k_{1}&=f(x_{i-1},t_{i-1})\\
k_{2}&=f(x_{i-1}+k_{1}/2,t_{i-1}+h/2)\\
k_{3}&=f(x_{i-1}+k_{2}/2,t_{i-1}+h/2)\\
k_{4}&=f(x_{i-1}+k_{3},t_{i-1}+h)
\end{align*}
である。

----
==== 例 ====

次のプログラム runge-kutta.R は先程の例の問題において、ルンゲ・クッタ法で数値計算を行う。
比較のために赤い線で解析解を表示する。

<file rsplus runge-kutta.R[enable_line_numbers="true", highlight_lines_extra="30,31,32,33,34"]>
# ルンゲ・クッタ法

# 関数
f <- function(x, t) {
  return(4.9 - 9.8 * t)
}

# tの範囲
a <- 0
b <- 1

# 初期値
x0 <- 0
t0 <- a

# 区間数
n <- 10

# 刻み幅
h <- (b - a) / n

# x,tのベクトル
x <- c(x0)
t <- c(t0)

# ルンゲ・クッタ法による繰り返し
xi <- x0
ti <- t0
for (i in 1 : n) {
  k1 <- f(xi, ti)
  k2 <- f(xi + k1 * 0.5, ti + h * 0.5)
  k3 <- f(xi + k2 * 0.5, ti + h * 0.5)
  k4 <- f(xi + k3, ti + h)
  xi <- xi + (k1 + 2 * k2 + 2 * k3 + k4) * h / 6
  ti <- a + i * h
  x <- c(x, xi)
  t <- c(t, ti)
}

# 結果の出力
print(sprintf("h = %f", h), quote = FALSE)
print(sprintf("t = %f, x(t) = %f", t[n + 1], x[n + 1]), quote = FALSE)

# 解との比較
plot(t, x, type = "o", pch = 18)
curve(4.9 * t - 0.5 * 9.8 * t ^ 2, xname = "t", add = TRUE, col = "red")
</file>

<code rsplus>
[1] h = 0.100000
[1] t = 1.000000, x(t) = -0.000000
</code>

{{ :numericalipa:rk_n10.png?nolink |ルンゲ・クッタ法（n=10）}}

$n$ は $10$ と小さいにもかかわらず、精度の良い解が得られていることが分かる。

ルンゲ・クッタ法はオイラー法よりも良い結果をもたらすことが多いが、ルンゲ・クッタ法が苦手な問題もある。

===== 小テスト② =====

[[KMS>|Moodle Server（非公式）]]で第6回の小テスト②を受験しなさい。

===== 定期試験 =====

  * 日時：2026年5月28日（木）1限
  * 場所：205教室
  * Moodle で行う
  * 持ち込み可（スマホは不可）

今までの小テストと課題に似た問題を出すので復習しておくように。

===== 授業アンケート =====

  * 科目名：数値情報処理a
  * 授業コード：？

「er-学籍番号8桁@k.koeki-u.ac.jp」に空メールを送信すると返事のメールが届くので、後はその指示に従うこと。