WANtaroHP (F90 MATH)

ページ概要

Contents
連立一次方程式
固有値解析
Gauss-Legendre 積分
高速フーリエ変換・逆変換
複素数の計算例
ベッセル関数，ガンマ関数 (Fortran 2008)
3次スプライン補間
代数方程式の解法 (DKA 法)
完全楕円積分

連立一次方程式

プログラム概要

連立一次方程式を解くプログラムです．
Gauss-Jordan の掃き出し法，バンドマトリックス用のCholesky 法は，ここに含まれています．これらは，ＦＥＭプログラムのサブルーチンで用いられています．

ソースコードと関連文書

Filename	Description
f90_GJ.txt	プログラムソースコード

固有値解析

プログラム概要

固有値解析を行うプログラムです．
Jacobi 法，一般化 Jacobi 法が含まれています．
Jacobi 法は，主成分分析のサブルーチンとして用いられています．
一般化 Jacobi 法は，FEMの座屈解析および固有振動数解析のサブルーチンとして用いられています．

ソースコードと関連文書

Filename	Description
f90_JACOBI.txt	プログラムソースコード

Gauss-Legendre 積分

プログラム概要

ここでは2つのプログラムを紹介しています．
ひとつは，ガウス積分点の座標と重みを計算するものです．
もうひとつは，Gauss-Legendre 積分のサンプルプログラムです．
この文書 'http://www7.ocn.ne.jp/~kawa1/numeric.pdf' をプログラム作成において参照しました．

ソースコードと関連文書

Filename	Description
f90_GAUSSIP.txt	ガウス積分点座標と重み計算プログラムソースコード
f90_CALPI.txt	Gauss-Legendre 積分サンプルプログラムソースコード

高速フーリエ変換・逆変換

プログラム概要

高速フーリエ変換・逆変換を行うプログラムです．
コノプログラムではFortranの複素数変数は用いられていません．実数部と虚数部の配列をそれぞれ準備して計算を行っています．
このプログラムはFourierスペクトルの計算プログラムなどのサブルーチンとして用いられています．

ソースコードと関連文書

Filename	Description
f90_FFT.txt	プログラムソースコード
inp_FFT_0.txt	入力データサンプル
out_FFT_0.txt	出力データサンプル

複素数の計算例

プログラム概要

このプログラムは，複素数を用いたいくつかの計算例を示すものです．

ソースコードと関連文書

Filename	Description
f90_COMPLEX.txt	プログラムソースコード

ベッセル関数，ガンマ関数 (Fortran 2008)

プログラム概要

このプログラムは，Fortran2008に含まれる組み込み関数であるベッセル関数とガンマ関数の使用事例を示すものです．
gnuplot によるpngファイルも示しています．このグラフの中にはgnuplotの関数でプロットしたものも参考のために含ませてあります．
実行用バッチコマンドは以下のとおりです．

gfortran -o f08_BESSEL.exe f08_BESSEL.f08
gfortran -o f08_GAMMA.exe f08_GAMMA.f08
gfortran -o f08_GPL_PLOT.exe f08_GPL_PLOT.f08

f08_BESSEL > test_b.txt
f08_GAMMA > test_g.txt
f08_GPL_PLOT test_b.txt test_g.txt

ソースコードと関連文書

Filename	Description
f08_BESSEL.txt	ベッセル関数使用例プログラム
f08_GAMMA.txt	ガンマ関数使用例プログラム
f08_GPL_PLOT.txt	gnuplotによる作図用プログラム
fig_png_bessel.png	ベッセル関数画像出力
fig_png_gamma.png	ガンマ関数画像出力

3次スプライン補間

プログラム概要

このプログラムは，3次スプライン補間を行うものです．
連立方程式の解法には，バンドマトリックス用のコレスキー法が用いられています．
結果を確認するための gnuplotのスクリプトも付けてみました．

実行用バッチコマンド

gfortran -o f90_SPL.exe f90_SPL.f90
f90_SPL inp_data.txt 5 > out_data.txt

コマンドライン引数の並びは以下のとおりです．

f90_SPL fnameR md > fnameW

f90_SPL	コンパイルされたF90プログラム
fnameR	入力ファイル名
md	あるプロットと次のプロット間の分割数
fnameW	出力ファイル名

ソースコードと関連文書

Filename	Description
f90_SPL.txt	F90プログラムソース
inp_data.txt	入力データサンプル
out_data.txt	出力データサンプル
gpl_spline.txt	gnuplot 用スクリプト
fig_png_spline.png	出力画像サンプル

代数方程式の解法 (DKA 法)

プログラム概要

このプログラムは，代数方程式の解を，複素数解も含めて全て求めることができます．
代数方程式の解法には，DKA法が用いられています．
計算には，4倍精度実数が用いられており，その有効桁数は約30桁です．
重解を含む場合は，精度が落ちます．
このプログラムでは複素数タイプの変数は用いられていません．実数部と虚数部は別々に扱われています．
このFortranプログラムは，以下の文書に示されているCプログラムを参考に作られました．

http://www.ecs.shimane-u.ac.jp/~kyoshida/c6(2002).pdf

実行用バッチコマンド

gfortran -o f90_DKA.exe f90_DKA.f90

f90_DKA inp_test10.txt > out_test10.txt
copy _data_circ.txt _data_circ10.txt
copy _gmt_inp.txt _gmt_inp10.txt

f90_DKA inp_test20.txt > out_test20.txt
copy _data_circ.txt _data_circ20.txt
copy _gmt_inp.txt _gmt_inp20.txt

call bat_gmt_draw

コマンドライン引数の並びは以下のとおりです．

f90_DKA fnameR > fnameW

f90_DKA	コンパイルされたF90プログラム
fnameR	入力ファイル名
fnameW	出力ファイル名

出力ファイルfnameWは結果の数値を確認するためのものです．この他にも収束過程をGMTで描画・確認するため，固定ファイル名で以下の出力を行っています．

_iteration.txt	作業用ファイル
_gmt_inp.txt	GMTで読み込むためのデータファイル(各係数の収束過程を示す座標データ)
_data_circ.txt	GMTで読み込むためのデータファイル(Aberthの初期値の分布域を示す円の座標ザータ)

ソースコードと関連文書

Filename	Description
f90_DKA.txt	Program for DKA method
inp_test10.txt	入力データサンプル (Case-1)
inp_test20.txt	入力データサンプル (Case-2)
inp_test30.txt	入力データサンプル (Case-3)
out_test10.txt	出力データサンプル (Case-1)
out_test20.txt	出力データサンプル (Case-2)
out_test30.txt	出力データサンプル (Case-3)
bat_gmt_draw.txt	GMT作図用バッチファイル
fig_gmt_out10_1.png	画像出力事例 (Case-1)
fig_gmt_out10_2.png	画像出力事例 (Case-1 detail)
fig_gmt_out20_1.png	画像出力事例 (Case-2)
fig_gmt_out20_2.png	画像出力事例 (Case-2 detail)
fig_gmt_out30_1.png	画像出力事例 (Case-3)
fig_gmt_out30_2.png	画像出力事例e (Case-3 detail)
TeX_DKA.pdf	説明文書

Case-1	f(x)=(x+1)^10=0
Case-2	f(x)=(x-20)(x-19)...(x-1)=0
Case-3	f(x)=(x-30)(x-29)...(x-1)=0

完全楕円積分

プログラム概要

このプログラムは，第一種完全楕円積分 K(p) および第二種完全楕円積分 E(p) の値を求めるためのものです．
数値積分法には，級数展開による方法と，Gauss-Legendre 積分による方法が用いられています．
Gauss-Legendre 積分による方法のほうが級数展開による方法より圧倒的に速いです．
計算事例を以下に示します．下表中のita. は，収束に要した繰り返し回数を示します．

級数展開法による結果

  Mehod by Series expansion
  ita.   p              K(p)           E(p)
    1  0.0000000E+00  0.1570796E+01  0.1570796E+01
    5  0.1000000E+00  0.1574746E+01  0.1566862E+01
    7  0.2000000E+00  0.1586868E+01  0.1554969E+01
    9  0.3000000E+00  0.1608049E+01  0.1534833E+01
   11  0.4000000E+00  0.1640000E+01  0.1505942E+01
   14  0.5000000E+00  0.1685750E+01  0.1467462E+01
   19  0.6000000E+00  0.1750754E+01  0.1418083E+01
   27  0.7000000E+00  0.1845694E+01  0.1355661E+01
   41  0.8000000E+00  0.1995303E+01  0.1276350E+01
   83  0.9000000E+00  0.2280549E+01  0.1171697E+01
54842  0.9999000E+00  0.5645147E+01  0.1000515E+01
 time=   47.875000     (sec)

Gauss-Legendre 積分による結果 (積分点 = 10)

  method by Gauss-Legendre rule (n=10)
  ita.   p              K(p)           E(p)
    2  0.0000000E+00  0.1570796E+01  0.1570796E+01
    2  0.1000000E+00  0.1574746E+01  0.1566862E+01
    2  0.2000000E+00  0.1586868E+01  0.1554969E+01
    2  0.3000000E+00  0.1608049E+01  0.1534833E+01
    2  0.4000000E+00  0.1640000E+01  0.1505942E+01
    2  0.5000000E+00  0.1685750E+01  0.1467462E+01
    2  0.6000000E+00  0.1750754E+01  0.1418083E+01
    2  0.7000000E+00  0.1845694E+01  0.1355661E+01
    3  0.8000000E+00  0.1995303E+01  0.1276350E+01
    3  0.9000000E+00  0.2280549E+01  0.1171697E+01
   53  0.9999000E+00  0.5645148E+01  0.1000515E+01
 time=   0.0000000     (sec)

Gauss-Legendre 積分による結果 (積分点 = 2)

  method by Gauss-Legendre rule (n=2)
  ita.   p              K(p)           E(p)
    2  0.0000000E+00  0.1570796E+01  0.1570796E+01
    3  0.1000000E+00  0.1574746E+01  0.1566862E+01
    4  0.2000000E+00  0.1586868E+01  0.1554969E+01
    4  0.3000000E+00  0.1608049E+01  0.1534833E+01
    5  0.4000000E+00  0.1640000E+01  0.1505942E+01
    6  0.5000000E+00  0.1685750E+01  0.1467462E+01
    6  0.6000000E+00  0.1750754E+01  0.1418083E+01
    7  0.7000000E+00  0.1845694E+01  0.1355661E+01
    9  0.8000000E+00  0.1995303E+01  0.1276350E+01
   13  0.9000000E+00  0.2280549E+01  0.1171697E+01
  323  0.9999000E+00  0.5645148E+01  0.1000515E+01
 time=  0.12500000     (sec)

ソースコード

Filename	Description
f90_CEIS.txt	級数展開法
f90_CEIG.txt	Gauss-Legendre 積分