積分方法の比較

FORTRANの環境作ったし、実際に計算させてみた
やっぱJavaより早い気がするｗ
今回は積分のプログラミング方法３パターンの比較を行って、収束具合、正確性を見ていく

　積分の計算種類

$$\begin{eqnarray} \int_a^b f(x) dx \end{eqnarray} =\lim_{n \to \infty} \sum^n_{i=0}f(x)\Delta x$$ 基本的に積分は上記式のように区分求積法で計算していくわけだけど、「$n \to \infty$」はプログラムでは実現できない。なので、様々な計算法で近似して数値解を求める必要がある。

1. 長方形近似
一番基本形。積分する領域を長方形に切り刻んで全部足す方法
・数式 $$\begin{eqnarray} \int_a^b f(x) dx \end{eqnarray} =\sum^n_{i=a}f(i)\Delta x$$ 刻み幅$\Delta x$はループ数$N$を使い、$\Delta x = (b - a) / N$となる
$x$座標の始点、終点、ループ回数を引数にした関数だとこんな感じ

	FUNCTION SEKIBUN1(UP, BOTTOM, N)
	F(x) = 4 / (1 + x ** 2)
	real SUM, BOTTOM, UP, dx
	INTEGER N
	SUM = 0.0
	dx = (UP - BOTTOM) / N
	do 10 i = 0, N
	x = dx * i + BOTTOM
	SUM = SUM + F(x) * dx
	10 continue
	SEKIBUN1 = SUM
	RETURN
	END

view raw gistfile1.f90 hosted with ❤ by GitHub
2. 台形近似
長方形近似の上位互換。積分する領域を台形に刻むことで、勾配を配慮した形
・数式 $$\begin{eqnarray} \int_a^b f(x) dx \end{eqnarray} =\sum^n_{i=a}\frac{f(i)+f(i+\Delta x)}{2}\Delta x$$ 刻み幅$\Delta x$はループ数$N$を使い、$\Delta x = (b - a) / N$となる
$x$座標の始点、終点、ループ回数を引数にした関数だとこんな感じ

	FUNCTION SEKIBUN2(UP,BOTTOM,N)
	F(x) = 4 / (1 + x**2)
	real SUM, BOTTOM, UP, dx
	INTEGER N
	SUM = 0.0
	dx = (UP - BOTTOM)/ N
	do 11 i = 0, N
	x = dx * i + BOTTOM
	SUM = SUM + (F(x) + F(x + dx))* 0.5 * dx
	11 continue
	SEKIBUN2 = SUM
	RETURN
	END

view raw gistfile1.f90 hosted with ❤ by GitHub
3. シンプソンの公式
f(x) を二次関数 P(x) で近似する数値解析方法
・数式 $$\begin{eqnarray} \int_a^b f(x) dx \end{eqnarray} =\frac{h}{3}\Bigl\{f(a)+4\sum^m_{i=1}f(x_{2i-1})+2\sum^{m-1}_{i=1}f(x_{2i})+f(b)\Bigr\}$$ $h=(b-a)/(2m)$　$x_i=a+ih$　$(i=1,2,3...2m-1)$
$x$座標の始点、終点、ループ回数を引数にした関数だとこんな感じ


	FUNCTION SEKIBUN3(UP,BOTTOM,N)
	F(x) = 4 / (1 + x**2)
	real SUM, BOTTOM, UP, dx
	INTEGER N
	dx = (UP - BOTTOM)/ N
	sum = F(BOTTOM) + F(top)
	do i = 1, N
	if(MOD(i, 2).EQ.0) then
	sum = sum + 4 * F(BOTTOM + i * dx)
	else
	sum = sum + 2 * F(BOTTOM + i * dx)
	end if
	end do
	sum = sum * dx / 3.0
	SEKIBUN3 = SUM
	RETURN

view raw gistfile1.f90 hosted with ❤ by GitHub

　積分の比較

次に３つの積分方法の正確性、収束性、計算スピードについてみていく

1. 正確性
下記式の左辺の数値解をそれぞれを求め、$\pi$にどれだけ近づくか見てみる $$\begin{eqnarray} \int_0^1 \frac{4}{1+x^2} dx \end{eqnarray} =\pi$$


計算方法	解析解	正確性[%]
長方形近似	3.17157602	100.9544002
台形近似	3.16147637	100.6329184
シンプソンの公式	3.15489316	100.4233683

刻み回数100回の場合の積分結果 どの結果も同じオーダーなので、違いが分かりにくい。。。


2. 収束性
解の収束度合はどうだろうか。
x軸に刻み回数、y軸に数値解とし、刻み回数と解の収束をグラフ化した
sekibun1→長方形近似、sekibun2→台形近似、sekibun3→シンプソンの公式、pi→$\pi$

収束度合の数値化がちょっとできてないけど、シンプソンの公式がかなり早い段階で収束することがわかる。


3. 計算スピード
最後に計算速度の差について
この程度なら差を感じないけど、計算に一週間かかる場合とかだと重要になってくる。それぞれのアルゴリズムの計算回数からだいたいの計算時間のオーダーを求める
$f(x)$の計算ステップが$F$、分割数$N$とすると

計算方法	オーダー
長方形近似	$N(4+F)$
台形近似	$N(6+2F)$
シンプソンの公式	$N(6+F)$



つまり、長方形近似＞シンプソンの公式＞台形近似の順で早いことになる。
特に$f(x)$の計算ステップが異常にかかる系だと、台形近似が大変

　結論

シンプソンの公式が収束率、計算速度が優れていてよさげ
でも、実装が一番楽な長方形近似で分割数を膨大に増やして計算するのが、実装ミス、手間が省けて効率いいかも思った
以上

WINDOWS7でFORTRAN77の開発環境を作る

やっぱ数式とかプログラムしていく上で必要になってきたので、
FORTRAN77の開発環境を作っていく。

　FORTRAN77とは

科学技術計算に向いた手続き型プログラミング言語。その長い歴史の間に開発された非常に多くの数学関数やサブルーチンを数値解析ソフトウェアの形で持っている

by Wiki

まぁ簡単言えば計算特化型のプログラム言語。むしろ計算しかできない

　環境構築

1. ダウンロード
http://kkourakis.tripod.com/g77.htm
ココから下記をダウンロード
• g77exe.zip
• g77lib.zip
この二つが「bin」と「lib」。あと一つの「g77doc.zip」ってのは説明書だからいらない


2. 環境変数の設定
ダウンロードしたzipを解凍する。
「g77exe.zip」の中に「bin」ファイル
「g77lib.zip」の中に「lib」ファイル
があるので、それぞれのパスを環境変数に登録
環境変数の登録はwin7の場合
「スタートボタン」⇒「コンピュータ」を右クリック「プロパティ」⇒「システムの詳細設定」⇒「環境変数」⇒「システムの環境変数」をスクロールして、「Path」編集
最後尾にbinとlibのパスを追記。

例）
Cドライブ直下に「g77」フォルダを作り、その中に「bin」フォルダ、「lib」フォルダを格納した場合
を環境変数の最後に追記する。
「;」はパスの切れ目を表すので、書き忘れ注意
環境変数の画面を全部「OK」終わらせると、設定が反映される

　コンパイルと実行

サンプルで下記のファイル「sample.f」を作る コマンドプロンプトで下記を実行して実行できる確認

　手軽なエディタ

FORTRANのIF文とかに色がついてくれる便利なエディタがあれば開発が楽になるので紹介
「Notepad++ 5.8.5」が簡単だった
インストールは下記のURL
http://filehippo.com/jp/download_notepad/8806/

• 開発言語をFORTRANに設定する
タブの「言語」⇒「F」⇒「Frotran」


• 日本語記述で文字化けしない様にする
「フォーマット」⇒「UTF-8(BOMなし)エンコード」

完成！！