sikino のすべての投稿

vmware playerの容量を減らす

VMware playerの容量が大きくなってきた時に大きく効果のあったことです。

ここでは、
ホストOS:windows8.1
上に
vmware playerにインストールされてるゲストOS:Linux Mmint 17.1 rebecca
がある、という状況で行った手順をまとめます。

ゲストOS(VMwareを立ち上げたlinux内)のコマンドラインで以下のことを行ってください。VMware toolsを使います。インストールされてない方はLinuxMint導入時にやることに書いたので、VMware toolsのインストールを行ってください。

sudo vmware-toolbox-cmd disk list
sudo vmware-toolbox-cmd disk shrink /

で後は1~2時間待てばいいです。
僕はこれによって、windows上のフォルダ”Virtual Machines”の占める容量を43GB→22GBにまで減らせました。


  1. 仮想マシンのシャットダウンと起動(サスペンドではなく、シャットダウンです。)
  2. ファイル整理
    (特にwindowsから仮想OSへドラッグ&ドロップする際に生成される

    ~/.cache/vmware/drag_and_drop/

    の削除)

  3.     sudo vmware-toolbox-cmd disk list
        sudo vmware-toolbox-cmd disk shrink /

    の実

参考文献


VMWareのゲストイメージ圧縮方法

三角波、のこぎり波、矩形波、その他の数式

正弦波ではない周期関数、三角波、矩形波等々の数式による表現です。
if文は使っていません。
床関数\(floor(x)\)を主に用いています。

これらが唯一の作り方ではなく、最善の作り方でもありません。

非連続になる関数の場合の特定な関数の場合、中間の値を取るようにしています(フーリエ級数の収束より)。


階段関数

——-
f1
z_f1_x_c
——-
f2
z_f2_x_c
——-
fl
z_fl_x_c


のこぎり波

s1
z_s1_x_c
s2
z_s2_x_c
saw
z_saw_x_c


矩形波

k1,k2,k,ka,kb,kk
z_k1_abx_c

z_k2_abx_c

z_k_abx_c

z_ka_abx_c

z_kb_abx_c

z_kk_abx_c


その他

sl,s,spk
tr,se,de,mo
well1,well2

sl
z_sl_abx_c

s
z_s_ax_c
三角波は、\({\rm acos}(\cos(x))\)とも書けます。

spk
z_spk_abx_c

tr
z_tr_abx_c

se
z_se_abx_c

de
z_de_abx_c

mo
z_mo_abx_c

well1
zwell1_abx_c

well2
z_well2_abx_c

well

追記
\(\displaystyle
f(a,b,x) = sgn(abs(2*(x-b)/a)-1e0)
\)
でokです。aは井戸の幅、bは井戸の中心を示します。


まとめ
periodic_functions_c

全てをまとめたgnuplotのコードは以下のものです。

set size ratio -1
set yr[-2:2]
set xr[-5:5]
set samples 101
set grid
set xtics 1
set ytics 1
set mxtics 2
set mytics 2

f1(x)=floor(x)
f2(x)=-floor(-x)-1e0
fl(x)=0.5e0*(floor(x)-floor(-x)-1e0)

s1(x)=x-floor(x)
s2(x)=x+floor(-x)+1e0
saw(x)=x-0.5e0*(floor(x)-floor(-x)-1e0)

k1(a,b,x)=f1((x+a)/(a+b))-f1(x/(a+b))
k2(a,b,x)=f2((x+a)/(a+b))-f2(x/(a+b))
k(a,b,x)=fl((x+a)/(a+b))-fl(x/(a+b))
kk(a,b,x)=2e0*(k(a,b,x)-0.5e0)

ka(a,b,x)=0.5e0*(sgn(k(a,b,x)-0.25e0))+0.5e0
kb(a,b,x)=0.5e0*(sgn(k(a,b,x)-0.75e0))+0.5e0

sl(a,b,x)=kk(a+b,a+b,x+0.5e0*b+a)*k(a,b,x)

s(a,x)=2e0*abs(s2(x/(a*2))-0.5e0)*kk(a*2e0,a*2e0,x+a)
tr(a,b,x)=k(a,b,x)*s(0.5e0*(a+b),x-0.5e0*b)*(1e0+b/a)+sl(b,a,x-b)

spk(a,b,x)=-abs(kk(a,b,x))+1e0

se(a,b,x)=(x-b-(a+b)*fl(x/(a+b)))/a*k(a,b,x)+0.5e0*(k(a,b,x)*(1e0-b/a)+1e0)*spk(a+b,a+b,x)

de(a,b,x)=-se(a,b,x)*2e0*(k2(a+b,a+b,x-0.5e0*b)-0.5e0)
mo(a,b,x)=de(a*0.5e0,b+a*0.5e0,x+a*0.5e0)+de(a*0.5e0,b+a*0.5e0,-(x+a*0.5e0))

well1(a,b,x)=0.5e0*(sgn((abs(2e0*(x-b)/a)-1e0)+ka(1e0,1e0,abs(2e0*(x-b)/a))-0.5e0))+0.5e0
well2(a,b,x)=0.5e0*(sgn((abs(2e0*(x-b)/a)-1e0)+kb(1e0,1e0,abs(2e0*(x-b)/a))-0.5e0))+0.5e0

exp(ikx)が進行波であることの証明

\(e^{i{\bf k\cdot r}}\)が\({\bf k}\)の方向に進行すること、\(e^{-i{\bf k\cdot r}}\)が\(-{\bf k}\)の方向に進行することの所以を証明します。

まとめ

波は、実部と虚部を並べた時、虚部がある方向に進行していきます。

証明


一次元自由粒子の時間発展を考えます。
自由粒子の時間依存シュレーディンガー方程式は
\(
\displaystyle i\hbar\frac{\partial}{\partial t}\Psi(x,t)=\left(-\frac{\hbar^2}{2m}\frac{\partial^2}{\partial x^2}\right)\Psi(x,t)
\)

と書けます。時間依存しないシュレーディンガー方程式を求めるために、
\(\Psi(x,t)=\psi(x)e^{-i\frac{E}{\hbar}t}\)
と置いて代入すると
\(
\displaystyle -\frac{\hbar^2}{2m}\frac{d^2}{dx^2}\psi(x)=E\psi(x)
\)

となります。
上の微分方程式の解は、
\(
\psi(x)=e^{\pm i{\sqrt{\frac{2mE}{\hbar^2}}x}}
\)

です。通常、波数\(k=\sqrt{\frac{2mE}{\hbar^2}}\)を定義して、
\(
\psi(x)=e^{\pm ikx}
\)

と書きます。連続状態を考えているので、エネルギー\(E\)は、
\(
\displaystyle E=\frac{\hbar^2k^2}{2m}
\)

です。

時間依存しないシュレーディンガー方程式が解けたので、時間依存するシュレーディンガー方程式の解\(\Psi(x,t)\)は、角周波数\(\omega=\frac{E}{\hbar}\)を用いて、
\(
\begin{align}
\Psi(x,t)&=e^{\pm ikx}e^{-i\frac{E}{\hbar}t} \\
&=e^{i(\pm kx-\omega t)}
\end{align}
\)

と記述されます。

\(e^{ikx}\)の解について


時間依存しない解が
\(e^{ikx}\)
に対する時間発展を考えます。
この時、時間依存する解は\(e^{i(kx-\omega t)}\)ですが、変形して、
\(
\displaystyle e^{ik(x-\frac{\omega}{k} t)}
\)

と書きます。
時刻\(t=0\)に対して、波動関数は\(e^{ikx}\)
時刻\(t=\Delta t\)に対して、波動関数は\(e^{ik(x-\frac{\omega}{k}\Delta t)}\)
です。よって、時刻\(t=0\to \Delta t\)の間に、関数が\(f(x)\to f(x-a)\)になった、つまり位置\(+a\)だけ平行移動された関数になった、とみることができます。
なので進行する波であると言えるのです。
また、その時の進行速度を考えますと、時間\(\Delta t\)の間に距離\(\frac{\omega}{k}\Delta t\)だけ進んだわけですから、速度\(v\)は、
\(
\displaystyle v=\frac{\frac{\omega}{k}\Delta t}{\Delta t}=\frac{\omega}{k}
\)

です。
グラフで時間経過を表せば、以下のようになります。
時刻t=0で\(e^{ikx}\)で表される波形(実線、実部:紫線、虚部:緑線)が時間経過して\(e^{ik(x-\frac{\omega}{k}\Delta t)}\)となると破線になります。
expikx_c
↑画像のgnuplotスクリプト

\(e^{-ikx}\)の解について


次に時間依存しない解が
\(e^{-ikx}\)
に対する時間発展を考えます。
この時、時間依存する解は\(e^{-i(kx+\omega t)}\)ですが、変形して、
\(
\displaystyle e^{-ik(x+\frac{\omega}{k} t)}
\)

と書きます。
時刻\(t=0\)に対して、波動関数は\(e^{-ikx}\)
時刻\(t=\Delta t\)に対して、波動関数は\(e^{-ik(x+\frac{\omega}{k}\Delta t)}\)
です。よって、時刻\(t=0\to \Delta t\)の間に、関数が\(f(x)\to f(x+a)\)になった、つまり位置\(-a\)だけ平行移動された関数になった、とみることができます。
なので後進する波であると言えるのです。
また、その時の進行速度を考えますと、時間\(\Delta t\)の間に距離\(-\frac{\omega}{k}\Delta t\)だけ進んだわけですから、速度\(v\)は、
\(
\displaystyle v=\frac{-\frac{\omega}{k}\Delta t}{\Delta t}=-\frac{\omega}{k}
\)

です。
グラフで時間経過を表せば、以下のようになります。
時刻t=0で\(e^{-ikx}\)で表される波形(実線、実部:紫線、虚部:緑線)が時間経過して\(e^{-ik(x+\frac{\omega}{k}\Delta t)}\)となると破線になります。
exp-ikx_c
↑画像のgnuplotスクリプト


波動関数は、
進行する波は虚部が先に来て後から実部(実部の位相が\(\pi/2\)遅れる)
後進する波は実部が先に来て後から虚部(虚部の位相が\(\pi/2\)遅れる)
のように見えます。
言い換えれば、波は虚部が先行している方向に進む、ということです。

3次元の場合


3次元でも同じことが言えます。
3次元自由粒子のハミルトニアンは
\(
\displaystyle \hat{H}=-\frac{\hbar^2}{2m}\left(\frac{\partial^2}{\partial x^2}+\frac{\partial^2}{\partial y^2}+\frac{\partial^2}{\partial z^2}\right)
\)

で変数分離可能なので、3次元の1つの方程式ではなく、1次元の3つの方程式に変換可能であり、x,y,z方向は独立に決められます。
なので、
\(
\displaystyle \Psi(x,y,z,t)=e^{ik_x(x-\frac{\omega_x}{k_x} t)+ik_y(y-\frac{\omega_y}{k_y} t)+ik_z(z-\frac{\omega_z}{k_z} t)}
\)

と表されます。
ここで\(\omega_x, \omega_y, \omega_z\)は、1次元の3つの方程式に分けた時のそれぞれのエネルギー固有値\(E_x,E_y,E_z\)に対応する角周波数を示しています(\(\omega_i=E_i/\hbar,\;\; (i=x,y,z)\))。(エネルギーは\(E=E_x+E_y+E_z\)です。)

上の場合は\({\bf k}=(k_x,k_y,k_z)\)方向に進行する波を表します。
通常は波数ベクトルとして書いて、
\(
\begin{align}
& \Psi({\bf r},t)= e^{i{\bf k \cdot r}-\omega t} \\
&\;\;\;\;\; {\small (\omega=\omega_x+\omega_y+\omega_z)}
\end{align}
\)

と書かれます。

1次元調和振動子の直接的解法

1次元調和振動子の良く知られた解法は3種類あります。

  1. 昇降演算子を利用し\(\hat{a}\psi(x)=0\)から求める方法
  2. 微分方程式を境界条件の下、直接解く方法
  3. ハイゼンベルグ方程式から解く方法

です。ここでは直接的に解く2番目の方法を紹介します。

対象にする微分方程式は一次元調和振動子の原子単位系での表現
\(
\displaystyle \left[-\frac{1}{2}\frac{d^2}{dx^2}+Ax^2\right]\psi(x)=E\psi(x)
\)

です。ここでAは正の実数です。
この微分方程式を境界条件\(\left. \psi(x)\right|_{x\to \pm\infty}=0\)の下で解きます。

まとめますと、エネルギー\(E_n\)と固有関数\(\psi_n(x)\)は、
\(
\begin{align}
E_n&=\sqrt{2A}\left(n+\frac{1}{2}\right)\\
\psi_n(x)&=\frac{\left(2A\right)^{1/8}}{\pi^{1/4}\sqrt{2^n n!}}H_n( {\scriptsize (2A)^{1/4}}x)\exp\left(-\sqrt{\frac{A}{2}} x^2\right)\;\;\;\; (n=0,1,2,\cdots)
\end{align}
\)

です。


導出


両辺を\(-2\)を掛けて式変形して
\(
\displaystyle \left[\frac{d^2}{dx^2}-2Ax^2+2E\right]\psi(x)=0
\)

常套手段ですが、この微分方程式を満たすべき関数\(\psi(x)\)の漸近形を考えます。
何とか解ける部分だけを解いていこうという方針です。

漸近での解


もしも\(x\to \infty\)だったならば、上式の\(2E\)の項は定数であるため、\(-2Ax^2\)と比べて無視できるようになるはずです。
ただし、微分の項\(\frac{d^2}{dx^2}\)は比較することはできないため無視はしないで置いておきます。
つまり、関数\(\psi(x)\)は\(x\to \infty\)の漸近で
\(
\displaystyle \left[\frac{d^2}{dx^2}-2Ax^2\right]\psi(x)=0
\)

の形の微分方程式を満足しなければなりません。
式変形して
\(
\displaystyle \frac{d^2}{dx^2}\psi(x)=2Ax^2\psi(x)
\)

関数を2回微分したら\(x^2\times\)(自分自身)が出てくる、という性質を持つ関数は
\(
\displaystyle \psi(x)= C \exp(Bx^2)
\)

という形で掛けそうです。実際にこの形の2階微分は
\(
\begin{align}
\frac{d\psi}{dx}&=C\cdot 2Bx \exp(Bx^2) \\
\frac{d^2\psi}{dx^2}&=C\cdot 2B(1+2Bx^2)\exp(Bx^2) \\
&\sim C\cdot 4B^2x^2\exp(Bx^2)
\end{align}
\)
です。今、2階微分で括弧内に定数”1”が出てきていますが、\(x\to \infty\)の漸近を考えているのでこの定数は\(2Bx^2\)と比較して無視されるはずです。よって最後の式変形をしています。

漸近での微分方程式に代入すると、
\(
\begin{align}
\frac{d^2}{dx^2}\psi(x)&=2Ax^2\psi(x) \\
C 4B^2x^2\exp(Bx^2)&=2Ax^2 C \exp(Bx^2) \\
2C(2B^2-A)x^2\exp(Bx^2)&=0
\end{align}
\)
上記の方程式は\(x\)の値に依らず常に満たされなければなりません。
よって、\(2C=0\)もしくは\(2B^2-A=0\)である必要があります。
\(C=0\)の解は波動関数が至る所でゼロである自明な解です。なので物理的には適しません。

※漸近でゼロであるだけで漸近じゃないところではokなんじゃないの?という疑問が出てきますが、\(C=0\)の場合、漸近で\(\psi(x)=0\)であり、その微分\(\frac{d\psi}{dx}=C 2Bx \exp(Bx^2)\)もゼロです。波動関数は全領域に渡って滑らかに続いていなければならないため、いつまでたっても波動関数はゼロです。ということは、存在確率がゼロ、有限にはならないつまらない解なのです。通常、\(x\to \infty\)の漸近で\(\psi(x)=0\)という条件は入れますが、その微分\(\frac{d\psi}{dx}\)もゼロという条件は入れません。

なので物理的に意味のある解は
\(
\begin{align}
2B^2-A=0 \\
\to B=\pm \sqrt{\frac{A}{2}}
\end{align}
\)
というBでなければなりません。よって波動関数\(\psi(x)\)は漸近で2つの解の線形結合として書かれ、
\(
\displaystyle \psi(x)= C_1 \exp\left(-\sqrt{\frac{A}{2}} x^2\right)+C_2 \exp\left(\sqrt{\frac{A}{2}}x^2\right)
\)

という形で記述できます。
しかし、境界条件\(\left. \psi(x)\right|_{x\to \pm\infty}=0\)を満たさなければならない場合、係数\(C_2=0\)である必要があります。よって波動関数の漸近形は
\(
\displaystyle \psi(x)= C\exp\left(-\sqrt{\frac{A}{2}} x^2\right)
\)

と求められます。

漸近ではないところの解


漸近では波動関数は
\(
\displaystyle \psi(x)= C\exp\left(-\sqrt{\frac{A}{2}} x^2\right)
\)

という形を満たさなければならないことが分かりました。残りは漸近ではない場所の解です。
定数変化法を使います。定数をxに依存する関数として全領域に渡る波動関数を
\(
\displaystyle \psi(x)= f(x)\exp\left(-\sqrt{\frac{A}{2}} x^2\right)
\)

と仮定します。この操作は\(\psi(x)\)に関する方程式から\(f(x)\)に関する微分方程式に焼直す操作です。
全領域に渡る波動関数の満たすべき微分方程式は
\(
\displaystyle \left[\frac{d^2}{dx^2}-2Ax^2+2E\right]\psi(x)=0
\)

なので、これに代入します。
計算過程で出てくる文字を減らすために、\(B=-\sqrt{\frac{A}{2}}\)として、\(B\)で記述していきます。
\(
\begin{align}
&\left[\frac{d^2}{dx^2}-2Ax^2+2E\right]\psi(x)= \\
&\;\;\;\;\;\;\;\left\{\frac{d^2f}{dx^2}+4Bx\frac{df}{dx}+(4B^2-2A)x^2+(2B+2E)f\right\}e^{Bx^2}=0
\end{align}
\)
と求める事ができ、\((4B^2-2A)=0\)を利用すると、\(f(x)\)の満たすべき微分方程式は、
\(
\displaystyle \frac{d^2f}{dx^2}+4Bx\frac{df}{dx}+(2B+2E)f=0
\)

と求まります。
この微分方程式の解はエルミートの微分方程式と呼ばれる形をしています。
エルミートの微分方程式は
\(
\displaystyle \frac{d^2 H(x)}{dx^2}-2x\frac{dH(x)}{dx}+2nH(x)=0
\)

の微分方程式を満たす直交多項式です[1]。
そのままエルミートの微分方程式を今回の問題に適応することはできません。なぜならば1階微分の係数が違っているからです。
なので変数変換により、\(x=\alpha y\)とおいて、一階微分の係数\(4B\)をちょうど\(-2\)にする定数\(\alpha\)を探しましょう。
\(
\begin{align}
\frac{d}{dx}&=\frac{dy}{dx}\frac{d}{dy}=\frac{1}{\alpha}\frac{d}{dy}\\
\frac{d^2}{dx^2}&=\frac{1}{\alpha^2}\frac{d^2}{dy^2}
\end{align}
\)
を代入して、
\(
\displaystyle \frac{1}{\alpha^2}\frac{d^2 f(x)}{dy^2}+4By\frac{df(x)}{dy}+(2B+2E)f(x)=0
\)

xとyが混在しているので気持ち悪いです。\(f(x)\to f(y)\)とおきましょう。そして最後に求められたyを用いたf(y)をxの式f(x)に焼直します。
両辺に\(\alpha^2\)を掛けると
\(
\displaystyle \frac{d^2 f(y)}{dy^2}+4B\alpha^2 y\frac{df(y)}{dy}+(2B+2E)\alpha^2 f(y)=0
\)

です。今、一階微分の係数を\(-2\)にする定数\(\alpha\)を探していました。
なので、
\(
\begin{align}
4B\alpha^2&=-2 \\
\alpha&=\pm\sqrt{-\frac{1}{2B}}=\pm\left(\frac{1}{2A}\right)^{1/4}
\end{align}
\)
となります。ここの\(\alpha\)のプラスマイナスはどちらでもいいです。
プラスに取ることにすれば、\(f(y)\)についての微分方程式は、
\(
\displaystyle \frac{d^2 f(y)}{dy^2}-2 y\frac{df(y)}{dy}-\left(\frac{B+E}{B}\right) f(y)=0
\)

であるため、エルミートの微分方程式の形になりました。
エルミートの微分方程式は\(f(y)\)に掛かる項、すなわち\(-\left(\frac{B+E}{B}\right)\)の部分が\(2n\), (\(n=0,1,2,\cdots\))である時だけ解を満たします。
よって、
\(
\begin{align}
-\frac{B+E}{B}&=2n\\
E=E_n&=-2B(n+\frac{1}{2})\\
&=\sqrt{2A}(n+\frac{1}{2})\;\;\; (n=0,1,2,\cdots)
\end{align}
\)
を満たすエネルギー\(E=E_n(n=0,1,2,\cdots)\)でなければ解は存在しません。
その時の解\(f(y)\)も\(n\)で番号付けされて、
\(
f(y)=H_n(y)\;\;\; (n=0,1,2,\cdots)
\)

と記述され、右辺はエルミート多項式と呼ばれます。
エルミート多項式\(H_n(y)\)は
\(
\begin{align}
H_0(y)&= 1\\
H_1(y)&= 2y\\
H_2(y)&= 4y^2-2\\
H_3(y)&= 8y^3-12y\\
H_4(y)&= 16y^4-48y^2+12\\
H_5(y)&= 32y^5-160y^3+120y\\
H_6(y)&= 64y^6-480y^4+720y^2-120\\
\ldots
\end{align}
\)
と言う関数です。三項漸化式
\(
H_n(y)=2yH_{n-1}(y)-2(n-1)H_{n-2}(y)
\)

を用いて\(n=2\)以上のエルミート多項式は記述されます。

以上より微分方程式を満たす解とエネルギーは、整数(\(n=0,1,2,\cdots\))で番号付けされた解を持つことが分かりました。
\(y\)を変数変換に従い\(x\)に戻します。
\(
\begin{align}
x&=\alpha y \\
&\alpha=\left(\frac{1}{2A}\right)^{1/4}
\end{align}
\)
なので、
\(H_n(y)=H_n(x/\alpha)\)
です。
漸近形での波動関数と合わせて、本当に求めたかった解\(\psi(x)\)は、
\(
\begin{align}
\displaystyle \psi(x) &= f(x)\exp\left(-\sqrt{\frac{A}{2}} x^2\right) \\
&=C H_n(x/\alpha)\exp\left(-\sqrt{\frac{A}{2}} x^2\right)
\end{align}
\)

であることが分かりました。定数\(C\)の自由度がありますが、これは物理的に合った条件を用いてなされます。

おさらい


解きたかった問題は
\(
\displaystyle \left[-\frac{1}{2}\frac{d^2}{dx^2}+Ax^2\right]\psi(x)=E\psi(x)
\)

を境界条件\(\left. \psi(x)\right|_{x\to \pm\infty}=0\)の下で解くことでした(A>0)。

この微分方程式を何とか解くために、まず\(x\to \pm\infty\)で解の満たすべき形を考えました。
その後、全領域で境界条件を満たす解を定数変化法によって求めました。
計算の結果、微分方程式の解は与えられた適当なエネルギー\(E\)に対していつも存在するわけではなく、特別な\(E=E_n\)で番号付けされた値でなければ解が存在しないことが分かりました。
\(n\)番目のエネルギー\(E_n\)と関数\(\psi_n(x)\)は、
\(
\begin{align}
E_n&=\sqrt{2A}\left(n+\frac{1}{2}\right) \\
&\psi_n(x)= C_n H_n(x/\alpha)\exp\left(-\sqrt{\frac{A}{2}} x^2\right) \\
&\;\;\;\;\alpha=\left(2A\right)^{-1/4}
\end{align}
\)

です。与えられた条件下では規格化までは追求していないのでここでお仕舞いですが、規格化を考えてみましょう。

規格化


波動関数\(\psi_n(x)\)を規格直交化(正規直交化とも言う)します。
規格化は、全空間での存在確率
\(
\displaystyle \int_{-\infty}^{\infty}|\psi(x)|^2dx
\)

を1にすること。
直交化は異なる状態間をかけ合わせて積分した時ゼロになること・・・すなわち内積をゼロにします。
\(
\displaystyle \int_{-\infty}^{\infty}\psi_m^*(x)\psi_n(x)dx=C\cdot\delta_{mn}
\)

ここで\(C\)は任意の定数で、\(\delta_{mn}\)はクロネッカーのデルタを表します。
両方を考慮すると、規格直交化とは、
\(
\displaystyle \int_{-\infty}^{\infty}\psi_m^*(x)\psi_n(x)dx=\delta_{mn}
\)

を満たすように係数を決定します。

実際に波動関数を代入すると、
\(
\begin{align}
\int_{-\infty}^{\infty}\psi_m^*(x)\psi_n(x)dx&= C_m^*C_n \int_{-\infty}^{\infty} H_m(x/\alpha)H_n(x/\alpha)\exp\left(-\sqrt{2A} x^2\right)dx
\end{align}
\)

変数変換を行います。
\(
\begin{align}
x&=\alpha y \\
&\alpha=\left(\frac{1}{2A}\right)^{1/4}
\end{align}
\)
より、
\(
\begin{align}
& C_m^*C_n \int_{-\infty}^{\infty} H_m(x/\alpha)H_n(x/\alpha)\exp\left(-\sqrt{2A} x^2\right)dx \\
&= C_m^*C_n \int_{-\infty}^{\infty} H_m(y)H_n(y)\exp\left(-\sqrt{2A}\alpha^2y^2\right)\cdot \alpha dy \\
&= C_m^*C_n\alpha \int_{-\infty}^{\infty} H_m(y)H_n(y)\exp\left(-y^2\right) dy \\
\end{align}
\)
ここで、エルミート多項式の特性を使います。エルミート多項式の関係式の中に、
\(
\displaystyle \int_{-\infty}^{\infty} H_m(x)H_n(x)\exp\left(-y^2\right) dx=2^n n!\sqrt{\pi}\delta_{mn}
\)

という関係式があります。
なので、これを使うと、
\(
\displaystyle \int_{-\infty}^{\infty}\psi_m^*(x)\psi_n(x)dx=C_m^*C_n\alpha 2^n n!\sqrt{\pi}\delta_{mn}
\)

と求められます。今、\(({\mbox 右辺})=\delta_{mn}\)にしたいので、係数について、
\(
\begin{align}
C_m^*C_n\alpha 2^n n!\sqrt{\pi}\delta_{mn}&=\delta_{mn} \\
(C_m^*C_n\alpha 2^n n!\sqrt{\pi}-1)\delta_{mn}&=0
\end{align}
\)
にするように\(C_n\)を決定すればよいことになります。
クロネッカーのデルタより、\(m=n\)以外の時は係数に依らずゼロになるため、自明です。
なので、\(m=n\)の時,
\(
\begin{align}
|C_n|^2\alpha 2^n n!\sqrt{\pi}&=1\\
\to |C_n|&=\sqrt{\frac{1}{\sqrt{\pi}2^n n!\alpha}} \\
C_n&=\frac{1}{\sqrt{\alpha}\pi^{1/4}\sqrt{2^n n!}}\cdot \exp(i\theta)
\end{align}
\)
と規格化定数が求められます。位相\(\exp(i\theta)\)の任意性がありますが、この位相に特に意味は無く、自由に選べます。よって、式が簡単になる\(\theta=0\)を取ることにします。
以上の結果から、規格直交化された波動関数は\(\psi_n(x)\)は、
\(
\begin{align}
\psi_n(x)&=\frac{1}{\sqrt{\alpha}\pi^{1/4}\sqrt{2^n n!}}H_n(x/\alpha)\exp\left(-\sqrt{\frac{A}{2}} x^2\right)\\
\alpha&=\left(2A\right)^{-1/4}
\end{align}
\)
となります。\(\alpha\)を使わず、\(A\)だけで表せば、
\(
\displaystyle \psi_n(x)=\frac{\left(2A\right)^{1/8}}{\pi^{1/4}\sqrt{2^n n!}}H_n( {\scriptsize (2A)^{1/4}}x)\exp\left(-\sqrt{\frac{A}{2}} x^2\right)
\)

です。

参考文献

[1]小野寺嘉孝著「物理のための応用数学」裳華房 p.71

gnuplotで配列を使う

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

とスクリプトを書き、プロットすると
arraygnuplot2
という画像が得られます。

配列の利用


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)

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

アニメーション
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)
}

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

参考文献


Operate on 1D arrays in gnuplot

[adsense2][adsense3]

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

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

Ylm_eq_Y10_c

これを角度\(\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です。実際に書いて、出力するとそれぞれ以下のような画像が出力されます。
y10_normaly10_thetay10_phi
※左から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")

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

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

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

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

本来は多価関数なのですが、それゆえに上記方法では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で極座標形式\((r,\theta)\)で、こんな感じの2次元のグラフ
GPB-Concept-Dgrm-lg
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)

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

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

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

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

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

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

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

でこのようになります。
lorentz_sph_c

仕組みは極座標とは
\(
\begin{align}
x&=r\cos\theta \\
y&=r\sin\theta
\end{align}
\)
なので、媒介変数urに、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

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

これを解決するには

set pm3d depthorder

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

参考
gnuplot demo script: surface2.dem

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

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

目的は、

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

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

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


電動ガンの場合です。
バレル内部1
空気抵抗は初速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弾の速度、ピストンの速度を意味します。

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

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

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弾とバレル内部に隙間があることが写真で確認できます。

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

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

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


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

例(音源について)

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

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

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

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

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

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

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

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

画像の使用


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

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


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

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

    保護期間の満了により公有化された著作物。
    継承者不在により公有化された著作物。
    権利放棄により公有化された著作物。

がそれに当たります。

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

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

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

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


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

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

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

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

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

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


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

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

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

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


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

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

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

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

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


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


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


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


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


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

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


自分の作品にクリエイティブ・コモンズ・ライセンスを付けたい場合は、どこかに申請する必要など無く、自分で言い張ればいいんです。なぜならばあなたが作品を作り出した時に著作権は発生しているのですから。その作品にCCライセンスを付けるか否かは貴方の自由です。
もしもCCライセンスを付けたい場合はこちらのページに行けばどのように付ければいいかわかります。

※後からクリエイティブ・コモンズ・ライセンスの変更は行ってはなりません。
なぜならば、例えばはじめCC0ライセンスで公開し、皆が使い始めたところライセンスを変えたりしたらまるで使い物にならなくなってしまうからです。
なので、重要な事はその権利者、例えばページを運営している人が信用、信頼できるかどうかを見極めることが重要です。

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

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

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

実例


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

表示のみの場合CC BY 2.0

りんご -wikipedia です。
Red_Applere_compressed
©Abhijit Tembhekar(2009)/Adapted/CC BY 4.0

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

改変したのでAdaptedを付けました。また、この段階でこの画像は2次創作物になり、その時点からCCライセンスは僕が決める事ができます。表示と同じ、もしくはそれよりも制限が強いものに指定することができます(FAQ よくある質問と回答)。従って、この上の改変した画像に対し、All rights reservedとすることも可能なのです。今回は元のものと同じCC BY 2.0 表示にしました。

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

球面調和関数 -wikipedia にある画像を利用します。
1200px-Harmoniki_compressed
©Sarxos(2007)/Adapted/CC BY-SA 3.0

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

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

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

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

自然 wikipediaにあるパブリックドメインの画像です。
Bachalpseeflowers
クレジットの表記はまるで必要ありません。営利目的であっても、許可を得ずに複製・改変・翻案・配布することが可能です。

これはCC0の例です。
https://commons.wikimedia.org/wiki/File:%22A_mixed_flock_of_ducks_and_geese_fly_from_a_wetland_area%22.jpg
-A_mixed_flock_of_ducks_and_geese_fly_from_a_wetland_area-
CC0もクレジットの表記はまるで必要ありません。営利目的であっても、許可を得ずに複製・改変・翻案・配布することが可能です。

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

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

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

補足

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

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

gnuplotとアニメーション

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

  1. 解析解の表示
    1. 1次元のアニメーション
    2. 2次元のアニメーション
    3. 点をグラフ上に沿って動かす
    4. グラフをだんだんと書く
  2. データの表示
    1. 複数のデータを次々とプロット
    2. 複数のデータの同時プロット
    3. グラフをだんだんと書く
  3. 応用
    1. グラフを回転させる
    2. 大きさの違う球体
    3. PSO2の転送装置
    4. 「見せてあげよう!ラピュタの雷を!!」
    5. 魔女の宅急便
    6. 第5使徒ラミエル1
    7. 第5使徒ラミエル2
    8. 魔貫光殺砲
    9. マリオ
    10. SCP-1968 「世界を包む逆因果の円環」っぽい数式
    11. FGO召喚
    12. パンジャンドラム
  4. 動画の作成について
    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とか変更してください。実行するとこんな感じ。
gauss_low

上の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次元だろうとこのように簡単にできます。
gauss2d_analytic
上の動画のスクリプトはこれです。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

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


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

#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点だけを表示することができます。

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


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

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プロットがあり、それを回転させてみるような動画がほしいとします。例えばこんな感じのアニメーションが作りたいとします。
Riemann
これはリーマン面と呼ばれる関数です。(カラーバーの消去は ‘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すればいいです。

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


gnulotofdata_c
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)

とすればいいです。

大きさの違う球体


sphere3

PSO2の転送装置


PSO2の転送装置gnuplot

「見せてあげよう!ラピュタの雷を!!」


bulse2

魔女の宅急便


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

kiki

第5使徒ラミエル1


ramiel

第5使徒ラミエル2


ramiel_trans

魔貫光殺砲


mk

マリオ


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

mario

[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”という動画が作られます。