%=====================================================================
% 秋田高専 5E 計算機応用 テキスト
%テーマ 練習(5) 関数
%    本テキストの内容
%        ・関数
%          ・関数とは
%          ・関数の作り方と使い方
%          ・データを渡す方法
%          ・戻り値を返す方法
%          ・グローバル変数によるデータの受け渡し
%          ・main関数への引数渡し
%          ・再帰コール
%
% last updated 2005.5.25
%    created by  Masashi Yamamoto
%     e-mail yamamoto@akita-nct.jp
%=====================================================================
\documentclass[10pt,a4paper]{jarticle}
\usepackage{graphicx,amsmath,amssymb,ascmac,theorem,float}
\usepackage{listings, jlisting, html}
\oddsidemargin 0mm  %左の余白 25.4mm-0mm奇数ページ
\evensidemargin 0mm %左の余白 25.4mm-0mm偶数ページ
\textwidth 160mm
\newcommand{\command}[1]{「\texttt{#1}」}
\newcommand{\tw}[1]{\texttt{#1}}
\newcommand{\vm}[1]{\boldsymbol{#1}}
\newcommand{\diff}[0]{\mathrm{d}}
\newcommand{\pdiff}[3]{
 \if 1#1 \frac{\partial #2}{\partial #3}
 \else \frac{\partial^{#1} #2}{\partial #3^{#1}}\fi
}
%
\renewcommand{\lstlistingname}{リスト}
\lstset{language=C,%
        basicstyle=\footnotesize,%
        commentstyle=\texttt,%
        classoffset=1,%
        keywordstyle=\bfseries,%
	frame=tRBl,framesep=5pt,%
	showstringspaces=false,%
        numbers=left,stepnumber=1,numberstyle=\footnotesize%
	}%
%
\begin{document}
\title{C言語の学習\\プリプロセッサー}
\author{山本昌志\thanks{独立行政法人　秋田工業高等専門学校　電気工学科}}
\date{2004年6月2日}
\maketitle
%
%=====================================================================
\section{本日の学習内容}
%=====================================================================
本日の内容は、教科書の15章である。主にデータ構造に関わる12〜14章は省く。数値計算
のデータ構造は単純なので、ほとんど配列で間に合う。興味のある者はそれらのデータ構
造に関係する章を自分で学習せよ。
%
\begin{list}{}{\leftmargin=5zw \rightmargin=2zw \labelsep=1zw
 \labelwidth=5zw \itemindent=0pt}
 \item[\textbf{15章}]\textbf{プリプロセッサー}
\end{list}
%
%
%=====================================================================
\section{プリプロセッサー(15章)}
%=====================================================================
%---------------------------------------------------------------------
\subsection{プリプロセッサの実行タイミング(p.282)}
%---------------------------------------------------------------------
\command{gcc}を使って、ソースプログラムをマシン語に直すとき、通常はコンパイルす
ると言うが、実際にはコンパイル以外の作業も行っている。教科書のp.282に示されてい
るように、コンパイルの前にプリプロセッサーがソースファイルを書き直している。そして、
コンパイルの後、リンカーがいろいろなファイルを寄せ合わせて実行可能なファイルを作
成する。

プリプロセッサーができることは、教科書のp.283の表に示されている。いろいろな機能
があるが、本講義で使うのは
\begin{itemize}
 \item　ファイルの挿入を行う\tw{\#include}
 \item　マクロ定義を行う\tw{\#define}
\end{itemize}
くらいである。
%---------------------------------------------------------------------
\subsection{プリプロセッサの使い方}
%---------------------------------------------------------------------
プリプロセッサーの書式は、以下の通りです。今まで、さんざん書いてきているので、大
体理解はできると思います。
\begin{quote}
 \setlength{\baselineskip}{12pt}
 \begin{verbatim}
#コマンド　パラメーターリスト
 \end{verbatim}
\end{quote}

このプリプロセッサーの書式には、以下の約束があります。
\begin{itemize}
\item \#の前後に空白が有っても良い。\#とコマンドの間に空白が有っても良いのですが、
      プログラムがわかりにくくなりますので、通常は\#とコマンドの間に空白は置きま
      せん。
\item プリプロセッサーコマンドは、その行で終わりです。C言語の文のように';'が来る
      まで有効ということは有りません。したがって、文の区切りを表す';'もありませ
      ん。
\item ただし、プリプロセッサーコマンドを複数行にわたって、記述したい場合、行の終
      わりに'\'をつけて、次行と接続することができます。
\end{itemize}
%
%----------------------------
\subsubsection{ファイルの挿入}
%----------------------------
\tw{\#include “ファイル名”}は、指定されたテキストファイルとこの行を置き換えま
す。ファイル名が\tw{<filename>} で書かれた場合には、システムに定義されたファイルです。
unix では、\tw{/usr/include} にあるファイルで置き換えられます。プログラマーが作
成したファイルを挿入したい場合は、　\tw{\#include "myfile.h"} のように記述します。
プログラマーが自分でヘッダーファイルを書くことがあります。例えば、大規模なプログ
ラムを作成する場合、ファイルは１つではなく、いろいろな部分を別々のファイルに書い
て、それをあとであわせて、１つのプログラムにすることがあります（分割コンパイル）。
その場合には、共有する構造体や大域変数を１つのファイルにして、\tw{myfile.h} とい
うファイルを作り、それぞれのファイルで\tw{\#include} することがよくあります。こ
のようにすると、同じ文をそれぞれのファイルに記述する必要が無くなり、タイピングが
楽になるとともに、ミスが減ります。

ここでの学習では、分割コンパイルをするほどのプログラムを書くことはないでしょう。
将来、比較的大きなプログラムを書くときに、勉強してください。　本来、次のサンプル
のような使い方はしませんが、\tw{\#include}の動作の理解のために、実行させてみてください。
\tw{\#include}により、そこの行が指定のファイルに置き換わっているのが理解できます。　実
行の結果、使い方によっては、便利そうであることが分かりましたか?。ただ、注意して
欲しいことは、この行の置き換えはコンパイルの前に行われることです。コンパイルによ
り実行ファイルが出来上がった後に、ヘッダーファイルを書き換えても、それは反映され
ません。言いたいことは、ここにデータを書いた場合、その変更を反映させるためには、
再度コンパイルが必要ですということです。

\tw{Hello World !!}を出力するおなじみのプログラムです。プログラムの一部が、
\tw{myfile.h}にかかれています。以下が、C言語のソースプログラムです。変なプログラ
ムですが、ヘッダーファイルがちゃんと書かれていれば、実行可能です。
%
\begin{lstlisting}[caption=マクロを使ったプログラム,label=prog:macro]
#include <stdio.h>
#include "myfile.h"
 return(0);
}
\end{lstlisting}

これをコンパイルして、実行ファイルを作るためには、以下のヘッダーファイル
(\tw{myfile.h})も必要です。当然、ファイル名はソースプログラムの指定通りにする必
要があります。いかがヘッダーファイルです。
%
\begin{lstlisting}[caption=ヘッダーの例,label=prog:my_head]
int main(){
 printf("Hello World !!\n");
\end{lstlisting}
%
%------------------------
\subsubsection{マクロ定義}
%------------------------
\tw{\#include}はその行を指定のファイルで置き換えましたが、\tw{\#define}はファイル
内の文字列を指定の通りに置き換えます。単純な文字列の置き換えと、引数を含む文字列
の置き換えがあります。この引数を含む文字列の置き換えをマクロ定義と言います。それ
ぞれについて、説明します。

\underline{文字列置換}\\
これは、以下のように記述します。
\begin{quote}
 \setlength{\baselineskip}{12pt}
 \begin{verbatim}
#コマンド　パラメーターリスト
 \end{verbatim}
\end{quote}
　
このプリプロセッサーコマンドがあると、この行以降の'文字列1'が'文字列2'に、
\begin{quote}
 \setlength{\baselineskip}{12pt}
 \begin{verbatim}
#undef 文字列1
 \end{verbatim}
\end{quote}
があるまで、置き換わります。もし、\tw{\#undef}が無いと、そのファイルの最後の行までコマ
ンドは有効になります。通常、文字列1はすべて大文字で記述します。別に小文字でも良
いのですが、ほとんどのプログラマーは大文字を使います。その方が、プログラムがわか
りやすくなります。

\underline{引数付き置換}\\
先に説明した単純な文字列の置き換えと似ています。ただ、関数のように引数が使えます。
記述は、以下の通りです。
\begin{quote}
 \setlength{\baselineskip}{12pt}
 \begin{verbatim}
#define マクロ名(引数)　引数を含む文字列
 \end{verbatim}
\end{quote}

このプリプロセッサーコマンドも、
\begin{quote}
 \setlength{\baselineskip}{12pt}
 \begin{verbatim}
#undef マクロ名(引数)
 \end{verbatim}
\end{quote}
があるまで、有効です。もし、\tw{\#undef}が無いと、そのファイルの最後の行までコマンドは
有効になります。
%---------------------------------------------------------------------
\subsection{プリプロセッサを使ったプログラム例}
%---------------------------------------------------------------------
%------------------------
\subsubsection{微分方程式}
%------------------------
\tw{\#define}を使ったサンプルとして、2階の常微分方程式の解を数値計算で求めます。2階の
 常微分方程式の一般形は、
%
\begin{align}
 \frac{\diff^2 y}{\diff x^2}=g(x,y)
 \label{eq:general_2_diff_eq}
\end{align}
%
です。いきなり、この微分方程式の解を求めるとなると、難しそうですが、いたって簡単
です。次ページに、プログラムを載せてあります。このプログラムのなかで、微分方程式
を計算しているのは、たった1行です。プログラムの後半にある
%
\begin{quote}
 \setlength{\baselineskip}{12pt}
 \begin{verbatim}
y[i]=2*y[i-1]+DY2DX2(x-dx,y[i-1])*dx*dx - y[i-2];
 \end{verbatim}
\end{quote}
%
だけです。たったこの1行で、皆さんが苦労した微分方程式が計算できます。

簡単なだけではなく、どんな形の微分方程式でも解けます。皆さんが数学の授業で苦労し
た微分方程式は、解析解(解が初等関数の組み合わせ)があるものばかりです。この方法を
使うと、解析解が無いものまで計算可能です。驚いたでしょう。このようなことから、コ
ンピューターによる数値計算では、実に有益な情報が得られることが理解できるでしょう。

それでは、重要なこの1行の内容を示します。やっと数値計算の授業らしくなってきまし
た。まず、微分方程式を数値計算で解く場合、テイラー展開が重要になります。テイラー展開は、
%
\begin{align}
 f(x_0+\Delta x)
 &=f(x_0)+f^\prime(x_0)\Delta x
 +\frac{f^{\prime\prime}(x_0)}{2!}\Delta x^2
 +\frac{f^{(3)}(x_0)}{3!}\Delta x^3
 +\frac{f^{(4)}(x_0)}{4!}\Delta x^4
 +\cdots \nonumber \\
 &=\sum_{n=0}^\infty\frac{f^{(n)}(x_0)}{n!}\Delta x^n
\label{eq:teylar_plus}
\end{align}
%
です。これは、$x=x_0$に関するすべての導関数の値が分かれば、他の場所の関数の値がわか
るというものです。不思議ですね。次に、$-\Delta x$を考えます。式
(\ref{eq:teylar_plus})と同じように、
%
\begin{align}
f(x_0-\Delta x)
 &=f(x_0)
 -f^\prime(x_0)\Delta x
 +\frac{f^{\prime\prime}(x_0)}{2!}\Delta x^2
 -\frac{f^{(3)}(x_0)}{3!}\Delta x^3
 +\frac{f^{(4)}(x_0)}{4!}\Delta x^4
 -\cdots \nonumber \\
\label{eq:teylar_minus}
\end{align}
%
となります。次に、(\ref{eq:teylar_plus})と(\ref{eq:teylar_minus})式の各辺どおし
を足し合わせます。すると、
%
\begin{align}
 f(x_0+\Delta x)+f(x_0-\Delta x)
 &=2f(x_0)
 +f^{\prime\prime}(x_0)\Delta x^2
 +O(\Delta x^4)
\end{align}
%
となります。この式の4次以上の項を無視して、ちょっとだけ変形すると、
%
\begin{align}
f(x_0+\Delta x)=2f(x_0)+f^{\prime\prime}(x_0)\Delta x^2-f(x_0-\Delta x)
\end{align}
%
となります。$y=f(x)$なので、
%
\begin{align}
 f(x_0+\Delta x)=2f(x_0)+\frac{\diff^2 y}{\diff x^2}\Delta x^2-f(x_0-\Delta x)
\end{align}
%
となります。この式の右辺第2項は、与えられた微分方程式から、計算できます。式
(\ref{eq:general_2_diff_eq})を代入すると
%
\begin{align}
 y(x_0+\Delta x)=2y(x_0)+g(x_0,\,y_0)\Delta x^2-y(x_0-\Delta x)
 \label{eq:y_dx}
\end{align}
%
となります。

この式の右辺は、$y(x_0)$と$y(x_0-\Delta x)$と$g(x_0, y_0)$の値がわかれば、
$y(x_0+\Delta x)$の値が計算できることを示しています。関数$g(x,y)$は問題により与えれ
ています。残りは、通常初期条件として、与えられます。すると、$y(x_0+\Delta x)$が計算できま
す。この値をまた、式(\ref{eq:y_dx})に代入して、今度は、$y(x_0+2\Delta x)$を計算します。これを順次繰り
返せば、いくらでも計算できます。この?xごとの計算結果が、まさに2階の常微分方程式
(\ref{eq:general_2_diff_eq})の解になっています。
%

リスト\ref{prog:diff_eq}に示すプログラムでは、
\begin{align}
 \frac{\diff^2 y}{\diff x^2}=-y
\end{align}
%
を計算している。初期条件は、
%
\begin{align}
 &y(0)=0 \\
 &\frac{\diff y}{\diff x}=1\qquad\text{$x=0$のとき}
\end{align}
%
です。この初期条件から、y(?x)を計算する必要があります。そのために、テーラー展開を使います。以下の通りです。
%
\begin{align}
y(\Delta x)=y(0)+y^\prime(0)\Delta x+\frac{y^{\prime\prime}}{2}\Delta x^2
\end{align}
%
これで、$y(0)$と$y(\Delta x)$が初期条件より決められたので、あとは順次、
$y(2\Delta x),\,y(3\Delta x),\,y(4\Delta x),\cdots$を計算するだけです。この方程
式の解は、まさに、$\sin x$です。リストのプログラムで確認してみましょう。

実は、常微分方程式を数値計算により解く、もっと強力な方法があります。4次のルンゲ・
クッタの方法です。初期条件もこちらのほうが、簡単に表現できます。ただし、2階の微
分方程式を解くときには、ほんのちょっとだけ、工夫が必要です。今の段階で、それを説
明するのは、早すぎます。ルンゲ・クッタの方法を学習するときに説明します。
%
%--------------------------
\subsubsection{プログラム例}
%--------------------------
\begin{list}{}{\leftmargin=6em \labelwidth=3em}
 \item[1行]解くべき微分方程式、$\diff^2y/\diff x^2=-y$を引数付き置換で定義している。
 \item[2行]微分方程式の初期条件、$y(0)=0$を引数付き置換で定義している。
 \item[3行]微分方程式の初期条件、$y^\prime(0)=1$を引数付き置換で定義している。
 \item[4行]計算ステップ数を文字型置換で設定している。プログラム中のNが全て1000
	   に置き換わる。
 \item[5行]計算の上限を文字型置換で設定している。プログラム中の\tw{MAX\_N}が全て、
	   10.0に置き換わる。
 \item[20-21行]$(x,\,y)$の値をファイルに書き込んでいる。
 \item[28行]微分方程式を計算している。
 \item[29行]計算結果$(x,\,y)$の値をファイルに書き込んでいる。
\end{list}
%
\begin{lstlisting}[caption=微分方程式を計算するプログラム,label=prog:diff_eq]
#define DY2DX2(x,y) (-y)                /* 微分方程式 */
#define Y0 0                            /* 初期条件 y=0      at x=0 */
#define DYDX0 1                         /* 初期条件 dy/dx=1  at x=0 */
#define N 1000                          /* 計算ステップ数 */
#define MAX_X 10.0                      /* 計算の上限 ここまで計算 */
#include <stdio.h>

int main() {
  double dx, x, y[N+1];
  int i; FILE *result;

  result = fopen("cal_result","w");

  /* ------ 計算条件の設定 -------- */

  dx = MAX_X/N;                         
  y[0] = Y0;
  y[1] = Y0+DYDX0*dx+1/2*DY2DX2(0,y[0])*dx*dx;

  fprintf(result,"%e\t%e\n",0.0,y[0]);
  fprintf(result,"%e\t%e\n",dx,y[1]);


  /* ------ 微分方程式の計算 --------- */
  
  for(i=2; i <= N; i++){
    x = i*dx;
    y[i]=2*y[i-1]+DY2DX2(x-dx,y[i-1])*dx*dx - y[i-2];
    fprintf(result,"%e\t%e\n",x,y[i]);
  }

  fclose(result);
  return 0;

}
\end{lstlisting}
%
%--------------------------
\subsubsection{実行と表示}
%--------------------------
作成したプログラムを実行すると、計算結果のファイル(\tw{cal\_result})が作成されま
す。これは、各行に$(x,y)$の値がかかれています。これを、gnuplotというグラフ作成の
プログラムを用いて、可視化しましょう。方法は、以下の通りです。
\begin{enumerate}
 \item ターミナルからgnuplotを起動させます。コマンドは、\command{gnuplot}です。
       \begin{quote}
	\setlength{\baselineskip}{12pt}
	\begin{verbatim}	
	[yamamoto]$ gnuplot
	\end{verbatim}
       \end{quote}
 \item gnuplotが起動したら、計算結果の描画です。コマンドは、plot "ファイル名"で
       す。コマンドを送ると、図\ref{fig:result_by_gnuplot}のようなグラフが書かれ
       るはずです。
       \begin{quote}
	\setlength{\baselineskip}{12pt}
	\begin{verbatim}	
       gnuplot> plot "cal_result"
	\end{verbatim}
       \end{quote}
\item gnuplotを終了するときにコマンドは、exitです。
       \begin{quote}
	\setlength{\baselineskip}{12pt}
	\begin{verbatim}	
       gnuplot> exit
	\end{verbatim}
       \end{quote}
\end{enumerate}
%
\begin{figure}[hbpt]
 \begin{center}
  \includegraphics[keepaspectratio, scale=0.7]{figure/sin.eps}
  \caption{gnuplotによる計算結果のグラフ}
  \label{fig:result_by_gnuplot}
 \end{center}
\end{figure}
%
\begin{quote}
 \begin{itemize}
  \item[\textbf{[練習1]}]リスト\ref{prog:diff_eq}のプログラムを実行せよ。
  \item[\textbf{[練習2]}]リスト\ref{prog:diff_eq}と式の関係を調べよ。
  \item[\textbf{[練習3]}]リスト\ref{prog:diff_eq}を書き換えて、他の微分方程式を
			 解いて見よ。
  \item[\textbf{[練習4]}]計算ステップ数と誤差の関係を調べよ。
 \end{itemize}
\end{quote}
%======= 参考文献 =====================================================
\bibliographystyle{junsrt}
\bibliography{reference} 
%=====================================================================
\end{document}
