\documentclass[12pt]{jarticle}

\oddsidemargin 0mm
\evensidemargin 0mm
\textwidth 165mm
\topmargin -15mm
\textheight 230mm

\author{}
\title{計算物理学: No.4}
\date{}

\begin{document}
\maketitle

物理学を含む自然現象一般はしばしば微分方程式の形で表される。
ここではその中でも基本的な常微分方程式の数値解法を扱う。

\section{常微分方程式の数値解法}

\subsection{常微分方程式の正規形}

常微分方程式の問題は、常に1階の微分方程式の組に置き換えることができる。
例えば、2階の微分方程式
\begin{equation}
 \frac{d^2 y}{d x^2} + q(x) \frac{d y}{d x} = r(x)
\end{equation}
は新しい変数 $z$ を補助的に用いて、連立の1階微分方程式
\begin{eqnarray}
 \frac{dy(x)}{dx} &=& z(x) \nonumber\\
 \frac{dz(x)}{dx} &=& r(x)- q(x) z(x)
\end{eqnarray}
に書き直すことができる。

一般に（多変数、高階の）常微分方程式の問題は、関数 $y^i,i=1,2,\ldots,m$ 
に対して、
\begin{equation}
 \frac{d y^i(x)}{dx} = f^i (x,y^1,y^2,\ldots,y^m), \qquad i=1,\ldots, m
\end{equation}
の正規形をした $m$ 元連立1階微分方程式に帰着する。
初期条件 $y^i(0)$ の組を与えると、これを数値的に解くことができる。

\subsection{Euler 法}

刻み幅 $h$ に対し、独立変数の離散的な値を
\begin{equation}
 x_n = a + n h \qquad (n=0,1,2,\cdots)
\end{equation}
とする。常微分方程式を差分に直すと
\begin{equation}
 \frac{d y^i}{d x} (= f^i (x_n))\approx \frac{y^i_{n+1}-y^i_{n}}{x_{n+1} - x_{n}}
  = \frac{y^i_{n+1}-y^i_{n}}{h}
\end{equation}
となる。これから $x_n$における $y^i(x_n)$ の近似値を $y^i_n$ とした時、
次のステップ $x_{n+1}$ での未知関数 $y^{i}_{n+1}$ を
\begin{equation}
 \fbox{$
 y^i_{n+1} = y^i_n + h f^i_n, \qquad f^i_n \equiv f^i(x_n, y^1_n,\cdots,y^m_n)
 $}
\end{equation}
によって次々と決めるのがオイラー法である。

\subsection{Euler 法の誤差}

$y(x)$ のテイラー展開式を考えると、
\begin{equation}
 y(x+h) = y(x) + h y'(x) + \frac{1}{2} h^2 y'' (x) + \cdots
\end{equation}
従って常微分方程式を差分に直す誤差は、
\begin{equation}
 \frac{d y^i}{d x} = \frac{y^i_{n+1}-y^i_{n}}{h} ＋O（h^2)
\end{equation}
計算する区間を $a$ から $b$ までとし、それを $n$ 等分して計算を進めると
$
 h = (b-a)/n
$
なので、$x=a$ から $x=b$ までの $n$ ステップの誤差は、
\begin{equation}
 n O(h^2) = \frac{(b-a)}{h} O(h^2) = O(h)
\end{equation}
になる。
つまり、刻み幅を半分にすれば誤差は半分になる。

\section{課題}

課題 1.

\begin{equation}
 \frac{d y}{d x} = 1- y^2, y(0) = 0
\end{equation}
の数値解をオイラー法によって求めよ。この問題の解析解 $\tanh x$ と比較せ
よ。刻み幅 $h$ に対する誤差の振舞いを確認せよ。

\bigskip

課題 2.

単振動の微分方程式
\begin{equation}
 \frac{d^2 y}{d x^2} + y = 0
\end{equation}
は2元1階常微分方程式系
\begin{eqnarray}
 \frac{d y}{d x} &=& z, \nonumber \\
 \frac{d z}{d x} &=& -y
\end{eqnarray}
に置き換えられる。

この数値解をオイラー法によって求めよ（初期条件は適当に選べ）。結果を解析解
と比較せよ。刻み幅 $h$ に対する誤差の振舞いを確認せよ。

\end{document}