代数方程式の解法 (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 | コンパイルされたF90プログラム |
| fnameR | 入力ファイル名 |
| fnameW | 出力ファイル名 |
出力ファイルfnameWは結果の数値を確認するためのものです．
この他にも収束過程をGMTで描画・確認するため，固定ファイル名で以下の出力を行っています．
| _iteration.txt | 作業用ファイル |
| _gmt_inp.txt | GMTで読み込むためのデータファイル(各係数の収束過程を示す座標データ) |
| _data_circ.txt | GMTで読み込むためのデータファイル(Aberthの初期値の分布域を示す円の座標ザータ) |
ソースコードと関連文書
| 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)
ソースコード
