sikino のすべての投稿

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

Fortran90で構造体を使う

2021年9月2日 sikino コメントする

Fortran90で自分で型を定義し、それを使用します。
本稿でいう”型”は、一般には構造体と呼ばれます。
Fortran使用者にとっては”型”という名称が分かりやすいでしょう。

整数型、倍精度型などと同列にあるもの、と考えて差し支えないです。

本稿では基本的なことしか書いていません。

二階偏微分の計算の時でも構造体の考えを使用していますので、もう少し例が欲しい方はこちらもどうぞ。
Hyper-dual numbersによる二階偏微分の計算 -シキノート

概要

簡単な実例


module CEOM_mod

  implicit none

  type Particle

     double precision::m,x,v

  end type Particle

end module CEOM_mod



program main

  use CEOM_mod

  implicit none

  type(Particle)::p1



  p1%m  = 1d0

  p1%x  = 2d0

  p1%v  = 3d0



  write(6,*)p1%m,p1%x,p1%v



  stop

end program main

実行結果

$ ./a.out
1.0000000000000000 2.0000000000000000 3.0000000000000000
$

説明

上は、粒子p1が質量m, 位置x, 速度vを持つ、ということを想定しています。

そこで、Particleという新しい”型”（複素数型みたいに、実部と虚部という２つの独立な変数を１つの複素数型変数が持つ、みたいな変数の型）を定義しました。
Particle型を持つ変数は、すべてm, x, vを内包するというイメージです。

mainプログラムではparticle型が定義されているモジュールを呼び出すために

use CEOM_mod

という文を入れ、実際にparticle型の変数p1を定義するために

type(Particle)::p1

と書いています。変数p1が定義されると自動的にp1が持つm,x,vの固有な値も定義されます。固有の値をいじる場合は

p1%m = 1d0

と書いて、変数に%をつけてm,x,vのいずれかを書いて代入します。

例

配列


module CEOM_mod

  implicit none

  type Particle

     double precision::m,x,v

  end type Particle

end module CEOM_mod



program main

  use CEOM_mod

  implicit none

  integer::i,N

  type(Particle),allocatable::p(:)



  N=5

  allocate(p(1:N))

  do i=1,N

     p(i)%m  = 1d0

     p(i)%x  = dble(i)

     p(i)%v  = 0.1d0*i

  enddo



  do i=1,N

     write(6,*)i,p(i)%m,p(i)%x,p(i)%v

  enddo



  stop

end program main

実行結果

$ ./a.out
1 1.0000000000000000 1.0000000000000000 0.10000000000000001 
2 1.0000000000000000 2.0000000000000000 0.20000000000000001 
3 1.0000000000000000 3.0000000000000000 0.30000000000000004 
4 1.0000000000000000 4.0000000000000000 0.40000000000000002 
5 1.0000000000000000 5.0000000000000000 0.50000000000000000
$

サブルーチンの利用


module CEOM_mod

  implicit none

  type Particle

     double precision::m,x,v

  end type Particle

end module CEOM_mod



subroutine twicemass(N,p)

  use CEOM_mod

  implicit none

  integer,intent(in)::N

  type(Particle),intent(inout)::p(1:N)



  integer::i

  

  do i=1,N

     p(i)%m  = 2 * p(i)%m

  enddo

  

  return

end subroutine twicemass



program main

  use CEOM_mod

  implicit none

  integer::i,N

  type(Particle),allocatable::p(:)



  N=3

  allocate(p(1:N))

  do i=1,N

     p(i)%m  = 1d0

     p(i)%x  = dble(i)

     p(i)%v  = 0.1d0*i

  enddo



  call twicemass(N,p)

  

  do i=1,N

     write(6,*)i,p(i)%m,p(i)%x,p(i)%v

  enddo



  stop

end program main

実行結果

$ ./a.out
1 2.0000000000000000 1.0000000000000000 0.10000000000000001 
2 2.0000000000000000 2.0000000000000000 0.20000000000000001 
3 2.0000000000000000 3.0000000000000000 0.30000000000000004 
$

補足

オブジェクト指向になれている方は、構造体の中身を外部から直接アクセスできる点で褒められた書き方ではないでしょう。
Fortran90にオブジェクト指向の考え方を持ってくるのは、計算速度というFortran90の長所を潰してしまうので、このくらいで良いでしょう。

数学, 物理学

古典的な一次元波動方程式のグリーン関数

2021年8月8日 sikino コメントする

本稿のpdfはこちらをどうぞ
https://slpr.sakura.ne.jp/qp/wp-content/uploads/2021/08/Green_1dclassical5.pdf

問題

古典的な一次元波動方程式で、非斉次項を持つ問題

を解くことを考えます。この形は、有限範囲内しか相互作用を引き起こさない物に十分遠方から波が入射してきた場合を考える際に現れたりします。

この形の方程式が出てくる状況を考えてみましょう。

の問題を考えます。移項して

ですので、グリーン関数を用いると

となります。ここで、

を満たす関数（一般解）であり、グリーン関数を
で与えられます。ここを満たす関数（特殊解）としました。
式(4)の右辺に解$f(x,t)$が含まれていますが、このままの形でとどめておきます。この形にしておくと、ほかの式に$f(x,t)$をダイレクトに代入できたり、近似を使うときに便利な形です。

（式（7）はありません）

本稿の目的は、$\hat{D}$が$\hat{D}=\frac{\partial^2}{\partial t^2}-a\frac{\partial^2}{\partial x^2}$で与えられている場合に、$G(x,t; x’,t’)$の具体的な形を導出することです。

導出

今、解きたい問題は、

です。デルタ関数をそのまま扱う場合、積分が絡んでこないと扱うのが難しいです。そうでなければ、フーリエ変換を用いて波数・周波数空間で解いていくのが良いでしょう。グリーン関数や、デルタ関数のフーリエ変換を考えると

となります。ここで、基底関数が$e^{ikx}, e^{-i\omega t}$と符号が異なることに注意しましょう。
また、積分区間はいつも負の無限大から正の無限大であり、いちいち書くのは面倒なので省略しています。
フーリエ変換として考えるのではなく、基底関数へ射影した時の係数と考えましょう。

実際に式(8)に代入すると

となります。

ここで$e^{i[k(x-x’)-\omega(t-t’)]}$はゼロにならないので、その前の項[ ]内がゼロにならなければなりません。よって

が導けます。よって式(9)に代入して

が計算できれば位置と時間のグリーン関数が得られます。

さて、この積分を計算してみましょう。
右辺の$e^{i\omega t}/\omega$の形を持つ積分は複素関数論において頻出する積分です。この積分の特徴として、$\omega$に渡る積分を実軸上で実行する際に、たまたま$\omega=\pm \sqrt{a}k$に等しい点を通ってしまうため、実軸上では被積分関数が発散する、という特徴があります。この発散する点は極と呼ばれており、積分を行う際に積分経路をうまく選んで極を迂回しなければなりません。しかし、厄介なことに極を上に回るか下に回るかで積分結果が変わってしまいます。

もし、天才でしたら当たり前のように適切な積分経路をすぐに思いつくでしょうが、我々凡人には具体的に試してみるくらいしか思いつきません。ですのでちょっとズルしながら因果律を満たす解を探しましょう。

今、我々はグリーン関数について考えています。ということは、計算した結果の境界条件は、因果律を満たすような結果が欲しいわけです。
つまり、グリーン関数が得られた場合、$\theta(t-t’)$の形が出てくるであろうことは予想できます。そのため、これを考慮して極を回る方向を考えてグリーン関数を求めていきましょう。

式(14b)の第一項

について考えます。ヘヴィサイド関数を積分表記で表すと、

と書けることを使います。ここで$i$は虚数単位です。もう嬉しいですね。かなり似通った形であることが分かるでしょう。
式(16)の形に持っていくためには変数変換を行えば良さそうです。つまり、極は上に回れば因果律を満たすグリーン関数が得られそうです。

では具体的に計算していきましょう。変数変換$y=\omega+\sqrt{a}k$を使って、

となります。分母に$+i\varepsilon$を加えるようにして

と求めることができます。

続いて式(14b)の第二項

についても同様に計算していきましょう。変数変換$y=-\omega+\sqrt{a}k$を使って、

ヘヴィサイド関数を考えると、

の形が使えそうです(式(16)の$\omega\to-\omega$の変数変換ですね)。ですので、

と求められます。

式(18c)と式(22c)より、式(14b)に代入すると

となります。

さて、$\frac{\sin x}{x}$の形はsinc関数と呼ばれているほど面白い性質を持つ関数です。そのフーリエ変換は矩形になることが知られています。その計算は実際に
\begin{align}
\int_{-\infty}^{\infty}\frac{\sin(k)}{k}e^{ikx}\frac{dk}{2\pi}=\frac{1}{2}\Bigl[\theta(x+1)-\theta(x-1)\Bigr]
\end{align}
となります。

立ち戻り、変数変換$\kappa=\sqrt{a}k(t-t’)$を考えれば、

と、グリーン関数が求められました。グリーン関数は$x-x’, t-t’$の形をしていることが分かります。そのため、

と書いても良いことが分かるでしょう。ここで、

です。

まとめと補足

まとめ

まとめますと、偏微分方程式

の因果律を満たす解は、

となります。以上から、波動方程式

の因果律を満たす一般解は

となります。

補足

補足します。因果律を満たすグリーン関数である、ということが暗に意味される場合、積分(29a)の時間に渡る積分はグリーン関数に含まれるヘヴィサイド関数を先に利用して、

と積分範囲を組み込んで表記する場合があります。実際の例が式(29c)であり、その時間積分の積分範囲のようになることを見越している、ということです。

少しグリーン関数の物理的な意味を考えてみましょう。
グリーン関数は、ある特定の時刻、ある特定の位置でデルタ関数という非常に特殊な衝撃が生じた結果を記述します。
言い換えれば、偏微分方程式で表されている系が、デルタ関数によってどのような応答を起こすか？という関数です。

式(27b)を見てみると、１つのヘヴィサイド関数$\theta(t-t’)$と１つの矩形関数$\text{rect}\left(\frac{x-x’}{2\sqrt{a}(t-t’)}\right)$が使われています。
$\theta(t-t’)$の部分は、まさに因果律を示しており、衝撃が加わる前の時刻$(t-t’ \lt 0)$においては系の応答は無いよ、ということを表しています。至極全うでしょう。

矩形関数の部分は、衝撃が生じた後にその衝撃は波及する範囲を表しています。具体的に、時刻$t=t’$, 位置$x=x’$において衝撃が生じた場合、矩形関数の中身が値を持つ範囲は$[x’-\sqrt{a}(t-t’), x’+\sqrt{a}(t-t’)]$の範囲です。つまり、衝撃が時々刻々と広がっていることを示しています。グリーン関数を図示すると、このような感じです。

衝撃が加わる$t=0$を境にグリーン関数が値を持ち始め、それが時間とともに広がっていくのが分かると思います。

ついでにグリーン関数の広がる速度を計算してみましょう。

矩形関数の中身が値を持つ範囲$x$は、時間$\Delta t=t-t’$の間に、距離$x-x’=\pm \sqrt{a}(t-t’)$だけ波及しますから、衝撃の結果が速度$\pm\sqrt{a}$で広がっていることが分かります。
これは、古典的な波動方程式でよく知られているように、系の速度が

の系の速度は$\pm\sqrt{a}$である、という事実と一致します。
つまり、衝撃が加わった時刻以外では、系そのものの性質しか存在しませんので、衝撃が波及していく速度が一致することは何ら不思議ではありません。むしろ、一致しなければなりません。

数学, 物理学

デルタ関数、ヘヴィサイド関数に関係する数式

2021年8月8日 sikino コメントする

本稿のpdfはこちら
https://slpr.sakura.ne.jp/qp/wp-content/uploads/2021/08/Green_related_functions1.pdf

定義

ヘヴィサイド関数

符号関数

矩形関数

デルタ関数

デルタ関数$\delta(x)$は、

を満たす関数です。この定義式を満たす場合、デルタ関数は以下の性質を持つことが分かります。

フーリエ関数

性質、関係式

積分表示

定積分が関係する場合

デルタ関数

ヘヴィサイド関数

sinc関数に関係する式

畳み込み積分

数学, 物理学

位相速度、群速度とは！？調べてみました！！(パロディ)

2021年6月13日 sikino コメントする

近頃SNSやインターネットでよく見る言葉、「位相速度」「群速度」。

元ネタや原理を分かった気になって使っている方、多いのではないでしょうか。
せっかくなので調べてみました！！！

（数式で殴りに殴ったものを以下のページに書きました。真面目に詳細を知りたい方はこちらをどうぞ。）
https://slpr.sakura.ne.jp/qp/group-phase-velocity-delay/

そもそも「位相速度」「群速度」って何？
学会騒然。ある波数近傍だけ値を持つ場合にテーラー展開!?
位相速度はどうやって!?驚きの仮定
簡単だと思いました？注意です！
いろいろなところに群速度に関する議論が！
最後に
Note

そもそも「位相速度」「群速度」って何？

なにはともあれwikipedia!!
関連するキーワードをササッと知るには強い味方です。
「位相速度」「群速度」を調べるとこのように説明されていました。

群速度（ぐんそくど、英: group velocity）とは、複数の波を重ね合わせた時にその全体（波束）が移動する速度のことである。

群速度　-wikipedia

位相速度（いそうそくど、英語: phase velocity）は、位相、すなわち波の山や谷の特定の位置が移動する速度のことである。

位相速度　-wikipedia

どうやら、波束や波の速度を表現する際に都合が良い言葉なようです。

波束とは何なのでしょう…？また変な言葉が出てきました。wikipediaの出番です！！

波束（はそく、英: wave packet, wave train）は、局所的に存在する波うち/波動であり、移動する1個の波動の塊のようにふるまう。

波束 -wikipedia

うーん、つまり一回だけ波打って、ある程度まとまった波みたいなものですかね。
水面に水が一滴だけ落ちたときに発生する波みたいなものでしょう。

具体的にどういう意味なのでしょうか。気になりますよね！

皆さんも群速度、位相速度の導出が気になったら「tweet」、お願いします！

学会騒然。ある波数近傍だけ値を持つ場合にテーラー展開!?

角周波数が波数の関数として書けている場合、分散関係があるといいます。
つまり、角周波数$\omega$が波数$k$に依存して$\omega=\omega(k)$と書けている場合、波形のほとんどは

とように書けます！！
本来、位置$x$と時間$t$で書ける関数は、$k$と$w$の二重積分で表現されなければなりませんが、分散関係のおかげで$\omega$が強制的に決まってしまうので、一重積分で良くなるのですね！
うーん、とても嬉しい！！！！

…でも、あまりに一般的過ぎて分かりませんね。
しかし！！
応用で重要なのは、ゆっくり振動する包絡線と早く振動する波の掛け算で表されていることが多いのです。

つまり、これから考える波が持つ波数は、$k= k_0$の周りしか値を持たないのです。ここで、早く振動する波の波数を$k_0$と置きました。

$k_0=5$の場合、これから考えていく$f(k)$はこんな感じの特徴を持っています！ババン！

こんな特徴を持つのならば、$k=k_0$以外ではどうせ$f(k)\to 0$となるので、考えてもしょうがないっていう近似が使えそうですよね。

そ　こ　で　！

波数$k=k_0$周りでいろいろテーラー展開して使用していきます。つまり、$k=k_0$の近傍さえ合っていれば、$k_0$以外でどんな振る舞いになろうとも$f(k)$がゼロになるので大丈夫でしょう！ということです。

早速、分散関係についてテーラー展開を行うと

となります。式を簡単に書くために$\omega'(k_0)\equiv \frac{d\omega}{dk}\bigr|_{k=k_0}$と置きました。

計算を進めると、

より、

のように書けます。つまり、式(4)が言っていることは、
時刻$t$の波形とは、初期状態$t = 0$の波形 $f(x,0)$を$\omega'(k_0)t$だけ平行移動させた波に、時間だけに依存する位相に関する変化分 $e^{−i[\omega(k_0)−k_0 \omega'(k_0)]t}$が掛け合わされたもの、と表すことができる。と言っています。

以上から、初期状態の波が時刻$t$に至るまでの速度は$\omega'(k_0)$で決まることが分かり、これは
群速度
と呼ばれます！！

位相速度はどうやって!?驚きの仮定

さて、ここまでで群速度が求まりましたが、位相速度が求められていません。
そこで、もう少し具体的な波の形を定義しましょう。
波が、

と書けている場合を想定します。ここで、振幅$A(x,t)$は実数値関数で波の包絡線を表しています。また$A(x,t)$の変化は $k_0$ に比べて非常にゆっくり変化しているものとします。
つまり、

今まで考えてきた仮定が全部使える

んですね！

式(5)に代入しますと

となります。つまり、振幅の変化はいつも初期状態の単なる平行移動として表現することができ、位相の変化も簡単に表せられます。

そこで、振幅が平行移動する速度を群速度$v_g$、位相が平行移動していく速度を位相速度 $v_p$ として定義すると、式(6c) よりそれぞれ

と定義してしまうのが良さそうです！

簡単だと思いました？注意です！

注意しなければならないのは、群速度と位相速度の概念は、

振幅と位相を記述する部分がはっきりと区別できるような状況でなければ、定義できない

ということです！振幅部が早く振動していたら、これまでの議論が使えなくなることに注意しましょう。
その場合、群速度と位相速度の概念が変わってしまうので、意味をなさなくなってしまいます。

いろいろなところに群速度に関する議論が！

様々なところにリンクがあります。やはり皆、群速度の概念を理解するのに苦労しているようですね。

「yam＠広島大」物理化学Monographシリーズ　”3. 物体の速度と物質波の速度　－E=hνの本質の理解－”

最後に

いかがでしたか？気に入ったなら是非、記事の下にある「like!」ボタンや「tweet」して皆さんに共有して下さいね!
それでは皆さん、良い物理ライフを!

Note

本記事の大枠の構成はいかがでしたか?　-ニコニコ大百科を参照しています。

Q氏とR博士

負の周波数とは何か

2021年5月30日 sikino コメントする

三角関数で負の周波数はいらない。
負の周波数という言葉が出てくるならば、それは指数関数で表現されていることを暗に意味する。

Q氏とR博士の会話１

Q「負の周波数というものがあると聞きました」
Q「周波数は正ではないのでしょうか」

R「それは”周波数”の意味が違うのだな」
R「詳しくは基底関数だ。正の時と負の時とで、基底関数が違う」

Q「意味が分かりません。もう少し詳しくお願いします」

R「そうだな」
R「別の言い方をすれば、”周波数”と言ったときに、大きく２つの周波数があるのだ」
R「１つはsin, cosの三角関数の周波数、もう１つはexpの指数関数の周波数だ」
R「この説明でどうかな」

Q「周波数は周波数でしょう。関数で意味が変わるのでしょうか。意図がつかめません。もう少しお願いします」

R「素直でよろしい。その通り、周波数の意味は関数によって変わるのだ」
R「具体的に考えよう。周波数を引数に持つ関数を考えてみようか」
R「わたしが、『時間tに依存する周波数f=1を持つ関数yを考える』と言ったら、君はどういう関数yを想像するかな」

Q「三角関数ですね。$y=\sin(2\pi t)$を想像します」

R「そうだ。そこが負の周波数が受け入れられない原因なのだ」
R「私は$y=\exp(i2\pi t)$という指数関数を想定して周波数f=1と伝えた」
R「しかし君は勝手に三角関数の周波数として理解してしまった、理解できてしまった」
R「三角関数の周波数ならば、その周波数は正の周波数しかありえない」（数式による補足）
R「なぜなら、三角関数の周波数が負になったところで、それは正の周波数を持つ波の定数倍で表現できてしまうからだ」
R「つまり、$\sin(-t)=-\sin(t)$だからな。負の周波数は不要だ。」

Q「なるほど」

R「負の周波数という言葉は指数関数を考えていることを暗に示している」
R「指数関数の場合、$\exp(-it)$を$\exp(it)$を使って表現することができないのだ」
R「だから負の周波数が絶対に必要になる」

Q「問題文の説明不足が原因でしたか。理解できました」

Q氏とR博士の会話２

Q「周波数の成分が少なく済む三角関数の方が便利ではないですか」
Q「負を考える必要がないのでしょう」

R「１つ、伝えていなかったな。どちらも同じ程度の便利さだ。」
R「先ほどsinの正の周波数だけで表現できるとしたが、実はこれでは不十分なのだ」
R「ある関数を三角関数で表すならば、必ずsin, cosの両方を用いなければならない」
R「つまりsin関数の周波数ゼロから正の無限大、cos関数の周波数ゼロから正の無限大の二つを用いる関数だ」
R「反対にexp関数を使う表現方法ならば、周波数は負の無限大から正の無限大の一つを用いる関数だ」
R「”周波数”にも”sin関数の周波数”、”cos関数の周波数”、”exp関数の周波数”の３つがあるわけだ」
R「指数関数ならば１つの関数系(exp)で済むから、その意味で簡単だ。しかし、負の周波数が出てきたり、複素関数として扱わなければならなくなる。」
R「三角関数ならば２つの関数系(sin, cos)を用いなければならない。どちらの周波数を考えているかいちいち説明しなければならないが、実部だけで話ができる利点がある。」

Q「なるほど。理解しました」
Q「どちらも場合によるのですね」
Q「フーリエ変換を考えたい場合は指数関数の方がよさそうですね」

R「その通り。まぁ、指数関数による表現がまさにフーリエ変換だがな」
R「ちなみに、基底関数という言葉は、ここでいう三角関数や指数関数のことだ」
R「どういう関数で展開したか？の”どういう関数”が基底関数と呼ばれる」

数式による補足

任意の関数$y(x)$を表現する場合、振動関数であるsin, cosまたはexpを考えます。
任意の関数を適当な関数系で表現するためには、その関数系が完全系であることが求められます。
つまり、三角関数を用いる場合、
$\displaystyle
y(x)=\int_0^{\infty} A(f) \sin(2\pi f x) df+\int_0^{\infty} B(f) \cos(2\pi f x) df
$
として表現できます。反対に、指数関数を用いるならば、
$\displaystyle
y(x)=\int_{-\infty}^{\infty} C(f) e^{i2\pi f x} df
$
と書かなければなりません。
つまり、負の周波数という言葉は指数関数が現れる場合にしかあり得ないのです。

もちろん、積分の変数変換によってsin, もしくはcos関数の周波数を負の無限大からゼロなどに変更することが可能ですが、こんなことをわざわざやる必要性、あり得ますかね。

私が昔見たことがある説明で
$\displaystyle
\begin{align}
y(x)&=\cos(\omega x) \\
&=\frac{1}{2}\cos(\omega x)+\frac{1}{2}\cos(-\omega x)
\end{align}
$
と変形できるから正の周波数と負の周波数になる、というものを見たことがありますが、これは完全に間違いを誘発させるための説明です。

なぜなら、$\cos$の係数に任意性があるためで
$\displaystyle
\begin{align}
y(x)&=\cos(\omega x) \\
&=100\cos(\omega x)-99\cos(-\omega x)
\end{align}
$
としても全く問題ないからです。三角関数の引数が負であるものと正であるものは独立ではないのです。

未分類

割り勘のお金のやりとりを表示する計算ツール

2021年5月22日 sikino コメントする

数人で旅行を行う際に、代表者がまとまって支払うことが多々あると思います。
旅行終了後に割り勘をする場合に計算が大変になりがちです。
その場合にスマホでも使える計算ツールを実装しました。

旅費計算ツール
https://slpr.sakura.ne.jp/qp/calc/splitcheck/split.html

発展版
https://slpr.sakura.ne.jp/qp/calc/splitcheck/split_v1.1.html
https://slpr.sakura.ne.jp/qp/calc/splitcheck/split_v1.2.html

例えば4人(A,B,C,Dさん)の旅行を考えます。グループとして支払った額を以下のように仮定します。
・旅館代としてAさんが4人分を3万円まとめて払ったとします。
・ご飯台としてBさんが4人分を1万円まとめて払ったとします。
・車台としてCさんが4人分を7000円まとめて払ったとします。
・Dさんは支払いをしていません。

旅行が終わった際に旅行全体で割り勘したいと考えたとき、計算とお金のやり取りをどうすればいいか悩みます。
少ないお金のやり取りで済ませたいですよね。そのためのツールがこれです。

人数を入れ、支払い額を入力し実行しますと結果が出ます。

上の例でお金のやり取りは、

・DさんがAさんに11750円支払う
・CさんがAさんに 4750円支払う
・BさんがAさんに 1750円支払う

ことで終わります。もし、100円単位でということであれば、”100円丸め”の計算を実行してください。

htmlとjavascriptで書いています。一度ページにアクセスしておきさえすれば、オフラインでも使えます。
計算実行時にサーバーにアクセスするわけではないので、私は入力した情報を参照することができません。安心してお使いください。
心配な方は、ソースを表示させて計算プログラムを確認してください。

お金のやり取りは最小になると思いますが、証明はしていません。最小に近いやりとり、です。

Blender

光の反射と質感

2021年2月17日 sikino コメントする

最近、Blenderを始めました。

3DCGでは見た目が重要視されます。そのため、物体がどのように光を反射しているのかが、見たときの質感に大きく影響を与えます。
Blenderではレイトレーシングでレンダリングを行うことができるので、どのように反射するのか知っていれば、それに近い質感を再現することができます。

反射の種類についてあまり知らなかったので物体の反射をまとめました。

https://slpr.sakura.ne.jp/qp/wp-content/uploads/2021/02/blender_relrection.pdf(14MB)

上のPDFは重いので、参考のjpg画像を置いておきます。

物理学

再学習用の電磁気学

2021年2月14日 sikino コメントする

電磁気学を可能な限り短くまとめました。電磁気学を一度学んだ方の再学習用を想定しています。
歴史を追うような考えをしていません。なので、証明には「なるからなる」を多く取り入れています。
経緯を知りたい方は各種参考書をご覧ください。

https://slpr.sakura.ne.jp/qp/wp-content/uploads/2021/02/electromagnetism.pdf

electromagnetism

間違いがありましたら、コメントなどいただけると助かります。

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

陰関数の等高線を数値計算で求める

2021年1月5日 sikino コメントする

まとめ

関数

を求める問題を考えます。
$f(x_0, y_0)=0$を満たす$(x_0, y_0)$が与えられている場合、$(x_1\equiv x_0+\Delta x, y_1\equiv y_0+\Delta y)$を求めることを考えていきます。
まず、

の変換を行い、初期値$(x_0, y_0)$を$(x’_0, y’_0)$に変換します。
ここで、$\theta$は定数で

です。$(x’_0, y’_0)$を新たな初期値として微分方程式

を解き、$x_0′$から微小距離$h$だけ発展させた$x_1’\equiv x_0′ + h$とその値$y_1’\equiv y_0′ + \Delta y$を求めます。
その後、逆変換

を経ることで$f(x_1,y_1)=0$を満たす$(x_1,y_1)$の値を得ることができます。

導出

陰関数

を満たす$(x,y)$を求めます。ある点$(x_0, y_0)$において$f(x_0, y_0)=0$を満たすことが分かっている場合、
陰関数の全微分より、$(x_0, y_0)$の近傍で

が満たされています。$f(x_0,y_0)=0$であり、$\Delta x\approx \Delta y \to 0$ならば、$O(\Delta^2)$の項を無視して

が成立する変化のみが許されます。両辺を$\Delta x$で割って、$\Delta x\to 0$の極限をとれば、

を得ます。すなわち、式(9)を何らかの手段、例えば数値計算で解けば等高線を得ることができます。

発散への対処

さて数値計算によって微分方程式(9)を解いて等高線を描く場合、$y$が$x$の一価関数になることはほとんどありません。
そのため、式(9)の右辺の分母がゼロになることがあり、微分が発散するため計算を行うことができなくなることがあります。言い換えれば、発散しない範囲しか計算できなくなるということです。
これを回避するためには様々な方法があると思いますが、例えば変数変換を適宜行いながら積分を実行していくことでこの問題を回避することができます。

発散を回避するためのアイデアは以下の通りです。
$y$を$x$の関数と見たとき、微分が発散するということは、分母がゼロになることを意味しています。
しかし、その点において90°回転させて$y$を$x$の関数としてみれば微分はゼロに等しいため、問題なく計算を進めることができます。
図示すれば以下の通りになります。

左図は座標を変化せずにそのままの関数ですが、右図のようにある点を計算したいときにその点近くの傾きと等しいように回転した座標から見るようにします。こうすることで発散を回避することができます。

これから変数変換について考えていきます。新たな変数$x’, y’$を用いて$x=x(x’, y’),~y=y(x’, y’)$（またはその逆関数）と書けているならば、

と書いても良いです。すなわち式(9)を解く代わりに$x’$と$y’$に関する方程式

を解いても良いのです。具体的に変数変換


を考えます。これは座標を$\theta$回転させた座標系への変換です。もちろん、$\theta=0$のとき、$x=x’, y=y’$が成立し、恒等変換となります。微分方程式は

であり、$(x’, y’)$の空間で解いていきます。

さて、回転角$\theta$は任意ですが、式(15)の右辺が発散しないような$\theta$を選ぶ必要があります。
また、数値計算を行う上では、傾きがほとんど0に近いほうが有利です。そのため、方針としては現在の傾きを計算し、その傾き分だけ逆方向に座標回転させることでこれを実現させます。
すなわち、

分だけ回転させます。現在の傾きと座標軸の傾きを一致させることで傾きがほとんどゼロの座標軸に回転させるのです。

プログラム

以下にプログラムを示します。
必要な情報は$f(x,y)$の$x, y$方向の偏微分$\frac{\partial f}{\partial x}, \frac{\partial f}{\partial y}$の解析解です。関数名”fx”, ”fy”にそれぞれを入力してください。
計算は刻み幅制御の陽的ルンゲクッタ法を用いています。

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

program main
implicit none
integer::i,N
double precision::x,h,tol,len,lenmax
double precision::xp,yp,xt,yt,xt0,yt0
double precision::theta,st,ct

double precision,allocatable::y(:)
double precision,external::fx,fy
double precision,parameter::pi=dacos(-1d0)
external::grk

N=1 ! Number of differential equations
allocate(y(1:N))

xt = -2d0
yt = 1.36465560765604d0 +1d-10
lenmax = 10d0 ! Limitation of calculation length
tol = 1d-8 ! Tolerance of Runge-Kutta

i=0
theta = 0d0
len = 0d0
h=1d-6
do while(len.le.lenmax)
xt0 = xt
yt0 = yt
theta = atan2(fx(xt0,yt0),fy(xt0,yt0))
st = sin(theta)
ct = cos(theta)

! From x to x-prime space
xp = xt0*ct - yt0*st
yp = xt0*st + yt0*ct

! Time propagation in x-prime space
x = xp
y(1) = yp
call drkf45_simple(grk,x,h,N,y,tol,theta)

! From x-prime to x space
xt = x*ct + y(1)*st
yt = - x*st + y(1)*ct

! Length in x-y space
len = len + sqrt((xt-xt0)**2+(yt-yt0)**2)
write(10,'(2e25.10e3)')xt,yt

i=i+1
enddo
write(6,*)"Number of referenced times -->",i

stop
end program main

double precision function fx(x,y)
implicit none
double precision,intent(in)::x,y

fx = 3d0*x**2 - 2d0*y

return
end function fx

double precision function fy(x,y)
implicit none
double precision,intent(in)::x,y

fy = 3d0*y**2 - 2d0*x

return
end function fy

subroutine grk(N,x,y,f,theta)
implicit none
integer,intent(in)::N
double precision,intent(in)::x,y(1:N),theta
double precision,intent(out)::f(1:N)

double precision::tfx,tfy
double precision::xt,yt
double precision::st,ct
double precision,external::fx,fy

st=sin(theta)
ct=cos(theta)
xt = x*ct + y(1)*st
yt = - x*st + y(1)*ct

tfx = fx(xt,yt)
tfy = fy(xt,yt)
f(1) = - ( ct*tfx - st*tfy )/( st*tfx + ct*tfy )

return
end subroutine grk

subroutine drkf45_simple(grk,x,h,N,y,tol,theta)
! if h < hmin, propagate forcibly with warning.
!
implicit none
integer,intent(in)::N
double precision,intent(in)::tol,theta
double precision,intent(inout)::x,h,y(1:N)

integer::i,j,FLAG,key
double precision::R,delta,tx,Sy,err
double precision,allocatable::ty(:),K(:,:),tf(:)
double precision,parameter::hmin=1d-10,hmax=0.5d0
integer,parameter::s=6
double precision::a(1:s,1:s),b1(1:s),b2(1:s),c(1:s),Rc(1:s)
external::grk

c(1:6)=(/0d0, 0.25d0, 0.375d0,&
0.9230769230769230769230769230769230769231d0, 1d0, 0.5d0/)
a(1:6,1:6)=0d0
a(1,1:6)=(/0d0, 0d0, 0d0, 0d0, 0d0, 0d0/)
a(2,1:6)=(/0.25d0, 0d0, 0d0, 0d0, 0d0, 0d0/)
a(3,1:6)=(/0.09375d0, 0.28125d0, 0d0, 0d0, 0d0, 0d0/)
a(4,1:6)=(/0.8793809740555302685480200273099681383705d0, &
-3.277196176604460628129267182521620391443d0, &
3.320892125625853436504324078288575329995d0, 0d0, 0d0, 0d0/)
a(5,1:6)=(/2.032407407407407407407407407407407407407d0,-8d0, &
7.173489278752436647173489278752436647173d0, &
-0.2058966861598440545808966861598440545809d0, 0d0, 0d0/)
a(6,1:6)=(/-0.2962962962962962962962962962962962962963d0,2d0, &
-1.381676413255360623781676413255360623782d0, &
0.4529727095516569200779727095516569200780d0,-0.275d0,0d0/)
b2(1:6)=(/0.1185185185185185185185185185185185185185d0, 0.d0,&
0.5189863547758284600389863547758284600390d0, &
0.5061314903420166578061314903420166578061d0, &
-0.18d0, 0.03636363636363636363636363636363636363636d0/)
b1(1:6)=(/0.1157407407407407407407407407407407407407d0, 0d0,&
0.5489278752436647173489278752436647173489d0, &
0.5353313840155945419103313840155945419103d0, -0.2d0, 0d0/)
Rc(1:6)=(/0.002777777777777777777777777777777777777778d0,0d0, &
-0.02994152046783625730994152046783625730994d0, &
-0.02919989367357788410419989367357788410420d0, 0.02d0, &
0.03636363636363636363636363636363636363636d0/)

key=0
allocate(ty(1:N),tf(1:N),K(1:s,1:N))
ty=0d0; tf=0d0; K=0d0

if(abs(h).ge.hmax)then
h=sign(1d0,h)*hmax
endif

FLAG=1

do while(FLAG.eq.1)
tx=x
do j=1,s
tx=x+c(j)*h
ty(1:N)=y(1:N)
do i=1,j-1
ty(1:N)=ty(1:N)+K(i,1:N)*a(j,i)
enddo
call grk(N,tx,ty,tf,theta)
K(j,1:N)=h*tf(1:N)
enddo

!step 4
R=0d0
do i=1,N
R=R+(Rc(1)*K(1,i)+Rc(3)*K(3,i)+Rc(4)*K(4,i) &
+Rc(5)*K(5,i)+Rc(6)*K(6,i))**2d0
enddo
R=abs(dsqrt(R)/h/dble(N))

Sy=0d0
do i=1,N
Sy=Sy+(y(i)*y(i))
enddo
Sy=dsqrt(Sy)
if(Sy.ge.1d0)then
err=tol*Sy
else
err=tol
endif

!step 5
if(R.le.err.or.key.eq.1)then
x=x+h
y(1:N)=y(1:N)+b1(1)*K(1,1:N)+b1(3)*K(3,1:N) &
+b1(4)*K(4,1:N)+b1(5)*K(5,1:N)
FLAG=0
endif

!step 6
! Avoid zero deviding.
if(R.ge.1d-20)then
delta=(err/(2d0*R))**0.25d0
else
delta=4d0
endif

!step 7
if(delta.le.0.1d0)then
!function changes dramatically.
h=0.1d0*h
elseif(delta.ge.4d0)then
!function changes loosely.
h=4d0*h
else
!function changes moderately.
h=delta*h
endif

!step 8
if(abs(h).ge.hmax)then
h=sign(1d0,h)*hmax
elseif(abs(h).lt.hmin)then
h=sign(1d0,h)*hmin
key=1
endif

enddo

if(key.eq.1)then
write(6,'(A,f10.5,A,f10.5)')&
"Strange point between ",x-h," and ",x
endif

deallocate(ty,tf,K)
return
end subroutine drkf45_simple

具体例

具体的に関数$f(x,y)=x^3-2xy+y^3$の等高線を考えていきましょう。
関数の値がそれぞれ、$0.01,1,2,5$に等しいときの等高線を考えます。
具体的な初期値は、
$\displaystyle
f(-2,1.3657)\approx 0.01,~f(-2,1.4646)\approx 1,~f(-2,1.5568)\approx 2,~f(-2,1.7977)\approx 5
$
とします。
具体的に解いたのが左の図です。

”正解”としたgnuplotの等高線表示の結果と同じ結果が得られました。座標回転がうまく働いているようですね。微分の発散を抑えるだけでなく、数値計算量も抑えられているようです。

刻み幅$h$が正である場合と負である場合で進む方向が変わります。
また計算の途中で分岐がある場合、どうなるかは保証できません。

また、閉じた等高線が複数ある場合、その全てを計算することは本稿のプログラムだけではできません。初期値を求める別のプログラムを作成する必要があります。
更に、等高線がある値に等しい量を計算することはできません。初期値の値に等しい等高線を引きます。

追記

調べ物をしていたら

照井章, 佐々木建昭著『27.代数関数の陰関数描画について』
http://www.kurims.kyoto-u.ac.jp/~kyodo/kokyuroku/contents/pdf/0941-27.pdf
という素晴らしい資料がありました。

この文献では２つの方法が紹介されており、
１．$f(x,y_0)=0$を解いて$x=x_1,x_2,\cdots,x_n$を解く
２．陰関数定理から求める
という方法です。資料では、１．の方法を採用していましたが、本稿の方法は２．の方法でした。
１．はもれなく等高線を見つけることができますので、考えられる結果を得たい場合によい方法です。
一方、２．は全てのゼロ点を見つけることができませんが、計算が早くできます。別の方法などで注目したい等高線が分かっている場合に良い方法でしょう。

効率的に解くためには、場合によって使い分けることが一番です。

弾道計算, 物理学

エアガンの集弾限界

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

理想的なエアガンとBB弾、それに風などが無い理想的な環境があったとします。

バレルとBB弾に直径差があることによって生じるばらつきは、

　30m先 | 直径3cm
　50m先 | 直径10cm

となります。この直径以下にばらつきを抑えることは構造的に不可能です（60cmのバレルで0.9Jで射出する場合）。
つまり、このばらつきよりも大きくばらつく場合、バレルとBB弾の直径差以外にばらつく原因が存在します。

バレルとBB弾に直径差があることによって生じるこのばらつきを軽減するには、
　・バレルの直径を小さくする
　・バレルを長くする
　・バレル内部の素材とBB弾の反発が起こりにくいものにする
にすることで軽減されますが、他の要因が場合はこの限りではありません。

ここでの理想的とは、

・球の直径のばらつきは無い
・いつも同じ初速で撃ち出せる
・無風

とします。この場合でも同じ点に集弾することはありません。
なぜなら、バレルの直径とBB弾の直径が異なる為です。この差によってどの位集弾性が悪くなるのか、見積もりましょう。

これは構造的な問題であり、ばらつきの原因の中で取り除く事ができない１つの原因です。

ばらつく原因として、以下の三つが考えられますが、まず本稿では理想状態のばらつきを考えます。

ばらつきの種類	理想状態のばらつき（本稿）	製品誤差によるばらつき	セッティングによるばらつき
BB弾重さ	–	あり	なし
BB弾の大きさ	なし	あり	なし
回転のばらつき	なし	なし	あり
手振れ	なし	なし	あり

ばらつきが生じる原因

何故ばらつきが生じるかを考えましょう。
全くの理想であれば、銃口から同じ初速、角度、回転量で射出されたBB弾にばらつきが生じることはありません。
しかし、現実にはバレルの大きさとBB弾の大きさに差があります。これによって銃口から射出する時に進行方向と垂直な平面に速度を持つことになります。
銃口から出てきたBB弾の速度$\mathbf{v}$を以下のように表現します。

ここで、$v_{z}, \mathbf{e}_{z}$はバレルの方向に沿う速度、単位ベクトルで、
$v_{x,y}, \mathbf{e}_{x,y}$はバレルの方向に垂直な面への速度、単位ベクトルです。
すなわち、壁面に垂直な方向に対する速度ベクトルの大きさ$v_\perp$は、$v_\perp=\sqrt{v_x^2+v_y^2}$と表されます。

ばらつきが無くまっすぐ飛ぶのであれば$v_{\perp}=0$であり、そうでなければ$v_{\perp}\ne 0$です。

BB弾のばらつきが生じる原因は、射出方向に垂直な方向（横方向）に有限の速度が生じている、と仮定します。
この横方向の速度が生じる原因の一つは、ピストンでBB弾を空気で押す際に圧力が一定ではないとか、回転を掛けるためのゴムで生じる、などいろいろ考えられます。
ここでは、パッキンを通過して、もうこれ以上横揺れを増やすような原因が生じないのだ、と仮定します。壁の反射によって変化はあるものの、速度は減衰する一方であるとします。

それでは、定式化をしていきましょう。

定式化

横方向の速度はそれほど大きくないだろう、と予想するので空気抵抗は考えません。この条件のもと考えていきます。

バレルの直径$d_s$とBB弾の直径$d_B$の差があることによって、バレルに垂直な方向にBB弾が自由に進める距離$l$は

と書けます。$n$回目の衝突から$n+1$回目の衝突までにかかる時間$t_n$は、その間のBB弾の垂直方向の速度$v_n$に依存して、

と書けます。よって、$N$回衝突するのにかかる時間$t$はそれぞれの衝突までに掛かる時間を足し合わせればよいので、

です。バレルの内壁との衝突が起こることによってBB弾の速度の変化するとします。するとBB弾とバレルとの反発係数$r$を用いれば、

と表現されます。初速度を$v_0$と書いてこれを代入すれば、

と書くことが出来ます。ただし、$r\ne 1$です。$r=1$の場合は$\displaystyle
t=\frac{l}{v_0}(N+1)
$
です。$r\ne 1$の場合に$N$について変形すれば、

と書くことが出来ます。
初速$v_0$、反発係数$r$, 銃口に到達するまでに掛かる時間$t$が分かれば、銃口から出た時のバレルに垂直な方向の速度$v_\perp$は式(5b)より、

と書けます。具体的に妥当な値を入れて射出時の、進行方向に垂直な方向の速度を計算してみましょう。

具体的な衝突回数の見積もり

60cmのバレル長を持ち、0.20gで0.9J程度の球が射出される場合、セットされた位置から射出までにかかる時間は約$t=0.015\mathrm{[s]}$と分かっていますのでこれを利用します（バレル内部の計算より）。
バレルBB弾の直径差は

とします。$r$はどうにかして求めることにして、$v_0$はとりあえず変数としましょう。

反発係数の見積もり

$r$を見積もります。
反発係数の見積もりは、実験して大雑把な値を見積もります。
$r$を見積もるために、家にある固い材質のものとBB弾とを衝突させて、高さを計測しました。以下の実測結果を得ました。

BB弾-フローリング 30cmから落として15cm上がる
BB弾-アクリル板 30cmから落として10cm上がる
BB弾-アルミ 10cmから落として6cmまで上がる

高さ$h_0$から静かに落として、反発した後に極大値をとるときの高さ$h_1$が分かっているとき、反発係数$r$は

と求められることが分かっているので、反発係数は

BB弾-フローリング　　$r\approx 0.71$
BB弾-アクリル板　　$r\approx 0.57$
BB弾-アルミ　　$r\approx 0.77$

となります。バレル内部はどちらかといえば金属に近いと思いますので、$r\approx 0.75$と仮定して計算を進めていきます。

横軸に進行方向に垂直な方向の初速度…つまり、ホップを掛けるためのゴムなどで、進行方向ではない方向の速度のことを意味します…、縦軸に射出時の垂直な方向の速度をプロットしました。横方向の速度が仮に50m/sになっていても、10m/sになっていようともあまり変わらないことが分かります。なので、10m/sと仮定しましょう。
ばらつきの上限を与えるにはよい指標です（※1）。

仮に$r=0.75, v_0=10\mathrm{[m/s]}$とすると、$N\approx 20.6$と計算されます。つまり、20.6回バレルに衝突してから飛び出していくわけです。実際には衝突回数は整数しかありえないので、小数点以下を切り捨てて20回衝突が起こって飛び出していきます。
ただし計算を行う上では衝突回数を切り捨てると不連続性が発生しますので、小数点の衝突回数を認めることにします。

この場合、銃口から飛び出す際の、進行方向に垂直な方向の速度$v_\perp$は

を得ます。この垂直方向の速度は上下左右に振れる可能性がありますが、最も触れる場合は横方向でしょう。

30, 50m先の広がり具合を考えます。
着弾までの時間$t_{30m}, t_{50m}$はそれぞれ$0.5, 2 [s]$分かっているので（弾道計算（BB弾）の結果）、着弾時の直径$d_{30m}, d_{50m}$は

と求められます。ここで垂直方向の速度が非常に遅いので、空気抵抗は無視して考えています。BB弾の重さは関係ありません。軽いものほど空気抵抗を受けて減衰するので、ぶれは小さくなります。

空気抵抗を考慮すると、これよりは小さくなるということです。

実測定との比較

さて、30mチャレンジというものがあります30mチャレンジ公式ランキング　2016年度。
この記録によると、30m先で直径15cm程度の範囲に収まるようです。
すなわち、バレル-BB弾の直径差による集弾性の悪化よりも大きな影響を及ぼす要因がある、ということです。

恐らくは回転量が一定ではないとか、BB弾の重さにばらつきがあったり、銃身の振動があったりという要因のほうが大きいということでしょう。

逆に言えば、スタンダードな大きさのバレルとBB弾を用いる限り、30m先で3cm以内に収めることは不可能、ということです。

注釈

※1
$v_0\to \infty$で射出時の速度$v_\perp$は発散しますが対数の発散です。非常にゆっくり発散するため、$v_0=10$であろうが$v_0=50$であろうがほとんど変化しません。どこかでカットオフ（これ以上はとらない上限）となるような値を取れればよいかもしれません。
豆知識ですが、対数の発散をほとんど無視する目的としてカットオフを設けるやり方は、量子力学の繰り込み論が有名ですね。だからと言って、エアガンの計算で量子力学が現れる！などは言わないでください。全く関係ありません。