sikino のすべての投稿

Fortran90で構造体を使う

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という新しい”型”(複素数型みたいに、実部と虚部という2つの独立な変数を1つの複素数型変数が持つ、みたいな変数の型)を定義しました。
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<br />
                    2   1.0000000000000000        2.0000000000000000       0.20000000000000001<br />
                    3   1.0000000000000000        3.0000000000000000       0.30000000000000004<br />
                    4   1.0000000000000000        4.0000000000000000       0.40000000000000002<br />
                    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<br />
                    2   2.0000000000000000        2.0000000000000000       0.20000000000000001<br />
                    3   2.0000000000000000        3.0000000000000000       0.30000000000000004<br />
$

補足


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

クーロンポテンシャルによる散乱(量子力学)

詳細


本稿では概要だけを説明し、導出過程や詳細は以下のpdfで説明します。
https://slpr.sakura.ne.jp/qp/wp-content/uploads/2021/08/scattering_coulomb.pdf

概要


電荷\(q_A\)を持つ原子核に、遠方から\(q_B\)を持つ電子が衝突する過程を考えます。
例えば、水素原子の原子核である陽子1個(\(q_B=+e\), eは電気素量)に対し、遠方から電子(\(q_A=-e\))が衝突する過程です。
図で表せばこのような感じです。

相対座標に対する時間依存しないシュレーディンガー方程式は、

と書くことができます。ここで、\(\hbar\)はプランク定数、\(\varepsilon_0\)は真空の誘電率、\(\mu\)はAとBの換算質量を表し, \(E\)は衝突のエネルギーを表します。
簡単にするために係数をまとめて

と書きます。ここで、

という量を定義しました。
これから、波数\(k\)で定義される平面波が入射したらどのように応答するか?を知りたいので、式(1)を境界条件

の下で解いていきます。

導出した結果は、

となります。この展開は\(|r-z|\to 0\)となる\(z\)軸上もしくはz軸近傍では使えないことを注記しておきます。元々の境界条件(4)は、\(z\to -\infty\)だけで決まっており\(x,y\)には条件はありませんが、結論として得られるものは\(z\)軸からは離れた領域でなければならない、とより制限された範囲でのみ導ける、となっています。

ここで

であり、\(\gamma\)はSommerfeld parameterと呼ばれる量であり、\(\nu=\frac{\hbar k}{\mu}\)は2粒子間の相対速度を表しています。また、

はクーロン散乱振幅と呼ばれる量であり、波数\(k\)を持って入射した平面波が、入射方向に対して角度\(\theta\)方向にどれだけ球面波が歪むか?を表す量です。

導出過程や詳細は以下のpdfで説明しています。
https://slpr.sakura.ne.jp/qp/wp-content/uploads/2021/08/scattering_coulomb.pdf

本稿では、散乱状態の導出にとどめております。
解釈、流束などについては、次回説明します。

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

本稿の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)を見てみると、1つのヘヴィサイド関数\(\theta(t-t’)\)と1つの矩形関数\(\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}\)である、という事実と一致します。
つまり、衝撃が加わった時刻以外では、系そのものの性質しか存在しませんので、衝撃が波及していく速度が一致することは何ら不思議ではありません。むしろ、一致しなければなりません。

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

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

定義


ヘヴィサイド関数

符号関数

矩形関数

デルタ関数

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

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

フーリエ関数

性質、関係式


積分表示

定積分が関係する場合


デルタ関数

ヘヴィサイド関数


sinc関数に関係する式

畳み込み積分


デジタル信号処理と偏微分方程式の関係

問題提起


デジタル信号で良く以下の図のように入力と出力が説明されます。

ここで疑問が生まれました。それは

  • 時間だけの関数であるならば、通過前後の位置はどこに記述されている?
  • 出力が、入力とシステムの畳み込み積分である根拠は?

です。通過位置については、どこにも位置が現れていないのです。また、畳み込みは確かに正しそうですが、どこから畳み込みが現れたのか根拠が分かりません。
これの解決を目指します。

本稿の結論は、散乱問題の考え方と同じで良い、です。つまり、

  • 位置と時間は、偏微分方程式の分散関係によって結び付いているので、時間又は位置について分かればもう片方は自動的に来まる。
  • 初期状態から終状態が生まれる際にグリーン関数を用いて書くことができるので、それが畳み込みの形となり、出力となる。

ということです。

序論


デジタル信号処理に現れる入力・出力を表す関数は何か?


波は偏微分方程式として記述されることが多く、基本的に波は位置と時間の関数として書かれます。
しかし、デジタル信号処理の分野では位置に関する表記を消去して時間のみの関数(時間信号)について解くだけで良し、とされます。
波はそもそも偏微分方程式ですので、何らかの仮定をおいて位置は考えなくて良し、としているはずですが、その仮定とは何なのでしょうか。
調べても『当たり前』とされている部分のようで、数式で示されている資料を見つけることができませんでした。

この疑問を詳しく説明します。
デジタル信号処理では、波と相互作用する何かを”システム”という言葉で表現し、その”システム”を特徴づける”インパルス応答”やそのフーリエ変換である”周波数特性”として扱います。
そして、”入力”と”システム”の相互作用を畳み込み積分した結果を”出力”として表します。
ここで私が疑問に思ったのは、

  • 「インパルス応答を用いて畳み込み、は妥当そうだが、それはなぜokなのか?」
  • 「現実には”システム”を通過する前後の位置という情報があるはずなのに、それはどこに行った?」

という2点です。

畳み込み積分が出てくる数学上の話はグリーン関数が想像されますが、やはり検索しても見当たらないので自分で導出を試みました。
いろいろ悩みましたが、以降の導出は私の導き、納得できた結論です。考え方は散乱問題と同じです。
つまり、畳み込みが出てくるということは非斉次方程式が出てきて、グリーン関数を用いて解が表される方程式があるということです。

導出


関係性の導出


波動方程式

を解くことを考えます。ここで、\(\hat{D}\)は位置\(x\)と時間\(t\)に関する微分演算子です。具体的に例えば
\begin{align}
\hat{D}=\frac{\partial^2}{\partial t^2}-a\frac{\partial^2}{\partial x^2}
\end{align}

みたいな感じです。相互作用を表現する\(s(x,t)\)の持つ性質として
\begin{align}\displaystyle
s(x,t)\Bigr|_{x\gt |a|}=0
\end{align}
を仮定します。

この仮定の下で、初期状態は\(x\lt -|a|\)の領域のみで波形が与えられると考えます。つまり、初期状態\(f_0(x,t)\)は

を満たす状態として与えられます。この条件下で、式(1)の解を

の形で探します。\(g(x,t)\)は未知の、これから求めたい関数です。
式(1)に代入すると\(g(x,t)\)について

を満たす関数として書けます。式(4)の解\(g(x,t)\)は

とグリーン関数を用いて書くことができます。ここで、\(g_0(x,t)\)は
\begin{align}
\Bigl[\hat{D}(x,t)+s(x,t)\Bigr]g_0(x,t)=0
\end{align}
を満たす解であり、\(G(x,t;x’,t’)\)は

を満たす解です。補足しますが、右辺を見ると\(x-x’, t-t’\)の形でしか存在しません。これはつまり、左辺のグリーン関数も\(x-x’, t-t’\)の関数になっていることを意味します。つまり、
\begin{align}
G(x,t;x’,t’)=G(x-x’,t-t’)
\end{align}
という形になっており、グリーン関数は4変数関数ではなく、2変数関数であることが分かります。どちらの表記でも構わないので適宜用いることにします。
もし初期状態が存在しない、つまり\(f_0(x,t)=0\)ならば\(g(x,t)=0\)が想定されるため、

でなければなりません。よって、

と書けます。ここで\(H(x,t;x’,t’)\)は

を満たす関数で、俗にいうインパルス応答を記述します。ここでもGと同様に
\begin{align}
H(x,t;x’,t’)=H(x-x’,t-t’)
\end{align}
と書くことができます。以上から、式(1)の解は

と書けます。ここで解の意味を考えます。右辺第一項は初期状態(入力)であり、右辺第二項は、初期状態が\(s(x,t)\)と相互作用した結果生じる相互作用項(出力)を意味します。

仮に初期状態が

とかけていた場合を考えます。式(10)の右辺第二項を計算すると

とまとめられます。ここで、

と定義しました。以上から、解は

となります。良く知られているように、波数空間(もしくは、初期状態を角周波数\(\omega\)の変数とするならば角周波数空間)で相互作用の結果は、相互作用(インパルス応答)のフーリエ変換 \(H(k,\omega(k))\)と初期状態のフーリエ変換との積となります。

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

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

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

  1. そもそも「位相速度」「群速度」って何?
  2. 学会騒然。ある波数近傍だけ値を持つ場合にテーラー展開!?
  3. 位相速度はどうやって!?驚きの仮定
  4. 簡単だと思いました?注意です!
  5. いろいろなところに群速度に関する議論が!
  6. 最後に
  7. 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) よりそれぞれ

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

簡単だと思いました?注意です!


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

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

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

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


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

群速度の謎 -平林 浩一, (C) 2006

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

最後に


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

Note


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

負の周波数とは何か

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

Q氏とR博士の会話1


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

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

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

R「そうだな」
R「別の言い方をすれば、”周波数”と言ったときに、大きく2つの周波数があるのだ」
R「1つはsin, cosの三角関数の周波数、もう1つは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博士の会話2


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

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

としても全く問題ないからです。三角関数の引数が負であるものと正であるものは独立ではないのです。

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

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

旅費計算ツール
https://slpr.sakura.ne.jp/qp/calc/splitcheck/split.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で書いています。一度ページにアクセスしておきさえすれば、オフラインでも使えます。
計算実行時にサーバーにアクセスするわけではないので、私は入力した情報を参照することができません。安心してお使いください。
心配な方は、ソースを表示させて計算プログラムを確認してください。

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

ガウス求積法の導出

ガウス求積法は、関数の数値積分を行う方法です。
N個の点で関数を参照すると、2N-1次以下の多項式で表される関数を厳密に数値積分することができる方法です。

ガウス求積法の凄さを少し述べます。
通常、3個の任意の点における値があれば、その3点を結ぶ2次関数を作ることができるので、うまい選び方をすればN個の点でN-1次の関数を厳密に表現できそうであり、それが厳密に積分できてもそこまで不思議ではありません。

ですが、ガウス求積法では分点の位置にも制限を加えることで積分の次数を倍にすることができます。
分点の位置を制限するにもかかわらず、高々二倍になるだけか…とも思われますが、これは同じ精度に達するまでに関数の評価点数を半分に減らせることと同じです。
評価点数を少なくできるというこの性質は、多次元の計算において顕著に現れます。なぜなら、1次元で半分になるので2次元では4分の1になります。

多少無理してもガウス求積法を使う価値があります。
科学技術計算ではガウス求積法を使うことをお勧めします。ガウス求積法を推奨します
ガウス求積法を使いましょう(圧)。

本稿のPDF版はこちらのリンクからどうぞ。

https://slpr.sakura.ne.jp/qp/wp-content/uploads/2021/03/Gauss_legendre_quadrature_derivation_20210904.pdf

導出手順


導出は以下の順番で行います。

  1. 点数と次数の関係
  2. 具体的な位置と重みの計算

『点数と次数の関係』では、M個の離散的な点を用いてN次多項式を厳密に積分できるとき、NとMの関係を導出します。
『具体的な位置と重みの計算』では、具体的に、M点の位置とその重みを導出します。

ガウス求積法を利用する方は、1. が分かればよいと思います。
ガウス求積法を研究する方は、2. まで分かればよいと思います。

また、本稿ではガウス求積法の中で最も基本的な
ガウス=ルジャンドル求積法
について記述します。他の方法も、ガウス=ルジャンドル求積法の導出を理解できれば、おのずと分かるでしょう。

先に結論をまとめます。
\(2N-1\)次以下の多項式\(f(x)\)の積分を、N個の点だけを用いて厳密に計算できます。
つまり、

を厳密に計算する点の位置\(x_i\)と\(w_i\)が存在し、N個の\(x_i\)とN個の\(\omega_i\)は

から求めることができます。ここで、\(P_n(x)\)はn次ルジャンドル多項式です。

点数と次数の関係


N次多項式を有限区間\([-1,1]\)で、M個の点を用いて数値的に厳密に積分することを考えます。つまり、関数\(f(x)\)がN個の既知の係数\(a_n\)を用いて

と書けている場合に、積分

を考えます。この積分を離散的な点だけを用いて計算します。つまり、適当なM個の点を用いて

となるような右辺の\(x_i, \omega_i\)の組み合わせを求めることが目標です。

式(3)をできるだけ計算していきましょう。

まず、式(3)の左辺を積分すると、

となります。

続いて、式(3)の右辺を計算すると、

と書けます。ここで、\(\omega_i\)は重みであり、適当な係数です。もしも任意の\(a_n\)について、式(3) の右辺と左辺が厳密に一致するならば、\(n=0, 1,\cdots, N\)について、

が成立します。この式は合計で\(N+1\)個の方程式の数があります。今、右辺の未知の係数は\(\omega_i, x_i, (i=1,\cdots, M)\)なので、右辺の自由度は\(2M\)個あります。
つまり、両辺を一致させるような係数の組み合わせがあるならば、

である必要があります。よって、もっとも小さくなるような\(M\)であるならば、上式の等号が成立します。できるだけ少ない点数で計算したいので、等しい場合を採用して

となります。つまり、M個の点を利用すれば\(2M-1\)次の多項式を厳密に計算できる点の位置\(x_i\)と\(w_i\)が存在します。

具体的な位置と重みの計算


具体的な分点の位置

積分

を求めるために、離散的なN個の点\(x=x_i\)における\(f(x)\)の値だけを利用します。つまり、

の点は既知とします。ここで、\(f_i\equiv f(x_i)\)と定義しました(この時点でまだ\(x_i\)は決まっていません)。
新たな関数\(f(x)\)を定義します。\(f(x)\)は、式(12)のN個の点を全て通る\(N-1\)次多項式であり、一意に決まります。具体的には、ラグランジュ多項式を用いて

と書けます。ここで、更に\(f(x)\)と\(f(x)\)の差の関数\(g(x)\)

を定義します。\(x\)の次数を確認しておくと、\(g(x), f(x)\)は\(2N-1\)次、\(f(x)\)は\(N-1\)次多項式です。
ここで、\(f(x_i)=F(x_i),~(\mbox{for all } i)\)なので、\(g(x)\)は

の性質を持ちます。積分を考えると

であり、\(f(x)\)を離散的に表すと、

です。今、\(f(x)\)は\(2N-1\)次であり、\(f(x)\)はN個の点で表現されているので、前節の議論から式(17)をゼロにするような\((x_i, \omega_i)\)の組み合わせが存在することが分かります。
今, \(g(x)\)の満たす性質として式(15)が分かっていますので、N個の\(x_i\)について、

と書き換えることができます。
ここで、\((x-x_1)(x-x_2)\cdots(x-x_N)\)はN次多項式で、\(G(x)\)は\(N-1\)次多項式です。つまり、

を満たすN個の\(x=x_i\)を探すことができれば嬉しいわけです。ここで、\(G(x)\)は\(N-1\)次多項式なので、\(N-1\)次のlegendre多項式を用いた基底関数で展開可能です。なので、

と書くことができます。また、

はN次なので、N次多項式の関数で書くことができます。

つまりN次多項式としてN次legendre多項式\(P_n(x)\)を選ぶことにします。
単にN次であるだけだと、N次以下の全てのlegendre多項式の和で書けなければなりませんが、点\(x_1, x_2, \cdots, x_N\)を\(P_N(x)=0\)を満たすN個の点を作為的に選ぶことにしましょう。

このように選ぶと、N次以下の和ではなく、N次のlegendre多項式だけで書くことができて、

とする定数\(c\)を見つけることができます。

これからの議論において定数\(c\)を具体的に決める必要はありませんが明確に決めたい場合

の関係から決まります。

式(19)に式(20), 式(22)を代入すると

となります。これはつまり、N次legendre多項式のゼロ点を求積法におけるゼロ点として選ぶと

を満たすように\(g(x)=0\)にできるということです。

ここまでの議論で、分点の位置\(x_i\)が求まりました。
残るは、まだ決定していない重み\(\omega_i\)を決めることです。

具体的な重み

上記までの議論から、\(2N-1\)次の多項式で記述される\(f(x)\)の積分は、N個の点を用いて

を厳密に積分する分点の位置と重み\(x_i, \omega_i\)の組があることを示し、分点の位置\(x_i\)は、N次legendre多項式のゼロ点、すなわち

を満たすN個のゼロ点であることが分かりました。次は重みを決めます。
式(16)より\(2N-1\)次の\(f(x)\)をN次のlegendre多項式のゼロ点を通る多項式の積分に置き換えできます。

つまり、式(13)を用いて

のように置き換えできます。つまり、重みは

ですので、これを求めることが問題となります。

積和は書き換えがやりにくいですので、意図的に消していきましょう。
式(22)から、

分かっています。意図的にk番目の式についてまとめると、

となります。また、式(35)の両辺を微分すると

となります。特に、\(x=x_k\)ならば、

式(36)と式(40)を組み合わせると、

なので、重み\(\omega_k\)は

となります。この積分を計算するために多項式の性質を用いて計算を進めましょう。
右辺の形は、クリストッフェル=ダルブーの定理から導けそうな形であることを当たりをつけて計算します(私のような凡人には思い付きません。天下り的に使います。導けるからそれで良いです)。

直交多項式(Legendre多項式やLaguerre多項式など)にはクリストッフェル=ダルブーの定理 (Christoffel–Darboux theorem) があります。それは、n次の直交多項式を\(P_n(x)\)と書くとき、\(P_n(x)\)には

という性質であります。ここで、\(k_n\)は、直交多項式\(P_n(x)\)の\(x^n\)の係数で、

で表されるときの係数です。この係数を具体的に求めてみましょう。そのために、
legendre多項式の性質の一つである三項漸化式

を用います。xの次数で両辺を比較すると

なので、次数を比較すれば最高次の係数について

の関係があります。よって、

を得ます。また、\(h_N\)は直交性から導かれる係数で、

を満たします。legendre多項式の場合、\(w(x)=1\)であり、\(h_n=\frac{2}{2n+1}\)です。

重みの計算に使いたいため、クリストッフェル=ダルブーの定理に対し、\(N\to N-1\)にして、legendre多項式の場合に適用、さらに\(y=x_k\)の場合を考えて、

より、

を得ます。ここで\(P_{N}(x_k)=0\)を利用しました。

式(43)に式(56)を代入して、

となります。\(P_n(x)\)の微分は次数が一つ下がるので、\(P_{n-1}(x)\)と同じ次数になります。つまり、両者に何か関係がありそうですので頑張って見つけてみます。

微分に関する三項漸化式

を用います。特に\(n=N,~x=x_k\)の場合、

の関係があります。代入すると重みは

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

まとめ


\(2N-1\)次以下の多項式\(f(x)\)の積分を、N個の点だけを用いて厳密に計算できます。
つまり、

を厳密に計算する点の位置\(x_i\)と\(w_i\)が存在し、\(x_i, \omega_i\)は

から求めることができます。ここで、\(P_n(x)\)はn次ルジャンドル多項式です。

参考文献など


恐らく、ガウスが発表した論文がこれではないかと思います。ドイツ語のwikipedia(https://de.wikipedia.org/wiki/Gau%C3%9F-Quadratur)にリンクがありました。
Methodus nova integralium valores per approximationem inveniendi. In: Comm. Soc. Sci. Göttingen Math. Band 3, 1815, S. 29–76

Abramowitz and Stegun, “Handbook of Mathematical Functions.”
https://www.math.ubc.ca/~cbm/aands/page_887.htm

L. Bos et al, “An Orthogonality Property of the Legendre Polynomials”, Constr Approx (2016)
https://math.indiana.edu/promotion/files/legendrepaperinprintCA.pdf

LECTURE 16GAUSS QUADRATURE
https://coast.nd.edu/jjwteach/www/www/30125/pdfnotes/lecture16_19v17.pdf

Identities and properties for associated Legendre functions
https://www.mat.univie.ac.at/~westra/associatedlegendrefunctions.pdf

Christoffel–Darboux formula -wikipedia
https://en.wikipedia.org/wiki/Christoffel%E2%80%93Darboux_formula

光の反射と質感

最近、Blenderを始めました。

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

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

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

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