関数の絶対値の実軸上での微分

関数の絶対値の微分です。

まとめ


実数値関数の場合
\(
\begin{align}
\frac{d}{dx}|f(x)|&=\frac{f}{|f|}\frac{df}{dx} \\
\frac{d^2}{dx^2}|f(x)|&=\frac{f}{|f|}\frac{d^2f}{dx^2}
\end{align}
\)

複素関数の、実軸上での微分の場合
\(
\begin{align}
\frac{d}{dx}|f(x)|&=\frac{1}{2|f|}\left(\frac{df^*}{dx}f+f^*\frac{df}{dx}\right) \\
\frac{d^2}{dx^2}|f(x)|&=\frac{1}{2|f|}\left[ \frac{d^2f^*}{dx^2}f + f^*\frac{d^2f}{dx^2}+\frac{1}{2}\left(\frac{1}{f^*}\frac{df^*}{dx}-\frac{1}{f}\frac{df}{dx}\right)\left(f^*\frac{df}{dx}-f\frac{df^*}{dx}\right)\right]
\end{align}
\)

導出


関数の絶対値\(|f(x)|\)は、
\(
\displaystyle |f(x)|=\sqrt{f^*(x)f(x)}
\)

と書かれます。これを利用して絶対値を書き換えます。
\(
\begin{align}
\frac{d}{dx}|f(x)|&=\frac{d}{dx}\sqrt{f^*f} \\
&=\frac{1}{2}(f^*f)^{-1/2}\cdot \left(\frac{df^*}{dx}f+f^*\frac{df}{dx}\right) \\
&=\frac{1}{2|f|}\left(\frac{df^*}{dx}f+f^*\frac{df}{dx}\right)
\end{align}
\)
と求められます。また、2階微分は1階微分を微分すればよいので、
\(
\begin{align}
\frac{d^2}{dx^2}|f(x)|&=-\frac{1}{2}\frac{1}{|f|^2}\cdot\frac{1}{2|f|} \left(\frac{\partial f^*}{dx}f+f^*\frac{\partial f}{dx}\right)^2+\frac{1}{2|f|}\left(\frac{d^2f^*}{dx^2}f+2\frac{df}{dx}\frac{df^*}{dx}+f^*\frac{d^2f}{dx^2}\right) \\
&=\frac{1}{2|f|}\left[-\frac{f}{2f^*}\left(\frac{df^*}{dx}\right)^2-\frac{f^*}{2f}\left(\frac{df}{dx}\right)^2-\frac{df}{dx}\frac{df^*}{dx}+\frac{d^2f^*}{dx^2}f+2\frac{df}{dx}\frac{df^*}{dx}+f^*\frac{d^2f}{dx^2}\right] \\
&=\frac{1}{2|f|}\left[ \frac{d^2f^*}{dx^2}f + f^*\frac{d^2f}{dx^2}+\frac{1}{2}\left(\frac{1}{f^*}\frac{df^*}{dx}-\frac{1}{f}\frac{df}{dx}\right)\left(f^*\frac{df}{dx}-f\frac{df^*}{dx}\right)\right]
\end{align}
\)

\(f(x)\)が実関数の場合は\(f^*=f\)なので、
\(
\begin{align}
\frac{d}{dx}|f(x)|&=\frac{1}{2|f|}\left(\frac{df}{dx}f+f\frac{df}{dx}\right) \\
&=\frac{f}{|f|}\frac{df}{dx}
\end{align}
\)
であり、
\(
\begin{align}
\frac{d^2}{dx^2}|f(x)|
&=\frac{1}{2|f|}\left[ \frac{d^2f}{dx^2}f + f\frac{d^2f}{dx^2}+\frac{1}{2}\left(\frac{1}{f}\frac{df}{dx}-\frac{1}{f}\frac{df}{dx}\right)\left(f\frac{df}{dx}-f\frac{df}{dx}\right)\right]\\
&=\frac{f}{|f|}\frac{d^2f}{dx^2}
\end{align}
\)

もしくは、符号関数を\({\rm sgn}(x)\)用いて
\(
\begin{align}
\frac{d}{dx}|f(x)|&={\rm sgn}(f(x))\frac{df}{dx} \\
\frac{d^2}{dx^2}|f(x)|&={\rm sgn}(f(x))\frac{d^2f}{dx^2}
\end{align}
\)
と表されます。符号関数を用いる定義のほうが良いと思います。なぜなら、符号関数は\({\rm sgn}(0)=0\)と、ちゃんと定義されているからです。

tips)

正則関数と非正則関数について
(複素平面上のいたる点で一階微分可能)=(複素平面上のいたる点で無限回微分可能)=(正則関数、コーシー・リーマンの関係式が成立)
それ以外の関数は非正則関数と呼ばれる。

参考文献

複素関数の絶対値の微分 -biochem_fan’s note

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

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


スポンサーリンク


疑似配列の四則演算


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

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

スポンサーリンク

 

角度に依存して色を付ける(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つの方法でプロットします。


スポンサーリンク


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θをとみなせば書くことができます。

スポンサーリンク

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)
バレル内部でのBB弾の方程式←今ここ


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

いくつか仮定があります。

  1. 気体の状態方程式\(PV=nRT\)が成立している(空気を理想気体として扱う, 補遺1)
  2. BB弾とバレルの間に空気の漏れは無い(補遺2)

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

ピストン-BB弾間の体積\(V\)はピストンの位置とBB弾の位置の2つに依存します。
なので、ピストンの位置による時間変化を\(x_0(t)\)、BB弾の位置を\(x(t)\)、バレルの断面積を\(S\)とすると
\(
V(t)=S\cdot (x(t)-x_0(t))
\)

と記述されます。
気体の状態方程式より、
\(
PV=nRT
\)

であり、右辺を定数\(c_0\)と置き、圧力\(P(t)\)をBB弾に掛かる力F(t)、すなわち\(P(t)=\frac{F(t)}{S}\)とすれば、
\(
\begin{align}
\frac{F(t)}{S}\cdot S\cdot (x(t)-x_0(t))=c_0 \\
F(t)=\frac{c_0}{x(t)-x_0(t)}
\end{align}
\)
と書けます。よって、空気抵抗と共に書けば、運動方程式は
\(
\displaystyle m\frac{d^2}{dt^2}x(t)=\left(\frac{c_0}{x(t)-x_0(t)}-P_0\cdot S\right)-6\pi \eta R \cdot v(t)-\frac{1}{2}C_d \rho \pi R^2 \cdot v^2(t)
\)

と導けます。ここで、\(P_0\)は大気圧、\(\eta\)は粘性率、\(R\)はBB弾の半径、\(C_d\)は抗力係数、\(\rho\)は空気の密度を意味します。
空気抵抗に関する詳細は球体の空気抵抗と係数をご覧ください。

あとはこの運動方程式を解いて、BB弾が出てくるバレル位置において適当な速度になるような物理的な意味を十分に持つ\(c_0, x_0(t)\)を入れて計算すればいいのです。

ここで1つ言えることは、外部と内部の圧力が一定になる最適なバレル長というものが存在する、ということです。
このバレル長を使えば、BB弾が射出されたときに、銃口で気圧差は無く、それによって生じる不確定な要素を最小にすることができるでしょう。

そのバレル長とは、式(4)の右辺が0になる時であり、計算によって求めることができます。
いくつになるかは実際求めてみないと想像できません。数cmであるのか数十mであるのか、はたまた30cmとかその位なのか・・・。

ちなみに実銃ではほぼ延々に火薬によって非常に多くのガスが生成され、加速されるためこのようなことは起こらないでしょう(あるかもしれませんが、とてつもなく長いバレルでなければ実現しないと思われます)。

補遺


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

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

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

補遺2
BB弾とバレル内側との間における空気の漏れはどの程度あるのかはわかりません。

バレルの内径:6.08mm (東京マルイのものはこれらしいですが、信頼できるページが見当たりませんでした。)
BB弾の直径:5.95mm [2]

完全に理想気体として、時間を考慮しないで扱うのであれば、少しでも空間が空いていればどれだけ早く押し込もうとも気圧は一定になろうとするため、BB弾に掛かる力はゼロです。
なので時間を考慮した方程式が必要でしょう。
連続の式
\(
\displaystyle \frac{\partial \rho}{\partial t}+\nabla \cdot \vec{j}=0
\)

を用いれば時間と共に漏れ出る空気量を知ることが出来そうです。
ただ、まずは簡単のため、空気が漏れない場合をまず考えましょう。

簡単に空気の漏れを表現するには、あたかもピストン-BB弾間の体積\(V\)が時間と共に増加していく、と考えることである程度の効果を入れることができると考えられます。

実際に式を解いてから、空気の漏れが起こる場合と起こらない場合でどれほどの差が出てくるのか?を考えることでこの効果を入れるべきか入れないべきかを考えるのがいいでしょう。

良く空気の漏れがーと言われているので本当かどうか理論的に検証するべきでしょう。

参考文献


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

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