「プログラミングと数値計算」カテゴリーアーカイブ

プログラミングと数値計算, 物理学

非線形Schrödinger方程式のソリトン解

コメントする

非線形シュレディンガー方程式
\(
\displaystyle i\frac{\partial \Psi}{\partial t}=-\frac{1}{2}\frac{\partial^2 \Psi}{\partial x^2}-g|\Psi|^2\Psi
\)
にはある解析解が存在します。それがソリトン(soliton)解と呼ばれるもので,上式のソリトン解は
\(
\displaystyle \Psi(x,t)=\sqrt{\Omega}\ {\rm sech}\left\{\sqrt{\Omega}\left(x\sqrt{g}-\frac{gV}{2}t\right)\right\}\cdot \exp\left\{i\frac{V\sqrt{g}}{2}x-i\frac{g}{2}\left(\frac{V^2}{4}-\Omega\right)t\right\}
\)
です。(\(g>0\),\(\Omega\):ソリトンの振幅、\(V\):ソリトンの速度に関するパラメータ、ソリトン自体の速度は\(V\sqrt{g}/2\))

[adsense1]

ソリトンの歴史的背景

「非線形」とは重ね合わせの原理が成り立たない系です。

1844年、スコットランドのJ.Scott-Russellによって孤立した波(solitary wave)を観測した事が報告されました J.Scott-Russellによる報告”Report on Waves”(リンク先のSR44.pdf, 16.3MB))。
当時の認識では、波は波動方程式で記述され、その波の速度\(v\)は\(v=f\lambda\)の元、一定である。だからパルス状の波は異なる波長の波の重ね合わせで書けているはずで、時間と共に分散していくはず。なのになぜ時間が経過しても孤立した波が存在できるのか?という事で大きな論争となりました。

60年後の1895年、オランダのKortewegとde Vriesによって”浅い水の波”を記述する非線形偏微分方程式(KdV方程式)が提出され、この方程式の特解として孤立波が与えられました。
孤立波は、

  1. 空間的に局在した波が、その性質(速さや形)を変えずに伝搬する
  2. 孤立波は互いの衝突に対して安定であり、各々の個別性を保持する

という性質を持つ非線形波動と定義されます[1]。
2番目の、粒子のような性質を持つことから、solitary に接頭語-on をつけ、soliton(ソリトン)と名づけられました。

その後、1981年に佐藤幹夫がソリトンの統一理論(佐藤理論やKP理論)を発表しました。
これによりソリトン方程式(ソリトンを記述し，かつ厳密に解ける方程式)に決着が付きました。
ソリトン方程式は非線形なのに厳密に解ける、可積分系である。

ソリトン方程式を解く方法は([4]を引用しますが)

上で指摘したように，logistic方程式が解けるからくりとソリトン方程式が解けるからくりはよく似ています．違いは，logistic方程式が変数変換一発で線形常微分方程式になってしまったのに対し，ソリトン方程式の場合は変数変換で双線形形式になり，双線形形式の解として行列式が現れ，行列式の中身に簡単な線形方程式が現れるというところです．しかし，離散化で保存するべき構造は明らかです．まず，解の中身の線形方程式を離散化し，行列式の構造をそのまま使って双線形形式を作る．最後に変数変換して非線形のレベルに戻ればよい．

となるそうです。

また、ソリトン方程式の特徴である、無限個の対称性(無限個の保存量)は、Gardner変換という変換をすることで証明できるそうです[5]。
これ以上はこの分野の専門家ではないので話せません。

ちなみに津波もソリトンの一つとみなせます。

ソリトン解が生まれるイメージ

なぜソリトン解が生まれるのでしょうか。
今、孤立した波(空間的に凸)を考えます。この時、

エネルギー的に安定になろうとして密度を均一にするために広がろうとする効果
と
粒子間を結び付ける引力相互作用(例えば水面だったら水と水との分子間力等)のため集まろうとする効果

のつり合いによって、丁度均衡が保たれるとき、このソリトン解が生まれます。

･･･実は、ソリトン解には2種類あります。
それは明るい(Bright)ソリトン解と暗い(Dark)ソリトン解です。
今まで話していたのは全て明るいソリトン解です。
暗いソリトン解とはどういったものでしょう。
暗いソリトン解とは、ある部分が空間的に凹んでいる、孤立した解です。

エネルギー的に安定になろうとして密度を均一にするためにその凹みを埋めようとする効果
と
粒子間の斥力相互作用のために粒子間を避けようとする効果

のつり合いによって、丁度均衡が保たれるとき、この暗いソリトン解が生まれます。
暗いソリトン解が生まれるのは斥力相互作用の時で、斥力相互作用を持つ系というのは、調べた限りでは量子力学のボーズ・アインシュタイン凝縮体で、暗いソリトンは渦ソリトンという形で現れるそうです。これ以上の具体例は分かりませんでした。もしも具体例を知っているという方は教えていただければ幸いです。
暗いソリトンの解析解は参考文献[1]の本に紹介されているので、それをご参考にしてください。

非線形シュレディンガー方程式におけるソリトン解

では本題の、非線形シュレディンガー方程式における(明るい)ソリトン解を考えましょう。
下の形の非線形シュレディンガー方程式を考えます。
\(
\displaystyle i\frac{\partial \Psi}{\partial t}=-\frac{1}{2}\frac{\partial^2 \Psi}{\partial x^2}-g|\Psi|^2\Psi
\)
ここで、\(\Psi=\Psi(x,t)\)で、\(g\)は正の値で相互作用の強さ(この場合、引力相互作用)を表します。

この非線形シュレディンガー方程式のソリトン解\(\Psi(x,t)\)は、
\(
\displaystyle \Psi(x,t)=\sqrt{\Omega}\ {\rm sech}\left\{\sqrt{\Omega}\left(x\sqrt{g}-\frac{gV}{2}t\right)\right\}\cdot \exp\left\{i\frac{V\sqrt{g}}{2}x-i\frac{g}{2}\left(\frac{V^2}{4}-\Omega\right)t\right\}
\)
であり、\(\Omega\)はソリトンの振幅の大きさ、\(V\)はソリトンの速度を決めるパラメータを表します。ソリトン自体の速度は\(V\sqrt{g}/2\)となります([2]を参考)。
また、\({\rm sech}(x)\)は双曲線関数の一種(双曲線正割と呼ばれる)であり、
\(
\displaystyle {\rm sech}(x)=\frac{1}{\rm cosh(x)}=\frac{2}{e^{x}+e^{-x}}
\)
を表します。

解析解のプロット

解析解をプロットします。gnuplotコードは下のほうに載せておきます。
\(g=2, V=1, \Omega=1\)とすると、以下の振る舞いが観測されます。
ここで紫はソリトン解の実部、緑は虚部、青は絶対値2乗を表します。
動画は1枚当たり原子単位系で0.1秒、合計で10秒間のシミュレーションです。
また、このソリトンの速度は\(V\sqrt{g}/2\sim 0.7071\)です。
soliton1

[adsense2]

gnuplotコード

非線形シュレディンガー方程式のソリトン解(解析解)を出力します。
gnuplot上で下のスクリプトを実行してください。
(ただし、gnuplot ver4.6以降に限ります。)

omega=1e0
V=1e0
g=2e0
x0=-5e0 # initial position

set xr[-10:10]
set yr[-1.5:1.5]
set samples 1000
set xl "x[a.u.]"

sech(x)=2e0/(exp(x)+exp(-x))
amp(x,t)=sqrt(omega)*sech(sqrt(omega)*((x-x0)*sqrt(g)-g*V*t*0.5e0))
phase(x,t)=V*sqrt(g)*x*0.5e0-g*t*0.5e0*(V*V*0.25e0-omega)
soliton(x,t)=amp(x,t)*exp({0e0,1e0}*phase(x,t))


#set term gif animate delay 10 optimize size 960,720
#set output 'movie.gif'
do for[i=0:100:2]{
   t=i*0.1e0
   plot abs(soliton(x,t))**2 lw 3 lc 1 lt 1 ti sprintf("|\psi|^2, t=%2.1f",t),\
        real(soliton(x,t)) lw 3 lc 2 lt 2 ti "Real",\
        imag(soliton(x,t)) lw 3 lc 3 lt 3 ti "imag"
}
#set out
#set terminal wxt enhanced

もっとソリトンについて知りたい方はまず参考文献[3]を読むことをお勧めします。
その後、[4]を読み、[5]を読み、[1]の本を読むのが良いと思われます。
[3],[4]は簡単な表現を用いてソリトンとその後の発展について記述されています。

参考文献

[1]和達三樹著 『非線形波動 (現代物理学叢書) 』岩波書店 (2000年) p.7
[2]和達三樹著 『非線形波動 (現代物理学叢書) 』岩波書店 (2000年) p.29

[3]ソリトンの数学 – Researchmap
[4]ソリトン ～ 不思議な波が運んできた，古くて新しい数学の物語 ～
[5]〔連載〕非線形波動―ソリトンを中心として―第５章　逆散乱法

↑画像のフォントはキユマヤ園様による数式フォント -びゅんびゅん→SSSです！

プログラミングと数値計算, 数学

離散フーリエ変換と高速フーリエ変換(fortran90)

2件のコメント

こちらのページは後程消去いたします。

以下のページに統合しますので、ご参照ください。
https://slpr.sakura.ne.jp/qp/dft-for-numerical-calculation/

本ページにはミスもあり、上の方が正しいですのでご参照ください。

数値計算で離散フーリエ変換を行う方法です。
言語は
fortran90
で
intelのMKL(マス・カーネル・ライブラリー)
を用いて離散フーリエ変換を行いたいと考えます。

高速フーリエ変換は離散フーリエ変換を高速に行う手法のことです。

通常の離散フーリエの計算量のオーダーは、\(O(N^2)\)であり、
基数2の高速フーリエの計算量のオーダーは、\(O(Nlog_2 N)\)です。

基数2は良く言われる通常の配列の個数\(2^m\)個のことを表しています。

離散フーリエ変換を行う際には周期境界条件が課されます。
分点の端で周期境界条件が満たされていない場合、それはその点でステップ関数が含まれることを意味します。
依って予期しない高周波成分が出現するので、注意してください。
一番良い判定方法は関数の対数微分が一致するかどうかを判定するのが良いでしょう。

MKLの離散フーリエ変換ルーチンの概要1

MKLの離散フーリエ変換ルーチンは非常に優秀です。
渡す配列サイズによってプログラムが自動的に判断し、最適な手法で離散フーリエ変換を行います。
MKLのマニュアルによると、

注 : DFT関数は任意の長さをサポートしている。
これらのルーチンは、基数 2 だけでなく、基数 3、5、7、11 に対しても高
性能で広範な機能性を提供する。対応する基数の一覧は、使用しているラ
イブラリー・バージョンの「インテルMKLテクニカル・ユーザー・ノー
ト」を参照のこと。

とあります。
高速フーリエ変換の考え方は\(2^m\)の時だけでなく、\(3^m, 5^m, 7^m, 11^m, \cdots\)の時にも同様に考えることができます。
この時、何の整数乗にするかを区別するために “基数” という言葉を使います。
プログラムを作る場合、通常は一番シンプルな基数”2″が選ばれます。

具体的に配列サイズ\(N\)が
\(N= 2^m\)、ただしmは正の整数(配列をAとするとサイズ\(A(0:2^m-1)\))の時に基数2の高速フーリエ変換が行われます。
それ以外の時は通常の離散フーリエ変換が行われているはずです。

余り自信が無いのではっきりとしたことは言えませんが、分点量(2~20000程度まで)に対するMKLのフーリエ変換の計算速度はほとんどわからないくらいでした。特定の基数の時だけ早いんだろうと予想していましたが、そんなことはなく、それだけ優秀であるという事でしょう。

もう6年前に書かれている記事ですので今はどうなっているかわかりませんが、FFTWと計算速度を比較したページがありましたので載せておきます。
一言でいえば、MKLによる離散フーリエ変換はFFTWよりも早いよ、ということです。
インテル謹製の数値演算ライブラリ「MKL」を使ってプログラムを高速化 3ページ -OSDN, (2009)

MKLの離散フーリエ変換ルーチンの概要2

MKLでの1次元、長さnの順方向離散フーリエ変換は、
\(
\displaystyle z_k=\sum_{j=0}^{n-1}w_je^{-i2\pi jk/n}
\)
に従い行われます(MKLリファレンスマニュアル、11-26)。
数値計算の世界では、余分な係数\(\frac{1}{2\pi}\)をなくすため、指数関数の肩に\(2\pi\)が掛けられたものが良く使用されます。

※ここでは、フーリエ変換前の空間を”位置空間“と呼び、フーリエ変換後の空間を”波数空間“と呼ぶことにします。
位置⇔波数ではなく、時間⇔周波数としても同じことなのでどちらでもいいのですが、ここでは、位置⇔波数の呼び方に統一いたします。

DFTexampleNeven
上図のように位置空間上で区間\(L\)の領域で\(N\)個に分割します。
上図はNが偶数の場合です。
MKLの離散フーリエ変換ルーチンはいかなるNに対してもフーリエ変換が可能ですが、奇数の場合、以下のように波数成分の格納が変わります。
周波数格納順

\(
f(x)=\sin(\alpha x)
\)
を離散フーリエ変換した場合、
波数空間でのフーリエ変換後のピークは\(k=\frac{\alpha}{2\pi}\)に現れます。
この時、f(x)のフーリエ変換後の関数\(g(k)=F_{[f(x)]}\)は複素数として格納され、上記の場合では
\(
\displaystyle g(k)=i\left\{\frac{N}{2}\delta(k-\frac{\alpha}{2\pi})+\frac{N}{2}\delta(k+\frac{\alpha}{2\pi})\right\}
\)
となります。ここで\(i\)は虚数単位を表し、\(\delta(k)\)はデルタ関数を表します。
波数空間での最大の周波数値\(f_c\)はナイキストの定理より、
\(
\begin{align}
f_c = \frac{1}{2\Delta}, \Delta=\frac{L}{N}
\end{align}
\)
と表されます。これは、波数空間では\(\pm f_c\)までサンプリング可能であることを示しています。
フーリエ変換後の値\(g(k)\)の分点上の和\(S\)、すなわち
\(
\displaystyle S=\sum_{j=0}^{N-1}g(k_j)
\)
は離散フーリエ変換の性質により、\(S\)は分点の総数\(N\)に等しくなります。

逆離散フーリエ変換によって元の空間に戻る場合、規格化因子\(N\)で割らなければなりません。

MKLの1次元離散フーリエ変換ルーチンのプログラム例

使い方はFFT Code Examples -Intel®Developer Zone を参考としました。

例として
関数
\(
f(x)=\sin(4x)
\)
をフーリエ変換し、元に戻すプログラムを書きましょう。
プログラムの流れとしては、

関数\(f(x)\)の位置空間での書き出し(ファイル”before.d”に書き出し、横軸:位置\(x\), 縦軸\(f(x)\))
↓ 順方向フーリエ変換(規格化因子で割られない)
関数\(g(k)=F_{[f(x)]}\)の書き出し(ファイル”k.d”に書き出し、横軸:\(k/(2\pi)\), 縦軸\(g(k)\))
↓ 逆方向フーリエ変換(規格化因子Nで割られる)
関数\(f(x)=F^{-1}_{[g(k)]}\)の書き出し(ファイル”after.d”に書き出し、横軸:x, 縦軸f(x))
になっています。

コンパイルを行うにあたり重要な事が2点あります。
1, MKLを使う事
2, fftサブルーチンの上にinclude文を付け加える事(かも？)

このincludeは、MKLにある離散フーリエ変換ルーチンを使いますよーというサインであり、パスはMKLのバージョンや、インストールした状況によって変わります。
うまく、ファイル”mkl_dfti.f90″を見つけてください。パスを上から辿っていくのがいいと思います。

※2016/02/27時点での最新バージョン(ifort –version で確認できます)
ifort (IFORT) 16.0.2 20160204
において、mkl_dfti.fのデフォルトでの場所は、
“/opt/intel/compilers_and_libraries_2016/linux/mkl/include/mkl_dfti.f90”
です。また、このバージョンだけかどうかわかりませんが、include文が無くてもどうやら動くようです。

プログラムのコンパイルと実行は

$ ifort -mkl main.f90
$ ./a.out

で良いでしょう。

サブルーチンのコンセプトは、
ただ順/逆方向離散フーリエ変換をしてくれるサブルーチンです。速度は追求していません。

速度を追及する書き方をしているのは下のベンチマークの項目です。プログラムの可読性が悪くなるので、中身を知らないと調整が出来ません。
ベンチマークに書いてあるように、サブルーチンdftの中身の

  Status = DftiCreateDescriptor(hand,DFTI_DOUBLE,DFTI_COMPLEX,1,n)
  Status = DftiCommitDescriptor(hand)
  Status = DftiFreeDescriptor(hand)

はFFTを使うための準備で、これは配列の大きさが変わらない限り1度実行するだけで良いのです。
これを排除すると3-4倍ほど速くなります。

上で用いているサブルーチン
dft,dfts,dftfはそれぞれフーリエ変換、フーリエ変換後の空間の値、ソート用のルーチンであることを意味し、その中身は以下のようなものです。例のごとく、MKLのルーチンはそのままでは使いにくいので1クッション挟みます。

ここ↑に書いてあるサブルーチン

  1. dft(N,z,FB)
  2. dftf(N,f,h)
  3. dfts(N,f,z,mf,mz)

を説明します。

  1. dft(N,z,FB)

    このルーチンによって順/逆方向フーリエ変換が行われます。
    引数はそれぞれ

    • (入力) N : integer, データ配列のサイズ
    • (入出力) z : complex(kind(0d0)), データ配列z(0:N-1)
    • (入力) FB : character, 順方向(“forward”)、逆方向(“backward”)の指定

    です。データ配列に上書きして値が返ってきます。
    サイズNのデータ配列z,順方向の離散フーリエ変換をしたい場合は、

    call dft(size(z,1),z,"forward")

    とするのがいいでしょう。
    規格化は順方向の時には行わず、逆方向の時にだけ\(1/N\)を掛けています。

  2. dftf(N,f,h)

    このルーチンは波数空間での横軸の値f(0:N-1)を与えます。
    引数はそれぞれ

    • (入力) N : integer, データ配列のサイズ
    • (出力) f : double precision, 波数空間の横軸の値f(0:N-1)
    • (入力) h : double precision, 位置空間でのサンプルレート(刻み幅)

    です。
    波数fの順番は少しややこしく、
    周波数格納順
    となっています。
    また、ここで出力されるfは波数を\(2\pi\)で割ったものが出力されます。
    具体的には、\(\sin(x)\)をMKLのフーリエ変換を行いグラフを書くとフーリエ変換後の空間(波数空間)で、ピークの位置は\(\pm 0.159115494309\cdots (\sim \frac{1}{2\pi})\)に表れることを意味します。
    また、変換先の空間をグラフで描画したい時に値を小さい順に並べ替えたい必要が出てきたときは並べ替え用のルーチンdftsを使用してください。

  3. dfts(N,f,z,mf,mz)

    このルーチンは変換先の波数空間でグラフを描画したい時に、順番を値を小さい順に並べ替えます。
    引数はそれぞれ

    • (入力) N : integer, データ配列のサイズ
    • (入力) f : double precision, 波数空間の横軸の値f(0:N-1)
    • (入力) z : complex(kind(0d0)), 波数空間でのデータ配列z(0:N-1)
    • (出力) mf : double precision, 並び替えられた波数空間の横軸の値mf(0:N-1)
    • (出力) mz : complex(kind(0d0)), 並び替えられた波数空間でのデータ配列z(0:N-1)

    です。上書きして引数を返す必要性がないだろう、と判断して別の配列に書き出させています。
    もしもメモリが足らないとか、気に食わない人は書き換えてください。

[adsense1]

出力ファイル”before.d”, “k.d”, “after.d”を書いてみます。
位置空間での区間\([-3\pi\sim3\pi]\)で、分点数\(N=288\)で行った場合、
位置空間での波形、この場合は\(\sin(4x)\)が描かれます。
位置空間

波数空間では
\(
\displaystyle g(k)=i\left\{\frac{288}{2}\delta(k-\frac{4}{2\pi})+\frac{288}{2}\delta(k+\frac{4}{2\pi})\right\}
\)
が得られるはずで、実際にプロットして確かめてみると
N288_c
となります。
※Nは2以上のどんな整数でもokです。偶数、奇数、素数であっても。

※実部に現れる小さなピークに関して
原因はよくわかりません。\(N=2^m\)ではない時になんかこれが出てきます。
不安な人は\(N=2^m\)に固定して使ってください。
N=1024で行うとこうなります。
N1024

ベンチマーク

同じ計算機上でのベンチマーク結果を載せます。
コンパイラ：ifort (IFORT) 16.0.2 20160204
MKL:/opt/intel/compilers_and_libraries_2016.2.181
計算スレッド数1
ベンチマークdft_mkl_c

ベンチマーク用プログラムはこちら↓

並列計算時の計算速度
call mkl_set_num_threads(Ncpu)

を用いてFFTを並列計算で行った時、計算時間のCPU個数(8コア16スレッド)と配列の要素数依存性を調べます。
フーリエ変換するのはガウス関数とします。
結果はこのようになりました。

複数の線は試行回数を表しています。プログラムは全く同じで、実行した時刻だけが違います。
”計算1回”とは（順方向→逆方向フーリエ変換）の1セットを”1回”としています。
8コア16スレッドの環境で実行したので、物理的なCPUの個数8個で計算速度はおおよそ打ち止めになっています。

計算に用いたコードはこちら

2次元のフーリエ変換について

ルーチン等々はほぼ同じです。
下のプログラムは
\(
f(x,y)=\sin(4x)+\sin(3y)
\)
をフーリエ変換するものです。
フーリエ変換実行前の,位置空間での関数の形(“before.d”, x軸:\(x\), y軸:\(y\), z軸\(f(x,y)\)),
フーリエ変換し、波数空間で見た関数の形(“k.d”, x軸:\(k_x/(2\pi)\), y軸:\(k_y/(2\pi)\), z軸\(Imag\{g(k_x,k_y)\}\)),
逆フーリエ変換した後の関数の形(“after.d”, x軸:\(x\), y軸:\(y\), z軸\(f(x,y)\))は
before_k_after_2D_c
となります。
波数空間で\(k_x\)だけを見た時,\(k_x=0\)にあたかもピークが見えます。
これはy成分があるために出てくるもので、証明は以下のようにできます。
関数\(f(x,y)\)は
\(
\begin{align}
f(x,y)&=\sin(4x)+\sin(3y)\\
&=i\left\{-\frac{1}{2}e^{i4x}+\frac{1}{2}e^{-i4x}-\frac{1}{2}e^{i3y}+\frac{1}{2}e^{-i3y}\right\}
\end{align}
\)
なので、連続のフーリエ変換を考えるとすると、
\(
\begin{align}
&\int f(x,y)e^{-ik_x x}e^{-ik_y y}dxdy\\
&=i\left\{-\frac{1}{2}\int e^{i(4-k_x)x} dx \int e^{-ik_y y} dy+\frac{1}{2}\int e^{-i(4+k_x)x} dx \int e^{-ik_y y} dy \right. \\
&\;\;\;\;\;\;\;\;\;\;\;\; \left.-\frac{1}{2}\int e^{-ik_x x} dx \int e^{i(3-k_y)y} dy+\frac{1}{2}\int e^{-ik_x x} dx \int e^{-i(3+k_y)y} dy\right\} \\
&=i\pi\left\{-\delta(k_x-4)\delta(k_y)+\delta(k_x+4)\delta(k_y)-\delta(k_x)\delta(k_y-3)+\delta(k_x)\delta(k_y+3)\right\}
\end{align}
\)
となります。なので、フーリエ変換によってピークが出てくる位置というのは、
\((k_x,k_y)=(4,0), (-4,0), (0,3), (0,-3)\)
の4点となります。なので、y成分がいるために\(k_x=0\)にあたかもピークがあるように見えてしまうんですね。(証明終わり)

2次元のフーリエ変換はそのまま2次元データの配列を渡すことは出来ないため、一度、1次元配列への変換(サブルーチンswap1d2d)を呼び出しています。本来は2次元データと1次元データの間の関係ではEquivalenceを使うべきなのですが、何をしているのかの(個人的な)見やすさを考慮してわざわざ別の配列を用意して格納しています。

include "/opt/intel/mkl/include/mkl_dfti.f90"
program main
  implicit none
  integer::i,j,Nx,Ny
  double precision,allocatable::x(:),kx(:),mkx(:)
  double precision,allocatable::y(:),ky(:),mky(:)
  double precision::hx,xmin,xmax
  double precision::hy,ymin,ymax
  complex(kind(0d0)),allocatable::z(:,:),mz(:,:)

double precision,parameter::pi=dacos(-1.d0)
  complex(kind(0d0))::func
  external::func

Nx=100; Ny=100
  allocate(x(0:Nx-1),kx(0:Nx-1),mkx(0:Nx-1))
  x=0d0; kx=0d0; mkx=0d0
  allocate(y(0:Ny-1),ky(0:Ny-1),mky(0:Ny-1))
  y=0d0; ky=0d0; mky=0d0
  allocate(z(0:Nx-1,0:Ny-1),mz(0:Nx-1,0:Ny-1))
  z=dcmplx(0d0,0d0); mz=dcmplx(0d0,0d0)

xmin=-3d0<em>pi
  xmax=3d0</em>pi
  hx=(xmax-xmin)/dble(Nx)
  ymin=-3d0<em>pi
  ymax=3d0</em>pi
  hy=(ymax-ymin)/dble(Ny)

do i=0,Nx-1
     x(i)=xmin+dble(i)<em>hx
  enddo
  do j=0,Ny-1
     y(j)=ymin+dble(j)</em>hy
  enddo
  do i=0,Nx-1
     do j=0,Ny-1
        z(i,j)=func(x(i),y(j))<br />
     enddo
  enddo

open(21,file="before.d")
  do i=0,Nx-1
     do j=0,Ny-1
        write(21,'(4e20.10e3)')x(i),y(j),dble(z(i,j)),dimag(z(i,j))
     enddo
     write(21,*)
  enddo
  close(21)

call dftf(Nx,kx,hx); call dftf(Ny,ky,hy) ! get frequency

call dft2d(Nx,Ny,z,"forward") ! forward dft

call dfts2d(Nx,Ny,kx,ky,z,mkx,mky,mz) !sort frequency
  open(21,file="k.d")
  do i=0,Nx-1
     do j=0,Ny-1
        write(21,'(4e20.10e3)')mkx(i),mky(j),dble(mz(i,j)),dimag(mz(i,j))
     enddo
     write(21,*)
  enddo
  close(21)

call dft2d(Nx,Ny,z,"backward") ! backward dft
  open(21,file="after.d")
  do i=0,Nx-1
     do j=0,Ny-1
        write(21,'(4e20.10e3)')x(i),y(j),dble(z(i,j)),dimag(z(i,j))
     enddo
     write(21,*)
  enddo
  close(21)

stop
end program main

function func(x,y)
  implicit none
  double precision::x,y
  complex(kind(0d0))::func

func=dcmplx(dsin(4d0<em>y)+dsin(3d0</em>x),0d0)

return
end function func

がメインのプログラムで、以下が上で呼び出しているサブルーチンの中身です。

subroutine dftf(N,f,h)
  implicit none
  integer,intent(in)::N
  double precision,intent(in)::h
  double precision,intent(out)::f(0:N-1)

integer::i
  double precision::mf(0:N-1)

if(mod(N,2).eq.0)then
     do i=0,N-1
        mf(i)=(dble(i+1)-dble(N)<em>0.5d0)/(dble(N)</em>h)
     enddo

<pre><code> do i=0,N-1
    if(i.le.N/2-2)then
       f(i+N/2+1)=mf(i)
    else
       f(i-N/2+1)=mf(i)
    endif
 enddo

else
do i=0,N-1
mf(i)=(dble(i)-dble(N-1)0.5d0)/(dble(N)h)
enddo



 do i=0,N-1

    if(i.le.(N-1)/2-1)then

       f(i+(N-1)/2+1)=mf(i)

    else

       f(i-(N-1)/2)=mf(i)

    endif

 enddo

endif

return
end subroutine dftf
!——————————–
subroutine dft2d(Nx,Ny,z,FB)
!sikinote
!date : 2015/12/29
!developer : sikino
!CC BY 4.0 (http://creativecommons.org/licenses/by/4.0/deed.ja)
use MKL_DFTI
implicit none
integer,intent(in)::Nx,Ny
complex(kind(0d0)),intent(inout)::z(0:Nx-1,0:Ny-1)
character(*),intent(in)::FB

complex(kind(0d0))::z1(0:Nx*Ny-1)
integer::Status
TYPE(DFTI_DESCRIPTOR),POINTER::hand
integer::id(1:2)

!DFT : Discrete Fourier Transform
!
!n –> number of data.
!z(i,j) –> value of data at x and y
!FB –> “forward” : Forward DFT
! –> “backward” : Backward DFT

call swap1d2d(nx,ny,z,z1,1)
id(1)=nx; id(2)=ny

Status = DftiCreateDescriptor(hand,DFTI_DOUBLE,DFTI_COMPLEX,2,id)
Status = DftiCommitDescriptor(hand)

if(trim(FB).eq.”forward”)then
Status = DftiComputeForward(hand,z1)
elseif(trim(FB).eq.”backward”)then
Status = DftiComputeBackward(hand,z1)
else
write(6,*)”DFT string different”
stop
endif
Status = DftiFreeDescriptor(hand)

call swap1d2d(nx,ny,z,z1,-1)

if(trim(FB).eq.”backward”)then
z=z/dble(nx*ny)
end if

return
end subroutine dft2d
!———————————————
subroutine swap1d2d(nx,ny,z2,z1,isign)
!change array matrix between 1D and 2D
! if isign == 1 –> from 2D to 1D
! if isign == -1 –> from 1D to 2D
implicit none
integer,intent(in)::nx,ny,isign
complex(kind(0d0)),dimension(0:nx-1,0:ny-1),intent(inout)::z2
complex(kind(0d0)),dimension(0:nx*ny-1),intent(inout)::z1
integer::j1,j2,k

if(isign.eq.1)then
do j2=0,ny-1
do j1=0,nx-1
k=j2nx+j1
z1(k)=z2(j1,j2)
enddo
enddo
elseif(isign.eq.-1)then
do j2=0,ny-1
do j1=0,nx-1
k=j2nx+j1
z2(j1,j2)=z1(k)
enddo
enddo
endif

return
end subroutine swap1d2d
!—————————-
subroutine dfts2d(Nx,Ny,fx,fy,z,mfx,mfy,mz)
implicit none
integer,intent(in)::Nx,Ny
double precision,intent(in)::fx(0:Nx-1),fy(0:Ny-1)
complex(kind(0d0))::z(0:Nx-1,0:Ny-1)
double precision,intent(out)::mfx(0:Nx-1),mfy(0:Ny-1)
complex(kind(0d0)),intent(out)::mz(0:Nx-1,0:Ny-1)
complex(kind(0d0))::mmz(0:Nx-1,0:Ny-1)
integer::i,j,k,l

if(mod(Ny,2).eq.0)then
do i=0,Ny-1
if(i.le.Ny/2)then
j=i+Ny/2-1
mfy(j)=fy(i)
mz(0:Nx-1,j)=z(0:Nx-1,i)
else
j=i-Ny/2-1
mfy(j)=fy(i)
mz(0:Nx-1,j)=z(0:Nx-1,i)
endif
enddo
else
do i=0,Ny-1
if(i.le.(Ny-1)/2)then
j=i+(Ny-1)/2
mfy(j)=fy(i)
mz(0:Nx-1,j)=z(0:Nx-1,i)
else
j=i-(Ny-1)/2-1
mfy(j)=fy(i)
mz(0:Nx-1,j)=z(0:Nx-1,i)
endif
enddo
endif

mmz=mz

if(mod(Nx,2).eq.0)then
do k=0,Nx-1
if(k.le.Nx/2)then
l=k+Nx/2-1
mfx(l)=fx(k)
mz(l,0:Ny-1)=mmz(k,0:Ny-1)
else
l=k-Nx/2-1
mfx(l)=fx(k)
mz(l,0:Ny-1)=mmz(k,0:Ny-1)
endif
enddo
else
do k=0,Nx-1
if(k.le.(Nx-1)/2)then
l=k+(Nx-1)/2
mfx(l)=fx(k)
mz(l,0:Ny-1)=mmz(k,0:Ny-1)
else
l=k-(Nx-1)/2-1
mfx(l)=fx(k)
mz(l,0:Ny-1)=mmz(k,0:Ny-1)
endif
enddo
endif

return
end subroutine dfts2d

昔のMKLの離散フーリエ変換ルーチンについて

MKLは汎用性を高くするために仕様変更が行われ、
以前は

call zfft1d(...)

とするだけで順/逆方向フーリエ変換は使えていたようですが、いつからか(2004~2006年辺り?)廃止され、使えなくなっています。

間違い、勘違いがあると思います。
いかなる問題が発生しても私は責任を一切負いません。
それを念頭に置いたうえ使用してください。

参考文献

インテル®マス・カーネル・ライブラリーリファレンスマニュアル(2006年)
ニューメリカルレシピinC

ニューメリカルレシピ・イン・シー 日本語版―C言語による数値計算のレシピ

[adsense2]

プログラミングと数値計算

fortran90でファイルの存在を確かめる

コメントする

fortran90です。
ファイルが存在するかしないかを確かめるためのプログラムです。
組み込み関数”access”を用いると簡単に確認できます。

program main
  implicit none
  integer::access
 
  write(6,*)access( "./temp.d", " ")
 
  stop
end program main

上の文は、もしもファイル “./temp.d” が今いるディレクトリに存在するかどうかを返すものです。

また、access関数はファイルの有無(引数は空白文字” “)のみならず、そのファイルの読み取り権限(引数は”r”)、書き込み権限(引数は”w”)、実行権限(引数は”x”)を調べることができます。

ゼロが返ってきたらその権限(ファイルの有無)があるということです。
実行例は

 $ ls
main.f90  temp.d
 $ gfortran main.f90
 $ ./a.out
           0
 $

となります。intel®fortranコンパイラでももちろんできます。

[adsense1]

ここより以下のものはaccess関数を用いない方法です。

ファイルをopenできたら存在する、errorが返ってきたら存在しない、と判断します。
そのプログラムとサブルーチンはこんな感じ。

program main
  implicit none
  integer::ox
 
  call checkfile("temp.d",ox)
 
  write(6,*)ox
 
  stop
end program main

subroutine checkfile(fname,ox)
  implicit none
  character(*),intent(in)::fname
  integer,intent(out)::ox
 
  open(100,file=trim(fname),status='old',err=999)
  close(100)
  write(6,'(3A)')"file '",fname,"' exist"
  ox=1
  return
 
999 continue
  close(100)
  write(6,'(3A)')"file '",fname,"' don't exist"
  ox=0

  return
end subroutine checkfile

”ox”は〇×のつもりで使っています。
ファイル名”temp.d”をinputパラメータとして渡し、そのファイルが存在したらox=1を、存在しなかったらox=0を返します。

実行例は、

 $ ls
   temp.d main.f90
 $ gfortran main.f90
 $ ./a.out
   file 'temp.d' exist
           1
 $ gfortran main.f90
 $ ./a.out
   file 'temp1.d' don't exist
           0

となります。２回目の実行ではtemp.dをtemp1.dに変えて存在しない場合を表示しています。

[adsense2]

参考
ファイルの有無を調べる -知識の箱

プログラミングと数値計算, 弾道計算

弾道計算(BB弾)のコード

17件のコメント

2016/03/13更新
2016/03/18更新
2016/03/20更新
2016/07/23更新
2016/07/25更新
弾道計算(BB弾)のfortran90での計算コードです。


2016/03/18 以降のコードに関して
コードは改変、複製、配布、営利、非営利目的に使用していただいて構いません。
ただし著作権を放棄している訳ではありません。
また、このコードを利用した場合、クレジットの表記をお願いします。

Linuxで、OSがUbuntuやLinuxMintの場合はfortranコンパイラは

sudo apt-get install gfortran

でインストールできます。
intel®fortranコンパイラとgfortranコンパイラで正常に動くことは確かめています。

もろもろのファイルをmain.f90に記述した後、
gfortranコンパイラなら

gfortran main.f90
./a.out

intel®fortranコンパイラなら

ifort main.f90
./a.out

で実行、ファイル”orbit.txt”が生成されます。

弾道計算に関するその他ページ
弾道計算(BB弾)の理論
BB弾の回転量について(実験との比較)
弾道計算(BB弾)の結果
弾道計算の結果2, 比較と詳細データ
弾道計算(BB弾)のコード(fortran90)←今ここ
弾道計算のコード(Excel)
バレル内部でのBB弾の方程式
水中下でのBB弾の弾道計算

できること

・BB弾の弾道計算(重力、慣性抵抗、粘性抵抗、ホップ有、回転の減衰を考慮)ができます。
式で表すと
BB弾の運動方程式(式1)
\(
\displaystyle m\frac{d^2 \vec{r}}{dt^2}=
-mg\vec{k}+\left\{-6\pi \eta R |\vec{V}|-\frac{1}{2}C_d \rho \pi R^2 |\vec{V}|^2\right\}\frac{\vec{V}}{|\vec{V}|}
-\frac{4}{3}\pi R^3 2\rho |\vec{\omega}| |\vec{V}|\frac{\vec{V}\times\vec{\omega}}{|\vec{V}\times\vec{\omega}|}
\)
とBB弾の回転角速度の減衰を記述する方程式(式2)
\(
\begin{align}
\frac{d\omega_z}{dt}&=N_z/I \\
N_z&=\frac{\rho C_f}{2}R^3\int_0^{2\pi}d\varphi\int_0^{\pi}d\theta |u\sin\varphi-R\omega\sin\theta|(u\sin\varphi-R\omega\sin\theta)\sin^2\theta
\end{align}
\)
の2つを考慮したものになります。

・ゼロインを決めた時の最適な軌道
・上下方向が最小の振れ幅になる軌道
も出力することができます。

※BB弾の回転の減衰を考えるか考えないかは、inputパラメータ”omgdecay”により制御できます。
回転の減衰に疑問がある場合はこれを利用してください。

注意事項

・とんでもないパラメータでは計算できない場合がありますので、そこら辺はご了承ください。
・僕のプログラミング知識ではこれが限界です。
僕はこのプログラムを使用したことによって生じた、いかなる場合でも一切の責任を負いません。

主なルーチンと説明

計算理論

式1の時間発展はルンゲ-クッタ-フェールベルグ法による適応刻み幅制御を行っています。デフォルトでは8桁の精度で(※1)一致するようになっています。モジュールGBL内のtol=1d-8を変更することで精度を制御することができます。
式2の右辺の積分は、デフォルトでは15次のガウス=クロンロッド求積法を2回用いています。今回の問題に対してはあまり良い方法とは言えません。ですが、5-6桁はロンバーグ積分法による結果と一致するので、実用上問題ないと判断し、計算速度を優先しました。
式2の回転減衰の効果はルンゲクッタ法を行う際の連立方程式に組み込んでいません。これだけはオイラー法で行っています。なので、1ステップ目の刻み幅をどうするかで、本来はずれないはずですが、変更するとそれが伝搬してゼロイン位置がずれます。デフォルトでは全ての1ステップ目の刻み幅は統一してあるので上記問題は起こりません。本当に厳密にやるならば、コードを変更して減衰を連立方程式内に組み込んでください。
2016/03/13更新
コードを変更して減衰を連立方程式内に組み込みました。計算は若干遅くなりますが、上記問題はもう起こり得ません。

※1 相対誤差もしくは絶対誤差で、です。この理由は、相対誤差に統一してしまうと値が非常に小さいとき、機械誤差によっていつまでたっても収束しないためです。なので、元になる値が1より小さいときは絶対誤差を、1より大きいときは相対誤差を取るようにしています。

2016/03/13更新
2016/03/18更新
2016/03/20更新
2016/07/23更新
2016/07/25更新
2017/10/07更新

inputファイルとコード

inputというファイル名のファイルを作り、中に

&input
m=0.25d-3,          ! Weight of bullet [kg]
a=3d-3,             ! Radious of bullet [m]
energy=0.90d0,      ! Energy of bullet [J]
g=9.80665d0,        ! Gravity [m s^-2]
temperature=20d0,   ! Temperature [degree of Celsius]
pressure=101325d0,  ! Atmospheric pressure [Pa]
moisture=0.6d0,     ! Value of moisture 0~1
omgdecay="yes",     ! Introduce omega decay. "yes" or "no"
search_zeroin="yes", ! Search good vz and omgy. "yes" or "no"
zeroinx=50d0,       !   +-- if yes, set zeroin length [m]
search_updown="yes",! Search good vz and omgy. "yes" or "no"
updownz=0.20d0,     !   +-- if yes, set up-down henght [m]
theta=-1.7d0,
rx=0.d0,
ry=0.d0,  
rz=1.0d0,
vy=0.d0,
ux=0.d0,
uy=0.d0,
uz=0.d0,
omgx=0.d0, !*2pi
omgy=-340d0, !*2pi
omgz=0.d0, !*2pi
outputfile="orbit", ! Output file name prefix
stept=0.01d0,       ! Output interval [s]
ik=0, ! ik=0 exact, ik=1 approximation
&end

と記述してください。このinputファイルの例では、

重さ(m)0.25g
エネルギー(energy)0.90J
重力加速度(g)9.80665m/(s^2)
温度(temperature)20度
大気圧(pressure)101325Pa
湿度(moisture)60%
回転角速度の減衰(omgdecay)有り
ゼロインによる最適な弾道の探索(search_zeroin)しない
↑がyesなら、ゼロインの位置(zeroinx)50m
上下幅による最適な弾道の探索(search)する
↑がyesなら、上下幅(updownz)0.1m
銃口の位置(rz)1m
y方向の速度(vy)0m/s
射出角度(theta)-0.5度
風速(ux=uy=uz=0)0m/s
ホップ(omgy)-210回転/s
出力ファイル名の接頭語(outputfile)orbit
出力の時間間隔(stept)0.01秒ごとに出力
回転の減衰を表す積分の近似方法(ik=0 : 積分の直接計算、ik=1 : 積分を近似)

となっています。
回転の減衰を考えたくない時は、omgdecay=”no”としてください。

そして下のコードを展開し、コピー&ペーストして一つのファイル(例えばmain.f90という名前)に入れてください。
(カーソルをmodule GBL…のmに合わせた後、一番下までスクロールしてshift+クリックで楽に選択できます。)
個人的にファイルをダウンロードして…は嫌いなのでこういう形で公開しています。

2016/03/13 更新
2016/03/18 更新
2016/03/20 更新
2016/07/08 更新
2016/07/23 更新
2016/07/25 更新
2016/08/06 更新
ソースコードは1523行あります。

使い方

・fortranコンパイラ(ifortやgfortran)
・gnuplot
を想定しています。

2017/02/04追)
Bernadotte66さんの報告より、windows上でgfortranコンパイラを用いる際、gfortranコンパイラがver5.1以前のものである場合エラーが出てしまうようです。ver6.2以上のコンパイラを推奨します。

Linuxでは、gfortranコンパイラver4.8.4で問題ないことは確かめています。

それでコンパイルして実行してください。
一連の流れとしては、以下の通りです。

$ gfortran main.f90
$ ./a.out
&INPUT
 M=  2.5000000000000001E-004,
 A=  3.0000000000000001E-003,
 ENERGY= 0.90000000000000002     ,
 G=  9.8066499999999994     ,
 TEMPERATURE=  20.000000000000000     ,
 PRESSURE=  101325.00000000000     ,
 MOISTURE= 0.59999999999999998     ,
 OMGDECAY="yes                                             ",
 SEARCH_ZEROIN="yes                                             ",
 ZEROINX=  50.000000000000000     ,
 SEARCH_UPDOWN="yes                                             ",
 UPDOWNZ= 0.20000000000000001     ,
 THETA=-0.50000000000000000     ,
 RX=  0.0000000000000000     ,
 RY=  0.0000000000000000     ,
 RZ=  1.0000000000000000     ,
 VY=  0.0000000000000000     ,
 UX=  0.0000000000000000     ,
 UY=  0.0000000000000000     ,
 UZ=  0.0000000000000000     ,
 OMGX=  0.0000000000000000     ,
 OMGY= -210.00000000000000     ,
 OMGZ=  0.0000000000000000     ,
 STEPT=  1.0000000000000000E-002,
 OUTPUTFILE="orbit                                           ",
 IK=          0,
 /
Reach ground at  0.980E+00
===============================
 search vz and wy by zeroin -->   50.00000
 ---- false_position
 Sequence 1 of 25 : omgy/2pi  -767.96000000 : zeroin   100.00000000
 Sequence 2 of 25 : omgy/2pi  -735.92000000 : zeroin   100.00000000
 Sequence 3 of 25 : omgy/2pi  -703.88000000 : zeroin   100.00000000
 Sequence 4 of 25 : omgy/2pi  -671.84000000 : zeroin   100.00000000
 Sequence 5 of 25 : omgy/2pi  -639.80000000 : zeroin   100.00000000
 Sequence 6 of 25 : omgy/2pi  -607.76000000 : zeroin   100.00000000
 Sequence 7 of 25 : omgy/2pi  -575.72000000 : zeroin    17.22935680
 Sequence 8 of 25 : omgy/2pi  -543.68000000 : zeroin    15.16813299
 Sequence 9 of 25 : omgy/2pi  -511.64000000 : zeroin    12.90696782
 Sequence 10 of 25 : omgy/2pi  -479.60000000 : zeroin    10.42356698
 Sequence 11 of 25 : omgy/2pi  -447.56000000 : zeroin     7.68923851
 Sequence 12 of 25 : omgy/2pi  -415.52000000 : zeroin     4.66717441
 Sequence 13 of 25 : omgy/2pi  -383.48000000 : zeroin     1.30910154
 Sequence 14 of 25 : omgy/2pi  -351.44000000 : zeroin    -2.44916818
  +--- Root find, go to conversion phase --+  
Conversion =>       0.00033361%
Conversion =>       0.00902931%
Conversion =>       0.24389215%
Conversion =>       6.77148580%
Conversion =>     104.78145034%
 ---- false_position
 Sequence 1 of 1 : omgy/2pi  -366.89192543 : zeroin    -0.57095843
  +--- Root find, go to conversion phase --+  
Conversion =>      0.00075051%
Conversion =>      0.08150935%
Conversion =>      8.77790138%
Conversion =>    579.19923386%
  -371.80132940800553       -2.5094934272766087    
Reach ground at  0.170E+01
===============================
 search vz and wy by up-down height -->    0.20000
 ---- false_position
 Sequence 1 of 25 : omgy/2pi  -959.96000000 : zeroin    19.21215678
 Sequence 2 of 25 : omgy/2pi  -919.92000000 : zeroin    17.82618788
 Sequence 3 of 25 : omgy/2pi  -879.88000000 : zeroin    16.40387340
 Sequence 4 of 25 : omgy/2pi  -839.84000000 : zeroin    14.95914779
 Sequence 5 of 25 : omgy/2pi  -799.80000000 : zeroin     9.80000000
 Sequence 6 of 25 : omgy/2pi  -759.76000000 : zeroin     9.80000000
 Sequence 7 of 25 : omgy/2pi  -719.72000000 : zeroin     9.80000000
 Sequence 8 of 25 : omgy/2pi  -679.68000000 : zeroin     9.80000000
 Sequence 9 of 25 : omgy/2pi  -639.64000000 : zeroin     9.80000000
 Sequence 10 of 25 : omgy/2pi  -599.60000000 : zeroin     9.80000000
 Sequence 11 of 25 : omgy/2pi  -559.56000000 : zeroin     1.45287636
 Sequence 12 of 25 : omgy/2pi  -519.52000000 : zeroin     1.10739833
 Sequence 13 of 25 : omgy/2pi  -479.48000000 : zeroin     0.80320434
 Sequence 14 of 25 : omgy/2pi  -439.44000000 : zeroin     0.54089117
 Sequence 15 of 25 : omgy/2pi  -399.40000000 : zeroin     0.32062896
 Sequence 16 of 25 : omgy/2pi  -359.36000000 : zeroin     0.14208536
 Sequence 17 of 25 : omgy/2pi  -319.32000000 : zeroin     0.00427472
 Sequence 18 of 25 : omgy/2pi  -279.28000000 : zeroin    -0.09464709
  +--- Root find, go to conversion phase --+  
Conversion =>       0.00183551%
Conversion =>       0.00989716%
Conversion =>       0.05297720%
Conversion =>       0.28399231%
Conversion =>       1.52644398%
Conversion =>       8.26417579%
Conversion =>      39.67268127%
Conversion =>     179.58494188%
 ---- false_position
 Sequence 1 of 1 : omgy/2pi  -312.86039405 : zeroin    -0.01405310
  +--- Root find, go to conversion phase --+  
Conversion =>      0.00065868%
Conversion =>      0.03115036%
Conversion =>      1.49063057%
Conversion =>    102.25144566%
  -317.78761751909860       -1.5473682498931858    
Reach ground at  0.149E+01
     2.681[CPU sec]
$ gnuplot

    G N U P L O T
    Version 4.6 patchlevel 4    last modified 2013-10-02
    Build System: Linux x86_64

    Copyright (C) 1986-1993, 1998, 2004, 2007-2013
    Thomas Williams, Colin Kelley and many others

    gnuplot home:     http://www.gnuplot.info
    faq, bugs, etc:   type "help FAQ"
    immediate help:   type "help"  (plot window: hit 'h')

Terminal type set to 'wxt'
gnuplot> plot "orbit.txt" u 2:4 w l
gnuplot> replot "orbit_h.txt" u 2:4 w l
gnuplot> replot "orbit_opt.txt" u 2:4 w l

gnuplot上での出力が
orbitGK
のようになったら成功です。

出力ファイルは最大で3つ生成(inputパラメータ “search_??????” )に依存します。
orbit.txtは、inputによって与えられたパラメータでそのまま計算したもの
orbit_opt.txtはそのエネルギー、重さで上下振れ幅が最小でゼロイン位置がinputのzeroinxのもの
orbit_h.txtはそのエネルギー、重さで上下振れ幅が最小で上下振れ幅がinputのupdownzのもの
となります。

出力ファイルはそれぞれ14列あり、
時刻t[秒],位置x,y,z[m]:x,y,z方向の速度[m/s]:BB弾の回転数x,y,z[回転/s]:x,y,z方向の風速:エネルギー[J]
になっています。
出力ファイルの中身は下のように書いてあるので、見てください。

# m:    0.25000E-3[kg]
# g:       9.80665[m s^{-2}]
# a:    0.30000E-2[m]
# temperature:      20.00000[degree Celsius]
# moisture:       0.60000[no-dimension]
# pressure:  101325.00000[Pa]
# eta:  0.18197E-4[kg m^{-1} s^{-1}]
# rho:     1.20123[kg m^{-3}]
# tol:  0.10000E-7
# eps:  0.10000E-8
#         t[s]         x[m]         y[m]         z[m]      vx[m/s]      vy[m/s]      vz[m/s]    wx[rot/s]    wy[rot/s]    wz[rot/s]   windx[m/s]   windy[m/s]   windz[m/s]    Energy[J]
     0.000000     0.000000     0.000000     1.000000    84.849583     0.000000    -0.740471     0.000000  -210.000000     0.000000     0.000000     0.000000     0.000000     0.900000
     0.010000     0.838863     0.000000     0.992911    82.937695     0.000000    -0.678265     0.000000  -209.261019     0.000000     0.000000     0.000000     0.000000     0.859890
    ...

実銃のパラメータでの計算

実銃は特に想定していませんが、計算ができると言えばできます。
実銃と比較して考慮されていない点は3つあり、

  1. 滑らかな球として考えているので抗力係数\(C_d\)が異なる
  2. 銃弾の歳差運動が考慮されていない(球ではなく細長いので)
  3. ライフリングによる回転が異なる
  4. コリオリ力を考慮に入れていない

1番目はかなり重大です。
私のプログラムを使うと1000メートル地点でのエネルギーの値が2倍ほど大きく出ます。
ライフリングによる横回転だけなら問題なく考慮できます。バックスピンで入っていたものを変えればいいだけです。
ただし、風が無い、またはコリオリ力を考慮しない状況では計算上は横回転があっても無くても変わりません。
風とコリオリ力がある時にライフリングによる弾道のずれが影響します。
2番目はどんな力が働いたりするか私には見当がつきません。専門の方じゃないと分かるまでに時間がかかると思います。
コリオリ力はプログラムで入れればいいです。これはそこまで辛くないでしょう。
ライフリングの効果は入れられるには入れられますが、ちょっと変な結果になるので使わないほうが良いでしょう。
なのでこの効果は入れていません。
パッと思いつく限り、違う点は以上です。

よって、本プログラムによる実銃を想定した計算とは、
風が無く、赤道上で実銃と同じ口径の滑らかな球が実銃のエネルギーで射出されたときの弾道
となります。

ちなみに、.338ラプアマグナム弾の場合を想定したinputパラメータはこちら
角度を使って計算しているのでそこのとこのプログラムは各自直してください。

&input
! For .338 lapua magnum( but treated as sphire)
m=16.2d-3,          ! Weight of bullet [kg]
a=4.29d-3,      ! Radious of bullet [m]
energy=6562.16d0,      ! Energy of bullet [J]
g=9.80665d0,        ! Acceleration of Gravity [m s^-2]
temperature=20d0,   ! Temperature [degree of Celsius]
pressure=101325d0,  ! Atmospheric pressure [Pa]
moisture=0.6d0,     ! Value of moisture 0~1
omgdecay="no",     ! Introduce omega decay. "yes" or "no"
search_zeroin="no", ! Search good vz and omgy. "yes" or "no"
zeroinx=50d0,       !   +-- if yes, set zeroin length [m]
search_updown="no",! Search good vz and omgy. "yes" or "no"
updownz=0.10d0,     !   +-- if yes, set up-down henght [m]
theta=0.6d0,        ! default, theta is neglected. change program main.
rx=0.d0,
ry=0.d0,  
rz=1.0d0,
vy=0.d0,
vz=0.d0,
ux=0.d0,
uy=0.d0,
uz=0.d0,
omgx=0.d0, !*2pi [omgx*2*pi rotate/s]
omgy=0.d0, !*2pi [omgy*2*pi rotate/s]
omgz=0.d0, !*2pi
outputfile="orbit", ! Output file name prefix
stept=0.01d0,       ! Output interval [s]
&end

上記パラメータで計算するとこんな軌道が得られます(点は0.05秒ごとの軌跡)
lapua338_c

.338ラプアマグナム弾のデータは↓
.338 ラプアマグナム弾 / .338 Lapua Magnum -MEDIAGUN DATABASE
を参照しました。

改良版

これから最適な回転量、角度を探すプログラムに更新したいと思いますが、とりあえず
改良したプログラムを載せておきます。
input
main.f90

プログラミングと数値計算

対角化

コメントする

fortran90でLapackを用います。

倍精度ルーチン

・複素一般行列の対角化, diag(N,A,Ev)
・複素エルミート行列の対角化, zdiag(N,A,Ev)
・実エルミート行列の対角化, herdiag(N,A,Ev), diag_slct(N,A,Nq,Ev)
・実対称三重対角行列の対角化, tridiag(N,A,Ev)

4倍精度ルーチン

・複素一般行列の対角化, qdiag(N,A,Ev)

複素一般行列　　　　　　　　　：行列要素は複素数で、値は適当に入っている。
複素エルミート行列　　　　　　：行列要素は複素数で、\(A_{i,j}=A^*_{j,i}\)の関係がある。
実エルミート行列（実対称行列）：行列要素は実数で、\(A_{i,j}=A_{j,i}\)の関係がある。
実三重対角行列　　　　　　　　：行列要素は実数で、\(A_{i,i}\ne 0, ~~A_{i,i-1}=A_{i,i+1}\)、それ以外の要素は0。

（エルミートな行列と言われたら、その固有値は必ず実数です。）

複素一般行列の場合

wolframのページで、
Eigenvectors[{-2,2,4},{-2,4,2},{-2,1,4}]
と入力することで固有値、固有ベクトルを出してくれます。これで確かめを行いました。

・この行列を入れた時、正しい固有値と固有ベクトルを出力してくれる事を確認しています。
・40×40程度の行列でも確認しています。
・複素数の固有値、固有ベクトルでも確認しています。

Lapackを使います。

ifort -mkl main.f90

か、g95, もしくはgfortranでlapackにリンク

gfortran -llapack main.f90

して、コンパイルしてください。

大きさN×Nの正方行列Aを入力すると, 固有値Ev(1:N)とそれぞれの固有値の要素に属する固有ベクトルを上書きして返します。(1列目に固有値Ev(1)に属する固有ベクトルが、2列目に固有値Ev(2)に属する固有ベクトルが･･･という感じに。)
また固有値は実部の大小を比較して小さい順に並べ替えられ、固有ベクトルもそれに対応した列に格納されます。
計算速度は\(O(N^3)\)です。

複素エルミート行列の対角化

対称行列に限る場合は計算を早くすることができます。
ここでは複素エルミート行列を対角化して固有値を出力するケースを考えましょう。

そのまま使うのはやはり面倒です。なぜなら、このルーチンを用いるためにはワーク配列を何個か宣言して用いる必要があるからです。
このワーク配列は計算に直接関係ありません。
僕らが欲しいのは行列を入れたら固有値と固有ベクトルを返してくれさえすればいいものです。
下はサブルーチンzdiagは僕が作ったもので、複素エルミート行列\(A(1:N,1:N)\)を入れると固有値の配列\(Ev(1:N)\)と\(A\)に上書きして固有ベクトル\(A(1:N,1:N)\)を入れて返します。

エルミート行列なので、実は全ての行列要素を求める必要ありません。
下の例では、
エルミート行列\(A(i,j)\)は,\(j\ge i\)の上三角の要素だけ分かっていればいいです。
なので、コメントアウトしている行列でも同じ固有値と固有ベクトルが得られます。
以下のプログラムをlapackと共にコンパイルし、実行すると、

>./a.out   
(-0.42419E+01)
( 0.21580E+01)
( 0.80839E+01)
(-0.88841E+00,-0.92008E-01)  (-0.15371E+00,-0.13692E+00)  (-0.39557E+00, 0.58357E-01)  
( 0.17660E+00, 0.14837E-01)  ( 0.55488E+00,-0.56209E+00)  (-0.35177E+00, 0.47012E+00)  
( 0.41334E+00, 0.00000E+00)  (-0.57775E+00, 0.00000E+00)  (-0.70382E+00, 0.00000E+00)  
>

という結果が得られます。
1番目の固有値(-0.42419E+01)に対応する固有ベクトルは1列目に、
2番目の固有値( 0.21580E+01)に対応する固有ベクトルは2列目…というように代入されていきます。

プログラム例)
確かめ用のリンクEigenvectors[{-2,2+i,4},{2-i,4,1-2i},{4,1+2i,4}] -wolfram alpha

行と列を入れ替えたければ、以下のルーチンを使用してください。

subroutine swap2d(N,A)
  implicit none
  integer::N,i,j
  double precision::temp,A(1:N,1:N)
 
  do i=1,N
     do j=i+1,N
        temp=A(i,j)
        A(i,j)=A(j,i)
        A(j,i)=temp
     enddo
  enddo
 
  return
end subroutine swap2d

実エルミート行列の場合

以下のサブルーチンを使うと良いでしょう。

\(N\times N\)の実対称行列で、固有値の小さい方から順に\(Nq\)個だけ欲しい場合

実エルミート行列の、
完全対角化(herdiag)と、
選択的に固有値を計算する(diag_slct)
の比較(実時間計測)を行いました。


横軸が選択的に何個固有値を計算するか、を表し、縦軸は実際に掛かった時間をあらわします。
色は行列サイズによって変化し、黄色,紫、緑の順に行列サイズが\(2^{16}\times 2^{16}, 2^{14}\times 2^{14}, 2^{12}\times 2^{12}\)を表します。
破線はherdiagによる対角化時間、実線はdiag_slctによる計算時間を表します。
また、図中の値2,4,8は計算に用いたCPUの数です。CPUは8コア16スレッドです。

ここから分かるのは、欲しい固有値の数が行列サイズ\(N\times N\)の\(N/4\)より小さい場合は選択的に計算した方が早い、ということです。
選択的に計算する方は計算時に要求されるメモリも少なくて済みます。


(2016/08/27)

4倍精度の一般行列の対角化

4倍精度用の一般行列(複素非エルミートを含む)を対角化し、その固有値と固有ベクトルを出力するfortranコードが以下のページ
４倍精度対角化 -ただの備忘録
で公開されていました。
ここで公開されているサブルーチン”qeispack.f90”
を使うことで対角化できます。
gfortranとifortで動くことは確かめました。
コンパイルは例えば

ifort qeispack.f90 main.f90

のようにしてください。
本ページと同じように書き加えて実装すると以下のようになります。

実対称三重対角行列の場合

実対称三角行列の場合も置いておきます。

プログラミングと数値計算, ルンゲ=クッタ法

ルンゲ=クッタ法の説明と刻み幅制御

5件のコメント

ルンゲ=クッタ法(Runge-Kutta method、RK法)とは?
僕の知る限りの知識で紹介します。

特に良く使われる陽的ルンゲ=クッタ法は、
・実装が簡単
・良いアルゴリズムではない
という手法です。

良いアルゴリズムである陰的ルンゲ=クッタ法は、
陰的ルンゲ=クッタ法
をご覧ください。

もくじ

  1. ルンゲ=クッタ法の系統
  2. 陽的ルンゲ=クッタ法の段数と次数について
  3. 誤差について
  4. Butcher tableによるルンゲ=クッタ法の記述
  5. 埋め込まれた陽的ルンゲ=クッタ法
  6. ルンゲ=クッタ=フェールベルグ法による刻み幅の自動制御(適応刻み幅制御)
  7. ルンゲ=クッタ=フェールベルグ法による刻み幅の自動制御のプログラム(実数/複素数 で １/２ 階微分方程式を解くプログラム)
  8. 不連続な点を含む場合
  9. 刻み幅制御のベンチマーク(振り子)
  10. （追記）ルンゲ=クッタ=ドルマンド=プリンス法
  11. 陽的ルンゲ=クッタ法の導出
  12. 参考文献

理論はいいから4次ルンゲ=クッタ法の数値計算ではどうやるかだけ教えてくれ！という方は4次ルンゲ=クッタ法へどうぞ。

ルンゲ=クッタ法の系統

ルンゲ=クッタ法は微分方程式の数値計算解を得るための手法のことを指します。

通常の参考書で微分方程式を解くために良く紹介されているのは、オイラー法と中点法、４次ルンゲ=クッタ法でしょうか。
オイラー法も中点法も4次ルンゲ=クッタ法も、”陽的ルンゲクッタ法” と呼ばれる枠組みの1つです。

オイラー法は正確には “陽的1段1次ルンゲ=クッタ法” と呼ばれ、
中点法は “陽的2段2次ルンゲ=クッタ法”、
4次ルンゲクッタ法(RK4)は、”陽的4段4次ルンゲ=クッタ法” と呼ばれています。

“段”と”次”とはなんなんでしょう？それは、
計算の大変さ(段)と、計算の正確さ(次)
です。”段“の値が小さければ小さいほど計算時間が少なくて済みますし、”次“の値が高ければ高いほど計算が正確です。

オイラー法は1という計算コストで正確さ1が得られますし、RK4は4という計算コストで正確さ4が得られます。

4次ルンゲ=クッタ法が使われる理由

理由は実装が簡単でそれなりの精度を持つから。です。

陽的ルンゲ=クッタ法において、pという計算コスト(p段)で、pより大きな正確さq(q次)を得ることはできません。
Derivation of Runge–Kutta methodsによれば、
\(q\)次の正確さ(q次のオーダー)を得たい場合、最低限必要な段数\(p_{\mbox{min}}(q)\)は

という関係にあります。

ここで注目するべきは4次の時までは計算コストに比例して計算精度が上がっていきます。
しかし、4次以上では、計算コストの増加と計算の正確さが見合わなくなっていきます。故に
計算効率が良いのは4次だろう、
と予想できます。
また、重要な理由として、4次ルンゲ=クッタ法に現れる係数が5次以降と比べて圧倒的にシンプルであることが挙げられます。4次では\(0,1/2,1/6\)程度の係数だけが使われ、プログラム作成時の入力ミスがほぼ生まれません。しかし、5次では\(28561/56430, -7200/2197\)といった係数が数多く出てきます。
これらの理由から,4次ルンゲ=クッタ法(RK4)が数値計算科学の世界でよく使われるのです。

陽的ルンゲ=クッタ法に限って言えばプログラムの実装が非常に簡単であることが挙げられます。陰的ルンゲ=クッタ法と呼ばれるアルゴリズムもあり、これは陽的ルンゲ=クッタよりも優れていますが、計算量が多くなり、若干複雑なアルゴリズムになります。陰的ルンゲ=クッタ法を詳しく知りたい方は陰的ルンゲ=クッタ法をご覧ください。

Q. オイラー法もものすごく細かい分点を取れば、その計算効率はRK4と同じなんじゃないの？
A. 刻み幅の乗数で効いてくるのでそうではありません。高次を使っても計算が信頼できるのであれば、大きなステップで進めるほうが早いです。例え、目標精度への計算時間が同じだとしても、計算機の有効桁数によって否定されてしまいます。
RK4で典型的にとられる時間ステップの間隔は、おおよそ\(10^{-2}\sim 10^{-4}\)程度であり、RK4のエラーのオーダーは\(O(h^5)\)です。
そして、科学計算で使う際の有効桁数は倍精度型で16桁です。
1ステップ当りの誤差は\(h\)の減少に伴い、解が\(h^4\)の早さで収束していく、と言えます。
だから16桁の計算では\(h=10^{-1}\to h=10^{-4}\)に変化させる時、誤差は\(O(h^5)=c 10^{-5}\to c 10^{-20}\)と変化します。
おおよそ\(c\approx 1\)と見積もれば、(有効桁数16桁を多少超えてしまいますが、)有効桁数いっぱいまで正しい値が出るであろうと期待できます。

これに対し、オイラー法で同じことをするには\(h\)を\(10^{-16}\)にしなくてはいけません。
\(t\)の値が\(10^{-16}\)変わった時に、桁落ちの問題を回避できるほど関数\(x\)の値に変化が生まれるか？
が問題になり、まぁそんな急激な変化は生まれないでしょう、と予想できます。これでは桁落ちの問題を回避するほどの変化は到底望めません。

よって計算の効率と有効桁数の限界から、RK4なのです。

また、あまりに高次の方法を使うとルンゲ現象に代表される不安定性といったことが起こるかもしれません。
高次は高精度という意味とイコールではないことに注意しましょう。この高次≠高精度については等間隔の分点における数値積分の時に書いたので気になる方はどうぞ。。

誤差について

4次ルンゲ=クッタ法の、1ステップ当りの誤差は\(h^5\)に比例,すなわち\(O(h^5)\)です。
しかし、通常は細かいルンゲ=クッタ法を何度も繰り返して計算します。
区間\([a,b]\)を刻み幅\(h\)の4次ルンゲ=クッタ法で\(N\)回のステップを繰り返し計算する場合、誤差は
\(
\displaystyle N\cdot O(h^5)=\frac{b-a}{h}\cdot O(h^5)=O(h^4)\)

となります。よって、\(N\)回繰り返すような計算では、オーダーが１つ落ちる事に注意しましょう。

[adsense1]

陽的ルンゲ=クッタ法の段数と次数について

さて、ここまで“段”は計算コスト、で“次”は計算の正確さ、という曖昧な表現でしたが、その表現をちゃんと知りましょう。
段と次を知るためにはルンゲ=クッタ法の計算方法を知る必要があります。
具体例を載せます。
\(
\displaystyle \frac{dx}{dt}=f(t,x)
\)
の、\(t_{n}\rightarrow t_{n}+h\ \ (=t_{n+1})\)における初期値問題に対する、
陽的1段1次ルンゲ=クッタ法(オイラー法)の計算スキームは、
\(
x_{n+1}=x_{n}+h\cdot f(t_{n},x_{n})
\)
です。

陽的4段4次ルンゲ=クッタ法(RK4)の計算スキームは、
\(
\begin{align}
k_1&=f(t_n, x_n) \\
k_2&=f(t_n+h/2, x_n+h k_1/2) \\
k_3&=f(t_n+h/2, x_n+h k_2/2) \\
k_4&=f(t_n+h, x_n+h k_3) \\
x_{n+1}&=x_{n}+{(k_1+2k_2+2k_3+k_4)}h/6
\end{align}
\)
として与えられます。

一般的に、陽的s段のルンゲ=クッタ法とは
\(
\begin{align}
g_i&=x_n+h\sum_{j=1}a_{i,j}k_j\ \ \ (j\lt i, \ i=1,2,…,s) \\
k_i&=f(t_n+c_ih, g_i) \\
x_{n+1}&=x_n+h\sum_{i=1}^s b_ik_i
\end{align}
\)
として書けます。
ここで行列形式で与えられる係数\(a_{i,j}, b_{i},c_{i}\)によって、そのs段ルンゲ=クッタ法が持つ次数が決められます。段数はここから由来します。

点\((t_n, x(t_n))\)周りで関数をテーラー展開し、その関数が点\((t_n+h\ \ (=t_{n+1}), x(t_{n+1}))\)で作る点を近似解とするのがルンゲ=クッタ法です。
故に、\(x(t_{n+1})\)は、
\(
\begin{align}
x(t_{n+1})=x(t_n)+\left.\frac{h}{1!}\frac{dx}{dt}\right|_{t=t_n}+\left.\frac{h^2}{2!}\frac{d^2x}{dt^2}\right|_{t=t_n}+\left.\frac{h^3}{3!}\frac{d^3x}{dt^3}\right|_{t=t_n}+\left.\frac{h^4}{4!}\frac{d^4x}{dt^4}\right|_{t=t_n}+…
\end{align}
\)
と書けます。
ここで、テイラー展開としてどの程度一致させて\(x(t_n+h)\)を決定するか？を表すのが次数に当たります。

言葉で書くなら、

局所的な精度は，高次項が解のテイラー(Taylor)展開とどの程度マッチしているかによって測定される．オイラー法は一次精度なので，誤差は一次高く\(h^2\)というベキ乗から始まる． ~NDSolveの”ExplicitRungeKutta”メソッドより

ともあります。

Butcher tableによるルンゲ=クッタ法の記述

行列形式で与えられるルンゲ=クッタ法での係数\(a_{i,j}, b_{i},c_{i}\)は何なのか？
具体的に記述してみましょう。
オイラー法(1段1次)はもっとも単純で、係数は
\(
\begin{align}
a_{1,1}&=0　\\
b_{1}&=1 \\
c_{1}&=0
\end{align}
\)
です。これを一般的な表記法の式に当てはめれば、
\(
\begin{align}
g_1&=x_n+h a_{1,1}k_1 \\
k_1&=f(t_n+c_1h, g_1) \\
x_{n+1}&=x_n+h b_1k_1
\end{align}
\)
となります。

中点法は、
\(
\begin{align}
a_{1,1}&=0 \\
a_{1,2}&=0 \\
a_{2,1}&=1/2 \\
a_{2,2}&=0 \\
b_{1}&=0 \\
b_{2}&=1 \\
c_{1}&=0 \\
c_{2}&=1/2
\end{align}
\)
という組で与えられます。

この係数行列の組はまとめてButcher tableと呼ばれる表記をするのが便利です。

これは、\(a_{i,j}, b_{i},c_{i}\)を
Buchertable
としてまとめて書く表記法です。

再び、オイラー法はButcher tableで書くと
euler_butcher
とまとめて書くことができます。
中点法は
midpoint_butcher
RK4は
rk4_butcher
です。

高次のルンゲ=クッタ法(10,12,14次)

4次、5次…とずっとあるわけです。
こんなページがありました[3]。
High-Order Explicit Runge-Kutta Methods
この上のページには

  • 17段10次(8次が埋め込まれてる)
  • 25段12次(10次が埋め込まれてる)
  • 35段14次(12次が埋め込まれてる)

といったButcher tableにおける係数の値が書かれています。埋め込まれてる、の意味は次の節で説明します。
ただし、上のページのbutcher tableは
highorderbutcher2
となっているので注意が必要です。

埋め込まれた陽的ルンゲ=クッタ法

“埋め込まれた“という表現が出てきたのでその説明を行いましょう。
日本語では『埋め込まれた陽的ルンゲ=クッタ法』、英語では『embedded explicit runge-kutta method』と呼ばれるものがあります。
これは、p段q次陽的ルンゲ=クッタ法を作ったら、別の次数の陽的ルンゲ=クッタ法も、係数行列\(a_{i,j}, c_{i}\)を使って作れるじゃありませんか!
というものです。

Butcher tableは、この場合extended Butcher tableと呼ばれ、こういう形式で書かれます。
embedded_rk
この埋め込まれたルンゲ=クッタ法のいいところは、

  1. 計算誤差の評価ができる
  2. 刻み幅を自動的に制御できる、適応刻み幅制御。(応用として。)

という点です。ルンゲ=クッタ法によって得られた解が真の解とどのくらい違っているのか？が評価できるんです。

例えば、4次のルンゲ=クッタ法を使って得られた解\(x^{(4)}(t)\)と5次のルンゲ=クッタ法を使って得られた解\(x^{(5)}(t)\)があったとします。
もしも、\(x^{(4)}(t)\)と解\(x^{(5)}(t)\)の解の差を調べ、その差が無かったらその数値計算解は真の解に限りなく近い、と判断することができ、差が大きかったらその解は真の解から離れていて、数値計算の精度が足らない、と判断することができます。どちらも1つだけの解では出来なかったことです。

精度が足らない場合、刻み幅を小さくすれば精度が上がります。また逆に、精度が十分に足りている場合、刻み幅を大きくし、計算時間を減らすことができます。
これが適応刻み幅制御なのです。

違った次数のルンゲ=クッタ法を、まるで別々に計算してもいいのですが、そうすると計算時間が単純に考えておおよそ2倍になります。
解を評価するために2倍の計算時間が必要というのは良くない計算効率です。
そこで考えられたのが埋め込まれたルンゲ=クッタ法なのです。

具体例を挙げましょう。
一番簡単な埋め込まれたルンゲ=クッタ法は、ホイン法と呼ばれています。
butcher_heun
1行目は2次のオーダーを持ち、2行目は1次のオーダーを持ちます。

また、4次と5次を持つ埋め込まれたルンゲ=クッタ法は、ルンゲ=クッタ=フェールベルグ(Runge-Kutta-Fehlberg)法と呼ばれています。
その埋め込まれたルンゲ=クッタ法は、
butcher_rkf45
と書かれます。1行目は5次のオーダー、2行目は4次のオーダーを持ちます。

ルンゲ=クッタ=フェールベルグ法による刻み幅の自動制御(適応刻み幅制御)

さて、次数の違う2つのルンゲ=クッタ法を用いて、適応刻み幅制御を行いたいと考えます。
刻み幅を制御するにあたって、適当に精度良かったから2倍にしてもまだ大丈夫だろ、とか差が大きすぎるから刻み幅半分にしよう、ということをやってはいけません。
適当にやったら計算時間が余計にかかり、精度が良くない変な結果が得られます。

[5~9]によれば、ルンゲ=クッタ=フェールベルグ法において区間\(i\)での最適な刻み幅\(h’\)は区間\(i\)の誤差評価の結果を使って、
\(
\displaystyle h’=\delta h=\left(\frac{\varepsilon h}{2|x^{(5)}_{i}-x^{(4)}_{i}|}\right)^{1/4} h
\)
と予想できます。ここで\(\varepsilon\)はエラーを制御する精度の目安で、おおよそ計算時に要求する相対誤差です。もちろん、この\(h’\)は区間\(i\)の最適な幅ですが、関数に劇的な変化は無いだろうとして、次の区間の計算の一番初めに用いる刻み幅を推定するのです。
なので、\(i+1\)番目の計算区間では、計算するときはこの\(h’\)の値を使えばいいんです。
(ちなみに、m次ルンゲ=クッタ法の場合では
\(
\displaystyle h’=\delta h=\left(\frac{\varepsilon h}{2|x^{(m+1)}_{i}-x^{(m)}_{i}|}\right)^{1/m} h
\)
と推測されます。)

詳しい理由は分かりませんが、5次オーダーではなく、4次です。5次のオーダーは誤差評価のためだけに用いられているようです。
ルンゲ=クッタ=フェールベルグ法の計算スキームは[7]に詳しく書かれています。
日本語訳して、その計算スキームを書けば下のようになります。


ルンゲ=クッタ=フェールベルグ法による刻み幅の自動制御(プログラム)

講義のレポート等の宿題で使うのは僕の意向と異なるので使用はご控えください。
研究目的、趣味、確かめの場合はミスがあるかもしれないことを念頭に置いたうえならば使用と改変をご自由にしてください。
このプログラムを使用して生じた責任は取りません。

fortran90によるプログラムです。ほぼ上の説明をそのままプログラミングしたものです。

  • 実数、1階微分方程式の場合

    • 微分方程式
    \(
    \displaystyle \frac{d}{dx}y=y\cos x,~~y(x=0)=1
    \)
    を倍精度実数、刻み幅制御で\(x=10\)まで解く事を考えます。
    解析解は
    \(
    \displaystyle y(x)=exp(sin(x))
    \)
    です。コードは以下の通りです。

  • 実数、2階微分方程式の場合

    • 微分方程式
    \(
    \displaystyle \frac{d^2}{dx^2}y=-\frac{1}{4}y,~~y(x=0)=1,~ y'(x=0)=0
    \)
    を倍精度実数、刻み幅制御で\(x=20\)まで解く事を考えます。コードは以下の通りです。

    ※ここで使われているサブルーチンdrkf45は実数一階微分のプログラム内にあるルーチンと一字一句同一です。

  • 複素数、1階微分方程式の場合

    • 微分方程式
    \(
    \displaystyle \frac{d}{dx}y=y\cos x,~~ y(x=0)=1+i\frac{1}{2}
    \)
    を倍精度実数、刻み幅制御で\(x=10\)まで解く事を考えます。コードは以下の通りです。

  • 複素数、2階微分方程式

    • 微分方程式
    \(
    \displaystyle \frac{d^2}{dx^2}y=-\frac{1}{4}y,~~ y(x=0)=1+i\frac{1}{2},~y'(x=0)=0+i0
    \)
    を倍精度実数、刻み幅制御で\(x=20\)まで解く事を考えます。コードは以下の通りです。

    ※ここで使われているサブルーチンcrkf45は複素数一階微分のプログラム内にあるルーチンと一字一句同一です。

等間隔の出力の場合は、以下の通りで実行できます。
サブルーチンはdrkf45は変わっていません。

[adsense2]

不連続な点を含む場合

不連続な点を含む場合、境界条件を指定しないと解くことはできません。

さて、ここで微分方程式
\(
\begin{eqnarray}
\frac{dy}{dx}=
\left\{
\begin{aligned}
0\;\;(x\le 0)\\
1\;\;(0\lt x)
\end{aligned}
\right.
\end{eqnarray}
\)
を初期条件\(y(-0.095)=0\)の下で考えます(意図的に境界条件は考えません)。
これを4次ルンゲ=クッタと適応刻み幅ルンゲ=クッタで解いてみましょう。
プログラム上ではそのまま解くことが出来ます。
実際に解かせてみますと、
rk4_rkf451_disc
となり、2つの結果(赤：4次ルンゲ=クッタ、緑：刻み幅制御ルンゲ=クッタ)は異なってしまいます。これは1階微分の不連続性のため発生します。
不連続点\(x=0\)で関数\(y(x)\)に境界条件を指定しない限り、どちらも正しい解なのです。

さて、なぜこんなことが発生するのでしょうか？以下のように問題を表すことにします。

不連続な点を含む1階の微分方程式を考えます。
ここで不連続、という意味は関数\(y(x)\)の一階微分が、点\(x’\)で
\(
\displaystyle \left. \frac{dy}{dx}\right|^{x=x’+0}_{x=x’-0}=a,\;\;\;(a\ne 0)
\)
であるような点を指しているとします。

上記の例題を考えてみましょう。
上記の例題では\(a=1\)です。微分方程式を解析的に解いてみますと、
\(
\begin{eqnarray}
y(x)=
\left\{
\begin{aligned}
C_0\;\;(x\le 0)\\
x+C_1\;\;(0\lt x)
\end{aligned}
\right.
\end{eqnarray}
\)
となります。ここで\(C_0, C_1\)は定数です。
\(C_0, C_1\)は\(y(x)\)が解きたい問題の境界条件によって決まります。

例えば、\(y(x)\)は全領域に対して繋がっている、という条件を課しましょう。この場合、不連続点\(x’\)で
\(
\displaystyle \left. y(x)\right|^{x=x’+0}_{x=x’-0}=0
\)
という境界条件を満たさなければなりません。この条件を課すと、\(C_1=C_0\)となり、初めて関数\(y(x)\)を一意に決めることが出来ます。

1階微分方程式を解く場合、適応刻み幅制御では関数\(y(x)\)は計算領域内で繋がっている事が課されています。しかし、4次ルンゲ=クッタではその条件は課されません。\(C_1\)の値は初期条件に依存し、一意に関数が決まりません。
どちらが悪いという話ではありません。

通常は適応刻み幅でも、4次ルンゲ=クッタでも\(y(x)\)にどこか連続ではない変な点がある場合、その点で区間を別々に分けて解きます。その後、境界条件に従って値を調節して全体の関数を構成します。

ベンチマーク用

微分方程式の解法がどれくらい正しそうかのベンチマーク問題として振り子(角度が大きい時)を考えましょう(振り子の詳しい解説はこちら)。
以下の\(\omega=1\)としたときの運動方程式
\(
\displaystyle \frac{d^2\theta}{dt^2}=-\sin\theta
\)
を
\(
t=0,\; \theta=0; \frac{d\theta}{dt}=1.9\cdots (k=0.95)
\)
の初期条件の下解いた場合、
1周期\(T\)は
\(
T=4K(0.95)=10.360044923498004876778\cdots
\)
となります。
この値はwolfram alphaから求めました。
4EllipticK[0.950.95] wolfram alpha

刻み幅制御を行い、45000周期目の値を考えます。45000周期目は時刻
\(
T_{45000}=466202.0215574102\cdots
\)
です。刻み幅制御による精度を\(10^{-12}\)に設定し、数値計算を行わせます。

すると実行結果として”fort.10″に

0.4662020113E+006  -0.2103356901E-001   0.1899883579E+001   0.4224109363E-002
0.4662020155E+006  -0.1300808922E-001   0.1899955473E+001   0.4223843994E-002
0.4662020198E+006  -0.4982881533E-002   0.1899993468E+001   0.4223658015E-002
0.4662020240E+006   0.3042061693E-002   0.1899997567E+001   0.4223520921E-002

というデータが出力されます。
1列目が時刻\(t\)、2列目が\(\theta(t)\),3列目が\(\frac{d\theta(t)}{dt}\),4列目が刻み幅\(h\)です。
1回のステップでの要求精度12桁に対し、最終的な結果は8桁まで正しい値を出しています。
この時、計算回数は77,852,488回\(~10^{8}\)回行われているので、最終的な結果として4桁は少なくとも正しいと考えられます。
4桁合っていればいい状況で8桁もあっているのは、被積分関数が変な振る舞いをしないからでしょう。

また、60000周期で確認してみると(\(T_{60000}=621602.695409880292\cdots\))

0.6216026888E+006  -0.1531918479E-001   0.1899938246E+001   0.4223959920E-002
0.6216026930E+006  -0.7293808996E-002   0.1899986003E+001   0.4223717084E-002
0.6216026973E+006   0.7312355417E-003   0.1899999862E+001   0.4223582630E-002
0.6216027015E+006   0.8756011575E-002   0.1899979827E+001   0.4223462029E-002

です。
1回のステップでの要求精度12桁に対し、最終的な結果は8桁まで正しい値を出しています。
この時、計算回数は103,803,513回\(~10^{8}\)回行われているので、最終的な結果として4桁は少なくとも正しいと考えられます。
･･･まだまだ大丈夫そうですね。

少し特殊な初期条件(\(k=1\))でやってみましょう。
この\(k=1\)は、振り子の質点がちょうど真上に来て静止する非常に不安定な状態です。
何秒間静止していられるか試してみましょう。刻み幅の制御等は上記条件と同じです。
横軸に時間\(t\)、縦軸に\(\theta(t)\)を取った時のグラフです。
ellk1
すぐに破綻しました。正しい値は\(t=20\)位までですかね。これは、機械的な誤差があることによって不安定な平衡点からほんの少しだけ外れてしまったからです。だからカオスとかを考えるときとかは注意しなければなりません。

ルンゲ=クッタ=ドルマンド=プリンス法

フェールベルグ法は初期の頃に考えられた埋め込まれた方法です。
研究が進み、実用上では別の係数の組み合わせが良いことが分かってきました。
その一つが1980年に発見されたドルマンド=プリンス(Dormand-Prince)法です。

これは、7段4(5)次の方法です。
フェールベルグ法は6段4(5)次の方法ですので、次数は同じです。

良く調べていませんが、この違いは、4次の結果を基準にして求めたブッチャー係数(フェールベルグ法)か、5次の結果を基準に求めたブッチャー係数(ドルマンド=プリンス法)か？の違いのようです。

単純に考えて、同じ次数なのにドルマンド=プリンス法の方が段数が増えていて効率が悪いです。
しかし、本来は7段なのですが、7段目に呼び出した結果を取って置けば、次のステップの1段目に同じ値が使えるように設計されているので、プログラム上は6段と(ほぼ)同じ関数の呼び出し回数になります。

プログラムはこんな感じになるでしょう。

適当な刻み幅で出力

等間隔(サブルーチンは上のものと同じなので省略)

program main
  implicit none
  integer::i,N,info,Nx
  double precision::h,tol
  double precision::xa,xb,tx
  double precision,allocatable::y(:),x(:),work0(:)
  double precision,external::grk

N=1
  allocate(y(1:N),work0(1:N))

Nx=101
  allocate(x(1:Nx))
  xa=0d0
  xb=10d0
  do i=1,Nx
     x(i)=(i-1)*(xb-xa)/dble(Nx-1)+xa
  enddo

!initial conditions
  y(1)=1d0  ! x (0)

tol=1d-8
  write(10,'(2e25.10e3)')x(1),y(1)
  do i=2,Nx
     info=-1
     h=x(i)-x(i-1)
     tx=x(i-1)
     do while(info.le.0)
        call dDP45(grk,tx,h,N,y,x(i),info,tol,work0)
     enddo
     write(10,'(2e25.10e3)')x(i),y(1)
  enddo

stop
end program main

function grk(N,x,y,s)
  implicit none
  integer,intent(in)::N,s
  double precision,intent(in)::x
  double precision,intent(in)::y(1:N)
  double precision::grk

grk=0d0
  if(s.eq.1)then
     grk=y(1)<em>cos(x)
  else
     write(6,</em>)"***Error grk"; stop
  endif

return
end function grk

4倍精度ルーチン

4倍精度のサブルーチンです。
計算速度は倍精度の30~50倍かかるので、必要なとき以外使わないようにしましょう。

陽的ルンゲ=クッタ法の導出

ルンゲ=クッタ法の導出は煩雑です。単に複雑なだけです。
導出過程について詳しく述べられているページは、早川様が記述された以下のpdfを参照すると良いと思います。
Runge-Kutta法についてのノート(早川尚男)
計算過程を含め記述されているので分かりやすいです。

参考文献

[1]Derivation of Runge–Kutta methods
[2]NDSolveの”ExplicitRungeKutta”メソッド
[3]High-Order Explicit Runge-Kutta Methods
[4]List of Runge–Kutta methods
[5]Runge-Kutta-Fehlberg Method (RKF45)
[6]Runge-Kutta-Fehlberg method
[7]Lecture:13Runge-Kutta-Fehlberg Method
[8]GPU acceleration of Runge Kutta-Fehlberg and its comparison with Dormand-Prince method
[9]William H. Pressら著『ニューメリカルレシピ・イン・シー 日本語版―C言語による数値計算のレシピ– 』(1993)

ニューメリカルレシピ・イン・シー 日本語版―C言語による数値計算のレシピ

↑この本が一番有益だと思います。Fortran版もありますが、英語しかありません。ちなみに、英語で,若干古いバージョンでいいのならば
公式のホームページ
Numerical Recipes in C (1992)
Numerical Recipes in Fortran 77 and Fortran 90 (1992, 1996)
で無料で公開されています。

プログラミングと数値計算

数値積分(等間隔)

コメントする

数値計算での積分方法、特に等間隔の分点の場合であるニュートン・コーツ積分(とロンバーグ積分)に関する理論とのプログラムを載せます。

  1. ニュートンコーツ型の公式
  2. ルンゲ現象
  3. ロンバーグ積分
  4. サブルーチン”integral”のコードと使い方
  5. サブルーチン”romberg”のコードと使い方

等間隔の分点で積分を行う場合、よく使われる方法は台形則、もしくはシンプソン積分です。

台形則

台形則は、分点間を台形近似して面積を求める方法であり、下図のようなイメージで積分を近似する方法です。
外形_台形_c
積分\(\displaystyle \int_a^b f(x) dx\)を、
\(
\displaystyle \int_a^b f(x) dx \sim \sum_{i=0}^{N-1} \frac{f(x_{i+1})+f(x_i)}{2}h
\)
として近似します。この時、誤差のオーダーは\(O(h^3)\)となります。
言い換えると、
台形積分とは、隣り合う分点間を1次関数で近似して求める積分
言えます。

シンプソン積分

シンプソン積分は、隣り合う分点間を2次関数で近似して求める積分です。イメージは
外形_シンプソン_c
な感じです。数式では
\(
\displaystyle \int_a^b f(x) dx \sim \sum_{i=0}^{N/2} \frac{f(x_{2i})+4f(x_{2i+1})+f(x_{2i+2})}{3}h
\)
となります。

高次≠高精度

これ以降もあります。分点間を3次、4次、…としていけば高次の積分をすることが可能になります。
この分点間の次数をどんどん上げていって積分をする方法をまとめてニュートンコーツの公式、と呼びます。
例えばの台形則はニュートンコーツの公式の1次に相当し、シンプソン則はニュートンコーツの公式の2次に相当しています。

じゃあ100次の公式作れば精度が凄い公式ができるんじゃないか？となりますが、これは違います。
高次であることと精度が高いことは違うのです。だから使うな、というわけではなくて知っておいてほしいことなんです。
ではなぜか？これはwikipediaのニュートンコーツの公式のページからとってきたものですが、この式に理由があります。
ニュートンコーツ型式_c

問題となるのが誤差項です。例えば台形近似を見ますと、誤差項は\(h^3 f^{(2)}(\xi)\)とあります。hは刻み幅なので、この項は刻み幅を小さくすれば小さくするほど精度がよくなることを言っています。
が、問題は\(f^{(2)}(\xi)\)の方です。この項は本当の関数が分からないと評価できません。
多項式だったら微分するごとに次数が1つづつ減少するのでこの項は高次になれば小さくなります。じゃあもしも微分するごとに値が増加していったら…？この場合、hを小さくしてもそれを上回るくらい\(f^{(2)}(\xi)\)が大きくなることがあります。これをルンゲ現象と呼びます。

ルンゲ現象

例えば関数
\(
\displaystyle f(x)=\frac{1}{1+25x^2}
\)
を考えます。wikipediaのページを参考にすれば、その導関数と値は、
ルンゲ現象_c
となります。高次になればなるほど値が増加する関数なのです。
実際には\(h^n f^{(m)}(\xi)\)の積の値が小さいか大きいか？なので、一概には高次は悪い！とは言えません。高次は危ないくらいでしょう。

これを可能な限り回避するために、低次の積分を組み合わせて使う方法がとられます。
こうすることで高次の値を増加を抑えるのです。

そのほかの手法もいろいろあります。分点の位置に縛られなければガウス求積法が最高でしょう。素晴らしい方法です。ぜひ調べてみてください。

[adsense1]

ロンバーグ積分

さて、高次を追い求めるのもいいですが、一線を駕す…とまではいきませんが、ロンバーグ(Romberg)積分というものがあります。
このロンバーグ積分とは、台形則による計算を基本とし、無限小区間での積分結果を補外によって求めよう、というものです。
ここでいう補外の概念は以下のような感じです。

刻み幅1で台形則の結果から積分値S1が分かる
↓
刻み幅0.5で台形則による結果から積分値S2がわかる
↓
刻み幅0.25で台形則による結果から積分値S3がわかる
↓
刻み幅0.175で台形則による結果から積分値S4がわかる

とある程度計算していきます。そうすると刻み幅hを変化させるとそれに応じて積分値がS1→S2→S3→S4…と推移していくわけです。
この推移の仕方を計算すれば、刻み幅無限小の場合の積分値Sが分かる、こういう仕組みになっています。収束加速法と言ったり、リチャードソン補外、なんて言ったりします。
この補外による方法はかなり優れています。先ほどの高次がどうこう、といった問題が発生しないのです。台形則しか使ってないんですから。有限桁計算に適した方法である、とかwikipediaに載ってたりします。

ロンバーグ積分の詳しい計算方法や理論は
ニューメリカルレシピ・イン・シー 日本語版―C言語による数値計算のレシピ

ニューメリカルレシピ・イン・シー 日本語版―C言語による数値計算のレシピ

を見ると詳しく書かれています。

日本語では第6章 数値積分と数値微分にロンバーグ積分の記述があります。
英語ではRomberg integralで調べると良いでしょう。参考までに、英語では
Romberg integration exampleや、Romberg Integrationなにかがいいと思います。

2016/03/01
上のRomberg Integrationを参考にして書いたプログラムはこちらです。これは、解析関数fをRomberg積分するプログラムです。
解析関数ではなく、サイズ\(2^{n}+1\)の配列に格納されている場合はこのページの下の方にあるプログラムを参照してください。

call romberg(a,b,s,pre)

で呼び出し、区間a~bの積分\(s=\int_a^b f(x) dx\)を精度preで求めるプログラムです。

サブルーチン”integral”のコードと使い方

下のサブルーチン”integral”は1,2,3次元の短冊近似、台形則、シンプソン積分、シンプソン3/8積分、ブール則(ニュートン・コーツ型積分)をカバーします。
ただし、分点の個数に注意してください。分点の個数を気にせず使える近似方法は、短冊近似、台形近似のみです。

ロンバーグ積分では引数の関係上、ニュートン・コーツ型積分との間をうまく処理できなかったため、別のサブルーチンとして書きます。もしもロンバーグ積分を使いたい場合はもう少し先に進んでください。

まずはニュートン・コーツ型積分です。
下のサブルーチンを使ってください。
使い方は、

 call integral(配列の大きさ,配列,刻み幅,積分結果,積分手法)

という感じです。
実際のプログラムでは、以下のように使用してください。配列yは複素数解列(ただし実軸上の積分)でもokです。

1次元

サブルーチンintegralの使い方
積分則 呼び方
短冊近似(長方形近似)
 call integral(size(w,1),w,h,s,"box")
 配列yの大きさに指定はない。
台形近似
 call integral(size(y,1),y,h,s,"trapezoid")
 配列yの大きさに指定はない。
シンプソン則
 call integral(size(y,1),y,h,s,"simpson")
 配列yの大きさが3,5,7,…個、\(2n+1\)個でないといけない
シンプソン3/8則
 call integral(size(y,1),y,h,s,"simpson38")
 配列yの大きさが4,7,10,…個、\(3n+1\)個でないといけない
ブール則
 call integral(size(y,1),y,h,s,"boole")
 配列yの大きさが5,9,13,…個、\(4n+1\)個でないといけない

2次元配列の場合は以下のように指定してください。

 call integral(size(w,1),size(w,2),w,hx,hy,s,"simpson")

3次元の場合はこう指定してください。

 call integral(size(w,1),size(w,2),size(w,3),w,hx,hy,hz,s,"simpson")

積分のためのモジュール

module integral_mod
  !developer --> sikino
  !date --> 2015/04/07
  implicit none
  interface integral
     module procedure &
          dintegral, &
          dintegral2d, &
          dintegral3d, &
          cintegral, &
          cintegral2d, &
          cintegral3d
  end interface integral
contains
  subroutine dintegral(N,y,h,s,method)
    integer,intent(in)::N
    double precision,intent(in)::h,y(1:N)
    character(*),intent(in)::method
    double precision,intent(out)::s
    integer::i
    double precision::y0,y1,y2,y3
       
    s=0.d0; y0=0.d0; y1=0.d0; y2=0.d0; y3=0.d0
   
    if(trim(method).eq."box")then
       s=h*sum(y(1:N-1))
    elseif(trim(method).eq."trapezoid")then
       y1=y(1)+y(N)
       do i=2,N-1
          y2=y2+y(i)
       enddo
       s=(y1+2.d0*y2)*h*0.5d0
    elseif(trim(method).eq."simpson")then
       if(mod(N,2).ne.1)then
          write(6,*)"=====cannot calculation with simpson"
          write(6,*)"=====program stop"
          stop
       endif

       y1=y(1)+y(N)
       do i=2,N-1,2
          y2=y2+y(i)
       enddo
       do i=3,N-2,2
          y3=y3+y(i)
       enddo

       s=(y1+4.d0*y2+2.d0*y3)*h/3.d0

    elseif(trim(method).eq."simpson38")then

       if(mod(N,3).ne.1)then
          write(6,*)"=====cannot calculation with simpson38"
          write(6,*)"=====program stop"
          stop
       endif

       y0=y(1)+y(N)
       do i=2,N-2,3
          y1=y1+y(i)
       enddo
       do i=3,N-1,3
          y2=y2+y(i)
       enddo
       do i=4,N-3,3
          y3=y3+y(i)
       enddo
       s=(y0+3.d0*(y1+y2)+2.d0*y3)*3.d0*h/8.d0        

    elseif(trim(method).eq."boole")then
       if(mod(N,4).ne.1)then
          write(6,*)"=====cannot calculation with boole"
          write(6,*)"=====program stop"
          stop
       endif

       y0=y(1)+y(N)
       do i=5,N-4,4
          y1=y1+y(i)
       enddo
       do i=2,N-1,2
          y2=y2+y(i)
       enddo
       do i=3,N-2,4
          y3=y3+y(i)
       enddo

       s=(14.d0*y0+28.d0*y1+64.d0*y2+24.d0*y3)*h/45.d0
    else
       write(6,*)"=====cannot calculation in integral"
       write(6,*)"=====program stop"
       stop
    end if

    return
  end subroutine dintegral
 
  subroutine dintegral2d(Nx,Ny,z,hx,hy,s,method)
    implicit none
    integer,intent(in)::Nx,Ny
    double precision,intent(in)::hx,hy,z(1:Nx,1:Ny)
    character(*),intent(in)::method
    double precision,intent(out)::s
    integer::i
    double precision::ty(1:Ny),r(1:Nx)

    s=0.d0
    ty(1:Ny)=0.d0
    r(1:Nx)=0.d0
    do i=1,Nx
       ty(1:Ny)=z(i,1:Ny)
       call integral(Ny,ty,hy,s,method)
       r(i)=s
    enddo
    call integral(Nx,r,hx,s,method)
       
    return
  end subroutine dintegral2d
 
  subroutine dintegral3d(Nx,Ny,Nz,w,hx,hy,hz,s,method)
    implicit none
    integer,intent(in)::Nx,Ny,Nz
    double precision,intent(in)::hx,hy,hz,w(1:Nx,1:Ny,1:Nz)
    character(*),intent(in)::method
    double precision,intent(out)::s
    integer::i
    double precision::tyz(1:Ny,1:Nz),r(1:Nx)
       
    s=0.d0
    tyz(1:Ny,1:Nz)=0.d0
    r(1:Nx)=0.d0
    do i=1,Nx
       tyz(1:Ny,1:Nz)=w(i,1:Ny,1:Nz)
       call integral(Ny,Nz,tyz,hy,hz,s,method)
       r(i)=s
    enddo
    call integral(Nx,r,hx,s,method)    
       
    return
  end subroutine dintegral3d
 
  subroutine cintegral(N,y,h,s,method)
    integer,intent(in)::N
    complex(kind(0d0)),intent(in)::y(1:N)
    double precision,intent(in)::h
    character(*),intent(in)::method
    complex(kind(0d0)),intent(out)::s

    double precision::res,ims

    s=dcmplx(0d0,0d0); res=0.d0; ims=0.d0

    call integral(N,dble(y),h,res,trim(method))
    call integral(N,dimag(y),h,ims,trim(method))
   
    s=dcmplx(res,ims)
   
    return
  end subroutine cintegral

  subroutine cintegral2d(Nx,Ny,z,hx,hy,s,method)
    integer,intent(in)::Nx,Ny
    complex(kind(0d0)),intent(in)::z(1:Nx,1:Ny)
    double precision,intent(in)::hx,hy
    character(*),intent(in)::method
    complex(kind(0d0)),intent(out)::s

    double precision::res,ims

    s=dcmplx(0d0,0d0); res=0.d0; ims=0.d0

    call integral(Nx,Ny,dble(z),hx,hy,res,trim(method))
    call integral(Nx,Ny,dimag(z),hx,hy,ims,trim(method))
   
    s=dcmplx(res,ims)
   
    return
  end subroutine cintegral2d

  subroutine cintegral3d(Nx,Ny,Nz,w,hx,hy,hz,s,method)
    integer,intent(in)::Nx,Ny,Nz
    complex(kind(0d0)),intent(in)::w(1:Nx,1:Ny,1:Nz)
    double precision,intent(in)::hx,hy,hz
    character(*),intent(in)::method
    complex(kind(0d0)),intent(out)::s

    double precision::res,ims

    s=dcmplx(0d0,0d0); res=0.d0; ims=0.d0

    call integral(Nx,Ny,Nz,dble(w),hx,hy,hz,res,trim(method))
    call integral(Nx,Ny,Nz,dimag(w),hx,hy,hz,ims,trim(method))
   
    s=dcmplx(res,ims)
   
    return
  end subroutine cintegral3d
end module integral_mod

サブルーチン”integral”を用いた例題

実際の使い方。参考にどうぞ。
下のコードは2次元の積分
\(
\displaystyle \int\int xe^{-x^2}e^{-y}dxdy=\frac{e^{-x^2}}{2}e^{-y}+\mbox{const}
\)
を数値積分します。積分範囲を\(x=0\sim 3, \; y=0 \sim 5\)にした場合、その解析解は
\(
\begin{align}
\int_0^3 dx\int_0^5 dy xe^{-x^2}e^{-y}&=\frac{1}{2}(1-e^{-5}-e^{-9}+e^{-14})\\
&=0.4965697373627734784608751894356320535936864348993604\cdots
\end{align}
\)
となります。実行すると

./a.out
simpson   :    0.496569737366759E+000
romberg,3 :    0.496569737362773E+000

となります。コードは、

program main
  use integral_mod
  implicit none
  integer::i,j,n,Nx,Ny
  double precision::hx,hy,ans
  double precision::xmin,xmax,ymin,ymax
  double precision,allocatable::x(:),y(:)
  double precision::s
  double precision,allocatable::w(:,:)
  double precision::f
  external::f
 
  n=10

  xmin=0d0
  xmax=3d0

  ymin=0d0
  ymax=5d0

  Nx=2**(n)
  Ny=2**(n)
  allocate(x(0:Nx),y(0:Ny),w(0:Nx,0:Ny))
  x=0d0; y=0d0; w=dcmplx(0d0,0d0)

  hx=(xmax-xmin)/dble(Nx)
  hy=(ymax-ymin)/dble(Ny)
 
  do i=0,Nx
     x(i)=xmin+hx*dble(i)
  enddo
  do i=0,Ny
     y(i)=ymin+hy*dble(i)
  enddo

  do i=0,Nx
     do j=0,Ny
        w(i,j)=f(x(i),y(j))
     end do
  end do

  call integral(size(w,1),size(w,2),w,hx,hy,s,"simpson")
  write(6,'(A,e25.15e3)')"simpson   : ",s
 
  call integral(size(w,1),size(w,2),w,hx,hy,s,"romberg",3)
  write(6,'(A,e25.15e3)')"romberg,3 : ",s
 
  return
end program main

function f(x,y)
  implicit none
  double precision,intent(in)::x,y
  double precision::f
 
  f=x*exp(-x*x)*exp(-y)
 
  return
end function f

3次元の計算の例題です。
下のコードは
\(
\displaystyle \int\int\int x^4y^4z^4 dxdydz=\left(\frac{1}{5}\{0.7^5-(-1)^5\}\right)^3
\)
を数値積分します。積分範囲を\(x,y,z=-1\sim 0.7\)にした場合、その値は
\(
\displaystyle \int_{-1}^{0.7}\int_{-1}^{0.7}\int_{-1}^{0.7} x^4y^4z^4 dxdydz = 0.127496010896795\cdots
\)
となります。実行すると、

$ gfortran integralmod.f90 main.f90
$ ./a.out
analysis  :   0.127496010896795E-01
simpson   :   0.10507613E-006  0.10507613E-006
romberg,3 :  -0.17347235E-017 -0.17347235E-017

と得られます。中身はこうです。分点の数はx,y,z軸各々で違っていて構いません。

program main
  use integral_mod
  implicit none
  integer::i,j,k,n,Nx,Ny,Nz
  double precision::hx,hy,hz,ans
  double precision,allocatable::x(:),y(:),z(:)
  complex(kind(0d0))::s
  complex(kind(0d0)),allocatable::w(:,:,:)
 
  n=6
  Nx=2**(n)
  Ny=2**(n-1)
  Nz=2**(n+1)
 
  allocate(x(0:Nx),y(0:Ny),z(0:Nz),w(0:Nx,0:Ny,0:Nz))
  x=0d0; y=0d0; z=0d0; w=dcmplx(0d0,0d0)
 
  hx=1.7d0/dble(size(w,1)-1)
  hy=1.7d0/dble(size(w,2)-1)
  hz=1.7d0/dble(size(w,3)-1)
 
  do i=0,Nx
     x(i)=hx*dble(i)-1d0
  enddo
  do i=0,Ny
     y(i)=hy*dble(i)-1d0
  enddo
  do i=0,Nz
     z(i)=hz*dble(i)-1d0
  enddo

  do i=0,Nx
     do j=0,Ny
        do k=0,Nz
           w(i,j,k)=dcmplx((x(i)*y(j)*z(k))**4,(x(i)*y(j)*z(k))**4)
        end do
     end do
  end do

  ans=(0.7d0**5.d0+1.d0)/5.d0
  ans=ans**3.d0
  write(6,'(A,e23.15e2)')"analysis  : ",ans
 
  call integral(size(w,1),size(w,2),size(w,3),w,hx,hy,hz,s,"simpson")
  write(6,'(A,2e17.8e3)')"simpson   : ",dble(s)-ans,dimag(s)-ans
 
  call integral(size(w,1),size(w,2),size(w,3),w,hx,hy,hz,s,"romberg",3)
  write(6,'(A,2e17.8e3)')"romberg,3 : ",dble(s)-ans,dimag(s)-ans
 
  return
end program main

ロンバーグ積分のプログラム

ロンバーグ積分
 call romberg(jx,x,y,s,pre)
 配列yの大きさは2,3,5,10,…,\(2^{jx}+1\)個でないといけない ロンバーグ積分では収束精度”pre”を指定すること。もしも収束精度に達しなかった場合、警告と共に、与えられた値での収束限界の積分結果を返す。積分精度を高めたければ刻み幅を小さくする操作(例えば分点数を増やす等)をしてください。”pre”は、積分結果の絶対値が1より大きくなる場合は相対誤差を、小さくなる場合は絶対誤差を取ります。この理由は以前の結果の補正を加え、収束をさせるためであり、非常に小さい誤差の場合いつまでたっても機械誤差のため収束しないからです。
module romberg_mod
  !developer --> sikino
  !date --> 2016/03/01
  !         2016/03/03
  !
  ! 1D case :
  !  romberg(jx,x,y,s,pre)
  !          |  | | | +- (in)  precision (e.g. 1d-8)
  !          |  | | +--- (out) integration result
  !          |  | +----- (in)  y(1:2**jx+1) f(x)
  !          |  +------- (in)  x(1:2**jx+1) x
  !          +---------- (in)  related to array size
  !
  ! 2D case :
  !  romberg(jx,jy,x,y,z,s,pre)
  !          |  |  | | | | +-- (in)  precision (e.g. 1d-8)
  !          |  |  | | | +---- (out) integration result
  !          |  |  | | +------ (in)  z(1:2**jx+1,1:2**jy+1) f(x,y)
  !          |  |  | +-------- (in)  y(1:2**jy+1) y
  !          |  |  +---------- (in)  x(1:2**jx+1) x
  !          |  +------------- (in)  related to array size y
  !          +---------------- (in)  related to array size x
  !
  ! 3D case :
  !  romberg(jx,jy,jz,x,y,z,w,s,pre)
  !          |  |  |  | | | | | +-- (in)  precision (e.g. 1d-8)
  !          |  |  |  | | | | +---- (out) integration result
  !          |  |  |  | | | +------ (in)  w(1:2**jx+1,1:2**jy+1,1:2**jz+1) f(x,y,z)
  !          |  |  |  | | +-------- (in)  z(1:2**jz+1) z
  !          |  |  |  | +---------- (in)  y(1:2**jy+1) y
  !          |  |  |  +------------ (in)  x(1:2**jx+1) x
  !          |  |  +--------------- (in)  related to array size z
  !          |  +------------------ (in)  related to array size y
  !          +--------------------- (in)  related to array size x
  !
  implicit none
  interface romberg
     module procedure &
          dromberg, &
          dromberg2d, &
          dromberg3d
  end interface romberg
contains
  subroutine dromberg(jx,x,y,s,pre)
    implicit none
    integer,intent(in)::jx
    double precision,intent(in)::x(1:2**jx+1),y(1:2**jx+1),pre
    double precision,intent(out)::s

    integer,parameter::jm=6 !--> precision: O(h^(2*jm))
    integer,parameter::nm=2**jm
    double precision,allocatable::tx(:),ty(:)
    double precision::ts
    integer::k

    s=0d0
    if(jx.ge.jm)then
       allocate(tx(1:nm+1),ty(1:nm+1)); tx=0d0; ty=0d0
       do k=0,2**(jx-jm)-1
          ts=0d0
          tx(1:nm+1)=x(k*nm+1:(k+1)*nm+1)
          ty(1:nm+1)=y(k*nm+1:(k+1)*nm+1)

          call romberg_sub(jm,tx,ty,ts,pre)
          s=s+ts
       enddo
       deallocate(tx,ty)    
    else
       call romberg_sub(jx,x,y,s,pre)
    endif

    return
  end subroutine dromberg

  subroutine romberg_sub(jx,x,y,s,pre)
    ! reference "http://www.math.usm.edu/lambers/mat460/fall09/lecture29.pdf"
    implicit none
    integer,intent(in)::jx
    double precision,intent(in)::x(1:2**jx+1),y(1:2**jx+1),pre
    double precision,intent(out)::s

    integer::i,j,k,n,dn
    double precision::h,ps,tmp
    double precision::T(1:jx+1,1:jx+1)

    n=2**jx+1

    h=x(n)-x(1)
    dn=(n-1)/2

    T(1,1)=0.5d0*h*(y(1)+y(n))
    s=T(1,1)
    ps=s
    h=0.5d0*h
    do j=2,jx+1

       ! trapezoidal rule
       tmp=0d0
       do i=1,2**(j-2)
          tmp=tmp+y(1+(2*i-1)*(dn))
       enddo
       T(j,1)=0.5d0*T(j-1,1)+h*tmp

       do k=2,j
          ! Richardson extrapolation
          T(j,k)=T(j,k-1)+(T(j,k-1)-T(j-1,k-1))/(dble(4**(k-1))-1d0)
       enddo
       s=T(j,j)

       ! precision check
       if(abs(s).ge.1d0)then
          if(abs((ps-s)/s).le.pre)exit
       else
          if(abs((ps-s)).le.pre)exit
       endif
       ps=s
       h=0.5d0*h
       dn=dn/2
    enddo

    if(j-1.eq.jx)then
       write(6,'(A)')" -+-+- didn't converge at romberg integral -+-+- "
       write(6,'(A)')"       Please change stepsize h of array x  "
    endif

    return
  end subroutine romberg_sub

  subroutine dromberg2d(jx,jy,x,y,z,s,pre)
    implicit none
    integer,intent(in)::jx,jy
    double precision,intent(in)::pre,x(1:2**jx+1),y(1:2**jy+1)
    double precision,intent(in)::z(1:2**jx+1,1:2**jy+1)
    double precision,intent(out)::s

    integer::i,nx,ny
    double precision::ty(1:2**jy+1),r(1:2**jx+1)
    nx=2**jx+1
    ny=2**jy+1

    s=0.d0
    ty(1:ny)=0.d0
    r(1:nx)=0.d0
    do i=1,nx
       ty(1:ny)=z(i,1:ny)
       call romberg(jy,y,ty,s,pre)
       r(i)=s
    enddo
    call romberg(jx,x,r,s,pre)
       
    return
  end subroutine dromberg2d
 
  subroutine dromberg3d(jx,jy,jz,x,y,z,w,s,pre)
    implicit none
    integer,intent(in)::jx,jy,jz
    double precision,intent(in)::pre,x(1:2**jx+1),y(1:2**jy+1),z(1:2**jz+1)
    double precision,intent(in)::w(1:2**jx+1,1:2**jy+1,1:2**jz+1)
    double precision,intent(out)::s
    integer::i,nx,ny,nz
    double precision::tyz(1:2**jy+1,1:2**jz+1),r(1:2**jx+1)
       
    nx=2**jx+1; ny=2**jy+1; nz=2**jz+1
    s=0.d0
    tyz(1:ny,1:nz)=0.d0
    r(1:nx)=0.d0
    do i=1,nx
       tyz(1:ny,1:nz)=w(i,1:ny,1:nz)
       call romberg(jy,jz,y,z,tyz,s,pre)
       r(i)=s
    enddo
    call romberg(jx,x,r,s,pre)  
       
    return
  end subroutine dromberg3d
end module romberg_mod

ロンバーグ積分(romberg_mod)の例

必要なのは、上で紹介したモジュール “romberg_mod” と下のメインプログラムです。
以下のプログラムは1次元の定積分
\(
\int_1^10 \frac{1}{x^2} dx = 0.9
\)
を分点数\(2^8+1(jx=8)\)個でロンバーグ積分するものです。
精度は\(O(h^{2\cdot min{jm,jx}})\)です。
ここで\(jm\)はモジュールromberg_modの中にパラメータとして宣言されています。

コンパイルは

ifort romberg_mod.f90 main.f90

とでもすればいいでしょう。

program main
  use romberg_mod
  implicit none
  integer::i,jx,n
  double precision::h,a,b,s
  double precision,allocatable::x(:),y(:)
  double precision::f
  external::f
 
  a=1d0
  b=10d0
  jx=8
  n=2**jx+1
  allocate(x(1:n),y(1:n))
  h=(b-a)/dble(n-1)
  do i=1,n
     x(i)=a+h*dble(i-1)
     y(i)=f(x(i))
  enddo
 
  call romberg(jx,x,y,s,1d-8)
  write(6,*)s,"romberg"
 
  stop
end program main

function f(x)
  implicit none
  double precision::f,x
 
  f=1d0/(x*x)
 
  return
end function f

実行例(要求精度は\(10^{-8}\))

>./a.out
  0.900000000062642  
>

2次元ではこう。

program main
  use romberg_mod
  implicit none
  integer::i,j,Nx,Ny,jx,jy
  double precision::hx,hy
  double precision::xmin,xmax,ymin,ymax
  double precision,allocatable::x(:),y(:)
  double precision::s
  double precision,allocatable::w(:,:)
  double precision::f
  external::f
 
  jx=10
  jy=8
  xmin=0d0
  xmax=3d0

  ymin=0d0
  ymax=5d0

  Nx=2**(jx)
  Ny=2**(jy)
  allocate(x(0:Nx),y(0:Ny),w(0:Nx,0:Ny))
  x=0d0; y=0d0; w=dcmplx(0d0,0d0)

  hx=(xmax-xmin)/dble(Nx)
  hy=(ymax-ymin)/dble(Ny)
 
  do i=0,Nx
     x(i)=xmin+hx*dble(i)
  enddo
  do i=0,Ny
     y(i)=ymin+hy*dble(i)
  enddo

  do i=0,Nx
     do j=0,Ny
        w(i,j)=f(x(i),y(j))
     end do
  end do

  call romberg(jx,jy,x,y,w,s,1d-8)
  write(6,*)s
 
  stop
end program main

function f(x,y)
  implicit none
  double precision,intent(in)::x,y
  double precision::f
 
  f=x*exp(-x*x)*exp(-y)
 
  return
end function f

実行例

>./a.out
  0.496569737374163
>

[adsense2]

プログラミングと数値計算

ファイル読み込み

コメントする

ファイルの読み込みに関するサブルーチンを記述します。

 do i=1,4
     do j=1,83
        write(30,*)i,j,i*j
     end do
     write(30,*)
  end do

によってファイル”fort.30″が作られたとします。
今、fort.30を読み込んで
“i”の数”4”
と
“j”の数”83″
を取り出したいとします。

read文を使って読み込みますが、そのまま読み込むと空白部分を読み込んだり、読み込まないだったりします。
ここでは上のdoループによって作られたファイルの場合に使えるものを想定します。

ここでは、
大きな塊を表す数である”4″をblock,
塊の中の要素数を表す”83″をelement
と呼ぶことにします。

この問題を解く戦略は、ファイルを読み込む2種類の方法です。
1つは本当の行数(空行含む)を数え、もう1つは空行を飛ばして読み込む方法です。
本当の行数を数えるサブルーチンは下の”linecount”
であり、
空行を読み飛ばして行数を数えるサブルーチンは下の”linecount_eff”
です。
この二つと、一番下の行に追加される余分な１行を対処するために”breaklinecheck”というルーチンを使います。

これによってblockとelementを出力するサブルーチン”blockelement”を作っています。

下のプログラムを実行すると

$ gfortran main.f90
$ ./a.out
 ===Apply for fortran file will done===
 Nblock     ::         4
 Nelement   ::        83

という結果が得られるかと思います。

program main
  !developer => sikinote
  !date      => 2015/03/31
  implicit none
  integer::NBlock,Nelement
  character(48)::filename
 
  filename='./fort.30'
  call blockelement(filename,Nblock,Nelement)

  write(6,'(A,i10)')" Nblock     ::",Nblock
  write(6,'(A,i10)')" Nelement   ::",Nelement

  stop
end program
!===================================
   
subroutine blockelement(filename,Nblock,Nelement)
  !developer => sikino
  !date      => 2015/03/31
  implicit none
  character(*),intent(in)::filename
  integer,intent(out)::Nblock,Nelement
  integer::c1,c2
 
  call linecount(c1,filename)
  call linecount_eff(c2,filename)
  call breaklinecheck(c1,c2)
  Nblock=c1-c2+1
  Nelement=c2/Nblock
 
  return
end subroutine blockelement
!------------------------------
subroutine linecount(c,filename)
  implicit none
  integer,intent(out)::c
  character(*),intent(in)::filename

  integer::ier
  character(len_trim(filename))::fname
 
  fname=trim(filename)
  c=0
  open(100,file=fname,status='old',iostat=ier,err=990)
  do while(.true.)
     read(100,*,end=999)
     c=c+1
  enddo
999 continue
  close(100)
   
  return
!+-----------------------+
990 write(6,'(A)')"!!!!==error when open file",trim(fname),"info-->",ier
  write(6,*)"======program stop at linecount"
  stop

end subroutine linecount
!--------------------------------
subroutine linecount_eff(c,filename)
  implicit none
  integer,intent(out)::c
  character(*),intent(in)::filename

  integer::ier
  character(100)::cc
  character(len_trim(filename))::fname
 
  fname=trim(filename)

  c=0
  open(100,file=fname,status='old',iostat=ier,err=990)
  do while(.true.)
     read(100,*,end=998)cc
     if(len_trim(cc).gt.0)c=c+1
  enddo

998 continue
  close(100)
  return

990 write(6,'(A)')"!!!!==error when open file",trim(fname),"info==>",ier
  write(6,*)"======program stop at linecount_eff"
  stop

end subroutine linecount_eff
!-------------------------------------
subroutine breaklinecheck(c1,c2)
  implicit none
  integer,intent(inout)::c1
  integer,intent(in)::c2
  integer::Nb
 
  Nb=c1-c2+1
  if(Nb.eq.2.or.mod(c2,Nb).ne.0)then
     write(6,*)"===Apply for fortran file will done==="
     c1=c1-1
     Nb=c1-c2+1
     if(mod(c2,Nb).ne.0)then
        write(6,*)"line is different(may be last break)"
        write(6,*)"program stop at subroutine __breaklinecheck__"
        stop
     end if
  endif
 
  return
end subroutine breaklinecheck

データを読み込むには？

さて、ブロックの数と要素の数が上のサブルーチンを使うことにより求められることがわかりました。
実際にデータを配列に代入するためにはどうすればいいんでしょう？
型に応じて使うサブルーチンを変えます。
その手続きは下のモジュールを記述することでokです。これを書いた上で、
メインプログラムを以下のように書きます。そうすれば配列x(,)とy(,)に値がはいります。

program main
  use read1
  implicit none
  integer::NBlock,Nelement,i,j
  character(48)::filename
 
  double precision,allocatable::x(:,:),y(:,:)
 
  filename='./fort.30'
  call blockelement(filename,Nblock,Nelement)

  write(6,'(A,i10)')" Nblock     ::",Nblock
  write(6,'(A,i10)')" Nelement   ::",Nelement
 
  allocate(x(1:Nblock,1:Nelement),y(1:Nblock,1:Nelement))
  call read_filedata(size(y,1),size(y,2),x,y,filename)
 
  do i=1,Nblock
     do j=1,Nelement
        write(11,*)x(i,j),y(i,j)
     enddo
     write(11,*)
  enddo
 
  stop
end program

総称名を用いる場合の手続き(read_filedata())

module read1
  implicit none
  interface read_filedata
     module procedure &
          ! dx1 -> double precision array, x(:)
          ! cy2 -> complex array, y(:,:)
          ! xyy -> coloum of file, | x y y |
          read_dx0_dy1_xy, &
          read_dx0_cy1_xyy, &
          read_dx1_dy1_xy, &
          read_dx1_cy1_xyy, &
          read_dx1_dy2_xy, &
          read_dx1_cy2_xyy, &
          read_dx2_dy2_xy, &
          read_dx2_cy2_xyy
  end interface read_filedata
contains
  subroutine read_dx0_dy1_xy(Ne,y,place,col12)
    integer,intent(in)::Ne,col12
    character(*),intent(in)::place
    double precision,intent(out)::y(1:Ne)

    integer::i,ier
    double precision::a,b
    character(len_trim(place))::fn

    y=0d0
    fn=trim(place)

    open(28,file=fn,status='old',iostat=ier,err=977)
    if(col12.eq.1)then
       do i=1,Ne
          read(28,*)a,b
          y(i)=a
       enddo
    elseif(col12.eq.2)then
       do i=1,Ne
          read(28,*)a,b
          y(i)=b
       enddo
    else
       go to 977
    endif
             
    close(28)
    return

977 write(6,'(3A,i5)')"!!!!==error when open file",trim(fn),"info==>",ier
    write(6,*)"======program stop at read_dx0_dy1_11"
    stop
  end subroutine read_dx0_dy1_xy

  subroutine read_dx0_cy1_xyy(Ne,y,place)
    integer,intent(in)::Ne
    character(*),intent(in)::place
    complex(kind(0d0)),intent(out)::y(1:Ne)

    integer::i,ier
    double precision::a,b,c
    character(len_trim(place))::fn

    y=0d0
    fn=trim(place)

    open(28,file=fn,status='old',iostat=ier,err=977)
    do i=1,Ne
       read(28,*)a,b,c
       y(i)=dcmplx(b,c)
    enddo
    close(28)
    return

977 write(6,'(3A,i5)')"!!!!==error when open file",trim(fn),"info==>",ier
    write(6,*)"======program stop at read_dx0_cy1_xyy"
    stop
  end subroutine read_dx0_cy1_xyy

  subroutine read_dx1_dy1_xy(Ne,x,y,place)
    integer,intent(in)::Ne
    character(*),intent(in)::place
    double precision,intent(out)::x(1:Ne),y(1:Ne)

    character(len_trim(place))::fn
    double precision::a,b
    integer::i,ier

    x=0d0; y=0d0
   
    fn=trim(place)
    open(28,file=trim(fn),status='old',iostat=ier,err=977)
    do i=1,Ne
       read(28,*)a,b
       x(i)=a
       y(i)=b
    enddo
    close(28)
    return

977 write(6,'(3A,i5)')"!!!!==error when open file",trim(fn),"info==>",ier
    write(6,*)"======program stop at read_dx1_dy1_xy"
    stop
  end subroutine read_dx1_dy1_xy
 
  subroutine read_dx1_cy1_xyy(Ne,x,y,place)
    integer,intent(in)::Ne
    character(*),intent(in)::place
    double precision,intent(out)::x(1:Ne)
    complex(kind(0d0)),intent(out)::y(1:Ne)

    integer::i,ier
    double precision::a,b,c
    character(len_trim(place))::fn

    y=0d0
    fn=trim(place)

    open(28,file=fn,status='old',iostat=ier,err=977)
    do i=1,Ne
       read(28,*)a,b,c
       x(i)=a
       y(i)=dcmplx(b,c)
    enddo
    close(28)
    return

977 write(6,'(3A,i5)')"!!!!==error when open file",trim(fn),"info==>",ier
    write(6,*)"======program stop at read_dx1_cy1_xyy"
    stop
  end subroutine read_dx1_cy1_xyy

  subroutine read_dx1_dy2_xy(Nb,Ne,x,y,place)
    integer,intent(in)::Nb,Ne
    double precision,intent(out)::x(1:Ne),y(1:Nb,1:Ne)
    character(*),intent(in)::place
   
    integer::i,j,ier
    character(len_trim(place))::fn
    double precision::a,b
   
    x=0d0; y=0d0
   
    fn=trim(place)
    open(28,file=trim(fn),status='old',iostat=ier,err=977)
    do i=1,Nb
       do j=1,Ne
          read(28,*)a,b
          if(i.eq.1)x(j)=a
          y(i,j)=b
       enddo
    enddo
    close(28)
    return

977 write(6,'(3A,i5)')"!!!!==error when open file",trim(fn),"info==>",ier
    write(6,*)"======program stop at read_dx1_dy2_xy"
    stop
  end subroutine read_dx1_dy2_xy
 
  subroutine read_dx1_cy2_xyy(Nb,Ne,x,y,place)
    integer,intent(in)::Nb,Ne
    double precision,intent(out)::x(1:Ne)
    complex(kind(0d0)),intent(out)::y(1:Nb,1:Ne)
    character(*),intent(in)::place

    integer::i,j,ier
    character(len_trim(place))::fn
    double precision::a,b,c

    x=0d0; y=0d0

    fn=trim(place)
    open(28,file=trim(fn),status='old',iostat=ier,err=977)
    do i=1,Nb
       do j=1,Ne
          read(28,*)a,b,c
          if(i.eq.1)x(j)=a
          y(i,j)=dcmplx(b,c)
       enddo
    enddo
    close(28)
    return

977 write(6,'(3A,i5)')"!!!!==error when open file",trim(fn),"info==>",ier
    write(6,*)"======program stop at read_dx1_cy2_12"
    stop
  end subroutine read_dx1_cy2_xyy

  subroutine read_dx2_dy2_xy(Nb,Ne,x,y,place)
    integer,intent(in)::Nb,Ne
    double precision,intent(out)::x(1:Nb,1:Ne),y(1:Nb,1:Ne)
    character(*),intent(in)::place
   
    integer::i,j,ier
    character(len_trim(place))::fn
    double precision::a,b
   
    x=0d0; y=0d0
   
    fn=trim(place)
    open(28,file=trim(fn),status='old',iostat=ier,err=977)
    do i=1,Nb
       do j=1,Ne
          read(28,*)a,b
          x(i,j)=a
          y(i,j)=b
       enddo
    enddo
    close(28)
    return

977 write(6,'(3A,i5)')"!!!!==error when open file",trim(fn),"info==>",ier
    write(6,*)"======program stop at read_dx2_dy2_11"
    stop
  end subroutine read_dx2_dy2_xy

  subroutine read_dx2_cy2_xyy(Nb,Ne,x,y,place)
    integer,intent(in)::Nb,Ne
    double precision,intent(out)::x(1:Nb,1:Ne)
    complex(kind(0d0))::y(1:Nb,1:Ne)
    character(*),intent(in)::place
   
    integer::i,j,ier
    character(len_trim(place))::fn
    double precision::a,b,c
   
    x=0d0; y=dcmplx(0d0,0d0)
   
    fn=trim(place)
    open(28,file=trim(fn),status='old',iostat=ier,err=977)
    do i=1,Nb
       do j=1,Ne
          read(28,*)a,b,c
          x(i,j)=a
          y(i,j)=dcmplx(b,c)
       enddo
    enddo
    close(28)
    return

977 write(6,'(3A,i5)')"!!!!==error when open file",trim(fn),"info==>",ier
    write(6,*)"======program stop at read_dx2_cy2_12"
    stop
  end subroutine read_dx2_cy2_xyy
end module read1
プログラミングと数値計算

fortranコンパイラのコマンド

コメントする

ここではgfortranとifortでのコンパイル時に指定できるコマンド(最適化とデバッグ、mklの使用とopenmp)について記述します。

ここで対象とするファイル名をmain.f90にします。

最適化

  • gfortran
    gfortran -O3 -fopenmp main.f90
  • ifort
    ifort -mkl -openmp -xHOST -ipo -O3 main.f90

としてコンパイルするといいかと。並列計算がいらない場合は-fopenmp, -openmpのオプションを消してください。

デバッグ

デバッグとはバグを取り除いて望み通りのものにする作業です。

  • gfortran
    gfortran -Wall -fbounds-check -O -Wuninitialized -ffpe-trap=invalid,zero,overflow -fbacktrace -g main.f90
  • ifort
    ifort -mkl -check all -warn all -gen_interfaces -fpe0 -ftrapuv -traceback -g main.f90

とするのがいいでしょう。

※gfortranでは、オプション -ffpe-trap=invalid,zero,overflowが悪さすることがあるようです。fortranでのエラーメモにも書いたように、エラーを起こすはずがないコードでもなぜかエラーと言われることがあります。その時は

gfortran -Wall -fbounds-check -O -Wuninitialized -fbacktrace -g main.f90

でデバッグをやると良いでしょう。

[adsense1]

gfortranでMKLを使う

例えばgfortranコンパイラでIntel® Math Kernel Library (Intel® MKL) を並列計算も含めて使いたいとします。

今、こういう状況にあるとします。

・MKLがあり、/opt/intel/の中にあるディレクトリ”mkl”
が存在する状況下で
・gfortranコンパイラ
がある。

とします。MKLのバージョンによってコンパイルオプションが異なるので注意してください。
Intel® Math Kernel Library Link Line Advisorを見ればどういうコマンドを打てばいいのかを教えてくれます。

ここではMKLのバージョンが
Ver. 11.1
と
Ver. 10.2.5.035
の場合の具体例を載せます。

MKLver11.1の場合
gfortran -fdefault-integer-8 -fopenmp -m64 -I/opt/intel/mkl/include -Wl,--start-group /opt/intel/mkl/lib/intel64/libmkl_gf_ilp64.a /opt/intel/mkl/lib/intel64/libmkl_core.a /opt/intel/mkl/lib/intel64/libmkl_gnu_thread.a -Wl,--end-group -ldl -lpthread -lm main.f90
MKLver10.2.5.035の場合
gfortran -fopenmp -m64 -I/opt/intel/mkl/include /opt/intel/mkl/lib/em64t/libmkl_solver_lp64.a -Wl,--start-group /opt/intel/mkl/lib/em64t/libmkl_gf_lp64.a /opt/intel/mkl/lib/em64t/libmkl_core.a /opt/intel/mkl/lib/em64t/libmkl_gnu_thread.a -Wl,--end-group -lpthread -lm main.f90

※もしかしたら/opt/intel/mkl/ではなくて、/opt/intel/mkl/10.2.5.035/じゃないとだめかもしれません。

バージョンが違ったらどうするかはIntel® Math Kernel Library Link Line Advisorを参考にしてください。上のものを得るには下のようにすれば大丈夫です。

 Select Intel® product: (MKLのバージョンは？)  Intel(R) MKL10.2
 Select OS: (OSは何？)  Linux*
 Select usage model of Intel® Xeon Phi™ Coprocessor: 
(走らせようとしているコンピュータのcpuの種類がXeon Phi っていうcpuが数百個あるような特別なもの？)		  –
 Select compiler: (何のコンパイラでMKLを使いたいの？, gfortranだったらGNU fortran.)  GNU Fortran
 Select architecture:(良く分かりません！Help me →I don’t use any flag と進んで、”getconf LONG_BIT” と端末で打って何て表示されたかで入力しました。多分cpuが何ビットか、かな？)  Intel(R) 64
 Select dynamic or static linking:

(動的リンク(dynamic)がいい？静的リンク(static)がいい？
　　　　動的→汎用性あり、容量は軽い、動作は遅い
　　　　静的→汎用性なし、容量は重い、動作は早い )		  Static
 Select interface layer:(プログラムの中に整数型で2^31越えてるものはある？ない？越えてなかったら32bit integer, 越えてたら64bit integer)  LP64 (32-bit integer)
 Select sequential or multi-threaded layer:
(mklのサブルーチンを1つのCPUだけ(sequential)で走らせたい？それとも並列(multi-threaded)にする？)		  multi-threaded
 Select OpenMP library:(openMPの種類は何使う？intel か GNUか→僕は両方試したら、intelでは動かず、GNUでは動きました。)  GNU (libgomp)
 Select cluster library:　ここから下は何を言っているのか分かりませんでした。空白でもいけたから多分直接は関係ないオプションなのだろうと。  

次に下の方で作られたオプションをgfortranコンパイル時にくっつける。
$(MKLROOT)/libem64t/libmkl_solver_lp64.a -Wl,–start-group $(MKLROOT)/libem64t/libmkl_gf_lp64.a $(MKLROOT)/libem64t/libmkl_core.a $(MKLROOT)/libem64t/libmkl_gnu_thread.a -Wl,–end-group -lpthread -lm
と
-fopenmp -m64 -I$(MKLROOT)/include
が作られるかと思います。この２つをくっつければokです。
$(MKLROOT)はディレクトリ”mkl”へのパスです。シェルスクリプトの場合、()は要りません。

ここで、この通りコンパイルしようとすると/libem64t/というディレクトリはないよ、っていわれて出来ませんでした。
実際にたどってみると/lib/em64t/ならあって、これに変えたらコンパイルが成功しました。ケースバイケースかな。

参考

Fortranデバッグ用オプション

[adsense2]

プログラミングと数値計算, 物理学

二重振り子

コメントする

座標の取り方は下図のように取ります。棒の伸び縮みは無いものとします。
二重振り子座標
どういう解き方でもいいですが、ここでは

  1. デカルト座標\(L(x,y)\)でラグランジアンを記述
  2. デカルト座標から座標変換し、\((r,\theta)\)でラグランジアンを記述
  3. 新たな座標系で運動方程式を導く

という順で解いていきます。

[adsense1]

1, デカルト座標でのラグランジアンLは(運動エネルギーK)-(位置エネルギーU)と書けるため、
\(
L(x_1,\dot{x}_1,y_1,\dot{y}_1,x_2,\dot{x}_2,y_2,\dot{y}_2)\\
\displaystyle =\frac{1}{2}m_1(\dot{x}_1^2+\dot{y}_1^2)+\frac{1}{2}m_1(\dot{x}_1^2+\dot{y}_2^2)-(-m_1gy_1-m_2gy_2)
\)
と書けます。

2, デカルト座標から座標変換
式を簡単にするために座標変換を行います。新しい座標\((r_1,\theta_1,r_2,\theta_2)\)とデカルト座標\((x_1,y_1,x_2,y_2)\)の関係式は
\(
\begin{align}
x_1&=r_1\sin{\theta_1}\\
y_1&=-r_1\cos{\theta_1}\\
x_2&=r_1\sin{\theta_1}+r_2\sin{\theta_2}\\
y_2&=-r_1\cos{\theta_1}-r_2\cos{\theta_2}
\end{align}
\)
という関係があります。各々を時間で微分すれば、
\(
\begin{align}
\dot{x}_1&=\dot{r}_1\sin{\theta_1}+r_1\dot{\theta}_1\cos{\theta_1}\\
\dot{y}_1&=-\dot{r}_1\cos{\theta_1}+r_1\dot{\theta}_1\sin{\theta_1}\\
\dot{x}_2&=\dot{r}_1\sin{\theta_1}+r_1\dot{\theta}_1\cos{\theta_1}+\dot{r}_2\sin{\theta_2}+r_2\dot{\theta}_2\cos{\theta_2}\\
\dot{y}_2&=-\dot{r}_1\cos{\theta_1}+r_1\dot{\theta}_1\sin{\theta_1}-\dot{r}_2\cos{\theta_2}+r_2\dot{\theta}_2\sin{\theta_2}
\end{align}
\)
これらをラグランジアン\(L(x_1,\dot{x}_1,y_1,\dot{y}_1,x_2,\dot{x}_2,y_2,\dot{y}_2)\)に代入します。すると、新たな座標系でのラグランジアン\(L(r_1,\dot{r}_1,\theta_1,\dot{\theta}_1,r_2,\dot{r}_2,\theta_2,\dot{\theta}_2)\)が得られます。
\(
\begin{align}
L(r_1,\dot{r}_1,\theta_1,\dot{\theta}_1,& r_2,\dot{r}_2,\theta_2,\dot{\theta}_2) \\
=&\frac{1}{2}m_1(\dot{r}_1^2+r_1^2\dot{\theta}_1^2)+\frac{1}{2}m_2\left[\dot{r}_1^2+\dot{r}_2^2+r_1^2\dot{\theta}_1^2+r_2^2\dot{\theta}_2^2 \right. \\
&\left.+2(\dot{r}_1 r_2 \dot{\theta}_2-r_1\dot{r}_2 \dot{\theta}_1)\sin{(\theta_1-\theta_2)}+2(\dot{r}_1 \dot{r}_2 +r_1 r_2 \dot{\theta}_1 \dot{\theta}_2)+\cos{(\theta_1-\theta_2)}\right] \\
&+m_1gr_1\cos{\theta_1}+m_2g(r_1\cos{\theta_1}+r_2\cos{\theta_2})
\end{align}
\)

僕は先ほど式を簡単にするために座標変換をする、といいました。しかし、新しい座標系におけるラグランジアンはどう見ても元のデカルト座標系のラグランジアンに比べて複雑です。この理由は物理的な意味から来ています。
振り子をつないでいる棒が伸び縮みしないとすると系の自由度は角度\(\theta_1,\theta_2\)の２つです。
となると運動方程式は最高で2本の独立した方程式になるはずです。
しかし、デカルト座標の場合うまく自由度を落とすことができず、運動方程式は4つになってしまいます。
そこで棒が伸び縮みを簡単に表すことができる座標系に移ることでうまく方程式の数を減らせます。

新しい座標系でのラグランジアンで棒の伸び縮みがないという条件を表すには
\(
\begin{align}
r_1&=l_1\ (l_1\mbox{は定数}) \\
r_2&=l_2\ (l_2\mbox{は定数})
\end{align}
\)
と書けるわけで、また、
\(
\begin{align}
\dot{r}_1&=0 \\
\dot{r}_2&=0
\end{align}
\)
となるわけです。

\(m_1=m_2=m,\ l_1=l_2=l\)という場合を特に考えると、ラグランジアンは
\(
\displaystyle L(\theta_1,\dot{\theta}_1,\theta_2,\dot{\theta}_2)=ml^2\left[\dot{\theta}_1^2+\frac{1}{2}\dot{\theta}_2^2+\dot{\theta}_1\dot{\theta}_2\cos{(\theta_1-\theta_2)}\right]+mgl(2\cos{\theta_1}+\theta_2)
\)
と書けます。あとはラグランジュの運動方程式を当てはめて計算します。

3, 新たな座標系で運動方程式を導く
保存力場中でのラグランジュの運動方程式は
\(
\begin{eqnarray}
\left\{
\begin{aligned}
\frac{d}{dt}\left(\frac{\partial L}{\partial\dot{\theta}_1}\right)-\frac{\partial L}{\partial\theta_1}&=0 \\
\frac{d}{dt}\left(\frac{\partial L}{\partial\dot{\theta}_2}\right)-\frac{\partial L}{\partial\theta_2}&=0
\end{aligned}
\right.
\end{eqnarray}
\)
なので、代入し、\(\ddot{\theta}_1,\ddot{\theta}_2\)に関する運動方程式にすれば
\(
\begin{eqnarray}
\left\{
\begin{aligned}
\ddot{\theta}_1&=\frac{1}{2}\left\{-\ddot{\theta}_2\cos{(\theta_1-\theta_2)}-\dot{\theta}_2^2\sin{(\theta_1-\theta_2)}-2\frac{g}{l}\sin{\theta_1} \right\} \\
\ddot{\theta}_2&=-\ddot{\theta}_1\cos{(\theta_1-\theta_2)}+\dot{\theta}_1^2\sin{(\theta_1-\theta_2)}-\frac{g}{l}\sin{\theta_2}
\end{aligned}
\right.
\end{eqnarray}
\)
となります。これは非線形の２階連立微分方程式です。カオスです。解けません。
数値計算で解きます。用いるのは4次ルンゲ・クッタ法です。

実際に解いてみると、こうなります。

２重振り子の実験もされています。

[adsense2]

カオスが現れる場合、初期値鋭敏性という性質があります。
初期値がほんのちょっと変わるだけでその後の時間発展の様子が大きく変わる性質です。
ちょっと値を変えると上の動画とはまるで違う動きになります。
この動画は上の動画よりも中心に近いほうの粒子の初期速度が4%違います

機械精度の違いでさえも十分な時間発展の後には大きな違いが現れます。複雑であると同時に面白い現象です。