ラベル 物理 の投稿を表示しています。 すべての投稿を表示
ラベル 物理 の投稿を表示しています。 すべての投稿を表示

2023年10月27日金曜日

鏡像法(3)

鏡像法(2)からの続き

前回は,導体球が接地されている状況を考えた。$z$軸上の電荷$+q$の鏡像電荷$-q'$に相当する電荷は導体球表面に分布しており,これを鏡像電荷が代表して表わしていることになる。

次に,接地されておらず,電荷を持たない導体球を考える。導体球の中心を通る$z$軸上の点A$(0,0,d)$に電荷$+q$を置くと,静電誘導によって導体球表面には偏った電荷分布が生ずるとともに,球表面の電位は一定になる。ただし,この電荷分布を球表面について寄せ集めるとゼロになっている。

この状況を表現するためには,前回のモデルに加えて,導体球の中心に$+q'$相当の電荷を置けば良い。これによって,導体球の合計電荷はゼロになると同時に,導体球表面での電位一定の条件が満足されることになる。実際には,これらの電荷は導体球表面に分布しているのである。

図1:接地しない導体球と鏡像電荷

原点を中心とする半径$R$の接地していない導体球に対して,電荷$+q$とこれによって生ずる鏡像電荷$-q'$,$+q'$がつくる導体球外の電位の式は次のようになる。

$V(\bm{r}) = \dfrac{1}{4\pi\varepsilon_0}\Bigl\{ \dfrac{+q}{\sqrt{r^2+d^2-2 r d \cos\theta}} + \dfrac{-q'}{\sqrt{r^2+d'^2-2 r d' \cos\theta}} + \dfrac{+q'}{r} \Bigr\}$

接地しない導体球表面の誘導電荷密度は,$\sigma(\theta) = -\varepsilon_0 \dfrac{\partial V(r)}{\partial r}\Biggr |_{r=R}$で与えられる。
したがって,$\sigma(\theta) = \dfrac{1}{4\pi} \Bigl\{\dfrac{+q (r-d \cos\theta)}{(r^2+d^2-2 r d \cos\theta)^{3/2}} + \dfrac{-q' (r - d' \cos\theta )}{(r^2+d'^2-2 r d' \cos\theta )^{3/2}} + \dfrac{+q'}{r^2} \Bigr\}\Biggr |_{r=R}$
また,これによる導体球面上の全電荷は,$\displaystyle \int_0^{2\pi} \int_0^\pi \sigma(\theta) R^2 \sin \theta d \theta  d\phi$,すなわち$\ t = \cos\theta\ $とおけば,$\displaystyle 2\pi R^2 \int_{-1}^{1} \sigma(t) \bm{dt} $で与えられる。各項を$\ q_1,\ q_2,\ q_3 \ $とすると,$d>R>d'$なので,

$\displaystyle q_1 = \dfrac{q R^2}{2} \int_{-1}^{1} \dfrac{R - d\ t}{(R^2+d^2 - 2 R d\ t)^{3/2}} \bm{dt}$
$\displaystyle =  \dfrac{q R^2}{2} \Bigl [ \dfrac{R - d\ t }{Rd (R^2+d^2-2R d\ t)^{1/2}} \Bigr ]_{-1}^{1}- \dfrac{q R^2}{2} \int_{-1}^1 \dfrac{ -d}{Rd (R^2+d^2-2R d\ t)^{1/2}} \bm{dt}$
$\displaystyle =  \dfrac{q R^2}{2} \Bigl [ \dfrac{1}{Rd}\Bigl\{ \dfrac{R-d}{d-R}-\dfrac{R+d}{d+R} \Bigr\}- \dfrac{q R^2}{2} \Bigl [ \dfrac{ d}{(Rd)^2} (R^2+d^2-2R d\ t)^{1/2} \Bigr ]_{-1}^{1}$
$\displaystyle =  -\dfrac{q R}{d} - \dfrac{q}{2d} \Bigl\{ (d-R)-(d+R) \Bigr \} = -\dfrac{q R}{d} +\dfrac{q R}{d} = 0$

$\displaystyle q_2 = \dfrac{-q' R^2}{2} \int_{-1}^{1} \dfrac{R - d'\ t}{(R^2+d'^2 - 2 R d'\ t)^{3/2}} \bm{dt}$
$\displaystyle =  \dfrac{-q' R^2}{2} \Bigl [ \dfrac{R - d'\ t }{Rd' (R^2+d'^2-2R d'\ t)^{1/2}} \Bigr ]_{-1}^{1} + \dfrac{q' R^2}{2} \int_{-1}^1 \dfrac{ -d'}{Rd (R^2+d'^2-2R d'\ t)^{1/2}} \bm{dt}$
$\displaystyle =  \dfrac{-q' R^2}{2} \Bigl [ \dfrac{1}{Rd'}\Bigl\{ \dfrac{R-d'}{R-d'}-\dfrac{R+d'}{R+d'} \Bigr\} + \dfrac{q' R^2}{2} \Bigl [ \dfrac{ d'}{(Rd')^2} (R^2+d'^2-2R d'\ t)^{1/2} \Bigr ]_{-1}^{1}$
$\displaystyle =  \dfrac{q'}{2d'} \Bigl\{ (R-d')-(R+d') \Bigr \} = -q'$

$\displaystyle q_3=\dfrac{q' R^2}{2}\int_{-1}^{1} \dfrac{1}{R^2} \bm{dt}= q'$


図2:z軸からの角度 t=cos θの関数としての球表面電荷密度









2023年10月26日木曜日

鏡像法(2)

鏡像法(1)からの続き

鏡像法で電場を求める際の典型的な例は,点電荷と電荷を持たない導体球の系である。
これは次のような条件から幾何学的に求めることができる。

(0) 導体球の表面は等電位面である。
(1) 2次元(3次元)のユークリッド空間で2点からの距離が一定の比$k(\neq1)$となる点の集合は円(球面)である。
(2) 点電荷のつくる電位は測定点までの距離に反比例するので,2つの逆符号で大きさが異なる点電荷のつくる電位の和がゼロになる点は,距離が一定(電荷の絶対値)の比になる点の集合となる。


図:接地した導体球と点電荷が作る電場

原点$(0,0,0)$に中心Oがある半径$R$の接地された導体球面上の点をP$(x,y,z)$とする。球面上の電位はゼロで,$x^2+y^2+z^2=R^2$である。電荷$q$が点A$(0,0,d)\ (d>R)$にあり,電荷$-q'$が点B$(0,0,d')\ \ (d'<R)$に置かれている。

球面上の任意のP点の電位が0になる条件式は,$\displaystyle V(\bm{r})=\dfrac{q}{4\pi\varepsilon_0}\dfrac{1}{\sqrt{x^2+y^2+(z-d)^2}} + \dfrac{-q'}{4\pi\varepsilon_0}\dfrac{1}{\sqrt{x^2+y^2+(z-d')^2}} = 0 $である。したがって,$ R^2+d'^2-2zd' = \bigl(\dfrac{q'}{q}\bigr)^2 (R^2+d^2-2zd) $で,これが$z$によらずに成立するので次の2つの条件式が得られる。

$d' = \bigl(\dfrac{q'}{q}\bigr)^2 d \ $,$R^2+d'^2 =  \bigl(\dfrac{q'}{q}\bigr)^2 (R^2+d^2) \ $,$\therefore R^2 = d d'$,$|\dfrac{q'}{q}| =R/d = d'/R$
つまり,$(q,\ q',\ d,\ d',\ ,R)$の5変数に対して2条件式があるので,3つの量を与えると残りの2つが決定される。例えば,導体球$R$の外に1つの電荷$(q,d)$を置けば,導体球が等電位面となるのに必要なもう一つの電荷$(-q',d')$が定まる。

2023年10月25日水曜日

鏡像法(1)

しばらく前に一度だけ電磁気学の授業を担当したとき,静電場についてはポアッソン方程式の導出までで完結していた。鏡像法などで境界値問題を解くのはあまり好きになれないのでスキップした。

$xy$平面に導体面があって,$z$軸上の点$(0,0,a)$に電荷$q$を置くと,導体平面に誘導電荷が生ずる。導体平面上では電位は一定になり,電場は面に垂直な方向を向く。この境界条件を満足させるため,$(0,0,-a)$に電荷$-q$をおいて,電位と電場を計算する。

電位は,$\displaystyle V(\bm{r}) = \dfrac{q}{4\pi \varepsilon_0} \Bigl\{ \dfrac{1}{\sqrt{r^2+a^2-2 a r  \cos\theta}} - \dfrac{1}{\sqrt{r^2+a^2+2 a r \cos\theta}} \Bigr\}$であり,電場$\bm{E}(\bm{r})$の$z$成分は,$E_z = -\dfrac{\partial V(\bm{r})}{\partial z} = -\dfrac{\partial V(\bm{r})}{\partial r}\dfrac{\partial r}{\partial z} - \dfrac{\partial V(\bm{r})}{\partial \theta}\dfrac{\partial \theta}{\partial z} $となる。

$\dfrac{\partial r}{\partial z} = \dfrac{z}{r} = \cos \theta$なので,後に$\theta =\frac{\pi}{2}$を代入すると0となって,第1項の寄与はない。$\dfrac{\partial \theta}{\partial z} = \dfrac{\partial}{\partial z} \tan^{-1} \dfrac{\sqrt{x^2+y^2}}{z} = \dfrac{-\sqrt{x^2+y^2}}{x^2+y^2+z^2} = -\dfrac{\sin\theta}{r} \rightarrow -\dfrac{1}{r} \ (\theta = \pi/2)$

第2項の前半は,$- \dfrac{\partial V(\bm{r})}{\partial \theta}= \dfrac{q}{4\pi\varepsilon_0} \Bigl\{ \dfrac{ a r  \sin \theta}{(r^2+a^2-2ar\cos\theta)^{3/2}} + \dfrac{a r  \sin \theta}{(r^2+a^2+2ar\cos\theta )^{3/2}}\Bigr\} $

したがって,導体平面上の$\theta = \frac{\pi}{2}$における電場の$z$成分は,
$E_z = - \dfrac{\partial V(\bm{r})}{\partial \theta}\dfrac{\partial \theta}{\partial z} = - \dfrac{q}{4 \pi \varepsilon_0}\dfrac{2 a r }{(r^2+a^2)^{3/2}}\dfrac{1}{r}$
原点を中心とした導体面上の半径$r$の円環面$dS$の電荷面密度を$\sigma(r)$とすると,ガウスの法則から,$E_z(r) dS = \dfrac{\sigma(r) dS}{\varepsilon_0}\quad \therefore \sigma(r)  = \varepsilon_0 E_z(r) = - \dfrac{q }{2 \pi }\dfrac{a}{(r^2+a^2)^{3/2}}$

この誘導電荷の面密度をすべて加えると
$\displaystyle Q= - \dfrac{2 \pi q }{2 \pi } \int_0^\infty \dfrac{a r dr}{(r^2+a^2)^{3/2}} = -q \Bigl [  \dfrac{-a}{(r^2+a^2)^{1/2}} \Bigr ]_0^\infty = -q $

$z$ 軸上においた電荷と絶対値が等しく逆符号の電荷が導体面に誘導されることが確かめられた。

2023年10月24日火曜日

因果関係(2)

因果関係(1)からの続き

物理法則が時間を含む微分方程式で表されている場合はどうなるだろう。

ある量$A(t)$の時間での一階微分が,$\dfrac{dA(t)}{dt}=B(t)$として与えられているとする。このとき,$A(t+dt)=A(t)+B(t) dt$とかけるので,次の時間ステップにおける$A(t+dt)$を決めているのは,その直前の自分自身の値$A(t)$と$B(t)$であり,$B$が$A$の原因を構成しているといえるかもしれない。

一方,$C(t) \equiv \dfrac{dA(t)}{dt}$と定義すると,$B(t)=C(t)$であり,$B(t+dt)=B(t)+\dfrac{dB(t)}{dt}dt$,すなわち,$B(t+dt)=B(t)+\dfrac{dC(t)}{dt}dt=B(t)+\dfrac{d^2A(t)}{dt^2}dt$ と書ける。これを$A(t)$が原因で,$B(t)$が結果として得られたということは可能なのだろうか?


さて,ファラデーの電磁誘導の法則の微分形は,
$\nabla \times \bm{E}(\bm{r},t) = -\dfrac{\partial \bm{B}(\bm{r},t)}{\partial t}$であり,磁石をコイルに近づけたり遠ざけたりすれば,コイルにつながったランプを点灯させることができる。多くの場合これは,磁石の移動(磁場の時間的変化)が原因で,その結果ランプが点灯する(電場の回転の発生)と解釈している。太田さんらを除いて。

しかし,変位電流を含むアンペール=マクスウェルの法則の微分形は,
$\nabla \times \bm{H}(\bm{r},t) = \bm{J}(\bm{r},t)) + \dfrac{\partial \bm{D}(\bm{r},t)}{\partial t}$であり,上記の式と同形であることから,電束密度の変化が原因で磁場の回転が結果として生ずるというふうに解釈したくなる。また実際,電磁波の伝搬はこれによって説明してきた。しかしながら,それがおかしいというのが最近のトレンドになっているのだった。

ところで,これらは$\ \bm{B}(\bm{r},t)$と$\bm{D}(\bm{r},t)\ $が一階の時間微分を含んだ式になっているので,そちら側が因果関係の結果側でなければならないということにはならない傍証としてもよいのだろうか。

謎はつきない。

2023年10月23日月曜日

因果関係(1)

因果と相関(2)からの続き

変位電流が磁場を作るか否かいう議論において,因果と相関の区別ということが強調されることがある。これがなかなかやっかいだ。それが,電場や磁場を遅延時間の電荷密度や電流密度で記述するというジェフィメンコ式の評価につながり,近接相互作用による電場と磁場のからみあいのイメージを放棄させようと迫ってくる。

谷村省吾さんのように,相対論的な因果律(光速度を超える情報の伝搬はない)以外の因果関係を物理法則の中に読み取るべきではないという立場はクリアなのだけれども,そう簡単に割り切ることもできないのが凡人の浅ましさである。


最初に出てくるのが,ニュートンの運動方程式である。$\displaystyle m \dfrac{d^2 \bm{r(t)}}{d t^2}=\bm{F(t)}$は,力が原因が加速度がその結果だというのは疑いのないことで,逆にして,加速度が原因で力が結果として生ずるとは読めないだろう。当然だよね。というのが大方の意見である。

水平面内に摩擦がない誘導路で軌道をつくって小球を走らせるとする。その軌道は自分で設計することができ,摩擦が非常に小さければエネルギー損失はないので,速さは一定になる。加速度はその軌道を設計することで自由に与えることができる。3次元で考えればジェットコースターの設計だ。このときに働く誘導路の壁に与える力(慣性力と等しい)は,軌道による加速度を与えた結果として定まるということはできないだろうか。


2つの時間に依存する物理量$A(t), B(t)$を法則的に結ぶ関係式(方程式)$f(A(t), B(t))=0$があれば,その間には常に同時刻の相関関係が存在し,2つの物理量を変数とした座標でグラフを書くことができる。これを因果関係として理解できるためにはどのような条件が必要だろうか。

一つは,一方の物理量を人為的に自由に設定・制御できる場合,これが原因で他方の物理量が結果といえるというものだ。この場合でも法則は常に満足されるように現象は進行している。ということは,両方の物理量が設定・制御可能であれば,因果関係は関係式(方程式)では決まらず,実際の現象を実現する状況に依存するということになる。

これに近いのが,関係式(方程式)では因果関係は決まらず,初期条件や境界条件などが因果関係を決めるというものだ。確かにニュートンの運動方程式では普通,初期条件を与えて初めて運動が定まる。先ほどの誘導路の設定も広義の初期条件や境界条件の設定といえなくはない。

2023年10月22日日曜日

変位電流

変位電流が磁場をつくるか」という問題は過去からしばしば話題にされてきた。日本ではこの10年ほど,一つは半直線電流+端点のモデルにおける球対称電場を作る磁場がないことを根拠とした議論,もう一つは,ジェフィメンコ式(もしくは相対論的な4次元反対称テンソルによるマクスウェル方程式の表現)から,磁場が源としての実電流だけで記述されることを根拠とした議論が,鈴木,兵頭,中村・須藤などから提起されてきた。


太田浩一(1944-, 1970年代のOhta-Wakamatsuの太田さん)の電磁気学の基礎にもこのモデルを用いて,変位電流が磁場を作らないということが詳しく書かれていたため,その信奉者は増えてしまった。しかしその後,石原,斎藤,北野らによってこのモデルの問題点が明らかにされている。


そこで,主な議論をたどってみた。下記の他に,菅野礼司先生(変位電流と磁場の関係について)や高橋憲明先生の議論もある。最新の兵頭さんの論文にはまだ手が付いていない。


(1) マクスウェルアンペールの法則と変位電流

鈴木亨(物理教育60-1, 2012, 38-43

点電荷の変位電流から求めた磁場と半直線電流からビオ-サバールの法則で求めた磁場は等しくなる。球対称性から変位電流からの磁場は存在しない。

 

(2) 変位電流は磁場を“ 作る” 

兵頭俊夫(物理教育 62-1, 2012, 44-51

点電荷を小導体球に置き換えたモデルだが,本質的に(1)と同じ。球対称性が維持されるので点電荷(クーロン電場)の変位電流は磁場を作らない。

 

(3) 「変位電流は磁場を創らない」を考察するモデルについて

斎藤吉彦(物理教育 60-3, 2012, 209-2012

(1) (2)では荷電粒子の運動を無視していることになる。このためモデルの妥当性が失われていて誤った結論を導いている。ビオ-サバールの法則は近似法則である。↓(7)で訂正。

 

(4) 変位電流は磁場を作るのか?

中村哲・須藤彰三(物理教育 60-4, 2012, 268-273

電流(変位電流)が磁場をつくるというとき,源(source)としてか,作る(presence)としてかを区別する必要がある。因果関係としての源となるのは電流であり,変位電流はそれにあてはまらない。

 

(5) 変位電流と重ね合わせの原理について

石原諭(物理教育 61-4, 2013, 187-189

マックスウェル方程式の重ね合わせの原理を適用するには,それぞれの部分系で電荷保存則が満たされていなければ,正しい結果を与えない。半直線電流を点電荷と電流部分に分割した(2)は電荷保存則を満たしておらず誤った結論を導いている。

 

(6) 変位電流密度の役割

中村哲・須藤彰三(物理教育 62-4, 2014, 23-29

電流密度だけが磁場の因果的源(source )であり変位電流密度は源ではない。電場の時間微分(変位竃流密度)と磁場の空間微分(回転)は時空の属性である電磁場四次元テンソルの微分として統一的に扱うべき。因果律はマックスウェル方程式に内在せず,遅延解の選択で導入される。

 

(7) 半直線電流による電磁場の厳密解

斎藤吉彦(物理教育 62-4, 2014, 155-162

無限遠方から端点までの加速度運動する電荷の集合体がつくる磁場を計算した結果,(1) (2) のビオ-サバール法則の結果と一致した。これは電荷保存則を満たしマックスウェル方程式と矛盾しない。また,磁場の回転が球対称な解であり,Diracの磁気単極子解と同じ構造を持っている。

 

(8) 変位電流をめぐる混乱について

北野正雄(大学の物理教育 27, 2021, 22-25

非定常電流では,「電流がつくる磁場」「変位電流がつくる磁場」は定義できない。

変位電流と伝導電流は不可分でありこれを無視した分割はできない。

ビオ-サバールの式で計算される磁場には変位電流の効果が含まれる。
まだ少し気になっていることなど:

・閉回路(無限直線)でなくても,定常状態や準定常状態ならビオ-サバールの公式は使えて,しかもそれは変位電流の効果を含むものなのか
・微分方程式には因果律は含まれず,境界条件に含まれるのか。

2023年10月21日土曜日

磁気単極子

磁気単極子(ディラックの量子化条件)の話を耳にしたのは大学時代のことだった。基礎工学部の図書館でディラックの原論文を眺めていたような気もするが,ちゃんと読んだことはなかった。その書架には3つのクォークがトポロジカルな結び目の絵で表現されている本も並んでいた。

「変位電流が磁場を作るか」という物理教育の問題に関連して,半直線電流とその先端に電荷がたまる系がしばしば取り扱われる。このとき,磁場の回転が変位電流の球対称場となるような解が必要になる。これが,ディラックの磁気単極子で登場するベクトルポテンシャルの回転が球対称な磁場になる解と同じ構造をもっているという指摘がされていた[1]。

そこで,磁気単極子が作る磁場とベクトルポテンシャルの計算を確かめてみた。

3次元極座標系($\bm{e}_\rho,\ \bm{e}_\theta, \bm{e}_\phi$)では,
$\bm{A} = g \dfrac{ 1 - \cos \theta }{r \sin \theta} \bm{e}_\phi  = g \dfrac{\sqrt{x^2+y^2+z^2}-z}{\sqrt{x^2+y^2+z^2}\ (x^2+y^2)} (-y \bm{e}_x + x \bm{e}_y )$
このとき,
$\nabla \times \bm{A} = g \dfrac{\bm{r}}{r^3}$

ところで,これは$z$軸上では特異的であるから,その部分の構造を極限操作で取り出す。
以下では,円筒座標系($\bm{e}_\rho, \bm{e}_\phi, \bm{e}_z$)で考える。
$\bm{A}_\epsilon = \dfrac{g\ \theta(\rho - \epsilon) }{\rho} \Bigl( 1 - \dfrac{z}{\sqrt{\rho^2+z^2+\epsilon^2}} \Bigr) \bm{e}_\phi = A_\epsilon \bm{e}_\phi$
$\nabla \times \bm{A}_\epsilon = -\dfrac{\partial A_\epsilon}{\partial z} \bm{e}_\rho + \dfrac{1}{\rho} \dfrac{\partial(\rho A_\epsilon)}{\partial \rho} \bm{e}_z$
$=\dfrac{g\ \theta(\rho - \epsilon)}{\rho}\Bigl\{ \dfrac{\rho^2+\epsilon^2}{(\rho^2+z^2+\epsilon^2)^{3/2}}\bm{e}_\rho + \dfrac{\rho z}{(\rho^2+z^2+\epsilon^2)^{3/2}}\bm{e}_z \Bigr\}$
$ + \dfrac{g\ \delta(\rho-\epsilon)}{\rho}\Bigl\{ 1 - \dfrac{z}{\sqrt{\rho^2+z^2+\epsilon^2}}\Bigr\}\bm{e}_z$
$=\dfrac{g\ \theta(\rho - \epsilon)}{\rho^2+z^2+\epsilon^2}\Bigl\{\dfrac{\rho \bm{e}_\rho + z \bm{e}_z}{\sqrt{\rho^2+z^2+\epsilon^2}}\Bigr\}+ \dfrac{g\ \delta(\rho-\epsilon)}{\rho} \Bigl\{2 - 1 - \dfrac{z}{\sqrt{\rho^2+z^2+\epsilon^2}}\Bigr\}\bm{e}_z$
ここでは $z < 0$ という条件を課しているので,最後の$\{\ \}$の中の第2項と第3項の和はゼロになる( $\lim_{\epsilon \rightarrow 0}$ で $\dfrac{\rho\delta(\rho)}{2z^2}$となるため )。
また2次元のデルタ関数について次の関係式が成り立つ。
$\delta(\bm{\rho}) = \dfrac{\delta(\rho)}{\rho} \delta(\phi) = \delta(x) \delta(y)\quad \therefore  \dfrac{\delta(\rho)}{\rho} = 2\pi \delta(x)\delta(y)$
したがって,この条件を保証する$\theta(-z)$を含め,$r^2=\rho^2+z^2$として,
$\lim_{\epsilon \rightarrow 0} \nabla \times \bm{A}_\epsilon =  \nabla \times \bm{A} = g \dfrac{\bm{\hat{r}}}{r^2} + 4\pi g\  \delta(x)\delta(y) \theta(-z) \bm{e}_z$

この式の発散を計算すると,
$\nabla \cdot \nabla \times \bm{A} = - g \nabla \cdot \nabla \dfrac{1}{r} + 4 \pi g \delta(x) \delta(y)  \dfrac{\partial \theta(-z)}{\partial z} = 4 \pi g \delta(\bm{r}) - 4\pi g \delta (\bm{r})= 0$


[1]半直線電流による電磁場の厳密解(斎藤吉彦,2014)

2023年10月16日月曜日

導体球(4)

導体球(3)からの続き

ついでに,電場を取り除いて,導体球に電荷を与えて導体球表面に球対称一様電荷分布が生ずる状況を考える。

先ほどと同様に,観測点の位置ベクトル$\bm{r}$方向に$z$軸をとる。球対称性から$x$軸は自由に設定することができる。この結果,電位は次式で与えられる。

$\displaystyle V(\bm{r}) = \dfrac{\sigma R^2}{4\pi\varepsilon_0} \int \dfrac{ \sin \theta' d\theta' d\phi'}{ \sqrt{r^2+R^2-2rR \cos \theta'}} =  \dfrac{\sigma R^2}{2\varepsilon_0} \int \dfrac{ \sin \theta' d\theta'}{ \sqrt{r^2+R^2-2rR \cos \theta'}}$

再び,$\alpha = r^2+R^2 $,$\beta = 2 r R\ $と置いて,$\sqrt{\alpha -\beta}=| r-R |,\ \sqrt{\alpha + \beta}= r + R\ $である。$t = \cos \theta'$と変数変換して,$ dt = -\sin \theta' d\theta' \ $ なので,
$\displaystyle V(\bm{r}) = \dfrac{\sigma R^2}{2 \varepsilon_0} \int_{-1}^1 \dfrac{dt}{\sqrt{\alpha - \beta t }} = \dfrac{\sigma R^2}{2 \varepsilon_0} \Bigl\lvert \dfrac{-2}{\beta} \sqrt{\alpha - \beta t}\Bigr\rvert_{-1}^1 = \dfrac{\sigma R}{2 \varepsilon_0 r}(\sqrt{\alpha+\beta}-\sqrt{\alpha-\beta})$
$\displaystyle =  \dfrac{\sigma R}{2 \varepsilon_0 r} (r+R -|r-R|)$

したがって,$Q=4\pi R^2 \sigma$と置くと,次のように正しい静電ポテンシャルが得られた。
$\displaystyle V(\bm{r}) = \dfrac{Q}{4\pi \varepsilon_0 R}\quad (r<R)$
$\displaystyle V(\bm{r}) = \dfrac{Q}{4\pi \varepsilon_0 r} \quad (r>R)$

2023年10月15日日曜日

導体球(3)

導体球(2)からの続き

物理科学概説の授業で,表面に一様な電荷が分布する球殻内部の電場や電位の問題を説明しようとした。積分にまで踏み込めないが,立体角を使えばなんとか説明できる。ところでこれを真面目に積分計算しようとすると,電場中の導体球と同じ問題(面倒な楕円関数の積分が必要)が生ずることに今さらながら気がついた。

力学の重力ポテンシャルの場合も同じ問題があったはずで,これまでどうやって回避していたか思い出してみると,観測点の位置ベクトルの方向をz軸にとっている。これにより球対称性から簡単に積分ができていた。この方法が,一様電場中の導体球による表面電荷分布に対しても使えそうな気がしたので再挑戦してみる。

(1) 導体球の中心に置いた原点から観測点Pへの位置ベクトル$\bm{r}$の方向を$z$軸にとる。
そこで,$\bm{r} = (0,\ 0,\ r)$

(2) 一様電場ベクトル方向の導体球面上の位置ベクトル$\bm{e}$の$x-y$平面への射影を$x$軸にとる。このとき,$\bm{e}=(R \sin\lambda,\ 0,\  R \cos\lambda )$,ここで導体球の半径を$R$としている。

(3) 導体球面上の点Qへの位置ベクトルを,$\bm{r'}=(R \sin\theta' \cos\phi',\ R\sin\theta' \sin\phi',\ R\cos\theta')$とする。Qにある電荷要素は,$\rho(\bm{r'}) dS = \sigma R^2 \cos \omega \sin \theta' d\theta' d\phi'$である。ここで,$\sigma$は電荷面密度,$\cos\omega$は,$\bm{e}$と$\bm{r'}$のなす角度であり,$\cos\omega = \frac{\bm{e}\cdot\bm{r'}}{R^2} =  \sin \lambda \sin\theta' \cos\phi' + \cos \lambda \cos\theta' $である。

(4) 観測点Pと電荷要素点Qを結ぶ距離は,$|\bm{r} - \bm{r'}| = \sqrt{r^2+R^2-2rR \cos \theta'}$となる。

そこで,この電荷密度分布$\rho(\bm{r'})$がつくる静電ポテンシャル$V(\bm{r})$は次のようになる。
$\displaystyle V(\bm{r}) = \dfrac{\sigma R^2}{4\pi\varepsilon_0} \int \dfrac{(\sin \lambda \sin\theta' \cos\phi' + \cos \lambda \cos\theta' ) \sin \theta' d\theta' d\phi'}{ \sqrt{r^2+R^2-2rR \cos \theta'}}$

ここで,積分のうち,$\int_0^{2\pi} d \phi'$を実行すると,分子の$\cos\phi'$ を含む項はゼロになり,残りの項は$2\pi$倍となるので,
$\displaystyle V(\bm{r}) = \dfrac{2\pi \sigma R^2}{4\pi\varepsilon_0} \int \dfrac{( \cos \lambda \cos\theta' ) \sin \theta' d\theta'}{ \sqrt{r^2+R^2-2rR \cos \theta'}}$

さらに,$\alpha = r^2+R^2 $,$\beta = 2 r R\ $と置くと,$\sqrt{\alpha -\beta}=| r-R |,\ \sqrt{\alpha + \beta}= r + R\ $である。$t = \cos \theta'$と変数変換して,$ dt = -\sin \theta' d\theta' \ $ なので,
$\displaystyle V(\bm{r}) = \dfrac{\sigma R^2 \cos\lambda}{2 \varepsilon_0} \int_{-1}^1 \dfrac{t dt}{\sqrt{\alpha - \beta t }}$

この積分$I$は部分積分によって実行され,次のような結果を得る。
$\displaystyle I = \int_{-1}^1 \dfrac{t dt}{\sqrt{\alpha - \beta t}} = -\dfrac{4 \alpha}{3 \beta^2}\bigl(  \sqrt{\alpha -\beta} - \sqrt{\alpha + \beta} \bigr) -\dfrac{2}{3\beta} \bigl( \sqrt{\alpha - \beta} + \sqrt{\alpha + \beta}\bigr)$
すなわち,$\displaystyle I= \dfrac{2r}{3 R^2}\ (r<R),\quad I=\dfrac{2R}{3r^2}\ (r>R)$
最終的に,導体球の中の電位は線形になり,電場は一定になる。
$\displaystyle V(\bm{r}) = \dfrac{\sigma \cos\lambda}{3 \varepsilon_0} r  \quad (r<R)$

2023年10月14日土曜日

物理学科同窓会(2)


先週の土曜日に新大阪のワシントンホテルプラザで,阪大物理学科の同窓会があった。1972年(昭和47年)入学なので,昨年が50年目だった。10年前からこの同窓会が概ね毎年開催されるようになった。参加者は順次定年を迎えていくので,だんだん変化に乏しい日々が続き,健康や病気の話題の割合が増えてくる。

今年は,藤原さんのリクエストに端を発して,元原研の佐藤さんに「核融合の最近の現状」というタイトルで話をしていただいた。慣れないzoom経由で,音声トラブルがあったけれど,わかったことは次のとおり。

・ITERもJT60SAも何だかトラブっている。
・ステラレーターは連続だけれど,トカマクは準パルス的な動作をする。
・レーザー核融合のQ=1は,実質は1/100程度になるので,トカマクが有利。
・国内ベンチャーには,研究者がかんでいるが,まとまった炉システムを目指していない。

尾崎さんの,太陽光発電パネルにつけて効率アップする集光シートの特許の顛末もおもしろかった。実用化しようとするとほんの小さな金型だけで2000万円くらいになって,なかなかうまくいかなかったらしい。

2023年10月12日木曜日

導体球(2)

導体球(1)からの続き

一般の電荷密度分布$\rho(\bm{r'})$がつくる静電ポテンシャル$V(\bm{r})$は次のようになる。
$\displaystyle V(\bm{r}) = \dfrac{1}{4\pi\varepsilon_0} \int \dfrac{\rho(\bm{r'}) d\bm{r}'}{|\bm{r} - \bm{r'}|}$
ここで,ポテンシャルの位置座標は,$\bm{r} = (r \sin \theta \cos \phi, r \sin \theta \sin \phi , r \cos \theta)$,電荷素片の位置座標は,$\bm{r'} = (R \sin \theta' \cos \phi', R \sin \theta' \sin \phi' , R \cos \theta')$であり,$R\ $は導体球の半径。
また,静電誘導で導体球表面に誘起される電荷は $\rho(\bm{r'}) d\bm{r}' = \sigma \cos \theta' \sin \theta' d \theta' d\phi'$である。

したがって,
$\displaystyle V(\bm{r}) = \dfrac{\sigma}{4\pi\varepsilon_0} \int \dfrac{\cos \theta' \sin \theta' d\theta' d\phi'}{\sqrt{r^2+R^2-2rR (\sin \theta \sin \theta' \cos \phi' + \cos\theta \cos \theta')}}$
ただし,$\bm{r}$を含む平面内に$x$座標をとって,$\phi = 0$となるようにした。

$\alpha = r^2+R^2 -2rR  \cos\theta \cos \theta' \ge 0$,$\beta = 2 r R \sin \theta \sin \theta' \ge 0$と置くと,
$\displaystyle V(\bm{r}) = \dfrac{\sigma}{4\pi\varepsilon_0} \int \dfrac{\cos \theta' \sin \theta' d\theta' d\phi'}{\sqrt{\alpha - \beta \cos \phi'}}$
この$\phi'$による積分のところ$I_{\phi'}$は楕円積分となる。$\cos \phi' = 1- 2 \sin^2 \frac{\phi'}{2}$とすれば,
$\displaystyle I_{\phi'} = \int_0^{\pi/2} \frac{d\phi '}{\sqrt{(\alpha - \beta) + 2\beta \sin^2 \frac{\phi'}{2}}} = \int_0^{\pi/2} \frac{d\phi '}{\sqrt{(\alpha + \beta) - 2\beta \cos^2 \frac{\phi'}{2}}} $

仮にここまでできたとしても,$\alpha, \beta$に$\theta'$が含まれているものをさらに積分するのはどうしましょうとなった。チーン。

2023年10月11日水曜日

導体球(1)

非常勤講師をとして勤めるのも最後のセメスターになった。物理学概説という科目を担当することになり,原康夫さんの物理学基礎第5版のテキストを使っている。

講義の内容は,電磁気学現代物理学の範囲だ。以前,電磁気学の授業を担当したとき,静電気については,電荷分布から電場と電位(静電ポテンシャル)を導くという展開だった。時間の関係もあって,導体の概念や関連部分は飛ばしてしまっていた。これはだめです。

外部から一様電場をかけた導体球内部には電場が存在せず,電位は一定になる。これを簡単に求めるためには,導体球の中心に電気双極子を置いて,外部電場と重ね合わせるのが普通の教科書の手順だ。砂川さん理論電磁気学では静電ポテンシャルのルジャンドル展開を使ってもっとスマートに導出していた。

このとき,導体表面には球の中心Oを原点とし,電場方向を結ぶ座標軸からの角度の余弦に比例する電荷が分布することになる。それではこのような電荷分布から一様な電場が直接計算できるはずだが,残念ながら探してもその計算を具体的にしている資料はみつからない。

この計算では表面電荷分布を表す2つの角度について積分する必要がある。1つの変数での積分は楕円積分になるが,これをさらに積分するのはちょっと無理そうだ。しかたがないので,数値積分してみると,外部電場方向に垂直な等電位面が現れた。

f[r_, u_] := NIntegrate[ Cos[t] Sin[t]/
     Sqrt[r^2 + 1 - 2 r*(Sin[u] Sin[t] Cos[s] + Cos[u] Cos[t])],
             {t,  0, Pi}, {s, 0, 2 Pi}] / (r*Cos[u])

ここで,導体球の半径をR=1として,内部の点を(r sin u, 0, r cos u ),導体球面上の電荷要素の位置を (R sin t cos s, R sin t sin s, R cos t ) として,変数tと変数sで積分している。電荷分布は σ cos t で,積分結果を内部点のz座標 r cos u で割った。この結果が,内部点の座標変数 r, u を変えても一定になったので,等電位面が出現したことになる。



図:一様電場中の導体球(前野昌弘さんのテキストから引用)



2023年9月6日水曜日

日下周一

クリストファー・ノーラン映画オッペンハイマーは,まだ日本未公開なので,早くみたいものだと調べているうちに,オッペンハイマー(1904-1967)の弟子だった日下周一(1915-1947)に行き当たった。

日下周一は,大阪に生まれ,5歳で移住したカナダで育ち,ブリティッシュ・コロンビア大学からMITの大学院に進んだ。カリフルニア大学バークレー校でオッペンハイマーの指導を受け,1942年に博士号を取得した。中間子や核力の研究を行い,1943年には中間子論でパウリ(1900-1958)との共著論文(On the theory of a mixed pseudoscalar and a vector meson field)がある。1946年にプリンストン大学でウィグナー(1902-1995)の助教授になったが,1947年31歳で海水浴中に溺死した。

1940年に日下が日本を訪れたとき,湯川秀樹(1907-1981)や小林稔(1908-2001)や内山龍雄(1916-1990)に会って議論している。内山先生の一歳上だったのか。

1960年にオッペンハイマー夫妻が来日したとき,バークレー時代の弟子である日下周一の両親に会って弔意を表した。

1992年,大阪市立科学館の加藤賢一のところに吉永剛幸および大山幹男が来訪し,日下周一を顕彰することができないかとの話があった。2005年10月には大阪市立科学館(高橋憲明館長)において,関係者が中心の日下周一シンポジウムが開催されている。



写真:日下周一の御両親を弔問するオッペンハイマー夫妻(Esquire記事から引用)

[1]日下周一伝(加藤賢一データーセンター)
[4]物理学者日下周一の生涯と業績(星学館ブログ)
[5]NECROLOGY Kusaka Shuichi, 1915-1947(David Bohm, Robert R. Bush)
[6]Shuichi Kusaka  Theoretical Physicist ( When An Old House Whispers )
[7]Great Physicist Shuichi Kusaka  (You Tube)
[8]Notes on electrodynamics (J.R. Oppenheimer, University of California, Physics 207B, 1939)
[9]電気力学(オッペンハイマー講義録=日下周一のノート,小林稔訳)
[10]Einstein His Life & Times (Frank Phillipe, George Rosen, Shuichi Kusaka)
[12]Shuichi Kusaka (Find a Grave)


2023年9月4日月曜日

トリチウム(1)

話題のトリチウムだ。

放射性崩壊のベクレルという単位は,1秒間に1崩壊を表わしている。放射性物質の半減期をT [s] とすると,t秒における放射性粒子数は,N(t) = N(0) exp( - 0.693 t /T ) にしたがって減少する。そこで,1秒間の崩壊で減少する粒子数 x は, N(0)-N(1) ≒ N(0) × 0.693 1/T = x [Bq] となる。

トリチウムの半減期,12.3年は 12.3 [y] × 3.16 10^7 [s/y] = 3.88 10^8 [s] である。したがって,1モル = 6.02 10^23個のトリチウムは,1.08 10^15 Bq(1080兆ベクレルに相当する。福島第一原発のタンクにあるトリチウム水( HTO = 分子量20)が,830兆ベクレルだとすれば,0.77モルにあたるので,134万トンのトリチウム汚染水には約15gのトリチウム水が含まれている。

放射性物質の量を表すとき,ベクレル単位にすると非常に大きな量として印象づけられるが,g単位であれば影響が小さいように見せることもできる。面倒な話である。

トリチウム原子は1PBq = 10^15 Bqで1モルなので,これを使うと次の数値が得られる。
地球上存在する全トリチウム〜140kg,宇宙線で自然に生成されるトリチウム 200 g/年,原子力発電所などで生成されるトリチウム 300 g/年(うち1/4が放出カナダがCANDU炉で人工的に生成しているトリチウム 2 kg/年( 300万円/g),韓国のCANDU炉でも数百g/年 とある。


さて話はかわり,最近リバイバルしている核融合だ。国際熱核融合実験炉(ITER)のトカマクでは,DT反応が用いられる。すなわち重陽子とトリチウムが主な燃料となる。100万kWの核融合炉を1年運転させる(3×10^16J)ためには,効率100%としても14MeV(2×10^-12J)のエネルギーを持つ中性子が放出されるDT反応を10^28回 生起する必要があり,10^4mol = 30kgのトリチウムが必要になる。どこからかき集めてくるのだ

Wikipediaの「トリチウム」によれば,トカマクの場合点火時に3kg 用意すればよいとある。あとはブランケットのLiに中性子を吸わせて,6Li+n → T + 4He + 4.8MeV,7Li+n → T + 4He + n -2.5MeV で,トリチウム(T)を再生産するらしい。それでも消費燃料(D+T)が500g/日という記述になっている,どういうこと。レーザー核融合だったらベレットを作る手間が発生するが,どのみちリチウムブランケットからのトリチウム取り出しサイクルは必要になってくる。

核融合では高レベル廃棄物が出ないといえども,134万トンに薄まった15gのトリチウム水でばたついている人類が,数kgのトリチウムを自由に扱える日がくるのだろうか。


[2]よくわかる核融合炉の仕組み(日本原子力学会核融合工学部会)
[4]環境トリチウム—その挙動と利用(高島良正,1991)

2023年8月14日月曜日

LK-99

7月22日にarxivで公開された常温常圧超伝導の話題が,8月に入るとすぐ盛り上がっていた。

LK-99 が,常温常圧超伝導を示すといわれる物質名である。話題の論文の筆頭著者(化学者)である,S. Lee (李石培 이석배) と J-H. Kim (金智勳 김지훈) が1999年に発見した。六方晶系の鉛アパタイト(Pb10(PO4)6O) の鉛のいくつかを銅で置換したものであり,論文には登録商標マーク(LK-99®)があって,特許も取得している

中央大学の田口善弘さんが,arxivに上がっているプレプリントのクオリティをディスっていた。もしかするとその一部は著者らの専門が化学であって物理分野とは違う文化であることに起因するのかもしれない。物質名(通称)の命名方法(著者名イニシャルを含める)とか登録商標についてもそのあたりなのだろうか。

追試過程の報告があちこちからでているが,ネガティブなものと気持ちポジティブなものが混在していて,いきなり全否定というわけでもないようだ。例の常温超高圧超伝導の件よりは少しマシかもしれない。なお,このグループは2020年にNatureに同趣旨の論文を投稿しているが不採用だった。さらに,今回の論文を巡っては著者グループ間には微妙な確執があるとかないとかいう話だ。

この実験を受けて,理論サイドでは,密度汎関数理論(DFT)などの第一原理計算によるシミュレーションがなんらかの可能性を示唆するという論文が続出している。LK-99がだめでも新しい物質の可能性があるのではないかという楽しげな雰囲気も漂っている。まあどんなケースでもそれらしい理論は作れてしまうというのが世の常なのだけれど。

1986年の高温超伝導フィーバーのときは,物性実験の人達がこぞって乳鉢で材料を調整し,論文を書いていたが,そのときの熱気に近いものが立ち上がりつつある。YouTubeの浮上実験の動画を見て,単なる反磁性だという説とか反磁性だとしてもすごいのではないかという説が入り交じり,DIY素人が実験に参戦しつつあるらしい

なお,最新の変化しつつある情報は,英語版WikipediaのLK-99に詳しい。



図:高温超伝導の歴史(Wikipediaから引用)

P. S.  8月第2週に入って,400Kにおける抵抗値の減少が不純物のCuSの1次相転移によるものであり,浮上は強磁性由来だということで決着しそうな気配がただよってきた。祭りは終了。


2023年7月31日月曜日

コーシー=シュワルツの不等式


数理統計学を真面目に勉強してこなかったのでいろいろ不都合が生じている。統計的因果推論とか深層機械学習とか量子測定理論とか,簡単に読み砕けない資料がたくさんたまる。

授業で扱った最小二乗法と実験誤差の話を整理しようとしても,背景には数理統計学が控えている。昔,阪大の南園グループによるベータ崩壊の実験と我々の理論を突き合わせたときに,χスクェアフィットの計算を散々繰り返したけれど,所与の公式を使うだけであってその理論的根拠をつきつめて考えたはしなかった。

そこで最初から勉強を始めようとすると,いきなり確率変数でつまづくのだった。コンピュータプログラムのサブルーチンや関数のようなものだと思えば納得できるといえばいえるのだけれど,自然言語と数学的記号を使って理解しようとするとなかなかその本質がつかみきれない。入門書は沢山あるけれど,どれも何だか気持ち悪い。

竹村彰道(1952-)さんの現代数理統計学の本(旧版)が手元にあって,読みやすいかなとページをめくってみると,記述統計の復習から始まった。これなら大丈夫かと思いきや,いきなり,標本相関関数の大きさが -1から 1の範囲に限定されることは,コーシー=シュワルツの不等式を用いて容易に示すことができると説明無しにあった。

n次元ユークリッド空間のベクトルの内積の話だと思えばそのとおりなのだけれど,証明したことはなかったかも。Wikipediaでは数学的帰納法で証明していた。$A_k=(a_1,\ a_2,\ \cdots,\ a_k),\ B_k=(b_1,\ b_2,\ \cdots,\ b_k),\ $として,$\displaystyle S^{aa}_k=\sum_{i=1}^k a_i^2,\ S^{bb}_k=\sum_{i=1}^k b_i^2,\ S^{ab}_k=\sum_{i=1}^k a_i b_i, \quad R^{ab}_k=\frac{S^{ab}_k}{\sqrt{S^{aa}_k S^{bb}_k}} $
つまり,$ \bigl( S^{ab}_k \bigr)^2  \le S^{aa}_k S^{bb}_k$を証明すれば良い。

$k=1$の場合は,$ \bigl( S^{ab}_1 \bigr)^2 -  S^{aa}_1 S^{bb}_1 = (a_1 b_1)^2- (a_1^2)(b_1^2) = 0 $

$k=2$の場合は,$ \bigl( S^{ab}_2 \bigr)^2 -  S^{aa}_2 S^{bb}_2 = (a_1 b_1+a_2 b_2)^2- (a_1^2+a_2^2)(b_1^2+b_2^2) =  -(a_1 b_2- a_2 b_1)^2  < 0 $

$k \ge 2$に対して,$ \bigl( S^{ab}_k \bigr)^2  \le S^{aa}_k S^{bb}_k$ が成り立つと仮定して,$k+1$の場合を考える。与式は,$ \bigl( S^{ab}_k + a_{k+1}b_{k+1} \bigr)^2 - \bigl( S^{aa}_k + a_{k+1}^2 \bigr) \bigl(  S^{bb}_k + b_{k+1}^2 \bigr) $
$= \bigl( S^{ab}_k  \bigr)^2 - S^{aa}_k S^{bb}_k - \Bigl( a_{k+1}^2 S^{bb}_k + b_{k+1}^2 S^{aa}_k -2 a_{k+1}b_{k+1} S^{ab}_k \Bigr)$
$= \bigl( S^{ab}_k  \bigr)^2 - S^{aa}_k S^{bb}_k - \sum_{i=1}^k \Bigl( a_{k+1}^2 b_i^2 + b_{k+1}^2 a_i^2 -2 a_{k+1}b_{k+1} a_i b_i \Bigr)$
$=\bigl( S^{ab}_k  \bigr)^2 - S^{aa}_k S^{bb}_k - \sum_{i=1}^k \Bigl( a_{k+1} b_i - b_{k+1} a_i \Bigr)^2 < 0$

Wikipediaの証明などでは,$a_i, b_i >0$の場合だけに妥当するものが多いのでちょっと困る。
まあ,$\displaystyle f_k(x) = \sum_{i=1}^k (a_i x - b_i)^2$ の判別式$D \le 0$から証明するのが最も簡単なのだけど。


[1]賢者に学ぶ統計学の智(西内啓×竹村彰通,ダイヤモンド社)

2023年7月30日日曜日

最小二乗法(6)

最小二乗法(5)からの続き

実験データを$y = a x + b$にフィットする場合,最小二乗法で$(a,\  b)$とその平均二乗誤差$(\sigma_a^2,\ \sigma_b^2)$を求めてきた。これを,$y = f(x) = a x^2 + b x + c\ $に拡張して,自由度3が登場するかどうかを確認してみる。吉澤康和さんの「新しい誤差論」には結果だけ書いてある。

(1) a, b, c を決定する正規方程式とその解

$ \begin{pmatrix}\overline{x^4} & \overline{x^3} & \overline{x^2} \\ \overline{x^3} & \overline{x^2} & \overline{x^1} \\ \overline{x^2} & \overline{x} & 1 \\ \end{pmatrix} \begin{pmatrix} a \\ b \\ c \\ \end{pmatrix}= \begin{pmatrix}\overline{x^2\ y}\\ \overline{x\ y} \\ \overline{y} \\ \end{pmatrix}$ 

$ \begin{pmatrix} a \\ b \\ c \\ \end{pmatrix}= \dfrac{1}{\Delta_3} \begin{pmatrix}\overline{x^2\ y}(\overline{x^2}-\overline{x}^2)+\overline{x\ y}(\overline{x^2}\overline{x}-\overline{x^3})+\overline{y}(\overline{x^3}\overline{x}-\overline{x^2}^2)  \\ \overline{x^2\ y}(\overline{x^2}\overline{x}-\overline{x^3}) + \overline{x\ y}(\overline{x^4}-\overline{x^2}^2) + \overline{y}(\overline{x^3}\overline{x^2}-\overline{x^4}\overline{x}) \\ \overline{x^2\ y}(\overline{x^3}\overline{x}-\overline{x^2}^2) + \overline{x\ y}(\overline{x^3}\overline{x^2}-\overline{x^4}\overline{x}) +\overline{y}(\overline{x^4}\overline{x^2}-\overline{x^3}^2) \\ \end{pmatrix}$ 

ただし,$\Delta_3 = \overline{x^4}\overline{x^2}+2\overline{x^3}\overline{x^2}\overline{x}-\overline{x^2}^3-\overline{x^3}^2-\overline{x^4}\overline{x}^2$

$y_i$を共通の平均二乗誤差$\sigma^2_y$を持つ独立変数として,誤差伝播の法則より,

$\displaystyle \sigma_a^2= \sum_{i=1}^n \Bigl( \frac{\partial a}{\partial y_i}\Bigr) ^2 \sigma_y^2, \quad \sigma_b^2= \sum_{i=1}^n \Bigl( \frac{\partial b}{\partial y_i}\Bigr) ^2 \sigma_y^2 , \quad \sigma_c^2= \sum_{i=1}^n \Bigl( \frac{\partial c}{\partial y_i}\Bigr) ^2 \sigma_y^2 $ 

さらに,真の値$f_0(x_i)=a_0 x_i^2 + b_0 x_i + c_0$に対して,$\varepsilon_i = y_i -f(x_i)+ f(x_i) -f_0(x_i) =  \delta_i + f(x_i) -f_0(x_i) $ として,$\displaystyle \sigma_y^2 = \frac{1}{n} \sum_{i=1}^n \varepsilon_i^2 = \frac{1}{n} \sum_{i=1}^n \Bigl\{ \delta_i^2 + \tilde{\sigma}^2_{f(x_i)} \Bigr\}$

ところで,$\displaystyle \tilde{\sigma}^2_{f(x_i)} = \sum_{j=1}^n \Bigl\{ \frac{\partial(a x_i^2 + b x_i + c)}{\partial y_j}\Bigr\}^2$ であり,この項を再度  $\sigma_y^2$ で表してもとの式に戻して計算すれば良い。

つまり,$\displaystyle \frac{\partial a}{\partial y_j}, \  \frac{\partial b}{\partial y_j},\  \frac{\partial c}{\partial y_j}$が計算できればよいことになる。
$\displaystyle \frac{\partial a}{\partial y_j}=\frac{1}{n \Delta_3}\Bigl\{ x_j^2 (\overline{x^2}-\overline{x}^2)+ x_j(\overline{x^2}\overline{x}-\overline{x^3})+(\overline{x^3}\overline{x}-\overline{x^2}^2) \Bigr\}$
$\displaystyle \frac{\partial b}{\partial y_j}=\frac{1}{n \Delta_3}\Bigl\{ x_j^2 (\overline{x^2}\overline{x}-\overline{x^3}) + x_j (\overline{x^4}-\overline{x^2}^2) + (\overline{x^3}\overline{x^2}-\overline{x^4}\overline{x}) \Bigr\}$
$\displaystyle \frac{\partial c}{\partial y_j}=\frac{1}{n \Delta_3}\Bigl\{ x_j^2 (\overline{x^3}\overline{x}-\overline{x^2}^2) + x_j (\overline{x^3}\overline{x^2}-\overline{x^4}\overline{x}) +(\overline{x^4}\overline{x^2}-\overline{x^3}^2) \Bigr\}$

Mathematicaの力を借りると,計算結果が因数分解できて分子から$\Delta_3$が出る。
$\displaystyle \sigma_a^2 = \sum_{j=1}^n \Bigl( \frac{\partial a}{\partial y_j}\Bigr) ^2 = \frac{1}{n \Delta_3} \bigl( \overline{x^2} -\overline{x}^2 \bigr) \sigma_y^2$
$\displaystyle \sigma_b^2 = \sum_{j=1}^n \Bigl( \frac{\partial b}{\partial y_j}\Bigr) ^2 = \frac{1}{n \Delta_3} \bigl( \overline{x^4} - \overline{x^2}^2 \bigr) \sigma_y^2$
$\displaystyle \sigma_c^2 = \sum_{j=1}^n \Bigl( \frac{\partial c}{\partial y_j}\Bigr) ^2 = \frac{1}{n \Delta_3} \bigl( \overline{x^4} \overline{x^2}-\overline{x^3}^2  \bigr) \sigma_y^2$

$\displaystyle \sum_{j=1}^n \Bigl( \frac{\partial a}{\partial y_j}\frac{\partial b}{\partial y_j}\Bigr)  = \frac{1}{n \Delta_3} \bigl( \overline{x}\overline{x^2} -\overline{x^3} \bigr) \sigma_y^2$
$\displaystyle \sum_{j=1}^n \Bigl( \frac{\partial b}{\partial y_j}\frac{\partial c}{\partial y_j}\Bigr)  = \frac{1}{n \Delta_3} \bigl( \overline{x^2}\overline{x^3} - \overline{x}\overline{x^4} \bigr) \sigma_y^2$
$\displaystyle \sum_{j=1}^n \Bigl( \frac{\partial c}{\partial y_j}\frac{\partial a}{\partial y_j}\Bigr)  = \frac{1}{n \Delta_3} \bigl( \overline{x^4} \overline{x^2}-\overline{x^3}^2  \bigr) \sigma_y^2$

このとき
$\displaystyle \tilde{\sigma}^2_{f(x_i)}= \sum_{j=1}^n \Bigl\{ \frac{\partial a}{\partial y_i} x_i^2 + \frac{\partial b}{\partial y_j} x_i + \frac{\partial c}{\partial y_j} \Bigr\} ^2 = \frac{\sigma_y^2}{n \Delta_3} $
$\Bigl\{\bigl( \overline{x^2}-\overline{x}^2 \bigr) x_i^4 + 2 \bigl( \overline{x}\overline{x^3} -\overline{x^3} \bigr) x_i^3 + \bigl( \overline{x^4}-\overline{x^2}^2 + 2( \overline{x^3}\overline{x} - \overline{x^2}^2) \bigr) x_i^2 $
$+ 2\bigl( \overline{x^2}\overline{x^3} - \overline{x} \overline{x^4} \bigr) x_i + \bigl( \overline{x^2}\overline{x^4}-\overline{x^3}^2 \bigr)  \Bigr\}$

$x_i$について平均操作するとMathematicaを使い分子から$\Delta_3$が出ると。
$\displaystyle \frac{1}{n}\sum_{i=1}^n  \tilde{\sigma}^2_{f(x_i)} = \frac{\sigma_y^2}{n \Delta_3}$
$\Bigl\{\bigl( \overline{x^2}-\overline{x}^2 \bigr) \overline{x^4}+ 2 \bigl( \overline{x}\overline{x^3} -\overline{x^3} \bigr) \overline{x^3} + \bigl( \overline{x^4}-\overline{x^2}^2 + 2( \overline{x^3}\overline{x} - \overline{x^2}^2) \bigr) \overline{x^2} $
$\displaystyle + 2\bigl( \overline{x^2}\overline{x^3} - \overline{x} \overline{x^4} \bigr) \overline{x}+ \bigl( \overline{x^2}\overline{x^4}-\overline{x^3}^2 \bigr)  \Bigr\} = \frac{3}{n} \sigma_y^2$

したがって,自由度n-3の場合の式が得られた。
$\displaystyle \sigma_y^2 = \frac{1}{n} \sum_{i=1}^n \Bigl\{ \delta_i^2 + \tilde{\sigma}^2_{f(x_i)} \Bigr\} =  \frac{1}{n} \sum_{i=1}^n  \delta_i^2 + \frac{3}{n} \sigma_y^2$
$\displaystyle \therefore \sigma_y^2 = \frac{1}{n-3} \sum_{i=1}^n \delta_i^2$

2023年7月25日火曜日

最小二乗法(5)

最小二乗法(4)からの続き

完全にスッキリしなくて何だか気持ち悪いのだけれど,いきなり自由度がとかいわれて$n-2$が出てくるのがいやなので,吉澤さんの本に従って話を進めてみる。

$\displaystyle \tilde{\sigma^2}_{y(x_i)} =  \frac{1}{n}\sum_{i=i}^n  \tilde{\varepsilon_i}^2 = \frac{1}{n}\sum_{i=i}^n  \Bigl\{ a x_i + b - a_0 x_i - b_0  \Bigr\}^2$
これから,$f(x_i) = y(x_i) =  a  \bm{x_i} + b$として,独立変数$y_j$について,
$\displaystyle \tilde{\sigma^2}_{y(x_i)} =  \sigma_y^2 \sum_{j=i}^n  \Bigl\{ \frac{\partial a}{\partial y_j}\bm{x_i} + \frac{\partial b}{\partial y_j} \Bigr\}^2 = \frac{\sigma_y^2}{n^2 \Delta^2} \sum_{j=i}^n  \Bigl\{ (x_j-\overline{x}) \bm{x_i} + ( \overline{x^2} -\overline{x} x_j ) \Bigr\}^2$
$\displaystyle = \frac{\sigma_y^2}{n^2 \Delta^2} \sum_{j=i}^n  \Bigl\{ ( \bm{x_i}-\overline{x} ) x_j + ( \overline{x^2} - \overline{x}  \bm{x_i} ) \Bigr\}^2$
$\displaystyle = \frac{\sigma_y^2}{n \Delta^2}  \Bigl\{ \overline{x^2} ( \bm{x_i}-\overline{x} )^2 + 2 \overline{x} (\bm{x_i} - \overline{x})(\overline{x^2} -\overline{x} \bm{x_i}) + ( \overline{x^2} - \overline{x}  \bm{x_i} )^2  \Bigr\}$
$\displaystyle = \frac{\sigma_y^2}{n \Delta^2}  \Bigl\{ \bm{x_i}^2 ( \overline{x^2} - \overline{x}^2) + 2 \bm{x_i} (\overline{x}^3 - \overline{x^2} \overline{x}) + ( \overline{x^2}^2 - \overline{x^2} \overline{x}^2 )  \Bigr\}$
$\displaystyle = \frac{\sigma_y^2}{n \Delta}  \Bigl\{ \bm{x_i}^2  - 2 \bm{x_i} \overline{x} + \overline{x^2}  \Bigr\}$

添え字 $i$について平均すると,$\displaystyle \frac{1}{n}\sum_{i=1}^n \tilde{\sigma^2}_{y(x_i)} =\frac{\sigma_y^2}{n \Delta}\Bigl\{ \overline{x^2}  - 2 \overline{x} \overline{x} + \overline{x^2}  \Bigr\} =  \frac{2 \sigma_y^2}{n}$
そこで,
$\displaystyle \sigma_y^2 =\frac{1}{n}\sum_{i=1}^n \Bigl\{ \delta_i^2 + \tilde{\varepsilon_i}^2 \Bigr\} = \frac{1}{n}\sum_{i=i}^n \delta_i^2 + \tilde{\sigma^2}_{y(x_i)} = \frac{1}{n} \sum_{i=1}^n \delta_i^2 + \frac{2 \sigma_y^2}{n}$
$\displaystyle \therefore \sigma_y^2 = \frac{1}{n-2}\sum_{i=1}^n \delta_i^2 = \frac{1}{n-2}\sum_{i=1}^n (y_i - a x_i -b )^2$

2023年7月24日月曜日

最小二乗法(4)

最小二乗法(3)からの続き

$(a,\  b)$  に対する平均二乗誤差,$(\sigma_a^2,\ \sigma_b^2)$を考える。$(a,\  b)$ は直接測定された$(x, \ y)$の関数であるが,このうち$x_i$の誤差は非常に小さく,$y_i$の誤差だけが$n$個の独立変数として伝搬して$(a,\  b)$ に反映すると仮定する。ただし,各$y_i$自身の平均二乗誤差は共通でありこれを$\sigma_y^2$とおく。

誤差伝播の法則より,
$\displaystyle \sigma_a^2= \sum_{i=1}^n \Bigl( \frac{\partial a}{\partial y_i}\Bigr) ^2 \sigma_y^2 = \frac{\sigma_y^2}{n^2 \Delta^2} \sum_{i=1}^n  \Bigl( x_i-\overline{x} \Bigr) ^2 = \frac{\sigma_y^2}{n \Delta^2} \Bigl( \overline{x^2}-\overline{x}^2 \Bigr) = \frac{\sigma_y^2}{n \Delta} $

$\displaystyle \sigma_b^2= \sum_{i=1}^n \Bigl( \frac{\partial b}{\partial y_i}\Bigr) ^2 \sigma_y^2 = \frac{\sigma_y^2}{n^2 \Delta^2} \sum_{i=1}^n  \Bigl( \overline{x^2}-\overline{x}x_i \Bigr) ^2 = \frac{\sigma_y^2 \ \overline{x^2}}{n \Delta^2} \Bigl( \overline{x^2}-\overline{x}^2 \Bigr) = \frac{ \sigma_y^2 \ \overline{x^2}}{n \Delta}$

残るは,$\displaystyle \sigma_y^2 = \frac{1}{n} \sum_{i=1}^n (\varepsilon_i)^2\  $を実験値から導くことになる。ここで,$ \varepsilon_i = y_i-(a_0 x_i + b_0)  = y_i - (a x_i + b) + (a x_i + b) -(a_0 x_i + b_0) = \delta_i + \tilde{\varepsilon_i}$ である。
ただし,$a_0 x_i + b_0$が未知の真値,$a x_i + b$が平均値に対応し,$ \delta_i$が残差, $\tilde{\varepsilon_i}$が平均値の誤差に相当する。

$\displaystyle \therefore \sigma_y^2 =\frac{1}{n}\sum_{i=i}^n \Bigl\{ \delta_i^2 + \tilde{\varepsilon_i}^2 \Bigr\}$ ここで,$\displaystyle \frac{2}{n} \sum_{i=1}^n \delta_i \tilde{\varepsilon}_i =0$ である。なぜならば$\tilde{\varepsilon}_i$は$x_i$の一次関数であり,正規方程式より, $\sum_{i=1}^n \delta_i = 0$ と $\sum_{i=1}^n \delta_i x_i=0$ が成り立つから。

そこで,$y(x_i)=a x_i+b$として,$\displaystyle \tilde{\sigma^2}_{y(x_i)} =  \frac{1}{n}\sum_{i=i}^n  \tilde{\varepsilon_i}^2 = \frac{1}{n}\sum_{i=i}^n  \Bigl\{ a x_i + b - a_0 x_i - b_0  \Bigr\}^2$を求めることになるが,ここで,$(a, b)$が $y_i$の関数として誤差伝搬の法則を再度使って,$\sigma_y^2$で表せばよい(と吉澤康和さんの「新しい誤差論(1989)」に書いてあった)。


2023年7月23日日曜日

最小二乗法(3)


物理量 $x$を設定したとき,$y$が測定される。$n$回測定では,$(x_1,\ y_1),\ (x_2,\ y_2),\ \cdots (x_n,\ y_n)$ が得られたとする。2つの物理量の間には,$y\ =\ a x + b$という1次関数の関係があって,$(a,\ b)$にも物理量としての意味がある。

この$(a, \ b)$を求めるため,$\displaystyle S(a,b)=\frac{1}{n}\sum_{i=1}^n (y_i-a x_i -b)^2$を最小化するという条件を課す。すなわち,$\frac{\partial S}{\partial a}=0, \frac{\partial S}{\partial b}=0, $これから次の$(a,\ b)$に関する連立方程式(正規方程式)が得られる。

$\displaystyle \frac{1}{n} \sum_{i=1}^n x_i \bigl( y_i - a x_i - b \bigr) = 0 \rightarrow \quad a \overline{x^2} + b \overline{x} = \overline{xy} $
$\displaystyle \frac{1}{n} \sum_{i=1}^n \bigl( y_i - a x_i - b \bigr) = 0  \quad \rightarrow \quad a \overline{x} + b = \overline{y} $

これを解くと次の解が得られる。ただし,$\Delta = \overline{x^2} - (\overline{x})^2$ である。
$a=\frac{1}{\Delta}\bigl(\overline{xy}-\overline{x} \cdot \overline{y} \bigr)$
$b=\frac{1}{\Delta}\bigl( (\overline{x^2}\cdot \overline{y}-\overline{x} \cdot \overline{xy} \bigr)$