\documentclass[12pt]{jarticle}
\usepackage[dvips]{graphics}

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

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

\begin{document}
\maketitle

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

\subsection{Runge-Kutta 法}

常微分方程式の数値解法で、オイラー法から修正オイラー法にすると
精度が劇的に改善した。より高精度な公式を使うと計算の能率はさらに上がる。
この中で古典的なものの１つは Runge-Kutta 法である。
4 次の Runge-Kutta 法
\begin{eqnarray}
k_1 &=& h f(x_n,y_n) \nonumber\\
k_2 &=& h f(x_n + \frac{h}{2}, y_n +\frac{1}{2} k_1 ) \nonumber\\
k_3 &=& h f(x_n + \frac{h}{2}, y_n +\frac{1}{2} k_2 ) \nonumber\\
k_4 &=& h f(x_n + h, y_n +  k_3 ) \nonumber\\
y_{n+1} &=& y_n +\frac{1}{6}(k_1 + 2 k_2 +2 k_3 + k_4) 
+O (h^5)
\end{eqnarray}
で、誤差は累積すると
$O(h^4)$
である。中間値として4回 $f(x,y)$ を計算し
なければならないが、精度の向上が著しいのでメリットがある。

Runge-Kutta 法の導出と、誤差の評価については
「工学のための応用数値計算法入門」 第8章3節 \cite{Shinozaki-Matsushita}
など参照。

多元連立の常微分方程式に対しては、
$k_1,k_2,k_3,k_4$ を配列に置き換えれば良い。手続きは修正オイラー法の連立
常微分方程式の項も参考のこと。

\subsection{その他}

常微分方程式の数値解法には、Runge-Kutta 法の他、予測子・修正子法、
Richardson 補外を応用したBulirsch-Stoer 法などある。

参考文献 \cite{Shinozaki-Matsushita,Numerical-Recipes}

\section{Fortran 文法について}

\subsection{配列}
Fortran では、配列宣言は基本的には 1 から始まる。0 または負の数から配列
宣言したい時には
\begin{verbatim}
	real*8 x(0:100,-10:20)
\end{verbatim}
と言った風にする。

よくあるミスは、配列宣言の範囲外の変数を読み込み・書き出しするものである。
この場合、コンパイルエラーなどにはならないのでしばしば気づき難いが、予想
とは異なる出力結果となる。

\section{レポート課題}

No.5 での課題の他、
\begin{equation}
 \frac{ d^2 y}{d x^2 } + \sin (y) =0
\end{equation}
（大振幅の振子の問題;初期値は各自与えること）
をオイラー法、修正オイラー法で扱ってみよ。刻み幅に対して誤差がどのように
収束するかも調べること。

単に数字を出力結果として出すだけでなく、gnuplot などでグラフの形にしてみ
よ。調和振動子や、振り子の課題では、$x,y(x)$ のグラフだけでなく、
$y,dy/dx$ (位相空間)でも表示して、物理的意味の考察もしてみよ。振り子では、
大振幅で回転運動になることに注意。

余裕があれば、
4次の Runge-Kutta 法にも取り組んでみよ。

この他、各自物理的に興味ある常微分方程式も扱ってみよ。

レポート提出期限は 6月20日まで。

\begin{thebibliography}{99}
\bibitem{Shinozaki-Matsushita}
	「工学のための応用数値計算法入門」 篠崎、松下 共著、コロナ社
\bibitem{Numerical-Recipes}
	``Numerical Recipes in C'',
	``Numerical Recipes in Fortran''(Cambridge University Press)
\end{thebibliography}

\end{document}
