sikino のすべての投稿

gnuplot

gnuplotで配列を使う

2015年12月27日 sikino コメントする

gnuplotで配列はver5.0の時点では実装されていません。
文字を利用することで疑似的な配列が使えます。

まとめ

配列の名前を ” a ” として考えます。

a = ''
do for [i = 1:5] {
a = a . sprintf('%f ', 0.2*i)
}

plot for [i=1:5] word(a,i)*x

とスクリプトを書き、プロットすると

という画像が得られます。

配列の利用

gnuplotでは配列が実装されていません。
doループがver4.6から実装されたことから、これから実装されるかもしれません。

今、
$
y=0.2x,\;\;0.4x,\;\;0.6x,\;\;0.8x,\;\;x
$
の合計5つの関数をプロットしたいとします。
もしも素直に考えるならば、スクリプトは

plot for [i=1:5] i*0.2*x

です。これくらいは簡単ですが、複雑になった場合はどうでしょう。そのために疑似配列の使い方を知っておきましょう。

gnuplotで配列は実装されていないので、配列に似たもの(疑似配列と呼ぶことにします)を作ります。
普通のプログラミングでの”配列”のイメージは、番号付けした変数をデータの数だけ用意し、そこに代入していくイメージです。すなわち、
$a(1)=0.2,\;\; a(2)=0.4,\;\; a(3)=0.6,\;\; a(4)=0.8,\;\; a(5)=1$
という具合にです。

が、gnuplotでの疑似配列は、空白によって値と値が分けられた1つの文字型変数を用いて作ります。すなわち、
$a= ” 0.2\;\; 0.4\;\; 0.6\;\; 0.8\;\; 1″ $
を用意し、呼び出し時に
word(配列名、番号)
というgnuplotに用意されている関数を用いてそれぞれの要素にアクセスします。

この疑似配列aの作り方はこうです。

a = ''
do for [i = 1:5] {
a = a . sprintf('%f ', 0.2*i)
}

1行目は配列”a”の初期化です。
2~4行目は配列”a”に実際に値を入れていきます。
この形を覚えておけば配列を作ることができます。

配列の呼び出しは、word(,)を用います。1つ目の引数に対象の疑似配列、2つ目の引数に何番目の要素を呼び出すかを指定します。以下のように、

a = ''
do for [i = 1:5] {
a = a . sprintf('%f ', 0.2*i)
}
plot for [i=1:5] word(a,i)*x

とすればグラフとして書き出されます。

[adsense1]

疑似配列の四則演算

疑似配列の要素に数字を加えたいとします。
その場合は

b = ''
do for [i=1:5] {
b = b . sprintf('%f ',word(a,i)+1)
}

とすれば疑似配列bに疑似配列aのそれぞれの要素に1加えた値が格納されます。
※同じ配列aにその結果を格納したい場合は最後にa=bを付け加えてください。

疑似配列同士の和を取りたければ、

c = ''
do for [i=1:5] {
c = c . sprintf('%f ',word(a,i)+word(b,i))
}

とすれば疑似配列$c(i)$に(疑似配列$a(i)$)+(疑似配列$b(i)$)の値が格納されます。
差ならば-, 積ならば*, 商ならば/に変えてください。

実際に、以下のスクリプトを動かしてみます。

a = ''
do for [i = 1:5] {
a = a . sprintf('%f ', 0.2*i)
}
print a

b = ''
do for [i=1:5] {
b = b . sprintf('%f ',word(a,i)+1)
}
print b

c = ''
do for [i=1:5] {
c = c . sprintf('%f ',word(a,i)*word(b,i))
}
print c

とすると実行結果

gnuplot> load "array.plt"
0.200000 0.400000 0.600000 0.800000 1.000000
1.200000 1.400000 1.600000 1.800000 2.000000
0.240000 0.560000 0.960000 1.440000 2.000000
gnuplot>

を得ます。

応用例

半径の違う複数の球

A = ''
do for [i = 1:5] {
A = A . sprintf('%f ', 1e0/i)
}

set parametric
set ur[0:pi]
set vr[0:2*pi]
set key below
set ticslevel 0

set view equal xyz

splot for [i=5:1:-1] word(A,i)*sin(u)*cos(v),word(A,i)*sin(u)*sin(v),word(A,i)*cos(u)

というスクリプトを用いますと、半径の違う球が掛けます。

アニメーション

num=5
A = ''
do for [i = 1:num] {
A = A . sprintf('%f ', 1e0/i)
}

do for [j=1:119] {
t=0.05*j
set view 68,j*1.5,1,1
A=''
do for [i = 1:num] {
A = A . sprintf('%f ', 1e0/sqrt(1+10*sin(t+i*(2e0*pi)/num)**2))
}
splot for [i= num:1:-1 ] word(A,i)*sin(u)*cos(v),word(A,i)*sin(u)*sin(v),word(A,i)*cos(u)
}

を表示に関する設定をして実行すると下のアニメーションが得られます。

参考文献

Operate on 1D arrays in gnuplot

[adsense2][adsense3]

gnuplot

角度に依存して色を付ける(gnuplot)

2015年12月13日 sikino コメントする

gnuplotです。
gnuplotで、関数、もしくはデータでsplotすると3列目のデータに依存してカラーリングがされます。新たに4列目にデータを加え、splotをする時に

splot "test.d" u 1:2:3:4 w pm3d

とすると、4列目に従ったカラーリングがされます。

データの場合

例えば、以下のデータ(球面調和関数$l=1,m=0$)を持っていたとします。
(ここにテキストデータ”y10.d”(500kB)を置きました。)
1列目:$\sin{\theta}\cos{\phi}\cdot Y_{l,m}$, 2列目:$\sin{\theta}\sin{\phi}\cdot Y_{l,m}$, 3列目:$\cos{\theta}\cdot Y_{l,m}$, 4列目：$\theta$, 5列目：$\phi$
そのまま出力する場合、以下のようなスクリプトで、以下の画像が取得できるはずです。

set pm3d
set pm3d depthorder
unset sur
set view equal xyz
set ticslevel 0
splot "y10.d" u 1:2:3

これを角度$\theta$に依存させて、形は変えずに色を付けましょう。
splotするときに、カラーを4列目に沿って付けるとすれば、

splot "y10.d" u 1:2:3:4 w pm3d

であり、もしも角度$\phi$に依存させて描くのであれば

splot "y10.d" u 1:2:3:5 w pm3d

でokです。実際に書いて、出力するとそれぞれ以下のような画像が出力されます。

※左から3列目($Y_{l,m}$の値)に従って付けたカラー、4列目$\theta$に従って付けたカラー、5列目$\phi$に従って付けたカラー
となります。

※球面調和関数の出力に関して
球面調和関数は球面座標で$\theta,\phi$の関数$Y_{l,m}(\theta,\phi)$です。
ただし、今出力したい座標はデカルト座標($x,y,z$)です。デカルト座標で球面座標で記述される画像を得たい場合は単位ベクトル$\; \vec{r}/|\vec{r}|\;$に球面座標の関数を掛ければよいのです。
よって
$
\displaystyle Y_{l,m}(\theta,\phi)\cdot \frac{\vec{r}}{|\vec{r}|}
$
を表示すればよい、すなわちgnuplot上では($r=1$の場合)、
1列目に $x\;\;\to \;\; \sin{\theta}\cos{\phi}\cdot Y_{l,m}(\theta,\phi)$
2列目に $y\;\;\to \;\; \sin{\theta}\sin{\phi}\cdot Y_{l,m}(\theta,\phi)$
3列目に $z\;\;\to \;\; \cos{\theta}\cdot Y_{l,m}(\theta,\phi)$
を表示させればokです。

多価関数の場合

続いて、多価関数として有名なリーマン面を出力します。
通常のリーマン面(3列目の高さがカラーの基準)を表示させるには

set pm3d depthorder
set parametric
set ur[0.01:5]
set vr[0:4*pi]
i={0,1}
splot u*cos(v),u*sin(v),real(sqrt(u)*exp(i*0.5*v))

とすれば下にある、左側の画像が得られます。角度に依存したカラーリングをしたい場合、以下のスクリプトでokです。

set isosamples 150
set pm3d depthorder

set table "param1.d"
set parametric
set ur[0.01:5]
set vr[0:4*pi]
i={0,1}
splot u,v,real(sqrt(u)*exp(i*0.5*v))
unset table

set table "param2.d"
plot "param1.d" u ($1*cos($2)):($1*sin($2)):3:($2) w vector
unset table

splot "param2.d" u 1:2:3:4 w pm3d

データでは無く、解析的な関数を直接、別の列を基準にして色付け、出力させることはできないようです。なので、関数を数値的なファイルに出力させるコマンド

set table "data.d"
...(ここに数値的に出力させたいグラフをプロット)
unset table

を書いて出力を行わせています。

また、角度を出力させるので周期的なカラーリングが良いと思います。色相環を元に下のカラーマップを使うといいでしょう。

set palette defined(0"#ff0000",1"#ff8000",2"#ffff00",3"#80ff00",4"#00ff00",5"#00ff80",6"#00ffff",7"#0080ff",8"#0000ff",9"#8000ff",10"#ff00ff",11"#ff0080",12"#ff0000")

上のスクリプトを実行すると右図のように、角度に依存したカラーリングがされます。

※Riemann面について
Riemann面は
$
f(z)=Re(\sqrt{z})
$
で表される多価関数です。
gnuplotで単純に考えて上のグラフをプロットしようとすると、

i={0,1}
splot real(sqrt(x+i*y))

で良いはずですが、これを実際にプロットするとこうなります。

本来は多価関数なのですが、それゆえに上記方法ではgnuplotでは記述されません。
極座標を用いて以下のように考えます。
$
\begin{align}
f(z)&=\sqrt{z} \\
&=\sqrt{r}e^{i\frac{1}{2}\theta} (0\le \theta\lt 4\pi)\\
\end{align}
$
であるため、gnuplotでは媒介変数表示u,vを用いて
$(x,y,z)=(u\cos{v},u\sin{v},\sqrt{u}e^{i\frac{1}{2}v})$
$(u:0\sim \infty, v:0\sim 4\pi)$
とすればいいのです。

参考文献

gnuplot：値に応じて線に色付け -TN’s blog
gnuplotで数値をファイルに書き出す（table）
How to output smooth cspline curve as a data file -Stack Overflow

gnuplot

gnuplotで極座標のグラフを書く

2015年12月6日 sikino コメントする

gnuplotで極座標形式$(r,\theta)$で、こんな感じの2次元のグラフ

 Spacetime is warped and twisted by the
mass and spin of the earth. -GRAVITY PROBE B Testing Einstein’s Universeより

の、緑色のようなグラフを書きたいと思ったのが始まりです。
グラフの式を以下のように適当に書きます。
$
\displaystyle \displaystyle f(r)=-\frac{1}{r^2+1}
$
とモデル化しましょう。これを2つの方法でプロットします。

[adsense1]

gnuplot上で

splot -1/(x*x+y*y+1)

とすると、

というように出力されます。

上の場合、これはグリッドがx,y方向に切られています。動径方向と角度方向に切りたい場合はどうすればいいでしょうか。

gnuplotでこれを実現する方法は少なくとも2つあります。

データとしてファイルを作った後、splot
媒介変数表示を使う

です。今回は2番目の媒介変数表示を使って表示させます。

とりあえず、これで目的のグラフを得ることができます。

set parametric
set ur[0:5]
set vr[0:2*pi]
splot u*sin(v),u*cos(v),-1/(u**2+1)

でこのようになります。

仕組みは極座標とは
$
\begin{align}
x&=r\cos\theta \\
y&=r\sin\theta
\end{align}
$
なので、媒介変数uをrに、vをθをとみなせば書くことができます。

[adsense2]

3dプロットでの破綻

表面をカラー表示にしたいと思った時に、カラーに破綻が起こりました。

set parametric
splot u*sin(v),u*cos(v),-1/(u**2+1)
set ur[0:5]
set vr[0:2*pi]
set pm3d at s
unset sur
rep

とするとこうなりました。

これを解決するには

set pm3d depthorder

を追加してreplotすれば破綻しなくなります。set hidden3d のカラー版ですね。

参考
gnuplot demo script: surface2.dem

サバイバルゲーム, 弾道計算, 物理学

バレル内部でのBB弾の方程式

2015年12月3日 sikino コメントする

バレル内部でのBB弾の運動方程式です。

目的は、

BB弾は何秒間バレル内部に存在しているのか？
バレルが長いと減速になりうるのか？

を知ることです。詳細は弾道計算本をご覧ください。

弾道計算(BB弾)の理論
 弾道計算(BB弾)の結果
 弾道計算の結果2, 比較と詳細データ
 弾道計算(BB弾)のコード(fortran90)
弾道計算のコード(Excel)
バレル内部でのBB弾の方程式←今ここ

電動ガンの場合です。

空気抵抗は初速0[m/s]から90[m/s]前後にまで加速されるわけですから、空気抵抗の粘性抵抗、慣性抵抗どちらの項も無視することはできないでしょう。

空気の漏れも考えます。

この条件下では、バレル内部の運動方程式は以下のように導くことができます。

内部圧力変化と力

ピストン-BB弾間の圧力変化による力$F_V(t)$は、
$
F_V(t)=S\cdot P(t)
$
と書けます。ここで$P(t)$はピストン-BB弾間の内部圧力です。
$P(t)$は、断熱変化を仮定すると、ピストン-BB弾間体積$V(t)$を用いて
$
P(t)V^{7/5}(t)=\text{const}
$
の関係があります。体積$V(t)$は
$
V(t)=x(t)S_b-x_0(t)S_s
$
と書けます。ここで、ピストンの位置による時間変化を$x_0(t)$、BB弾の位置を$x(t)$、バレルの断面積を$S_b$、シリンダーの断面積を$S_s$としました。
この式は時刻に依存せずに決まるので、時刻$t=0$の時に大気圧$P_0$で、体積$V_0$であれば、任意の時刻での圧力$P(t)$は
$
\displaystyle P(t)=P_0 \left(\frac{V_0}{V(t)}\right)^{7/5}
$
と書けます。

空気の漏れについて

十分短い時間の間、時刻$t$において内部圧力$P(t)$, 空気の密度$\rho$とすると、漏れ出る空気の速度$v(t)$はベルヌーイの定理から
$
\begin{align}
P(t)&=\frac{1}{2}\rho v(t)^2+P_0 \nonumber \\
&\to v(t)=\pm \sqrt{\frac{2}{\rho}(P(t)-P_0)}
\end{align}
$
が成立するとします。厳密には、ピストンの速度は漏れ出る空気の速度に対して無視できるほど小さい、という仮定の下で成立します。漏れ出る流量$Q$は、バレル-BB弾の隙間の断面積$A$、実験とのズレを調節する無次元の流量係数$c’$を用いて$Q=c’ A v(t)$と書けます。
また、断熱過程を圧力-体積の関係を用いたいので、漏れ出る空気はBB弾-ピストン間の体積の増加として扱います。

フルシリンダーの場合の運動方程式

その他、空気抵抗による力を入れると、BB弾の位置$x(t)$, ピストンの位置$x_0(t)$, 体積$V(t)$の運動方程式は
$
\begin{align}
m\frac{d^2}{dt^2}x(t)&=\left[P(t)-P_0\right]S_b-\frac{1}{2}C_d \rho \pi R^2 |v(t)|^2\cdot\frac{v(t)}{|v(t)|} \label{bbin1}\\
m_s\frac{d^2}{dt^2}x_0(t)&=-k\left[x_0(t)-x_B-l\right]-[P(t)-P_0]S_s -F_f\label{bbin2} \\
\frac{d}{dt}V(t)&=v(t)S_b-v_0(t)S_s+c'(S_b-\pi R^2)\mathrm{sgn}(P(t)-P_0)\sqrt{\frac{2}{\rho}|P(t)-P_0|} \label{bbin3}
\end{align}
$
と導くことが出来ます。

ここで、$P_0$は大気圧、$\eta$は粘性率、$R$はBB弾の半径、$C_d$は抗力係数、$\rho$は空気の密度, $S_b$はバレルの断面積、$S_b$はシリンダーの断面積、$k$はピストンのばね定数、$l$はばねの自然長、$F_f$はピストン-シリンダーの摩擦、$v(t),v_0(t)$はそれぞれBB弾の速度、ピストンの速度を意味します。

空気抵抗に関する詳細は球体の空気抵抗と係数をご覧ください。

ここで１つ言えることは、外部と内部の圧力が一定になる最適なバレル長というものが存在する、ということです。

2019/01/20 追記)
本の頒布日より1年以上たちました。
計算結果は載せないつもりでしたが、載せたくなりましたので少しだけ載せておきます。
まず、実測データと上で立てたモデルの運動方程式の比較をします。
比較対象はインターネット上で公開されていた二つのデータです。

実測データとの比較

1つ目(図の←)はさばなび様で公開されていた記事
https://www.saba-navi.com/2015/10/29/laboratory_work_barrel_cut/で実測されていたデータで、バレル長の長さと初速の関係です。ただし、既にリンク切れのようです。
ピストンの重さの記述が無かったので、典型的とされている重さ25gを仮定して計算しています。

2つ目(図の→)は石岡様のホームページ
○電動ガンバレル、シリンダの組み合わせによる初速実験
の結果との比較です。

どちらの計算結果もまぁまぁ合っていることから、私が仮定したモデルは見当はずれなものではないということが分かるかと思います。

バレル内部の様子

さて、上のモデルが正しいとした時、運動方程式を解いてバレルの動きや内圧を考えてみましょう。
解いてみますと、こういった図が得られます。

図は、東京マルイM16を想定した時のBB弾の位置、速度、ピストンの位置、内圧の時間変化です。
計算のパラメータは、
BB弾の重さ0.20g
ピストンの重さ24g
ばね定数431N/m
ばねの自然長150mm
押し切られた時のばねの長さ120mm
BB弾の半径5.95mm
バレルの直径6.05mm
とした時の結果です。
この場合、生み出せる最高初速は96m/s（バレル長63cmの時）。
例えばバレル長が50cmの時、ピストンが動いてからBB弾が射出されるまでにかかる時間は0.013秒ということが分かります。
すなわち、上のピストンの重さやBB弾の重さ、ばねの強さ、シリンダー容量の場合は最適なバレル長は63cmであるということです。
また、内圧の上昇によるピストンのブレーキなども計算できていることが分かるかと思います。

また、漏れ出る空気の量は$3000[\mathrm{mm^3}]$だということが見積もられました。この量はM16の場合は本来のシリンダー容量の15%程度です。

典型的なバレル長である50cmの場合、BB弾は射出されるまでに約0.013秒かかります。内圧の上昇によるブレーキがかからず、同じ過程でばねが引かれると仮定すると、ピストンは約0.010秒で戻ると考えられます。
理論上の効率的な射出サイクルの限界は

（1発）÷（0.013秒＋0.010秒）≒ 43発/秒

という事です。現在販売されている東京マルイのハイサイクルは25発/秒なので、理論上あと1.7倍早くすることが可能です。

BB弾の重さを0.25gにした結果や、ばねの強さを変えた時の結果は
弾道計算本として計算結果をまとめた本に載せているので、そちらをご参照ください。本の詳細については弾道計算本の自家通販をご覧ください。
また時間が経ったらここに載せる…かもしれません。

補遺

補遺1
理想気体として扱える条件は、
低い気圧(分子の数が少なく、衝突等が無視できる)かつ高い温度(分子間力が分子の運動エネルギーに比べて無視できる)
です。
どうやら実在気体では10気圧以下、室温以上でこの条件は良く満たされ、理想気体とのずれは1%以内のようです[1]。
[1]を引用すると、

一般に，沸点の低い酸素･窒素･水素･ヘリウム等は，室温またはそれ以上の温度で10atm以下の圧力の場合，理想気体の値の1％以内で理想気体に近い性質を示す。

とあります。ピストンで空気が圧縮されたとき、BB弾とピストン間の圧力が10気圧以上にならなければ良い近似だと言えるでしょう。

参考文献

[1] 実在気体 -第2節　気体の状態方程式
[2]PERFECT HIT -TOKYO MARUI

『集弾性アップへの道』　BB弾とバレル内部に隙間があることが写真で確認できます。

未分類

著作物とその使用について

2015年11月23日 sikino コメントする

著作物についてです。
私はこういう事に関する専門家ではありません。インターネットで調べた結果なので、過信しないでください。

著作物が自由に使える場合

著作物が自由に使える条件とは何でしょうか。それは、
著作物が自由に使える場合は？ -公益社団法人著作権情報センター
にまとめられています。
要件を満たす場合、著作権にかかわらず自由に使うことができます。

例(音源について)

例えば学園祭で音楽にのせてダンスを披露したい。この場合の音源に関する著作権はどうかというと、

非営利目的の演奏など（著作権法第38条）
営利を目的とせず、観客から料金をとらない場合は、著作物の上演・演奏・上映・口述（朗読）などができる。ただし、出演者などは無報酬である必要がある。

とあります。
つまり、無料かつ無報酬、非営利目的であれば著作物であっても使用可能なのです。

では学園祭でコンテストが行われ、優勝、準優勝したら賞金や商品などがある場合はどうなるのでしょう。
この場合は営利目的かそうでないか、また、それは利益なのかが問題となるようです。

参加団体が継続的に活動を行い、収入の一部になり、その団体メンバーに得た利益を分配していれば利益目的です。
しかし、「営利」、「非営利」を調べてみると、

「営利」とは、構成員（株主など）の経済的利益を追求し、団体の利益を構成員が分配することを意味します。
　営利組織である会社は、株主が出資して会社を運転し、あがった利益を株主に配当するしくみになっています。
　それに対して、「非営利」とは、団体が利益を上げてもその利益を構成員（会員など）に分配しないという「非分配」を意味します。つまり、「非営利」とは、利益を上げてはいけないという意味ではなく、「利益があがっても構成員（社員など）に分配しないで、団体の活動目的を達成するための費用に充てること」と説明することができます。-NPOWEBより引用

にあるように、団体の活動目的のために使用する状態であれば非営利活動であるとみなせるようです。
ただし、最終的な判断は裁判所にゆだねられています。「営利目的」と「非営利目的」の間に明確な定義が無いのです。

すなわち、参加するために掛かった活動費用等と照らし合わせ、受け取った額(物)が利益に当たらないと判断できる、もしくは活動維持のために必要であるくらいならば、コンテストで賞金受け取ったとしても非営利目的なのです(補足有り)。

画像の使用

Webで画像を使用するときの利用はどうなのでしょうか。

表記がいらなく、自由に使える場合

著作者の表記等がいらなく、営利目的であっても、許可を得ずに複製・改変・翻案・配布・上演・演奏することが可能であるものというのは、パブリックドメイン状態にあるもの、もしくはクリエイティブコモンズライセンスがCC0にあるものです。

パブリックドメインとは、保護期間が終了したり、権利が放棄されている状態のことを言います。1．著作権全般　-日本著作権教育研究会によると、

がそれに当たります。

クリエイティブコモンズライセンスCC0とは、著作権者が利益を放棄してパブリックドメインに置くことに同意している画像のことです。ただし、いくら著作権者が権利を放棄したくて完全にパブリックドメインに置きたいと考えても、著作権を完全に放棄することは上記条件が満たされない限り法的には困難であるため、パブリックドメインという言葉ではなく、そう呼んでいるのです。

パブリックドメインとクリエイティブコモンズライセンスCC0は、ほぼ等しいものなのです。
パブリックドメインは一般的な言葉、クリエイティブコモンズライセンスCC0は民間で作られた言葉で、パブリックドメインのことを指す。という考えで良いでしょう。

例えば、この記事↓
無料で商用利用できる写真を13のサイトから一気に串刺し検索してゲット可能な「LibreStock」-GIGAZINEより
で紹介されているサイト「LibreStock」の画像は全てクリエイティブコモンズライセンスCC0に属しています。
なので、勝手にとってきて保存して、営利目的であろうが改変しようが再配布しようが構わないのです。ただし、パブリックドメインもそうですが、著作権は消えても著作者人格権という権利は残されるため、自分が作ったものだと言い張って損害を与えたり、作者をけなすような悪質な改変行為等は許されません。

公表された著作物の利用の場合

パブリックドメインまたはクリエイティブコモンズライセンスCC0以外の画像は引用という形で利用することができます。

公表された著作物に対しては、基本的には改変等をせず、引用にあたるのであれば利用可能です。
引用の詳細は、

引用と言えるためには、[1]引用する資料等は既に公表されているものであること、[2]「公正な慣行」に合致すること、[3]報道、批評、研究などのための「正当な範囲内」であること、[4]引用部分とそれ以外の部分の「主従関係」が明確であること、[5]カギ括弧などにより「引用部分」が明確になっていること、[6]引用を行う必然性があること、[7]出所の明示が必要なこと(複製以外はその慣行があるとき）（第４８条)の要件を満たすことが必要です(第３２条第１項)。-著作権なるほど質問箱より

であり、上の条件[1]~[7]が満たされていれば公表されている著作物も引用として利用可能なのです。

しかし、権利者の中には
『自分が元ネタであることは明記してほしいけど、改変とかは自由にしてもらいたい。』
『営利目的でなければ、この作品を、改変や再配布をしてもok。』
と考える場合もあるでしょう。
この細かい権利、利用条件を分かりやすく表示させようとしたシステムがクリエイティブ・コモンズ・ライセンスなのです。

クリエイティブ・コモンズ・ライセンスに関して

クリエイティブ・コモンズ・ライセンスに関する定義は、クリエイティブ・コモンズ・ライセンスとはに詳しく書かれています。
引用すると、

CCライセンスとはインターネット時代のための新しい著作権ルールで、作品を公開する作者が「この条件を守れば私の作品を自由に使って構いません。」という意思表示をするためのツールです。-creative commons JAPANより引用

というものです。
このクリエイティブコモンズライセンスに法的な根拠はありません。民間で作られた表記方法が世間一般に広く浸透したものです。現在では様々な機関、論文でも使用されることがあり、一般的な扱いを受けています。
wikipediaの画像にはクリエイティブコモンズライセンスが採用されています。

さて、クリエイティブ・コモンズ・ライセンスの表記はどのようにされるのでしょうか。
クリエイティブ・コモンズ・ライセンスとはによると、以下のような表示のいずれかがなされます。使用条件の緩い順に並べるとこうです。

表示 CC BY 4.0

(表示クリエイティブ・コモンズ・ライセンスとは -creative commons JAPAN より)
「表示」のみの場合、

◆元の作品の「Ⓒ 著作権者の名前公表年」の３点セット（これを「著作権表示」又は「クレジット」と呼ぶことがあります）
→元の作品に記載されている場合には、必ず記載してください。

◆元の作品の作者名、スポンサー、タイトル
→元の作品に表示があれば記載してください。

◆元の作品の著作権表示かライセンス情報に関するページへの指定されたＵＲＬ
→元の作品に表示があれば記載してください。
CCライセンスの作品を利用する際、何を記載すればよいのでしょうか。また、クリエイティブ・コモンズ・ライセンスのURLやアイコン、バナー等も表示しなければならないのでしょうか。FAQ　よくある質問と回答

を記載すれば使うことができます。これを守っていれば、営利目的、改変、再配布したりすることができます。

表示-非営利 CC BY 2.1

(表示—非営利クリエイティブ・コモンズ・ライセンスとは -creative commons JAPAN より)
「表示」+「非営利」であれば使用可能であることを示しています。
「表示」は上に同じ。
「非営利」であることは結局は裁判所の判断となります(「NC（非営利）」アイコンのついている作品を使用しても良いですか？-CCJ Q&A)。もしかしたら、些細な広告収入があるだけでも営利目的であると判断されるかもしれません。迷う場合は、非営利がついたものは使わないほうがいいでしょう。

表示-継承 CC BY-SA 3.0

(表示—継承クリエイティブ・コモンズ・ライセンスとは -creative commons JAPAN より)
「表示」+「継承」するのであれば使用可能であることを示しています。
「表示」は上に同じ。
「継承」は、改変は許されるが、もし改変した場合、元の作品と同じCCライセンスで公開しなければなりません。ただし、その条件下で有れば営利目的での2次利用は許可されます。

表示-非営利-継承 CC BY-NC-SA 2.1

(表示—非営利—継承クリエイティブ・コモンズ・ライセンスとは -creative commons JAPAN より)
「表示」、「非営利」、「継承」は上に出てきているのでここでは割愛します。

表示-改変禁止CC BY-ND 2.1

(表示—改変禁止クリエイティブ・コモンズ・ライセンスとは -creative commons JAPAN より)
「表示」は上に同じ。
「改変禁止」は文字通り改変禁止です。トリミングも、拡張子の変換も、サイズの変更も許されません。ただし、それを守れば営利目的に利用、再配布もできるのです。

表示-非営利-改変禁止CC BY-NC-ND 2.1

(表示—非営利—改変禁止クリエイティブ・コモンズ・ライセンスとは -creative commons JAPAN より)
「表示」、「非営利」、「改変禁止」は上に出てきているのでここでは割愛します。

自分の作品にクリエイティブ・コモンズ・ライセンスを付ける

自分の作品にクリエイティブ・コモンズ・ライセンスを付けたい場合は、どこかに申請する必要など無く、自分で言い張ればいいんです。なぜならばあなたが作品を作り出した時に著作権は発生しているのですから。その作品にCCライセンスを付けるか否かは貴方の自由です。
もしもCCライセンスを付けたい場合はこちらのページに行けばどのように付ければいいかわかります。
、
※後からクリエイティブ・コモンズ・ライセンスの変更は行ってはなりません。
なぜならば、例えばはじめCC0ライセンスで公開し、皆が使い始めたところライセンスを変えたりしたらまるで使い物にならなくなってしまうからです。
なので、重要な事はその権利者、例えばページを運営している人が信用、信頼できるかどうかを見極めることが重要です。

公式のQ&Aにもありました。

・ライセンスを変更したい場合はどうしたらよいでしょうか。

クリエイティブ・コモンズ・ライセンスは取消ができません。つまり、クリエイティブ・コモンズ・ライセンスのついた作品を入手した誰かに対して、そのライセンスで認められている利用を止めるということはできません。
–ライセンスを変更したい場合はどうしたらよいでしょうか。日本法準拠版クリエイティブ・コモンズ・ライセンスのFAQより引用

実例

では実際にwikipediaから画像を取ってきて、表記をしてみましょう。

表示のみの場合CC BY 2.0

りんご -wikipedia です。

©Abhijit Tembhekar(2009)/Adapted/CC BY 4.0

必ず表記しなければならないのは、元の作品の「Ⓒ 著作権者の名前公表年」の３点セットです。上画像の詳細ページに進むと、Date,Authorの欄を見つけることができますので、クレジットの表記は「©Abhijit Tembhekar(2009)」でokです。また、必須ではありませんが、念のためにその画像に対するリンクを張っています。またその他タイトル等はページには見つけることができませんでした。なのでこれでCCライセンスの「表示」はokです。

表示・継承の場合CC BY-SA 3.0

球面調和関数 -wikipedia にある画像を利用します。

©Sarxos(2007)/Adapted/CC BY-SA 3.0

必ず表記しなければならないのは、元の作品の「Ⓒ 著作権者の名前公表年」の３点セットです。上画像の詳細ページに進むと、Date,Authorの欄を見つけることができますので、クレジットの表記は「©Sarxos(2007)」でokです。また、必須ではありませんが、念のためにそのオリジナル画像に対するリンクを張っています。またその他タイトル等はページには見つけることができませんでした。なのでこれでCCライセンスの「表示」はokです。

また、サイズの変更を行ったのでこれは改変に当たります。なので、それを示すAdaptedを記述しています。改変を行った時点で2次創作物であり、その時点からCCライセンスは僕が決める事ができますが(※)、今回は継承の表記があります。なのでCCライセンスは元作品と同じライセンスを提示しなければなりません。CCライセンスの提示は上のように何のライセンスであるかということと、リンク先を張ればokです。

※⒝ 元の作品に改変を加えて二次的著作物を創作した場合によると、2次創作物のCCライセンスは、元の作品よりも制限を課す場合でしか認められていません。今回は継承なので気にする余地はありませんが、勝手に改変したものを「表示」だけにして公表してはなりません。

パブリックドメイン、もしくはCC0の場合

自然 wikipediaにあるパブリックドメインの画像です。

クレジットの表記はまるで必要ありません。営利目的であっても、許可を得ずに複製・改変・翻案・配布することが可能です。

これはCC0の例です。
https://commons.wikimedia.org/wiki/File:%22A_mixed_flock_of_ducks_and_geese_fly_from_a_wetland_area%22.jpg

CC0もクレジットの表記はまるで必要ありません。営利目的であっても、許可を得ずに複製・改変・翻案・配布することが可能です。

wikipediaにはCC0の画像一覧が存在します。

その他CCライセンスの場合

非営利の場合はここでは載せません。というのも、広告収入がどう判断されるかわからないため、実例を挙げるわけにはいかないからです。また、非営利は難しい言葉ではないため、載せる必要がないとも考えています。
また、改変禁止も載せません。改変せずに載せればいいだけなので…。

補足

分かりやすく解説された著作権に関する事柄は著作者にはどんな権利がある？ -公益社団法人著作権情報センターをご覧ください。

そのクレジット表記で大丈夫？クリエイティブ・コモンズ（CC）画像の使い方。

gnuplot

gnuplotとアニメーション

2015年11月7日 sikino 8件のコメント

gnuplot上の様々なアニメーションを作ります。
dofor構文を使うのでgnuplot ver4.6以降じゃないとこの方法は出来ません。

解析解の表示
データの表示
応用
動画の作成について
1. gifアニメの作成方法
2. ffmpegを用いてavi動画を作る

[adsense1]

1次元の場合の具体例

試しに量子力学で現れる、ガウシアン波束の時間発展を考えましょう。
こんな形の方程式を考えます。
$
\displaystyle y(x,t)=\left(\frac{1}{1+t^2}\right)^{1/2}\exp\left[-\frac{(x-t)^2}{1+t^2}\right]
$

これをgnuplot上で時間発展させるスクリプトはこちら。

do for [i = 0:400 ] {
t=i*0.02
plot sqrt(1/(1+t*t))*exp(-(x-t)**2/(1+t*t))
}

適宜xrangeとかyrangeとか変更してください。実行するとこんな感じ。

上のgifアニメを得るには、

set term gif animate optimize delay 2 size 480,360
set output 'movie.gif'

do for [i = 0:400 ] {
t=i*0.02
plot sqrt(1/(1+t*t))*exp(-(x-t)**2/(1+t*t)) lw 2
}

set out
set terminal wxt enhanced

というスクリプトを使ってください。
同ディレクトリにmovie.gifが作成されます。
ターミナルとかは各々の環境に合わせてください。

2次元の場合の具体例

2次元だろうとこのように簡単にできます。

上の動画のスクリプトはこれです。gifアニメが欲しい場合はコメントアウトを外してください。

#set term gif animate optimize delay 10 size 480,480
#set output 'movie.gif'

set pm3d at b
set xr[-5:10]
set yr[-5:10]
set zr[0:1]
set cbr[0:1]
set isosamples 50

do for [i = 0:50 ] {
t=i*0.05
splot sqrt(1/(1+t*t))*exp(-(x-t)**2/(1+t*t))*sqrt(1/(1+t*t))*exp(-(y-2*t)**2/(1+t*t))
}

#set out
#set terminal wxt enhanced

点をグラフ上に沿って動かす

のような動画を作りたいとします。この時は、媒介変数を用いて以下のスクリプトで再現されます。

#set term gif animate optimize delay 5 size 960,720
#set output 'movie.gif'

set param
set size ratio -1
set samples 10000

e = 1
omega=0.1

set tr[1:600]
do for [i = 1:200 ] {
plot e*cos(omega*t)/sqrt(t), sin(omega*t)/sqrt(t)
set label 1 point pt 7 ps 3 at e*cos(omega*i*3)/sqrt(i*3),sin(omega*i*3)/sqrt(i*3)
}

#set out
#set terminal wxt enhanced

点の表示は、

set label 1 point pt 7 ps 3 at e*cos(omega*i*3)/sqrt(i*3),sin(omega*i*3)/sqrt(i*3)

の文である1点だけを表示することができます。

グラフをだんだんと書いていく

グラフをちょっとずつ書いていきます。
解析的な解の場合

set param
set samples 10000
set tr[0.01:1]
imax=100
tmax=20e0*pi
ht=tmax/real(imax)

#set term gif animate optimize delay 6 size 600,600
#set output 'orbit.gif'

do for [i=1:imax] {
th(t,i)=t*real(i)*ht
plot 10e0*sin(th(t,i))/th(t,i),10e0*cos(th(t,i))/th(t,i) , \
10e0*sin(th(t,i)-pi)/th(t,i),10e0*cos(th(t,i)-pi)/th(t,i) lt 1 lc 2
}

#set out
#set terminal wxt enhanced

データの場合（200行のデータ”outgraph.d”を一行目から順番に出力したい場合）

do for [i=1:200]{
plot "outgraph.d" every ::0:0:i:0 using 1:2 w l
}

グラフを回転させる

ある3Dプロットがあり、それを回転させてみるような動画がほしいとします。例えばこんな感じのアニメーションが作りたいとします。

これはリーマン面と呼ばれる関数です。(カラーバーの消去は ‘unset colorbox’ とすれば消せます。)
詳しい記述は角度に依存して色を付ける(gnuplot)に書いてあるので気になる人は参照してください。

set pm3d depthorder
set parametric
set ur[0.01:5]
set vr[0:4*pi]
i={0,1}
splot u*cos(v),u*sin(v),real(sqrt(u)*exp(i*0.5*v))

#set term gif animate optimize size 480,360
#set output 'movie.gif'

do for [j = 0:90 ] {
set view 60,4*j,1,1
replot
}

#set out
#set terminal wxt enhanced

参考先はamesyabodyの記述したモンテカルロ法のページからです。

複数のデータを次々とプロット

複数のデータを時間経過を表示させる場合はデータfort.100, fort.101, fort.102,…を用意し、dofor構文の中を変えて、

#set term gif animate optimize size 480,360
#set output 'movie.gif'
do for [j = 100:200 ] {
splot sprintf("./fort.%d", j) u 1:2:3 w l t sprintf("./fort.%d", j)
}

#set out
#set terminal wxt enhanced

してloadすればいいです。

複数グラフの同時プロット

data1.d, data2.d, data3.d,data4.d,data5.d
を同時に出力したい時は、

plot for [i=1:5] sprintf("data%d.d",i) w l ti sprintf("data%d.d",i)

とすればいいです。

魔女の宅急便

魔女の数式はこちらです。
witch like curve -wolfram alpha

▼gnuplotコード(この場に展開)

第５使徒ラミエル1

▼gnuplotコード(この場に展開)

第５使徒ラミエル2

▼gnuplotコード(この場に展開)

魔貫光殺砲

▼ここクリックでの場に展開

マリオ

マリオのグラフは↓から。
Mario like curve -wolfram alpha

▼ここクリックでの場に展開

[adsense2]

SCP-1968っぽい数式

SCP-1968 世界を包む逆因果の円環
に現れている数式です。
数式はTorus helix radius change equationより。

▼ここクリックでこの場に展開

FGO召喚

特定の面と3D視点で見た場合

このリンク先(https://slpr.sakura.ne.jp/qp/supplement_data/FGO_gnuplot.zip)のデータも必要です。
（なんか良く分かりませんがメジェド様がバグるんですよね…）
▼ここクリックでこの場に展開

パンジャンドラム

▼ここクリックでこの場に展開

gifアニメについて

gifアニメの描画時間を早くしたり遅くしたりするには、

set term gif animate delay 10 optimize size 480,360

などと「delay 10」という文を付け加えましょう。delayは遅延、後ろの数字は遅延時間で、単位は0.01sになっています。
例として、
delay 5の場合は1秒間に20枚、
delay 10の場合は1秒間に10枚、
delay 50の場合は1秒間に2枚
表示させるスピードだということです。
delayはどんなに早くしたくても4程度までがいいです。delayを1とかにするとそんなに早いgifアニメの描写が出来ず、ブラウザ上でうまく表示できず、結果として遅く表示されます。

Gif -gnuplot 4.2

ffmpegを用いてavi動画を作る

ffmpegのインストール方法はCompile FFmpeg on Ubuntu, Debian, or Mintを参考にしてください。
このページにもインストール方法を書きましたが、更新が早いので公式のページを見るほうが確実です。
ここではgnuplotの視点を変更した画像をたくさん用意して動画を作成する、という方法をとります。
ディレクトリ”dat/”と視点変更と連続した画像を出力するスクリプトは”image.plt”というファイルを用意して

do for [j = 0:180 ] {
set view 60,2*j,0.5,3
replot

call "out.plt" "temp"

if(0<=j && j<10){
cv=sprintf("mv temp.jpg ./dat/image0000%d.jpg",j)
}
if(10<=j && j<100){
cv=sprintf("mv temp.jpg ./dat/image000%d.jpg",j)
}
if(100<=j && j<1000){
cv=sprintf("mv temp.jpg ./dat/image00%d.jpg",j)
}
if(1000<=j && j<10000){
cv=sprintf("mv temp.jpg ./dat/image0%d.jpg",j)
}
if(10000<=j && j<100000){
cv=sprintf("mv temp.jpg ./dat/image%d.jpg",j)
}
system cv
}

#Using Imagemagick, convert
#cv=sprintf("convert -delay 20 -resize 100% -loop 0 ./dat/image*.jpg ./anime.gif")

#Using ffmpeg
cv=sprintf("ffmpeg -r 12 -i './dat/image%05d.jpg' -vcodec mjpeg -qscale 0 -s 1300x888 movieout.avi")
system cv

とすればディレクトリ./dat/の中に180枚の.jpg画像が作られ、その場所にファイル名”movieout.avi”という動画が作られます。

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

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

2015年11月3日 sikino コメントする

非線形シュレディンガー方程式
$
\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番目の、粒子のような性質を持つことから、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$です。

[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です！

海外

モスクワでSIMカードの購入

2015年10月30日 sikino コメントする

まとめ

2015年10月、モスクワにあるМегаФонでプリペイドsimカードを買った時のお話
店員にまずは英語で”I want a prepaid sim card.”と話す。
ダメなら”Предоплаченные Cим-карты, Пожалуйста.“と言う。
パスポートを渡し、サインして、お金払ってお仕舞い。
※標準, micro, nanoの違いはありません。simを切り離すことでどれでも対応できます。

モスクワ市内であればそこそこ英語は通じます。
可能な限り空港で買うのが良いでしょう。
シェレメーチエヴォ国際空港(SVO)の場合、一度空港の出発ロビーへ行くのが良いでしょう。確か到着ロビーは1階で、出発ロビーは3階か4階です。そこまで行けば歩いていれば確実に見つかります。

ロシアのモスクワでsimカードを購入してインターネットに繋げるまでの体験です。
この情報は2015年10月のお話で、あくまで僕が体験したことを思い出しながら書いています。
当たり前ですがsimフリーじゃなきゃできません。

英語は申し訳程度の保険的な扱いです。
モスクワでは英語は通じません。英語の案内版すらありません。

僕は
МегаФон
というお店で買いました(英語名Megafon,メガフォン)。

まずお店に入ります。
まずは店員にダメ元で英語で話しかけるのがいいでしょう。
通じれば何でもいいですが、例えば
“I want a prepaid sim card.”
とかでいいでしょう。
これでダメでしたらロシア語で、
プリペイドsimカードを下さい。

“Предоплаченные Cим-карты, Пожалуйста.”

と言うか、この文を見せればいいでしょう。
※スマートフォン等にgoogle翻訳のオフライン版を事前に入れておくことをお勧めします。
僕の場合はたまたま英語が通じる店員が店内にいたので困らずに変えました。

これで通じたら、どのプランがいいか、を聞かれると同時にパスポートの提示が求められます。
店員にパスポートを渡します。
その後、simカードが入ったカードの束が渡され、
「好きな電話番号を選んでくれ」
と言われます。
このプリペイドsimに付属する電話番号は(多分)ロシア国内でМегаФонの電話へ掛ける場合に使えるものです。
まぁ、掛けるつもりが無いのならば何でもいいでしょう。

2015年10月時点ではМегаФонのプリペイドsimのデータプランは5つ存在し、
199pub, 0.5GBまで, ロシアの番号でМегаФонの番号なら通じる, 通話可能時間300分
390pub, 3GBまで, ロシアの番号でМегаФонの番号と家への電話なら通じる, 通話可能時間400分
590pub, 4GBまで, ロシアの番号でМегаФонの番号と家への電話なら通じる, 通話可能時間600分
1290pub, 8GBまで, ロシアの番号全てに通じる, 通話可能時間1500分
2700pub, 10GBまで, ロシアの番号全てに通じる, 通話可能時間5000分
です。

で、サインを求められるので、サインしてお金を払って終了です。
このサインは漢字でも平仮名でもアルファベットでも何でも構いません。

最後にこんなsimカードがもらえます。

時間としては10~15分程度だったでしょうか。早かった印象があります。

数学

– 1 = 1 ?

2015年10月3日 sikino コメントする

-1=i×i =√(-1)×√(-1) =√((-1)×(-1)) =√1 =1 つまり、-1=1かにゃっ？

— 物理学科に入学した凛ちゃんbot (@RinPhysBot) 2015, 9月 19

が発端です。間違いであることは直観で分かります。
問題はどこに間違いがあるか？です。

改めて問題を書けば、
$
\begin{align}
-1&=i\times i \\
&=\sqrt{(-1)}\times\sqrt{(-1)} \\
&=\sqrt{(-1)\times(-1)} \\
&=\sqrt{1} \\
&=1
\end{align}
$
です。

結論を先に言えば、2行目から3行目にかけての変形がダメです。
定義域の考慮が抜けています。
なぜダメなんでしょうか。詳しく見ていきましょう。

間違いが生じるのはなぜ?

$i$を極座標形式で考えます。
すると、
$
\displaystyle i=0+i=e^{i\frac{\pi}{2}}
$
ですよね。これは正しいです。
そして次に2行目から3行目を極形式でゆっくりと書いていけば、
$
\begin{align}
-1&=i\times i \\
&=e^{i\frac{\pi}{2}}\times e^{i\frac{\pi}{2}} \\
&=\left(e^{i\pi}\right)^{1/2} \times\left(e^{i\pi}\right)^{1/2} \cdots \mbox{2行目に対応} \\
&=\left(e^{i2\pi}\right)^{1/2}\cdots \mbox{3行目に対応}
\end{align}
$
です。

もしも、$e^{i\theta}$の$\theta$の範囲を勝手に$0 \le \theta \lt 2\pi$にしていたら。
この場合、角度$2\pi$はこの定義域内に含まれていません。
$e^{i\theta}$は$2\pi$の周期性を持つはずだ、だから
$
e^{i2\pi}=e^{i0}=1
$
であるはずだ、と思ってしまうんです。
これを行ってしまうと、最後の式は、
$
\begin{align}
\left(e^{i2\pi}\right)^{1/2}&=\left(e^{i0}\right)^{1/2} \\
&=1^{1/2} \\
&=1
\end{align}
$
となります。

よって間違った答え$-1=1$が生まれます。

間違いを起こさないためには

この問題を解決するためには定義域をきちっと示してあげればいいのです。
別の複素数を掛ける場合、通常、$\theta$の定義域も考慮しなければなりません。
故に、厳密に書けば、
$
\displaystyle i=e^{i\frac{\pi}{2}}\ \ \ (0\le\theta \lt 2\pi)
$
であり、
$
\begin{align}
-1&=i\times i \\
&=e^{i\frac{\pi}{2}}\times e^{i\frac{\pi}{2}}\ \ \ \ (0\le\theta \lt 4\pi)\\
&=\left(e^{i2\pi}\right)^{1/2}\ \ \ (0\le\theta \lt 4\pi)
\end{align}
$
なのです。$e^{i\theta}$は本来, 周期関数であって多価関数であることを忘れてはいけません。
故に、定義域内に収まっていないからと言って$2\pi$を加える必要はないのです。
だから、形はそのままであり、計算は
$
\begin{align}
&\left(e^{i2\pi}\right)^{1/2}\ \ \ (0\le \theta\lt 4\pi)\\
&=e^{i\pi}\ \ \ (0\le \theta\lt 2\pi)\\
&=-1
\end{align}
$
だから、やはり$-1=-1$なのです。

※上の式$-1=i\times i\cdots$は長いので、
$
\sqrt{(-1)^2}
$
の値は何ですか？という問いかけでもいいと思います。
この場合、この答えは$\pm 1$です。どちらでも正解です。
なぜならば、$X=\sqrt{(-1)^2}$とおいて、
$
\begin{align}
X&=\sqrt{(-1)^2} \\
X^2&=1 \\
&(X+1)(X-1)=0 \\
\end{align}
$
だからです。

この問題が顕著に表れるのは数値計算です。コンピュータは定義域の考慮はできません。人為的に入れるしか手段は無く、必ず$2\pi$の範囲に収まっていなければならないのです。なので、数値計算で複素数を扱う場合、この問題が発生していないか確かめましょう。符号が知らない間に変わっている、軽視いているととんでもない問題になります。

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

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

2015年9月20日 sikino 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$が掛けられたものが良く使用されます。

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

上図のように位置空間上で区間$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倍ほど速くなります。

▼ここクリックでこの場に展開

include "/opt/intel/mkl/include/mkl_dfti.f90"
program main
implicit none
integer::i,n
double precision,allocatable::x(:),k(:),mk(:)
complex(kind(0d0)),allocatable::z(:),mz(:)
double precision::hx,xmin,xmax

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

n=130
allocate(x(0:n-1),k(0:n-1),z(0:n-1))
x=0d0; k=0d0; z=dcmplx(0d0,0d0)

allocate(mk(0:n-1),mz(0:n-1))
mk=0d0; mz=dcmplx(0d0,0d0)

xmin=-3d0*pi
xmax=3d0*pi
hx=(xmax-xmin)/dble(n)

do i=0,n-1
x(i)=xmin+dble(i)*hx
z(i)=func(x(i))
enddo

open(21,file="before.d")
do i=0,n-1
write(21,'(3e20.10e3)')x(i),dble(z(i)),dimag(z(i))
enddo
close(21)

call dftf(n,k,hx) ! Get frequency k

call dft(n,z,"forward") ! forward dft

call dfts(n,k,z,mk,mz) !sort frequency
open(21,file="k.d")
do i=0,n-1
write(21,'(3e20.10e3)')mk(i),dble(mz(i)),dimag(mz(i))
enddo
close(21)

call dft(n,z,"backward") ! backward dft
open(21,file="after.d")
do i=0,n-1
write(21,'(3e20.10e3)')x(i),dble(z(i)),dimag(z(i))
enddo
close(21)

stop
end program main

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

func=dcmplx(dsin(4d0*x),0d0)

return
end function func

!-------------------------------

subroutine dft(n,z,FB)
!sikinote
! developer : sikino
! date : 2015/09/20
use MKL_DFTI
implicit none
integer,intent(in)::n
complex(kind(0d0)),intent(inout)::z(0:n-1)
character(*)::FB

integer::Status
TYPE(DFTI_DESCRIPTOR),POINTER::hand

!DFT : Discrete Fourier Transform
!
!n --> number of data.
!z(i) --> value of data at x(i)
!FB --> "forward" : Forward DFT
! --> "backward" : Backward DFT

Status = DftiCreateDescriptor(hand,DFTI_DOUBLE,DFTI_COMPLEX,1,n)
Status = DftiCommitDescriptor(hand)
if(trim(FB).eq."forward")then
Status = DftiComputeForward(hand,z)
elseif(trim(FB).eq."backward")then
Status = DftiComputeBackward(hand,z)
else
write(6,*)"DFT string different"
stop
endif
Status = DftiFreeDescriptor(hand)

if(trim(FB).eq."backward")then
z=z/dble(n)
end if

return
end subroutine dft

subroutine dfts(N,f,z,mf,mz)
!sikinote
! developer : sikino
! date : 2015/09/20
implicit none
integer,intent(in)::N
double precision,intent(in)::f(0:N-1)
double precision,intent(out)::mf(0:N-1)
complex(kind(0d0)),intent(in)::z(0:N-1)
complex(kind(0d0)),intent(out)::mz(0:N-1)

integer::i,j

if(mod(n,2).eq.0)then
do i=0,N-1
if(i.le.N/2)then
j=i+N/2-1
mf(j)=f(i)
mz(j)=z(i)
else
j=i-N/2-1
mf(j)=f(i)
mz(j)=z(i)
endif
enddo
else
do i=0,N-1
if(i.le.(N-1)/2)then
j=i+(N-1)/2
mf(j)=f(i)
mz(j)=z(i)
else
j=i-(N-1)/2-1
mf(j)=f(i)
mz(j)=z(i)
endif
enddo
endif

return
end subroutine dfts

subroutine dftf(N,f,h)
!sikinote
! developer : sikino
! date : 2015/09/20
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)*0.5d0)/(dble(N)*h)
enddo

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

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

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

dft(N,z,FB)
dftf(N,f,h)
dfts(N,f,z,mf,mz)

を説明します。

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$を掛けています。
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を使用してください。
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\}
$
が得られるはずで、実際にプロットして確かめてみると

となります。
※Nは2以上のどんな整数でもokです。偶数、奇数、素数であっても。

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

ベンチマーク

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

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

▼クリックでこの場に展開

program main
implicit none
integer::i,n
double precision,allocatable::x(:)
complex(kind(0d0)),allocatable::z(:)
double precision::hx,xmin,xmax

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

call mkl_set_num_threads(1)

do n=2,1100
if(mod(n,50).eq.0)write(6,*)”n –> “,n
allocate(x(0:n-1),z(0:n-1))
z=dcmplx(0d0,0d0)

xmin=-3d0*pi
xmax=3d0*pi
hx=(xmax-xmin)/dble(n)

do i=0,n-1
x(i)=xmin+dble(i)*hx
z(i)=func(x(i))
enddo

call benchm(n,z)

deallocate(x,z)
enddo

stop
end program main

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

func=dcmplx(dsin(4d0*x),0d0)

return
end function func

subroutine benchm(n,z)
implicit none
integer::n,k,kmax
double precision::t0,t1
complex(kind(0d0))::z(0:n-1)

kmax=1250
call cpu_time(t0)

do k=1,kmax
call dft(n,z,”forward”) ! forward dft
call dft(n,z,”backward”) ! backward dft
enddo

call cpu_time(t1)

write(10,'(i5,e20.10e4,A)’)n,(t1-t0)/dble(kmax)

return
end subroutine benchm

subroutine dft(n,z,FB)
!sikinote
! developer : sikino
! date : 2015/09/20
use MKL_DFTI
implicit none
integer,intent(in)::n
complex(kind(0d0)),intent(inout)::z(0:n-1)
character(*)::FB

integer::Status
TYPE(DFTI_DESCRIPTOR),POINTER::hand

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

Status = DftiCreateDescriptor(hand,DFTI_DOUBLE,DFTI_COMPLEX,1,n)
Status = DftiCommitDescriptor(hand)
if(trim(FB).eq.”forward”)then
Status = DftiComputeForward(hand,z)
elseif(trim(FB).eq.”backward”)then
Status = DftiComputeBackward(hand,z)
else
write(6,*)”DFT string different”
stop
endif
Status = DftiFreeDescriptor(hand)

if(trim(FB).eq.”backward”)then
z=z/dble(n)
end if

return
end subroutine dft

subroutine dfts(N,f,z,mf,mz)
!sikinote
! developer : sikino
! date : 2015/09/20
implicit none
integer,intent(in)::N
double precision,intent(in)::f(0:N-1)
double precision,intent(out)::mf(0:N-1)
complex(kind(0d0)),intent(in)::z(0:N-1)
complex(kind(0d0)),intent(out)::mz(0:N-1)

integer::i,j

if(mod(n,2).eq.0)then
do i=0,N-1
if(i.le.N/2)then
j=i+N/2-1
mf(j)=f(i)
mz(j)=z(i)
else
j=i-N/2-1
mf(j)=f(i)
mz(j)=z(i)
endif
enddo
else
do i=0,N-1
if(i.le.(N-1)/2)then
j=i+(N-1)/2
mf(j)=f(i)
mz(j)=z(i)
else
j=i-(N-1)/2-1
mf(j)=f(i)
mz(j)=z(i)
endif
enddo
endif

return
end subroutine dfts

subroutine dftf(N,f,h)
!sikinote
! developer : sikino
! date : 2015/09/20
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)*0.5d0)/(dble(N)*h)
enddo

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

並列計算時の計算速度

call mkl_set_num_threads(Ncpu)

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

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

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

▼ここクリックでこの場に展開

program main
implicit none
integer::i,n,j,Nr,l,Ntry,jx
double precision,allocatable::x(:),k(:),mk(:)
complex(kind(0d0)),allocatable::z(:),mz(:)
double precision::hx,xmin,xmax
character(48)::fname

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

xmin=-3d0*pi
xmax=3d0*pi

Ntry=10
do l=1,Ntry
write(fname,'(A,i0,A)')"ffttime_",l,".d"
open(11,file=trim(fname))
do jx=5,25
n=2**jx
allocate(x(0:n-1),k(0:n-1),z(0:n-1))
x=0d0; k=0d0; z=dcmplx(0d0,0d0)

hx=(xmax-xmin)/dble(n)

do i=0,n-1
x(i)=xmin+dble(i)*hx
z(i)=func(x(i))
enddo

Nr=1000
do j=1,16
call mkl_set_num_threads(j)
call system_clock(ti)
call fft(n,z,0,Nr)
call system_clock(tf,tr)
write(11,'(i6,i8,f10.3,2i6)')j,jx,(tf-ti)/dble(tr),Nr
write(6,'(i6,i8,f10.3,2i6)')j,jx,(tf-ti)/dble(tr),Nr
enddo
write(11,*)
deallocate(x,k,z)
enddo
close(11)
enddo

stop
end program main

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

func=exp(-0.5d0*x*x)

return
end function func

include "/opt/intel/mkl/include/mkl_dfti.f90"
subroutine fft(N,z,ifb,Ncycle)
use MKL_DFTI
implicit none
integer,intent(in)::N,ifb,Ncycle
complex(kind(0d0)),intent(inout)::z(0:N-1)

integer::Status,i
TYPE(DFTI_DESCRIPTOR),POINTER::hand

Status = DftiCreateDescriptor(hand,DFTI_DOUBLE,DFTI_COMPLEX,1,N)
Status = DftiCommitDescriptor(hand)

!if(ifb.eq.1)Status = DftiComputeForward(hand,z)
!if(ifb.eq.-1)Status = DftiComputeBackward(hand,z)
do i=1,Ncycle
Status = DftiComputeForward(hand,z)
Status = DftiComputeBackward(hand,z)
z=z/n
enddo

Status = DftiFreeDescriptor(hand)

return
end subroutine fft

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)$)は

となります。
波数空間で$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=-3d0pi
xmax=3d0pi
hx=(xmax-xmin)/dble(Nx)
ymin=-3d0pi
ymax=3d0pi
hy=(ymax-ymin)/dble(Ny)

do i=0,Nx-1
x(i)=xmin+dble(i)hx
enddo
do j=0,Ny-1
y(j)=ymin+dble(j)hy
enddo
do i=0,Nx-1
do j=0,Ny-1
z(i,j)=func(x(i),y(j)) 
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(4d0y)+dsin(3d0x),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)0.5d0)/(dble(N)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