2024年2月7日水曜日

Apple Vision Pro(2)

Apple Vision Pro(1)からの続き

2月2日にAppleのVision Proが発売になって,数日が経過した。アメリカ国内だけの限定発売であり,日本語にもまだ対応していないが,YouTubeでは日本人によるファーストインプレッションも沢山あがってきた。

最低の256GBモデルでも3499ドル(52万円),多少オプションを付ければ60万円のシステムなので,おいそれとは手が出ない(昔の為替レートなら30万円なのだけど)。その上現地までの航空費や宿泊費を含めると,仕事関係の人か,収益化につながるYouTuberかよほどのアップルフリークでないと参戦していない(米国の連絡先/配送先住所と米国のAppleIDが必要)。

なおかつ,アップル製品への関心が高くてテクノロジーに関する経験の深い人々による印象なので,割り引いておかないと本当のところはなかなかわからないかもしれない。それでも,その評判はなかなかよいものだった。

視野角は左右90度で上下も狭いが,外部カメラからとりこまれる環境映像が自然にパススルーされている。そこに非常に高精細なVision OSのUI画面が重ねられ,丁寧に個人向けに調整された視覚環境は申し分ないようだ。ユーザインタフェースも視線と指による(音声は英語がひつようなのであまりレビューされていない)もので,まったく違和感なしに使えている。自分が2008年8月に最初に手にしたiPhoneのタッチ画面を楽々と使いこなせたのと同じようだ。

手元にMacbookがあればその画面を表示できて,マウスコントロールやデータのやりとりがVision OS固有アプリとの間で簡単にできるところも素晴らしい。

外付けになるバッテリは妥協できそうだが,問題は,650gの重さのようだ。あと,人によっては眼がつかれるとのこと。ブラウザのYouTubeは問題なく視聴できるらしく,Apple のVisionOS ネイティブアプリもあって,iOSアプリも利用できるらしいからほとんど問題はなさそうだ。

ただ,空間コンピューティングを実現する新しいディスプレイとコンピュータだと考えた場合は,macOSのように,Xcode環境はターミナルが動いて,JuliaとMathematicaまで動かなければちょっとためらわれる。もしそれが可能になれば,いまのMacBook Airとディスプレイを完全におきかえるものとなる(スペック上はM1 MacBook Airを越えられる)。ただし,G5とGPSが使えなければ,モバイル端末としてのiPhone/iPadを代替するまでにはならない。



写真:店頭のApple Vision Pro (Wikipediaから引用)

[4]Apple Vision Pro - できること全て(大川優介)

2024年2月6日火曜日

中国の小学生の数学

中国の小学生が解いている数学の問題という触れ込みで次の面積を求める問いがあった。

図1:4分円と半円の交わる部分の面積を求める

小学生にも出来るはずだということで,いろいろ考えたけれど,どうしても解けない。いや,$ \tan^{-1} \alpha = 1/2, \tan^{-1}\beta = 2$によって,図の2つの角度さえ求めてよいならば,扇型AOEGの面積が $S_1=\alpha a^2$,扇型DOEGの面積が$S_2=\beta (a/2)^2$,そして四角形OAGDの面積が,$S_3=a^2/2$であることを用いて,求める面積は $S = S_1+S_2-S_3$となる。

あるいは,解析幾何学を使ってよいのならば,2つの円の式の交点からG$=(4a/5, 2a/5)$となり,面積は積分を使って,$S=\displaystyle \int_0^{4a/5} \Bigl( \sqrt{(a/2)^2-(x-a/2)^2}-a+\sqrt{a^2-x^2} \Bigr)\  dx$となる。

いずれにせよ,答えは,正方形OABCの一辺を$a=4$として,$S \approx 3.847$ である。
結局,中国の小学生はどうやってこの問題を解いているのだろうか


小学生のとき,似たような問題で長いこと未解決でクラスのみんなであれこれ議論したものがあった。それは図2右のようなもので,正方形の中の四つの四分円の交わる領域の面積を求めるものだ。図2左は授業でもよく出てくる問題であり,これならみんな解ける。

図2:小学校のときの未解決問題(右図)

あるとき,塾に通っていた友人たちが,塾の先生から答えを聞いてきて披露したことがあった。それはだめでしょう。せっかくみんなで自分たちで答えをだそうとがんばっていたのに。その解法には正三角形の面積を求める過程がふくまれていて平方根が登場する。小学生には無理な問題だったのだ。

いや,じつはそれほど無理でもない。小学校5,6年のときだろうか,学校で一番頭の良いことで有名だった大杉君というのが,平方根の筆算による計算法(開平法をどこかで学んできて,みんなに教えてくれたことがあった。なるほど,そういうことかと計算できるようになった友達は多い。たぶん,ピタゴラスの定理もどこかで聞きかじっていたかもしれないので,実はもう少しで解けるあたりまでの知識は蓄積していたはずなのだ。

2024年2月5日月曜日

三角分布と変数変換

一様分布と変数変換からの続き

2つの確率変数の三角分布があったとき,これを変数変換したときの確率分布を考える。

与える分布は,$p(x)=2x\ \theta(x) \theta(1-x),\ q(y) = 2y\ \theta(y) \theta(1-y)\ $とする。したがって,$0 \le x \le 1,\ 0 \le y \le 1 \ $を満足する。このとき,$\displaystyle \int_0^1 \int_0^1 p(x) q(y)\ dx dy =1$
これは,$\displaystyle \int_0^1 dx \int_0^x dy\ p(x) q(y) + \int_0^1 dx \int_x^1 dy\ p(x) q(y) =1$とも書ける。


(1)$\underline{X = x+y,\ Y = x-y \quad (0 \le X \le 2,\ -1 \le Y \le 1)\ }$の場合
このとき,$x = (X+Y)/2, \ y=(X-Y)/2, \ J(X,Y)=\frac{1}{2}\ $
$p(x) q(y) = 4 p q = (X+Y)(X-Y) = X^2-Y^2$

積分領域は, $-X \le Y \le 2-X$ かつ $X \le Y \le X-2$
$f(X,Y)\ $の期待値は,$\langle f \rangle = \int_0^1 dX \int_{-X}^{X} f(X,Y) \frac{X^2-Y^2}{2} dY + \int_1^2 dX \int_{X-2}^{2-X} f(X,Y)  \frac{X^2-Y^2}{2} dY$

(1-1) $\langle 1 \rangle =  \int_0^1 dX \Big\lbrack X^2 Y - \frac{Y^3}{3} \Big\rbrack_{0}^{X} + \int_1^2 dX \Big\lbrack X^2 Y - \frac{Y^3}{3} \Big\rbrack_{0}^{2-X} $
$\quad = \int_0^1 \frac{2}{3}X^3 dX + \int_1^2 \frac{2}{3}(2-X)(X^2+2X-2) dX $
$\quad =  \Big\lbrack  \frac{1}{6}X^4 \Big\rbrack_0^1 +\Big\lbrack  -\frac{X^4}{6} + 2X^2-\frac{8}{3}X \Big\rbrack_1^2 = 1$

(1-2) $\langle |Y| \rangle =  \int_0^1 dX \int_0^X (X^2Y-Y^3) dY + \int_1^2 dX \int_0^{2-X}(X^2Y-Y^3) dY$
$\quad = \int_0^1 \Bigl(\frac{X^4}{2}-\frac{X^4}{4} \Bigr) dX +  \int_1^2 \Bigl\{ \frac{X^2(2-X)^2}{2} - \frac{(2-X)^4}{4} \Bigr\} dX $
$\quad =  \Big\lbrack \frac{1}{20}X^5 \Big\rbrack_0^1 +   \Big\lbrack  \frac{1}{20}X^5 -\frac{4}{3}X^3 + 4X^2 -4X  \Big\rbrack_1^2 = \dfrac{4}{15} \ $


(2)$\underline{X = x+y,\ Y = xy \quad (0 \le X \le 2,\ 0 \le Y \le 1)\ }$の場合
このとき,$x = (X \pm \sqrt{X^2-4Y})/2, \ y=(X \mp \sqrt{X^2-4Y} )/2, \ J(X,Y)=\frac{1}{\sqrt{X^2-4Y}}$
$p(x) q(y) = 4 p q = (X \pm \sqrt{X^2-4Y})(X \mp \sqrt{X^2-4Y} ) = 4Y$

積分領域は, $0 \le Y$ かつ $X-1 \le Y \le X^2/4$
$f(X,Y)\ $の期待値は,$x>y$と$y>x$の場合をそれぞれ加えることで,
$\langle f \rangle = 2 \int_0^1 dX \int_0^{X^2/4} f(X,Y) \frac{4Y}{\sqrt{X^2-4Y}} dY + 2 \int_1^2 dX \int_{X-1}^{X^2/4} f(X,Y) \frac{4Y}{\sqrt{X^2-4Y}} dY$

(2-1) $\langle 1 \rangle = 2 \int_0^1 dX \int_0^{X^2/4} \frac{4Y}{\sqrt{X^2-4Y}} dY + 2 \int_1^2 dX \int_{X-1}^{X^2/4}  \frac{4Y}{\sqrt{X^2-4Y}} dY$
$\quad = 2 \int_0^1 dX \Big\lbrack -\frac{X^2+2Y}{3} \sqrt{X^2-4Y}  \Big\rbrack_0^{X^2/4} + 2 \int_1^2 dX  \Big\lbrack -\frac{X^2+2Y}{3} \sqrt{X^2-4Y} \Big\rbrack_{X-1}^{X^2/4}$
$\quad = \int_0^1 \frac{2 }{3}X^3 dX + \int_1^2 \frac{2}{3}(2-X)(x^2+2X-2) dX$
$\quad =  \Big\lbrack \frac{1}{6}X^4 \Big\rbrack_0^1 +  \Big\lbrack  -\frac{1}{6}X^4+ 2 X^2-\frac{8}{3}X  \Big\rbrack_1^2= 1$

(2-2) $\langle \sqrt{X^2-4Y} \rangle = 2 \int_0^1 dX \int_0^{X^2/4}  4Y dY + 2 \int_1^2 dX \int_{X-1}^{X^2/4} 4Y dY$
$\quad = 2 \int_0^1 dX \Big\lbrack 2Y^2 \Big\rbrack_0^{X^2/4} + 2 \int_1^2 dX  \Big\lbrack 2Y^2 \Big\rbrack_{X-1}^{X^2/4}$
$\quad = \int_0^1 \frac{1}{4}X^4 dX + \int_1^2 \Bigl\{ \frac{1}{4}X^4-4(X-1)^2 \Bigr\} dX  $
$\quad = \Big\lbrack \frac{1}{20}X^5 \Big\rbrack_0^2 +  \Big\lbrack -\frac{4}{3}(X-1)^3  \Big\rbrack_1^2 = \dfrac{4}{15}$ 

2024年2月4日日曜日

一様分布と変数変換

確率変数の積からの続き

2つの確率変数の一様分布があったとき,これを変数変換したときの確率分布を考える。

与える分布は,$p(x)=\theta(x) \theta(1-x),\ q(y) = \theta(y) \theta(1-y)\ $とする。したがって,$0 \le x \le 1,\ 0 \le y \le 1 \ $を満足する。このとき,$\displaystyle \int_0^1 \int_0^1 p(x) q(y)\ dx dy =1$
これは,$\displaystyle \int_0^1 dx \int_0^x dy\ p(x) q(y) + \int_0^1 dx \int_x^1 dy\ p(x) q(y) =1$とも書ける。


(1)$\underline{X = x+y,\ Y = x-y \quad (0 \le X \le 2,\ -1 \le Y \le 1)\ }$の場合
このとき,$x = (X+Y)/2, \ y=(X-Y)/2, \ J(X,Y)=\frac{1}{2}\ $

積分領域は, $-X \le Y \le 2-X$ かつ $X \le Y \le X-2$
$f(X,Y)\ $の期待値は,$\langle f \rangle = \int_0^1 dX \int_{-X}^{X} f(X,Y) J(X,Y) dY + \int_1^2 dX \int_{X-2}^{2-X} f(X,Y) J(X,Y) dY$

(1-1) $\langle 1 \rangle =  \int_0^1 dX \Big\lbrack \frac{Y}{2}\Big\rbrack_{-X}^{X} + \int_1^2 dX \Big\lbrack \frac{Y}{2}\Big\rbrack_{X-2}^{2-X} = \int_0^1 X dX + \int_1^2 (2-X) dX = 1$

(1-2) $\langle |Y| \rangle =  \int_0^1 dX \Big\lbrack \frac{Y^2}{2} \Big\rbrack_0^{X} + \int_1^2 dX \Big\lbrack \frac{Y^2}{2} \Big\rbrack_0^{2-X} = \int_0^1 \frac{X^2}{2}dX +  \int_1^2 \frac{(2-X)^2}{2} dX $
$\quad =  \Big\lbrack \frac{X^3}{6} \Big\rbrack_0^1 +   \Big\lbrack \frac{(X-2)^3}{6} \Big\rbrack_1^2 = \dfrac{1}{3} \ $


(2)$\underline{X = x+y,\ Y = xy \quad (0 \le X \le 2,\ 0 \le Y \le 1)\ }$の場合
このとき,$x = (X \pm \sqrt{X^2-4Y})/2, \ y=(X \mp \sqrt{X^2-4Y} )/2, \ J(X,Y)=\frac{1}{\sqrt{X^2-4Y}}$

積分領域は, $0 \le Y$ かつ $X-1 \le Y \le X^2/4$
$f(X,Y)\ $の期待値は,$x>y$と$y>x$の場合をそれぞれ加えることで,
$\langle f \rangle = 2 \int_0^1 dX \int_0^{X^2/4} f(X,Y) J(X,Y) dY + 2 \int_1^2 dX \int_{X-1}^{X^2/4} f(X,Y) J(X,Y) dY$

(2-1) $\langle 1 \rangle = 2 \int_0^1 dX \int_0^{X^2/4} \frac{1}{\sqrt{X^2-4Y}} dY + 2 \int_1^2 dX \int_{X-1}^{X^2/4}  \frac{1}{\sqrt{X^2-4Y}} dY$
$\quad = 2 \int_0^1 dX \Big\lbrack -\frac{1}{2} \sqrt{X^2-4Y}  \Big\rbrack_0^{X^2/4} + 2 \int_1^2 dX  \Big\lbrack -\frac{1}{2} \sqrt{X^2-4Y} \Big\rbrack_{X-1}^{X^2/4}$
$\quad = \int_0^1 X dX + \int_1^2 (X-1) dX  = 1$

(2-2) $\langle \sqrt{X^2-4Y} \rangle = 2 \int_0^1 dX \int_0^{X^2/4}  dY + 2 \int_1^2 dX \int_{X-1}^{X^2/4}  dY$
$\quad = 2 \int_0^1 dX \Big\lbrack Y \Big\rbrack_0^{X^2/4} + 2 \int_1^2 dX  \Big\lbrack Y \Big\rbrack_{X-1}^{X^2/4}$
$\quad = \int_0^1 \frac{X^2}{2} dX + \int_1^2 (\frac{X^2}{2}-2X + 2) dX  = \Big\lbrack \frac{X^3}{6} \Big\rbrack_0^2 +  \Big\lbrack -X^2 + 2X \Big\rbrack_1^2 = \dfrac{1}{3}$ 

2024年2月3日土曜日

円の長さ

正方形の長さからの続き

ある図形の大きさの指標となる長さを,図形内に一様分布する2点の距離の期待値として定義することで,都道府県の形や大きさを,面積や周長だけでなく"長さ"で特徴づけるという話をしている。

練習として,正方形内のランダムな2点の距離の期待値$\ d \ $を計算できることを確認した。次にトライするのが円であるが,ネットで検索しても生成AIにきいてもあまり適切な解答が得られない。一番近いのが,Yahoo知恵袋の「半径1の円内の任意の2点間の距離の期待値は?」だ。これも結局解析的な答えがでなくて,数値計算で $\ d=0.9054\ $という値を出している。

$d = \int_0^{2\pi} (\frac{1}{2\pi}) d\theta_1 \int_0^{2\pi} ( \frac{1}{2\pi}) d\theta_2 \int_0^1 (2 r_1) dr_1 \int_0^1 (2 r_2) dr_2 \sqrt{r_1^2+r_2^2-2r_1 r_2 \cos(\theta_1-\theta_2) }$
ただし,各積分の()内がそれぞれの変数に対応する確率密度関数,$p(\theta_1)$,$ p(\theta_2)$,$q(r_1)$, $q(r_2)\ $であり,それぞれの変数で積分すると1になるように規格化されている。

これを計算するためには,与えられた4変数の確率分布関数から変数変換によって,積分可能な形に持ち込む必要があるが,なかなか難渋する。しかたがないので,とりあえずJuliaとMathematicaで数値計算してみる。
a=zeros(Float64,1000001,2)

function ju(a,n)
  k = 0
  for i in 1:n
      x = 2*rand()-1
      y = 2*rand()-1
      if x^2+y^2 < 1.0
        k = k + 1
        a[k,1] = x
        a[k,2] = y
      end
  end 
  return k
end

function su(a,n)
  m = div(n,4)*3
  sum = 0
  for i = 1:m
    for j = i:m
      sum = sum + sqrt((a[i,1]-a[j,1])^2+(a[i,2]-a[j,2])^2)
    end
  end
  return sum/binomial(m,2)
end

n=300000

@time su(a,n)
235821
24.223278 seconds (7.18 k allocations: 500.875 KiB, 0.07% compilation time)
0.9055523

生成AIの2つが答えた,$\dfrac{4}{\pi}=1.27324$はたぶん誤っていたということだろう。

解析的に計算できないかと思うのだが,角度積分が完全楕円積分の形になるので,これをさらに積分するのはちょっと難しそうだった。角度積分を後回しにしてもさらに面倒か。Mathematicaに投げてみたが,忍耐可能時間内には答えが出なかった。

2024年2月2日金曜日

確率変数の積

将来必要になりそうな,確率変数の積の確率分布関数を求める。

2つの確率変数$X$と$Y$が確率密度分布関数$p(x),\ q(y)$に対応している。このとき,確率変数$Z=X*Y$はどのような確率分布をするか,再び,緑川章一さんのノートで勉強する。

確率変数 $Z=X*Y$の確率分布関数を $r(z)$とすると,$r(z) = \int_0^1  \int_0^1  p(x) q(y) \delta(z- x*y) \ dx\ dy =  \int_{0}^{1} \dfrac{1}{|y|} p(z/y) q(y) \ dy $となる。ここでデルタ関数の性質,$\delta(a x) = \delta(x)/|a|$を用いた。この$\ z \ $の範囲は,$ 0 < z < \infty$ である

(1)$X$と$Y$が,それぞれ一様分布,$p(x)  =  1 \ (0 \le x \le 1)$ ,$q(y)  =  1 \ (0 \le y \le 1)$を満足している場合。ここで,$0< z/y<1\ $より,$z<y<1$である。したがって,

$r(z) = \int_z^1 \frac{1}{y} 1*1 \ dy= -\log z$

(2)$X$と$Y$が,それぞれ三角分布,$p(x) = 2x \ (0 \le x \le 1)$,$q(y) = 2y \ (0 \le y \le 1)$をしている場合(単位円内の点の一様分布の動径変数)。

$r(z) = \int_z^1 \dfrac{1}{y} \dfrac{2z}{y} (2y)\ dy = \int_z^1 \dfrac{4z}{y} \ dy = - 4z \log z$

うーん,あんまりうれしくないかもしれない。後々$\log$の計算が残るので。

2024年2月1日木曜日

曖昧な弱者

1月30日の日経朝刊の経済教室の伊藤昌亮(1961-)の記事が目を引いた。

「弱さ」を競い合う社会 「曖昧な弱者」存在認識をという表題である。



図:今日の左右対立の構造(伊藤昌亮 日本経済新聞から引用して改変)

日本経済新聞に掲載された伊藤の図を引用するが,自分の理解を深めるために若干修正している。一番気になっているのは,マスメディアは政治経済エリート側に包摂されてしまっているのではないかということ,リベラル・保守,左派・右派の従来の定義とスコープが機能しているのかということであり,それぞれ?を付けている。

オカケンさん[1]の助けを借りて,伊藤昌亮の論説を解読すると次のようになる。

左派やリベラル派にとって明白な弱者とは,搾取された労働者や貧困化の女性・若者であり,アイデンティティポリティックスの対象とされる,在日外国人,被差別部落,沖縄・アイヌ,障害者,LGBTQなどである。文化エリートはこれらとの連帯を強く主張する。

一方で,OECD諸国の中でも著しく「小さい」日本政府(OECD諸国最低レベルの社会福祉費と教育費)は,その福祉・教育機能を,企業や家庭に投げてきたが,グローバリズムの嵐の中でそのシステムは崩壊し,いわゆる中流階級は消滅して,激しい二極分化が生じた。

この結果,大量に生じているのが,従来の明白な弱者カテゴリーでは十分にすくい上げられない,曖昧な弱者である。社会的にはっきりと認知・共感されない彼らは,そのフラストレーションを,「あいつらだけ認知されるのはずるい(在日特権言説,生活保護・高齢者バッシング)」と明白な弱者に対して牙をむく。

それは,アメリカのトランプ現象やヨーロッパの移民排斥右翼の台頭とまさに軌を一にする動きになっている。こうして,ネットワーク上には,政治的な意図を持ってDAPPIなどが着火すれば容易に燃え上がるネトウヨ的な素地が醸成されてきたのだ。

ただ,これらに保守・右派というレッテルを貼ってよいかどうかは疑問だ。たしかに,リベラル勢力に対抗するためだけに,宗教右派は明白な弱者たたきを繰り返しているが,日本維新の会に代表されるようなネオリベラリズムは保守とはいえない。むしろ,既存秩序を崩壊させる中で,新しい権益を掠め取ろうという作戦に立っているので話は複雑だ。



2024年1月31日水曜日

三体

昨年12月に録画してあったWOWOWのテレビSFドラマ「三体」(全30話)をようやく視聴した。霊河影視制作(上海)有限公司の作品だ。

原作の三体(第一部)の著者は劉慈欽(1963-)であり,2015年にアジア作家の作品で初めてヒューゴー賞(1953-)の長編小説部門に輝いた。第三部まで出版翻訳されていて,読みたい本リストのトップレベルに置いてある。

最初の10回の前半は,それほどでもなかったけれど,三体のVRゲームのイメージや文化大革命後の紅岸基地のあたりから急に面白くなる。ドラマ三体は,中国の配信プラットホームであるテンセントの作品なのだが,開放改革以前の中華人民共和国の様子をあれくらい描写していてもOKなのか。紅岸基地での物理学的な謎解きやサスペンスの部分がよかった。SFXも素晴らしい。今の日本だとせいぜいがゴジラであって,これほどの骨太の作品はちょっと無理だ。

Netflix版の三体(三部作を予定か)も近々公開されるはずだけれど,どんなものなのだろうか。予告編を見たが,これはこれでいいけれど,やはりいつものアメリカナイズされすぎた映像と世界観が少し鼻につく。中国版の方が主人公もいいし新鮮な感じがする。




2024年1月30日火曜日

Apple ID

朝起きていつものようにMacbook Airを立ち上げると,何だかエラーが出ている。

iCloudにアクセスできないとかなんとか。パスワードを入力してもその先に進めない。パスワードの変更もできない。困った。とても困った。Appleのサポートページには,「iCloud に接続またはサインインできない場合」には丁寧な説明がある。が,そこからIDを入れて,CAPTCHAを通ったのに,そこではねられてしまう。どうやらApple IDがロックされていることに気がついた。

思い当たる節がある。いよいよ非常勤講師も最後なので,金曜日に大学のMacbookAirの掃除をしていた。iCloudにログインしたままだったので,ログアウトしようとした。OSが古くてきびきび動かないのでパスワードも何度か間違えてしまう。授業が始まりそうになったので,途中で作業を中断したまま放置してきた。どうもこれがあやしい。

仕方がないので,大学でもう一度状況を確認してから,心斎橋のアップルストアに行こうと考えた。大学のMacbookAirの方は再起動して簡単にiCloudからログアウトできた。もちろんこれだけでは,AppleIDのロックは解除されない。このため,自分のMacBook Air 2000だけでなく,iPhoneもiPadも,メールは届かないし,アップルストアの予約も出来ないし,ICOCAのチャージも出来ない

心斎橋アップルストアの予約のために電話をしてみた。これまた繋がるまでに15分以上待たされた。忍耐力あるもの達だけが通過できるシステムだ。最初に,アンケート協力用の電話番号を入れさせられ,さらに待ち時間用音楽のジャンルを選択するのだが,そんなサービスはいいからとにかく速く対応してほしい。

結局,心斎橋まで行くことはなくて電話だけでロック解除してもらえた。ただし作業は24時間以内なのでしばらく待たなければならない。月にいるSLIMの電源が復活しているかどうかは,地球から電波を送ってそのレスポンスを見る必要があるのだが,気分はこれと同じだ。30分,1時間,3時間では復活していなかったが,6時間でようやくApple IDが復活してほっとした。めでたし,めでたし。

写真:ようやくここまでたどり着いて原因がわかった地点。

2024年1月29日月曜日

SLIM(2)

SLIM(1)からの続き

SLIMは高度5mあたりで,当初予定(高度1.8m)されたように,LEV-1LEV-2(SORA-Q)という無人探査ロボットシステムを放出した。

LEV-1(2.1kg,26 × 40 × 60 cm)は,2.4m/sで月面に落ちる予定だったが,もし5mからの自由落下ならば,鉛直方向の速度は4m/sになる。自律的に跳躍移動しながら方位制御して地球との直接通信を確立する。

LEV-2(0.25kg,直径8cmの球状から変形可)は,月面のレゴリス上を移動して動作ログを保存し,着陸機SLIMの周辺を撮影して,画像データと動作ログをLEV-1経由で地球に送信する。その結果がSLIM(1)で示した画像である。途中にブロックノイズが入っていることや,解像度はもっと出ているのだが,通信上の制約で落として送信したらしい。

LEV-2別名SORA-Qは,タカラトミーとの共同開発である。この月面に行ったSORA-Qと同じサイズ,同じ変形機能,同じ走行機能・撮影カメラを持った 1/1スケールモデルを2万7千円で販売している。わぉ!思わず注文してしまいそうになる。


写真:SORA-Q実寸モデル(タカラトミーから引用)

2024年1月28日日曜日

SLIM(1)

SLIMは,JAXAの小型月着陸実証機だ。月周回衛星かぐやのデータと照らし合わせながら,月面の目的地にピンポイントで着陸し,無人小型ロボットシステムで月を探査しようというものだ。

2023年9月7日に種子島宇宙センターからH-IIAロケットで打ち上げられ,12月25日には近月点600km,遠月点2000kmの月周回軌道に投入された。1月20日に高度15kmから降下を開始し,世界で5番目の月着陸に成功した。太陽電池からの電源供給ができていないというニュースを聞いたとき,あーこれはちょっと残念かと思った。

1月25日にJAXAによる記者会見が行われた。そこで紹介されていた写真が次のものだ。


写真:LEV-2が撮影したSLIM本体が転倒している様子(JAXAから引用)

SLIMは航法カメラで月面を撮影しながら,自律的な航法誘導制御を行っている。50m上空に至ったところで,障害物を避けるモードに移行する。目標地点は,経度25.2°,緯度-13.3°の2つのクレーターの境界上の斜面である。

月面上50mまでは予定通り順調に飛行していて,この段階でのピンポイント着陸精度は3-4m程度と考えられる。従来の数キロメートルに比べて1/1000の精度である。ところが,2基搭載している500Nのメインエンジ


の1基が脱落して,推力が半分になってしまった。それでも,最終的に秒速1.4mと想定範囲よりゆっくり着陸することになる。

この異常により,横方向の速度が発生してしまい,結果的に55m東にズレた点に接地する。このため,2段階の受け身型の着陸条件を満足できずに機体はそのまま斜面に着陸した。この結果,メインエンジが上向きの転倒状態で静止した。太陽電池は正常な上向きではなく西向きになり,太陽光が当たらないため電源供給ができなくなった。


月の一日は約30日であり,今月の上弦が1月18日,満月が1月26日,下弦は2月3日である。着陸したのは1月20日で経度20°(1.5日相当)だから,下図において地球と緑の線で結ぶ位置の月面上に着地している。これは,月の一日でいえば午前9時に相当する。1月24日ごろにSLIM着地点は正午を迎え,2月1日には日没となる。月面上のSLIMに西日が差して,太陽電池が復活する可能性があるのは,2月1日までということか。



図:SLIMと太陽の位置関係(実際には経度20°の緑線まで回転)

追伸:1月29日に,ようやく西日で電源が復活したようだ。新しい画像も撮影している。

2024年1月27日土曜日

月の一日

2024年1月26日,今日は満月だ。非常勤の授業が終って平端の駅の1番ホームへ向かう午後3時半すぎ,乗客のかたまりがどんどんホームから階段を下りてきて引きも切らない。そうか,今日は天理教の春季大祭の日だった。夜は冷え込んでいるけれど空は曇っていて月は見えない。

SLIMの記事を書くためには,月の一日について調べておかなければならない。

(1)月の公転周期
ケプラーの第三法則というのか,ニュートンの運動方程式を解けば,公転周期は$T=\dfrac{2\pi}{\sqrt{GM}} a ^{3/2}$である。$G$は万有引力定数,$M$は地球の質量,$a$は月の軌道の長半径である。$\sqrt{GM}=g R$であり,重力加速度$\ g=9.82 {\rm m/s}^2$,地球半径$\ R= 6.37 \times 10^6 {\rm m}$を使えばよい。$a\ $の値は遠地点と近地点の平均値であり,$a=3.83\times 10^8 {\rm m}$。これらから,$T=27.3$日となる。

(2)月の一日の長さ
潮汐作用の結果,月の自転周期と公転周期$\ T\ $は一致し,地球から見える月は常に同じ面になる。月の周期の間に地球が太陽の回りを公転するため,月の一日,例えば日の出から次の日の出のまでの時間$\ t\ $は$\ T\ $ではなく,それよりも長くなる。地球の公転角速度を$\ \Omega$,月の公転角速度を$ \omega$ とすると,$\Omega t = \omega (t -T) $が成り立つ。これから $t=T/(1-\Omega/\omega) = 27.3/(1-27.3/365) = 29.5 $日が得られる。これを朔望月という。



図:月の一日(朔望月)

2024年1月26日金曜日

正方形の長さ

都道府県の長さからの続き

正方形の領域$\ (x, y), \ 0 \le x \le 1$,$0 \le y \le 1\ $を考えて,この中の2点を$\ (x_1,y_1),\ (x_2, y_2)\ $とする。これらの座標が$\ p_0(z)=1\ (0 \le z \le 1),\  =0\ (z <0,\  1 < z)\ $で一様分布している。

このとき,確率変数の和と差の説明により,$x=x_1-x_2\ $と$\ y=y_1-y_2\ $は,$p(z)=1+z\ (-1 \le z \le 0),\  =1-z\ (0 \le z \le 1)\ $という確率分布になる。また,$X=(x_1+x_2)/2\ $と$Y=(y_1+y_2)/2\ $の確率分布は,$q(z)=z\ (0 \le z \le 1),\  =2-z\ (1 \le z \le 2)\ $となる。

そこで,2点の期待値は,$d=\int \int \int \int \sqrt{(x_1-x_2)^2+(y_1-y_2)^2} \ p_0(x_1) p_0(x_2) p_0(y_1) p_0(y_2) \ dx_1 dx_2 dy_1 dy_2 $
$\quad = \int \int \int \int \sqrt{x^2+y^2}\  p(x) p(y) q(X) q(Y) \ dx dy dX dY$
$\quad = \int \int \sqrt{x^2+y^2} \ p(x) p(y) \ dx dy = 4 \int_0^1 \int_0^1 (1-x)(1-y) \sqrt{x^2+y^2} \ dx dy$

ここで,$y = x \sinh z \ $と変数変換して,$y\ $の積分すなわち$z\ $での積分を先に行う。このとき,$y: 0\rightarrow 1\ $より,$z:0 \rightarrow \sinh^{-1}(1/x) = z_x\ $ $\bigl( \cosh z_x = \sqrt{1 + (1/x)^2} \ \bigr)$ であり,$\sqrt{x^2+y^2}= x \cosh x\ $と$\ dy = x\ \cosh z\ dz\ $が成り立つ。

$f(x) = \int_0^1 (1-y) \sqrt{x^2+y^2} dy = \int_0^{z_x} (1 - x \sinh z ) \cdot x \cosh z \cdot x \cosh z\  dz$
$\displaystyle \quad = \frac{x^2}{2} \int_0^{z_x} (1 + \cosh 2z )\ dz -\frac{x^3}{3} \Bigl[ \cosh^3 z \Bigr]_0^{z_x}$
$\displaystyle \quad =  \frac{x^2}{2}  \Bigl( \sinh^{-1}(1/x) + \sinh z_x \cosh z_x \Bigr) -\frac{x^3}{3} \Bigl( \cosh^3 z _x -1 \Bigr)$
$\displaystyle \quad = \frac{x^2}{2} \sinh^{-1}(1/x) + \frac{1}{2} \sqrt{1+x^2} +\frac{1}{3} x^3 - \frac{1}{3} (1+x^2)^{3/2}$

次に,これに$(1-x)$をかけて,$x$で積分してから4倍すれば$d$が求まる。
$\displaystyle d= 4\int_0^1 (1-x) \Bigl\{ \frac{x^2}{2} \sinh^{-1}(1/x) + \frac{1}{2} \sqrt{1+x^2} +\frac{1}{3} x^3 - \frac{1}{3} (1+x^2)^{3/2} \Bigr\} dx$

$g_1(x)=4 \int (1-x) \frac{x^2}{2} \sinh^{-1}(1/x) \ dx $
$\quad = \frac{1}{6}(2+2x-x^2)\sqrt{1+x^2} +\frac{1}{6}(-2+4x^4-3x^4)\sinh^{-1}x$
$g_2(x)=4 \int (1-x) \frac{1}{2} \sqrt{1+x^2} \ dx = -\frac{1}{3} (2-3x+2x^2) \sqrt{1+x^2} + \sinh^{-1}x$
$g_3(x)=4 \int (1-x) \frac{1}{3} x^3 \ dx = \frac{1}{3} x^4 -\frac{4}{15}x^5$
$g_4(x)=4 \int (1-x) \frac{-1}{3} (1+x^2)\ ^{3/2} \ dx $
$\quad = \frac{1}{30} (8-25x+16x^2-10x^3+8x^4) \sqrt{1+x^2} -\frac{1}{2} \sinh^{-1}x$

$\therefore g(x) = g_1(x)+g_2(x)+g_3(x)+g_4(x) = \frac{1}{15}(5 x^4 -4x^5) +$
$\quad  \frac{1}{30}(8x^4-10x^3-9x^2+15x-2)\sqrt{1+x^2} +\frac{1}{6}(-3x^4+4x^3+1)\sinh^{-1}x$

これから,$d=g(1)-g(0)=\frac{1}{15}\Bigl\{2+\sqrt{2}+5 \log(1+\sqrt{2}) \Bigr\}= 0.521405\ $が得られた。


2024年1月25日木曜日

確率変数の和と差

都道府県の長さからの続き

2つの確率変数$X$と$Y$がある。それぞれはある確率密度分布関数$p(x),\ q(y)$に対応している。このとき,確率変数$X \pm Y$はどのような確率分布をするかという問題を考えたい。

これについては,緑川章一さんのノートが参考になった。やはり専門の近い物理屋さんが書いたものは読みやくて助かる。これをまとめてみる。

$X$と$Y$が,それぞれ一様分布,$p(x)  =  1 \ (0 \le x \le 1)$ ,$q(y)  =  1 \ (0 \le y \le 1)$をしている。このとき,確率変数$Z$を $Z=X \pm Y$として,その確率分布関数の $r_{\pm}(z)$を求める。これは,$r_{\pm}(z) = \int  \int  p(x) q(y ) \delta(z-(x \pm y)) \ dx\ dy =  \int_{0}^{1} p(z \mp y) q(y) \ dy $となる。
なお,この$\ z \ $の範囲は,$r_{+}(z) \rightarrow 0 < z < 2$,$r_{-} \rightarrow -1 < z < 1 $ である

$\therefore \ r_{+}(z) \ \rightarrow \ ( 0 \le y \le 1 \ \&\& \  z-1 \le y \le z )$,つまり,
$r_{+}(z) = z\ (0 < z < 1),r_{+}(z) = 2-z \ (1 < z < 2)$
$\therefore \ r_{-}(z) \ \rightarrow \ ( 0 \le y \le 1 \ \&\& \  -z \le y \le 1-z )$,つまり,
$r_{-}(z) = 1+z \ (-1< z <0),r_{-}(z) = 1-z \ (0 < z < 1)$

図:確率分布関数の範囲

2024年1月24日水曜日

積分漸化式

積分(3)からの続き

三角関数の積分の漸化式が教科書に載っていたことをいまごろ思い出した。これならば双曲線関数にも簡単にあてはめられるはずだ。

$I_n = \int \sin^n x \ dx =  (-\cos x) \sin^{n-1} x - \int  (-\cos x) (n-1) \sin^{n-2} x \cos x\ dx$
$\quad = (-\cos x) \sin^{n-1} x + (n-1) \int (1-\sin^2 x) \sin^{n-2} x \ dx$
$\quad = (-\cos x) \sin^{n-1} x + (n-1) (I_{n-2} - I_n )$
$\therefore I_n = -\frac{1}{n} \cos x \sin^{n-1} x + \frac{n-1}{n} I_{n-2} \quad (n \ge 2) $

$I_n = \int \cos^n x \ dx =  (\sin x) \cos^{n-1} x - \int  (\sin x) (n-1) \cos^{n-2} x (-\sin x) \ dx$
$\quad = (\sin x) \cos^{n-1} x + (n-1) \int (1-\cos^2 x) \cos^{n-2} x \ dx$
$\quad = (\sin x) \cos^{n-1} x + (n-1) (I_{n-2} - I_n )$
$\therefore I_n = \frac{1}{n} \sin x \cos^{n-1} x + \frac{n-1}{n} I_{n-2} \quad (n \ge 2) $

$I_n = \int \tan^n x \ dx =   \int  (\frac{1}{\cos^2 x} -1) \tan^{n-2} x \ dx$
$\quad = \int (\tan x)' \tan^{n-2} x\ dx -  \int \tan^{n-2} x \ dx$
$\therefore I_n = \frac{1}{n-1}\tan^{n-1} x - I_{n-2}  \quad (n \ge 2) $

$I_n = \int \sinh^n x \ dx =  (\cosh x) \sinh^{n-1} x - \int  (\cosh x) (n-1) \sinh^{n-2} x \cosh x\ dx$
$\quad = \cosh x\ \sinh^{n-1} x - (n-1) \int (1+\sinh^2 x) \sinh^{n-2} x \ dx$
$\quad = \cosh x\ \sinh^{n-1} x - (n-1) (I_{n-2} + I_n )$
$\therefore I_n = \frac{1}{n} \cosh x\ \sinh^{n-1} x - \frac{n-1}{n} I_{n-2} \quad (n \ge 2) $

$I_n = \int \cosh^n x \ dx =  (\sinh x) \cosh^{n-1} x - \int  (\sinh x) (n-1) \cosh^{n-2} x \sinh x \ dx$
$\quad = \sinh x \cosh^{n-1} x - (n-1) \int (\cosh^2 x -1) \cosh^{n-2} x \ dx$
$\quad = \sinh x \cosh^{n-1} x + (n-1) (I_{n-2} - I_n )$
$\therefore I_n = \frac{1}{n} \sinh x \cosh^{n-1} x + \frac{n-1}{n} I_{n-2} \quad (n \ge 2) $

$I_n = \int \tanh^n x \ dx =   \int  (1-\frac{1}{\cosh^2 x} ) \tanh^{n-2} x \ dx$
$\quad = -\int (\tanh x)' \tanh^{n-2} x\ dx + \int \tanh^{n-2} x \ dx$
$\therefore I_n = -\frac{1}{n-1}\tanh^{n-1} x + I_{n-2}  \quad (n \ge 2) $

2024年1月23日火曜日

藤岡作太郎

就寝中にトイレに行きたくなるとき,眠りが浅くなって夢を見る。いや,半分覚醒してまどろんでいる状態なので夢ではないのかもしれない。こうした夢と覚醒がシュレーディンガーの猫のようになって区別しにくい時間がしばしば訪れる。

昨晩のその時間は,「鈴木大拙」についての説明を誰かに一生懸命しようとしていた。ただ,名前が思い出せないのである。えーっと,金沢出身で,西田幾多郎と友達で,いるでしょう,禅の研究で(善の研究ではない)海外に名を馳せた,誰だったか,ほらあの(静かな水面のある落ち着いた記念館のイメージを想起しつつ),えーっと,三太郎とよばれていたから,本名は○太郎のはずだけれど,それではわからないし・・・

そうこうしているうちに目が覚めてトイレに行ったが名前の記憶はオフのまま。再び布団に潜ってもまだ思い出せない。そのまま眠りに入ると,明け方近くの夢の中でようやく思い出すことができた。あ,鈴木大拙だ!朝起きても思い出した量子状態は崩壊することなく維持されていた。


三太郎というのは誰だったろうかと,Wikipediaで鈴木大拙=貞太郎(1870-1966)を調べてみると「同郷の西田幾多郎(1870-1945)、藤岡作太郎(1870-1910)とは石川県専門学校( 1881- 第四高等中学校 1887-)以来の友人であり、鈴木、西田、藤岡の三人は加賀の三太郎と称された」とあった。

藤岡作太郎はどんな人かとさらに調べると,日本で最初の文学博士,国文学(国文学全史平安朝篇)の人だった。その長男が物理学者で物理教育学会の会長も務めた藤岡由夫(1903-1976),孫がレーザ工学の藤岡知夫(1935-2022),ひ孫がテレビでおなじみの指揮者の藤岡幸夫(1962-)だった。

藤岡作太郎の長女の綾が,長男の藤岡由夫の友人の中谷宇吉郎(1900-1962)と結婚しているが若くして亡くなっている。孫の藤岡知夫の妻は原子物理学の菊池正士(1902-1974)の長女であり,これをたどると箕作家(みつくりけ)を通じて初代阪大総長の長岡半太郎(1865-1950)までつながる。なお,長岡半太郎と本多光太郎(1870-1954)と鈴木梅太郎(1874-1943)は理研の三太郎だ。

2024年1月22日月曜日

積分(3)

都道府県の長さからの続き

次の積分 $I_n=\int x^n \sqrt{x^2 + y^2} \ dx $  が必要なのであった。
そこで,$x=y\ \sinh z$と変数変換して,$dx = y \cosh z\ dz $と$\sqrt{x^2+y^2} = y \cosh z$
から,$I_n = y^{n+2} \int \sinh^n z\ \cosh^2 z\ dz$となる。後で必要になるものとして,$J_n = \int \sinh^n z \ dz $を定義しておく。

$I_0 = y^2 \int \cosh^2 z\  dz = \frac{1}{2} y^2  \int (1 + \cosh 2z) \ dz $
$=  \frac{1}{2} y^2 (z + \frac{1}{2} \sinh 2 z) =  \frac{1}{2} y^2  \sinh^{-1}(x/y) + \frac{1}{2}x \sqrt{x^2+y^2} $

$I_1 = y^3 \int  \sinh z\ \cosh^2 z\ dz =  y^3 \int t^2 dt = \frac{1}{3}\bigl( \sqrt{x^2+y^2}\bigr)^3$

$I_2 = y^4 \int \sinh^2 z\ \cosh^2 z \ dz = y^4 \int (\sinh^2 z + \sinh^4 z )\ dz = y^4 (J_2 + J_4)$

$I_3 = y^5 \int \sinh^3 z\ \cosh^2 z\ dz = y^5 \int (\sinh^3 z + \sinh^5 z )\ dz = y^5 (J_3 + J_5)$

などとなる。

$J_2 = \int \sinh^2 z \ dz = \frac{1}{2} \int (\cosh 2z -1) \ dz =  \frac{1}{4} \sinh 2z - \frac{1}{2} z$

$J_3 = \int \sinh^3 z \ dz = \int (\cosh^2 z -1) \sinh z \ dz = \frac{1}{3} \cosh^3 z - \cosh z$

$J_4 =  \int \sinh^4 z \ dz = \frac{1}{4} \int (\cosh 2z -1)^2 \ dz = \int \bigl( \frac{3}{8}-\frac{1}{2}\cosh 2z + \frac{1}{2} \cosh 4z  \bigr) \ dz$
$\quad = \frac{1}{8} \sinh 4z -\frac{1}{4} \sinh 2 z + \frac{3}{8}z $

結局$J_n$がシステマティックに計算できればよいということか。続く。

2024年1月21日日曜日

双曲線関数

都道府県の長さからの続き

一様分布の確率密度関数で正方形の内部のランダムな2点の平均距離を求める際に,面倒な積分が必要になる。このとき双曲線関数への変数変換を行うのだが,久しぶりに使うと勘が鈍っていてなかなか計算が進まない。ので,復習する。

$\sinh x = \dfrac{e^x - e^{-x}}{2},\ \  \cosh x = \dfrac{e^x + e^{-x}}{2},\ \  \tanh x = \dfrac{\sinh x}{\cosh x} = \dfrac{e^x - e^{-x}}{e^x + e^{-x}}$
$\cosh^2 x - \sinh^2 x = 1, \ \ \tanh^2 x = 1 - \dfrac{1}{\cosh^2 x},\ \ \dfrac{1}{\tanh^2 x} = 1 +  \dfrac{1}{\sinh^2 x}$
$\frac{d}{dx}\sinh x = \cosh x,\ \ \frac{d}{dx}\cosh x = \sinh x, \ \  \frac{d}{dx} \tanh x = \dfrac{1}{\cosh^2 x}$
$\int \sinh x \ dx= \cosh x,\ \ \int \cosh x \ dx = \sinh x, \ \  \int \tanh x\ dx = \log( \cosh x)$


$\sinh ( x \pm y )= \sinh x \cosh y \pm \cosh x \sinh y$
$\cosh ( x \pm y )= \cosh x \cosh y \pm \sinh x \sinh y$
$\tanh ( x \pm y )= \dfrac{\tanh x \pm \tanh y}{1 \pm \tanh x \tanh y}$

$\sinh 2x  = 2 \sinh x \cosh x = 2 \sinh x \sqrt{1 + \sinh^2 x}$
$\cosh 2x = 2 \cosh^2 x - 1 = 2 \sinh^2 x + 1$

$\sinh 3x  = \sinh^3 x + 3 \sinh x \cosh^2 x$
$\cosh 3x = \cosh^3 x + 3 \cosh x \sinh^2 x$

$\sinh 4x  = 4 \sinh^3 x \cosh x + 4 \sinh x \cosh^3 x$
$\cosh 4x =  \sinh^4 x + 6  \sinh^2 x +\cosh^2 x + \cosh^4 x$

$\sinh^{-1}x = \log ( x + \sqrt{x^2+1} ) = -\log(\sqrt{x^2+1} - x)$
$\cosh^{-1} x = \log (x + \sqrt{x^2-1}) = \log(x - \sqrt{x^2-1})$
$\tanh^{-1} x = \dfrac{1}{2} \log{\dfrac{x+1}{x-1}}$

$\frac{d}{dx}\sinh^{-1} x = \dfrac{1}{\sqrt{x^2+1}},\ \ \frac{d}{dx}\cosh^{-1} x = \dfrac{1}{\sqrt{x^2-1}}, \ \  \frac{d}{dx} \tanh^{-1} x = \dfrac{1}{1-x^2}$

$\int \sinh^{-1} x \ dx= x \sinh^{-1} x - \sqrt{x^2+1}$
$\int \cosh^{-1} x \ dx = x \cosh^{-1} x  - \sqrt{x^2-1}$
$\int \tanh^{-1} x\ dx = x \tanh^{-1} x + \frac{1}{2}\log(1-x^2)$



図:双曲線関数の定義

2024年1月20日土曜日

都道府県の長さ

ACTIVE GALATTICさんが,次のように書いていた。
中核都市の金沢と被害が激しい輪島は東京と水戸くらい離れていると石川県の広さを説明する投稿を先日見かけたが、石川県の広さは34~35位とランキング下位なので、むしろ東京や大阪が都道府県としては異様に狭いと認識すべきなのだろう
面積と距離が微妙に絡まっていて,ややモニョる。ある図形が与えられたときに,その面積ではなくて,長さを図形の大きさの指標として与えることはできるだろうか。地図上の領域だと東西南北の端点を考えて,南北間距離,東西間距離を抽出できるが,これでは2つの数字の組になる。あるいは図形の周長も考えられるが,フラクタル的な図形だと面積が等しくても周長はいくらでも長くできる。

そこで考えたのが次のような指標だ。図形を格子状の正方形に分割する。これらの正方形の全ての対を考え,各対の正方形の中心間の距離の平均値を求める。分割する正方形のサイズを小さくした極限として,一つの長さが定まるのではないか。(方法1)

検索してみると,もっとスマートな方法があった。図形内の点をランダムに選ぶ。その際に縦横方向にそれぞれ座標値が一様分布であるとした確率分布を仮定する。こうして選んだ2点間の距離の期待値を求めればよいというものだ。(方法2)

まず,方法1のプログラムをJuliaで作成して計算してみた。
a=zeros(Int64,1000001,2)

function ju(a,m,n)
  k = 0
  for i in 1:m
    for j in 1:n
      k = k + 1
      a[k,1] = i
      a[k,2] = j
    end
  end
end

function su(a,m,n)
  mn = m*n
  sum = 0
  for i = 1:mn
    for j = i:mn
      sum = sum + sqrt((a[i,1]-a[j,1])^2+(a[i,2]-a[j,2])^2)
    end
  end
  return sum/binomial(mn,2)
end

m=800
n=800
ju(a,m,n)
@time su(a,m,n)/m
n=   25    time= 0.016720 0.521829
n=   50 time= 0.018188 0.521513
n= 100 time= 0.065998 0.521433
n= 200 time= 0.788281 0.521412
n= 400 time= 12.1924   0.521407
n= 800 time= 198.615   0.521405

nは正方形の1辺の格子点の数,timeは計算時間(秒)であり,最後が計算結果の平均長である。
方法2による理論値が解析的に得られており,( 2+√2+5 log(√2 + 1) )/15 = 0.521405 なので,nを増やすとともに理論値に収束し,n=800では理論値に一致している。

m = 50; n = 50;
a = Flatten[Table[{i, j}, {i, 1, m}, {j, 1, n}], 1];
Timing[Sum[
     Sqrt[(a[[i, 1]] - a[[j, 1]])^2 + (a[[i, 2]] - a[[j, 2]])^2], 
   {i, 1, n*n}, {j, i, n*n}]/Binomial[n*n, 2]/n // N]

Mathematicaでは,若干コードが簡単になるけれど,時間は n=50 で 12.4801秒と,600倍以上かかってしまった。

(付)このコードはそのまま乱数による計算に置き換えることもできる。Juliaでは,a[k,1] = rand(),a[k,2] = rand() として最後の結果の/nを除く。Mathematicaでは,Tableの中身を{i,j}から{RandomReal[1], RandomReal[1]}に置き換えて最後の結果の/nを除く。

[1]地域内距離(腰塚武志)

2024年1月19日金曜日

トリチウム(3)

トリチウム(2)からの続き

非常勤で担当している物理科学概説の授業も,後3回を残すばかりになった。最後の授業日の準備をしているが,テーマは原康夫さんの教科書である第5版 物理学基礎の第25章「原子核と素粒子」だ。

トリチウムのベータ崩壊で話を終らそうと思ったら,宇宙線と上層大気の衝突によって年間に生成されるトリチウム量のところでつまづいてしまった。茨城大学の鳥養さんの資料では,年間72 PBq/y(PBq=ベタベクレル=10^15ベクレル)生成されるとなっている。

そもそも,ベクレルは単位時間当たりの崩壊数なので時間の逆数になっている。これをさらに時間で割った量が生成量であるというのはどういうことかと,かつて理解していたところで再度引っかかってしまった。

70歳を過ぎるとこんなことが増えていくのだろう。一日中家の中で失せ物を探している時間がどんどん増加していくのと同様に,頭の中の失せ物を探す時間が増えていくのだ。こんなときに,生成AIが頼りになれば有難いのだけれど,これが現時点ではあまりあてにはならない。

さて,時間$\Delta t$の間に崩壊する原子核の数は,$\Delta N = \lambda N \Delta t$である。$\lambda = \frac{0.693}{T_{1/2}}$は崩壊定数であり,時間の逆数の次元を持っている。それは不安定な原子が崩壊する確率を表している。言い換えれば,放射性同位元素の物質量$N$に$\lambda$を掛けたものがその物質のベクレル数に等しいことから,放射性同位元素の物質量を,単位が異なるベクレル数で表現しても差し支えないだろうという考えだ。

あるいは,本質的に時間とともに変化する存在である放射性同位元素の量を表現するのに,時間的に不変な状態を想定しているモルやkgで表すのは適当ではなく,むしろその時点でのベクレル数で表わした上で,今後はこの割合で減少していくということに注意を喚起するという習慣があると善意に理解しよう。

まあ,トリチウムが放射平衡している場合は,時間とともに変化しないけれども,いつ何時,核施設の事故があるかもしれない。

さて,2000年のUNSCEARの資料[1] に,Table 4 Production rates and concentrations of cosmogenic radionuclides in the atmosphereという表がある。これによると,宇宙線によるトリチウムの単位面積,単位時間当たりの生成数は,$ 2500 /({\rm m^2 s}) $ であり,地球表面積,$ 5.1 \times 10^{14} {\rm m^2} $との積から,1秒間に,$\mu = 1.28 \times10^{18} $個/sのトリチウム原子が生成される。1年間($y =3.15 \times 10^{7} {\rm s} $)では $\mu y = 4.0 \times 10^{25}$個となる。一方,トリチウムの崩壊定数は,$\lambda = 0.693/ ( 12.3 × 3.15 × 10^7) /{\rm s} $ なので,$\mu y \lambda $によってベクレルに換算すれば,$72 \times 10^{15} {\rm Bq}$が得られる。

また,この自然の機構によって地球上に存在するトリチウムの総量$N(t)$は,次の微分方程式$\frac{dN(t)}{dt}=-\lambda N(t) + \mu $の平衡解 $N(\infty)$で与えられ,$N(\infty) = \frac{\mu}{\lambda} = 7.2 \times 10^{26}$個= $1.28 \times 10^{18}$Bqである。

これを使って,大気中の平均トリチウム濃度を計算してみる。資料[1]では対流圏の体積が,$3.62 \times 10^{18} {\rm m^3}$と与えられ,$0.35 {\rm Bq / m^3}$となる。ところが資料[1]では,$1.4 {\rm mBq / m^3}$となっていて,何だか250倍大きくなってしまうのだ。なんで?

あら,表にはfractional amount in atomosphereというのがあって,その係数が1/250=0.004になっていた。トリチウムはほとんどHTOの形態で存在しているので,ほとんどが雨水/海水に溶けてしまうということなのかもしれない。

図:トリチウムの概念図(東京電力から引用

[2]環境トリチウムについて(鳥養祐二)
[3]トリチウムの環境動態(百島則幸)
[4]大気中トリチウム濃度の変遷と化学形態別測定(宇田達彦・田中将裕)