惑星直列

 太陽系の地球を除く惑星を一直線上に隙間なく並べると，ちょうど地球−月の距離と等しくなるという話を見かけたので確かめてみた。

地球を除く7惑星の直径 [km] は，{水星，金星，火星，木星，土星，天王星，海王星}={4879.4, 12103.6, 6794.4, 142984.0, 120536.0, 51118.0, 49528.0}となる。また，地球−月の平均公転半径 [km]は，388440.0 である。

Mathematicaで次のコードを実行する。

d = {4879.4, 12103.6, 6794.4, 142984.0, 120536.0, 51118.0, 49528.0}
r = Table[Sum[d[[i]], {i, 1, k}] - d[[k]]/2, {k, 1, 7}]
g1 = Graphics[{Circle[{384400/2, 0}, 384400/2, {11 Pi/12, 13 Pi/12}], 
   Circle[{384400/2, 0}, 384400/2, {-Pi/12, Pi/12}]}]
g2 = Graphics[Table[Circle[{r[[k]], 0}, d[[k]]/2], {k, 1, 7}]]
g3 = Graphics[Table[Circle[{r[[k]] - d[[1]], 0}, d[[k]]/2], {k, 2, 7}]]
Show[g1,g2]
Show[g1,g3]

水星をふくめると公転半径から少しはみ出すが，水星を除けばすっぽりとおさまった。

図：地球−月の公転軌道半径と７惑星の直列配置

クローニッヒ・ペニーモデル（４）

クローニッヒ・ペニーモデル（３）からの続き

これまでの結果を用いて，具体的な数値を代入したグラフを作成する。与えられた$k_n$の値に依存して，$-1 \le \cos(k_n a) \le 1$の条件から，1電子のエネルギー$E$がとりうる範囲についての条件が定まる。このとき，電子が取り得るエネルギー領域を許容帯（allowed band），取りえないエネルギー領域を禁制帯（forbidden band）とよぶ。相互作用のない電子の多粒子系では，これらの離散化された許容帯に電子が充填されることになる。

Mathematicaによって，この様子を確認してみれば次のようになる。

a = 2; b = 0.04; hc = 1973; m = 1.022*10^6; p = 2*1.05017;
v = 2*hc^2*p/(1.022*10^6*a*b)
(hc)^2/(m a^2)
200.001
0.952233

\[Alpha][e_] := Sqrt[m*e]*(a - b)/hc;
\[Beta][e_] := Sqrt[m*(v - e)]*b/hc;
\[Gamma][e_] := Sqrt[m*e]*a/hc;
X[e_] := Cos[\[Alpha][e]] Cosh[\[Beta][ e]] + ((\[Beta][e]/b)^2 - (\[Alpha][e]/(a - b))^2)/(2*\[Alpha][ e]/(a - b)*\[Beta][e]/b) Sin[\[Alpha][e]] Sinh[\[Beta][e]]

Plot[{ X[e], 1, (Cos[\[Gamma][e]] + p*Sin[\[Gamma][e]]/\[Gamma][e]), -1, Cos[\[Gamma][e]]}, {e, 0, 200}, PlotRange -> {-1.5, 3.5}, PlotStyle -> {Red, Gray, Blue, Gray, Orange}]
Table[FindRoot[X[e] == 1, {e, (hc a n )^2/(2 m) }], {n, 1, 5}]
{{e -> 2.95059}, {e -> 37.6031}, {e -> 45.1073}, {e -> 150.412}, {e -> 158.267}}
Table[ FindRoot[X[e] == -1, {e, 0.9*(hc a n )^2/(2 m) }], {n, 1, 6}]
{{e -> 9.40077}, {e -> 16.0036}, {e -> 84.6069}, {e -> 92.3771}, {e -> 242.886}, {e -> 242.886 + 0. I}}

d0 = Plot[{10^6 (e - 2.95059), 10^6 (e - 9.40077), 10^6 (e - 16.0036), 10^6 (e - 37.6031), 10^6 (e - 45.1073), 10^6 (e - 84.6069), 10^6 (e - 92.3771), 10^6 (e - 150.421), 10^6 (e - 158.267)}, {e, 0, (4 Pi)^2}, PlotStyle -> Table[{Gray, Dotted}, 9], PlotRange -> {0, 13}];
f0 = Plot[{0, Pi, 2 Pi, 3 Pi, 4 Pi}, {e, 0, (4 Pi)^2}, PlotStyle -> Table[{Green, Dotted}, 5]];
g0 = Plot[\[Gamma][e], {e, 0, (4 Pi)^2}, PlotStyle -> {Orange, Dashed}];
g1 = Plot[ArcCos[X[e]], {e, 0, Pi^2}];
g2 = Plot[2 Pi - ArcCos[X[e]], {e, Pi^2, (2 Pi)^2}];
g3 = Plot[2 Pi + ArcCos[X[e]], {e, (2 Pi)^2, .98 (3 Pi)^2}];
g4 = Plot[4 Pi - ArcCos[X[e]], {e, .98 (3 Pi)^2, .95 (4 Pi)^2}];

Show[{d0, f0, g0, g1, g2, g3, g4}, PlotRange -> {-1, 15}]

図１：境界条件から決まる（energy-cos(ka)）のグラフ

図２：図１から決まる分散関係（energy-ka）のグラフ

ポテンシャルの幅$b$を格子周期$a$の2%に設定したため，$b \rightarrow 0$の近似がよく当てはまる。そこで，$X(pa)$と$\cos \gamma(e) + P \sin \gamma(e)/\gamma(e)\ $の誤差は1%以下となり，グラフ上では区別できなくなっている。

なお，$e=\frac{(\hbar k)^2}{2m}=\frac{(\hbar c)^2}{2mc^2 a^2}\cdot (ka)^2 = 0.952233 (ka)^2\ $であることに注意する。

フェルミ分布

月曜の授業の予習シリーズ。

位相空間（$\mu$空間）の細胞に含まれる状態数は，プランク定数を$h$として，$dn=\frac{1}{h^3} dx dy dz dp_x dp_y dp_z$である。電子のようなスピン1/2のフェルミ粒子を考えると，位相空間の各状態にスピンアップとダウンの2状態がともなう。そこで，単位体積をとって，運動量空間細胞に含まれる状態数は $dn' = \frac{2}{V_0} \int_V dn = \frac{2}{h^3}dp_x dp_y dp_z  = \frac{8 \pi}{h^3} p^2 dp$となる。ただし，$p^2=p_x^2+p_y^2+p_z^2$。粒子の質量を$m$，エネルギーを$w=\frac{p^2}{2m},\ p=\sqrt{2mw}$とすると，$dw=\frac{p}{m}dp\ $より，$dn'=\frac{8 \pi}{h^3} \sqrt{2m^3 w}\ dw\ $となる。

この系のフェルミ分布関数は，$f(w_i)=[\exp(\frac{w_i - w_F}{kT}) + 1]^{-1}$である。ただし，$k$はボルツマン定数，$T$は系の絶対温度，$w_F$はフェルミ準位を表わす。これから，エネルギーの分布関数は，$N(w)dw= \frac{8 \pi}{h^3} \sqrt{2m^3 w} [\exp(\frac{w - w_F}{kT}) + 1]^{-1} dw$となる。これを速度空間の表式にひきもどして，$v_x,v_y$で積分すると$\ v_z$の分布関数が得られる。

そこで，次の関係式に留意する。$\int_{-\infty}^{\infty}dv_x\int_{-\infty}^{\infty}dv_y\int_{-\infty}^{\infty}dv_z f(v_x,v_y,v_z) = 4\pi \int_{0}^{\infty} f(v^2)\ v^2 dv $

$w=\frac{m}{2}v^2$，$dw = m v\ dv$なので

$4\pi \int_{0}^{\infty}f(v^2)  v^2 dv = 4\pi \int_{0}^{\infty} f(v^2)  \frac{v}{m} dw =  \int_{0}^{\infty} 4\pi f(v^2)  \sqrt{\frac{2w}{m^3}} dw$

そこで，$N(w)dw =  4\pi f(v^2)  \sqrt{\frac{2w}{m^3}} dw\ $とおけば，$f(v^2)=N(w) \frac{1}{4\pi} \sqrt{\frac{m^3}{2w}}\ $となるので，$\ v_z$の分布関数 $N(v_z) dv_z$は次式で与えられる。

$N(v_z) dv_z = \int_{-\infty}^{\infty}dv_x\int_{-\infty}^{\infty}dv_y f(v_x,v_y,v_z) dv_z =  \int_{-\infty}^{\infty}dv_x\int_{-\infty}^{\infty}dv_y N(w) \frac{1}{4\pi} \sqrt{\frac{m^3}{2w}} dv_z$

$= \int_{-\infty}^{\infty}dv_x\int_{-\infty}^{\infty}dv_y \frac{2m^3}{h^3} [\exp(\frac{w - w_F}{kT}) + 1]^{-1} dv_z$

$v_x, v_y$平面での積分を2次元の極座標によって実行するため，$u^2=v_x^2+v_y^2$とおいて，$dv_x dv_y = 2\pi u du$となる。

$\therefore \quad N(v_z) dv_z = \frac{4\pi m^3}{h^3} \int_0^\infty du u  [\exp(\frac{\frac{m}{2}(u^2+v_z^2) - w_F}{kT}) + 1]^{-1} dv_z$

$= \frac{4\pi m^3}{h^3} \int_0^\infty \frac{1}{2} dt  [\exp(\frac{\frac{m}{2}(t+v_z^2) - w_F}{kT}) + 1]^{-1} dv_z$

ここで，$a=\exp(\frac{\frac{m}{2}v_z^2 - w_F}{kT})$，$b=\frac{m}{2kT}$とおけば，必要な積分は$\int_0^\infty \frac{1}{a \exp(bt) + 1}dt$となり，その値は$\ \frac{1}{b}\log(1+1/a)\ $である。これより，$N(v_z) dv_z = \frac{4\pi m^2 kT}{h^3} \log (1 + \exp(\frac{w_F - \frac{m}{2}v_z^2}{kT})) dv_z$

これらの分布関数をMathematicaでプロットすると次のようになる。

f[w_, kT_] := Sqrt[w]/(Exp[(w - 1)/kT] + 1)
Plot[Table[f[w, 0.01*k], {k, 1, 10, 2}], {w, 0, 2},PlotRange -> {0, 1}]
g[v_, kT_] := kT/2 Log[(Exp[(1 - .5*v^2)/kT] + 1)]
Plot[Table[g[v, 0.01*k], {k, 1, 10, 2}], {v, 0, 2},PlotRange -> {0, 0.5}]

図１：エネルギー分布関数　N(w)

図２：速度分布関数　N(v_z)

ファン・デル・ワールスの状態方程式

統計物理学の授業がはじまった。専門科目にしては受講者が多かったので，講義室は共通講義棟の１Ｆに割り当てられていた。授業の第１法則：受講者数は時間とともに指数関数的に減少する。授業の第２法則：受講者数の空間密度分布は教卓からの距離の逆数に比例して減少する。

最後にファン・デル・ワールスの状態方程式を紹介して授業を終えたところ，早速質問があった。理想気体では，$pV=n\ R\ T\ $だった状態方程式が，実在気体の場合$\bigl( p + \frac{a}{V^2} \bigl) (V - b) = R\ T\ $と書いてあるけれど，右辺にモル数の$\ n\ $が抜けているのかというものだ。

「あ，ごめんごめん，忘れてたわ」と返事して終ったものの，帰り道に380段の階段を下りながら考えてみると何だかおかしい。圧力の補正項$\ \frac{a}{V^2}\ $は分子間力によるものであり，気体密度の二乗に比例する項だと説明した。ところが，圧力は示強変数なのに，補正項の分母は示量変数の二乗になっている。つまり，$a\ $は定数ではなく示量変数の二乗に比例しなければならない。

演習課題の参考にしていた「熱学入門」（藤原邦男・兵頭俊夫）の式も同様の問題点を含んだままだった。そこで，いくつかの本などでファン・デル・ワールス方程式を調べてみると，気体のモル数を1モルに限定していたり，体積として$\ V\ $ではなく，モル当たり体積$\ V_m\ $を用いている。

ということで，正しくは，$n$モルの実在気体に対しては，$\bigl( p + \frac{a n^2}{V^2} \bigl) (V - n b) = n\ R\ T\ $としなければならない。$a, b\ $の値もこれまで考えたこともなかったが，お茶の水女子大学の理学部編入試験問題にも採用されているくらいであり，水（H2O）の場合，$a=5.54\times 10^{-1} {\rm [Pa \cdot m^6 \cdot mol^{-2}]}$，$b=3.05 \times 10^{-5} {\rm [m^3 \cdot mol^{-1}]}\ $だった。

これを使ってMathematicaで状態図をプロットしてみる。

図：水蒸気のファン・デル・ワールス状態図

a = 0.554; b = 3.05*10^-5; R = 8.314; 8 a/(27 R b)

647.331　・・・　（臨界温度）

p[V_, T_] := R T /(V - b) - a/(V^2)

q[V_, T_] := Abs[R T /(V - b) - a/(V^2)]

LogLogPlot[{q[V, 673], q[V, 647.33], q[V, 573], q[V, 473], q[V, 373], 

  q[V, 273]}, {V, 3*10^-5, 3*10^-4}, PlotRange -> {10^5, 10^9}]

Plot[{p[V, 673], p[V, 647.33], p[V, 573], p[V, 473], p[V, 373], 

  p[V, 273]}, {V, 3*10^-5, 3*10^-4}, PlotRange -> {-10^8, 10^8}]

abc予想（３）

abc予想（２）からの続き

$c > rad(abc)\ $が成り立つabcトリプルを例外的abcトリプルとよぶことにする。そうはいっても前回は，これが無限個存在することも示している。例外的トリプルの個数を計算してみると，$c<1,000\ $で31個（全トリプル151,895個），$c<10,000\ $で120個（全トリプル個15,196,785），$c<100,000\ $で418個（全トリプル1,519,805,376）となる（最大10^15のオーダーまでの数の素因数分解を10億回しなければならないので，めちゃ時間かかった）。

これを計算したMathematicaコードは次の通り。

Rad[n_] := Times @@ (First /@ FactorInteger[n]);
f[n_] := {k = 0; Do[If[Rad[a*b*(a + b)] < a + b && GCD[a, b, a + b] == 1, k = k + 1], {a, 1, n/2 - 1}, {b, a + 1, n - a - 1}]; Print[k]};
g[n_] := {k = 0; Do[If[GCD[a, b, a + b] == 1, k = k + 1], {a, 1, n/2 - 1}, {b, a + 1, n - a - 1}]; Print[k]};
f[1000]
31
g[1000]
151 895

そこで，少しだけ条件を変更することで，例外的トリプルが有限個になるようにできないかを考えたところ，次のような予想が立てられた。(1)と(2)は同等であり，これがabc予想といわれるものである。

(1) 任意の$\varepsilon > 0$ に対して $c > {\rm rad}(abc)^{1+\varepsilon}\ $を満足するabcトリプルは高々有限個しか存在しない。

(2) 任意の$\varepsilon > 0$ に対してある $K(ε) > 0$ が存在し、全ての abcトリプルについて$c < K(\varepsilon) {\rm rad}(abc)^{1+\varepsilon}\ $が成り立つ。

ロフテッド軌道（３）

 ロフテッド軌道（２）からの続き

初速度を第１宇宙速度$v=\sqrt{\frac{GM}{R}}$に固定した場合，Mathematicaによるプロットを試みる。このとき，$u(\varphi)=A(\theta) \cos \varphi + B(\theta) \sin \varphi +\frac{GM}{h^2}$，$A(\theta)= - \frac{1}{R \tan^2 \theta }$，$B(\theta) = -\frac{1}{R \tan \theta}$，$\frac{GM}{h^2}=\frac{1}{R \sin^2 \theta}$ であり，Mathematicaのコードは以下のとおり。

In[1]:= {G, R, M} = {6.67*10^-11, 6.37*10^6, 5.97*10^24}
Out[1]= {6.67*10^-11, 6.37*10^6, 5.97*10^24}
In[2]:= {g, v} = {G M /R^2, Sqrt[G M/R]}
Out[2]= {9.81344, 7906.43}

In[3]:= A[t_] := 1/R*(1 - 1/Sin[t]^2)
In[4]:= B[t_] := -1/(R*Tan[t])
In[5]:= GM[t_] := 1/(R Sin[t]^2)
In[6]:= r[s_, t_] := 1/(A[t] Cos[s] + B[t] Sin[s] + GM[t])
In[7]:=
Table[f1[i] = Plot[{r[s, i]/R, 1}, {s, 0, Pi},
GridLines -> Automatic, PlotRange -> {0, 2.0}];, {i, 0.1, 1.5, 0.2}];
In[8]:= Show[Table[f1[i], {i, 0.1, 1.5, 0.2}]]
In[9]:=
Table[f2[i] = PolarPlot[{r[s, i]/R, 1}, {s, 0, Pi},
GridLines -> Automatic, PlotRange -> {0, 2.0}];, {i, 0.05, 1.5, 0.2}];
In[10]:= Show[Table[f2[i], {i, 0.05, 1.5, 0.2}]]

In[11]:= f[t_] :=
NIntegrate[ 1/Sqrt[G M R Sin[t]^2] *1/(A[t] Cos[s] + B[t] Sin[s] + GM[t])^2, {s, 0.001, Pi/30}]

In[12]:= Plot[f[t], {t, 0, Pi/8}]

図１：ロフテッド軌道（横軸$\varphi$ラジアン，縦軸$r/R$）

図２：ロフテッド軌道（上記と同じものを極座標表示）

この近似のもとでは，6500kmの高度に達する弾道ミサイルはほぼ地球を半周するところまで到達できるので，北朝鮮から米国本土全体が射程に入ることになる。鉛直上方から60度の角度で射出すれば，地球を1/4周するので1万kmまで到達できる。

実際には，空気抵抗があることや，初速度はロケットエンジンによる一定時間の加速によって獲得されることなどから，もう少し真面目な計算をする必要がある。

［１］ロケットの運動と人工衛星の打ち上げ（冨田信之）

平方完成

 平方完成は，入試問題を解くときなど，条件設定の場面でたいへん重宝する技法だ。ちょっと手計算が面倒な式がでてきたので，Mathematicaに任せようと思った。

ところが，探してみてもMathematicaで平方完成する関数が組み込まれていないようなのだ。もしかしたら調べ方が足りないのかもしれないが，普通に考えるとイの一番に出てきても良さそうな機能なのだが。

それらしいユーザ定義関数がいくつか見つかったけれど，2変数の整式を代入しても思ったような変形ができず，望みのものではなかった。しかたがないので，自分で関数パーツを考えることにした。これを一般化するには，Mathematicaプログラミングにおける文法の知識が足りなさすぎる。

ここで考えたのは，ある変数の二次式を与えたときに，平方完成された部分と残余部分のリストを返すユーザ定義関数 sq[式, 変数]だ。変数がn個ある場合は，n回繰り返して使う必要があるという残念なコード素片だ。

sq[f_, v_] :=
　Module[{a, b},
　　a = Coefficient[f, v^2];
　　b = Coefficient[f, v];
　　{a (v + b /(2 a))^2,
　　c = f - a (v + b /(2 a))^2}] // Simplify

これを使って次のような計算ができる。

In[1]:= sq[2 x^2 - 4 x y + 2 y^2 + 24 x - 24 y + 288 + 3 x y, x] 

Out[1]= {1/8 (-24 - 4 x + y)^2, 216 - 18 y + (15 y^2)/8}

In[2]:= sq[%[[2]], y] 

Out[2]= {3/40 (24 - 5 y)^2, 864/5}

四項間漸化式

 三項間漸化式からの続き

前回の問題を次のように修正する。2次方程式$\ x^2-p x + q = 0\ $の解を$\alpha, \beta$として，$a_n = \alpha^n +\beta^n$と定義する。ただし，$a_0=2,\  a_1 =\alpha + \beta = p,\  a_{-1}=\frac{1}{\alpha}+\frac{1}{\beta} = \frac{\alpha + \beta}{\alpha \beta} = \frac{p}{q}$である。

このとき，次の漸化式，$a_{n+1}=p\ a_{n}-q\ a_{n-1}, \ (n=0,1,2 \cdots)$が成り立つ。

また，数列の一般項は $a_n = \alpha^n + \beta^n = \Bigl( \dfrac{p + \sqrt{p^2-4q}}{2} \Bigr)^n +  \Bigl(\dfrac{p - \sqrt{p^2-4q}}{2} \Bigr)^n$ で与えられる。

これを3次方程式に拡張すると次のようになる。3次方程式$\ x^3-p x^2 + q x -r = 0\ $の解を$\alpha, \beta, \gamma$として，$b_n = \alpha^n + \beta^n + \gamma^n$と定義する。ただし，$b_0=3,\  b_1 =\alpha + \beta + \gamma = p,\  b_{-1}=\frac{1}{\alpha}+\frac{1}{\beta} +\frac{1}{\gamma} = \frac{\alpha \beta + \beta \gamma + \gamma \alpha}{\alpha \beta \gamma} = \frac{q}{r}$である。

このとき，漸化式は，$b_{n+2}=p\ b_{n+1}-q\ b_{n}+r\ b_{n-1}, \ (n=0,1,2 \cdots)$となる。

さらに，数列の一般項は $b_n = \alpha^n + \beta^n + \gamma^n$なので，基本対称式の値，$p=\alpha+\beta+\gamma,\ q=\alpha \beta + \beta \gamma + \gamma \alpha,\ r=\alpha \beta \gamma$の多項式になる。2次方程式の場合のように，3次方程式の一般解を代入してもよいが，非常に煩雑な表現になってしまい，Mathematicaを持ってしてもExpandとSimplifyを組み合わせて n=8で行き詰まってしまった。もちろん，漸化式を使えば大丈夫なのだが。

a = x /. Solve[x^3 - p x^2 + q x - r == 0, x];

b[n_] := a[[1]]^n + a[[2]]^n + a[[3]]^n

Table[Expand[b[i]] // Simplify, {i, 1, 8}]

{p, p^2 - 2 q, p^3 - 3 p q + 3 r, p^4 - 4 p^2 q + 2 q^2 + 4 p r, p^5 - 5 p^3 q + 5 p q^2 + 5 p^2 r - 5 q r, p^6 - 6 p^4 q + 9 p^2 q^2 - 2 q^3 + 6 p^3 r - 12 p q r + 3 r^2, p^7 - 7 p^5 q + 14 p^3 q^2 - 7 p q^3 + 7 p^4 r - 21 p^2 q r + 7 q^2 r + 7 p r^2, p^8 - 8 p^6 q + 20 p^4 q^2 + 8 p^5 r - 32 p^3 q r + 24 p q^2 r + 2 q (q^3 - 4 r^2) - 4 p^2 (4 q^3 - 3 r^2)}

三角関数の積

その鈴木貫太郎の問題で次の式の値を求めよというのがあった。

$\cos{\dfrac{\pi}{33}}\cos{\dfrac{2\pi}{33}}\cos{\dfrac{4\pi}{33}}\cos{\dfrac{8\pi}{33}}\cos{\dfrac{16\pi}{33}}$

三角関数の積和の式を繰り返し使うのか，それにしても面倒だろう，どうするのかなあと解答を見ると，$\sin{\dfrac{\pi}{33}}$をかけて倍角の公式からドミノ倒しのようにしてあっという間に解けてしまった。なるほどね。

そこでMathematicaで一般化してみた。

f1[n_, m_] := Product[Cos[2^k Pi/(2^n + 1)], {k, 0, m}] // Simplify
g1[j_] := Table[f1[j, i], {i, j - 1, 20, j}]
g1[5]
{1/32, -(1/1024), 1/32768, -(1/1048576)}

f2[n_, m_] := Product[Cos[2^k Pi/(2^n - 1)], {k, 0, m}] // Simplify
g2[j_] := Table[f2[j, i], {i, j - 1, 20, j}]
g2[5]
{-(1/32), -(1/1024), -(1/32768), -(1/1048576)}

フェルマーの小定理

フェルマーの小定理は，素数を$p$とし，$p$とは互いに素な整数を$a$として，$a^{p-1} \equiv 1 \quad (\mod p \ ) $というものである。

頭の体操のためによく見ている鈴木貫太郎の YouTubeチャンネルでは，ちょっと目には解きにくそうな整数問題に，縦横無尽にmodを使って条件を絞り込んでいくというパターンがよく見られる。自分が，高校生のときには，こんなふうにしてmodを活用するということがなかったので，新鮮な感じがしている。ここまで自由に使えれば便利には違いない。

そんなわけで，フェルマーの小定理あるいは，その他のmod問題の勘を養成するための１行コードをMathematicaで書いてみた。

f[n_] := Table[ Table[Mod[i^k, Prime[n]], {i, 1, Prime[n]}], {k, 1, Prime[n] - 1}]
g[n_] := Table[Table[Mod[i^k, n], {i, 1, n}], {k, 1, n - 1}]

Do[{Print[g[i], " "]}, {i, 2, 7}]
{{1,0}}
{{1,2,0},{1,1,0}}
{{1,2,3,0},{1,0,1,0},{1,0,3,0}}
{{1,2,3,4,0},{1,4,4,1,0},{1,3,2,4,0},{1,1,1,1,0}}
{{1,2,3,4,5,0},{1,4,3,4,1,0},{1,2,3,4,5,0},{1,4,3,4,1,0},{1,2,3,4,5,0}}
{{1,2,3,4,5,6,0},{1,4,2,2,4,1,0},{1,1,6,1,6,6,0},{1,2,4,4,2,1,0},{1,4,5,2,3,6,0},{1,1,1,1,1,1,0}} 

定偏角分散プリズム

ペラン・ブロッカなのかペリン・ブロッカなのかペロン・ブロッカなのかよくわからないが，定偏角分散プリズムというものがあることがわかった。いやあまりわかっていない。幾何光学をまともに勉強していないので。

とりあえず，イメージを掴むために，Mathematicaでプリズムの光路シミュレーションプログラムを書いてみた。わかったような，わからないような・・・それでも高等学校で学ぶ幾何光学の法則さえ知っていればコードは書ける。

プリズムは四角形であり原点(0,0)の内角は90度（直角），y軸上の点（0,2√3-2）は75度，(√3,1)にある点は135度，x軸上の点（√3+1/√3,0）は60度となっていて，その屈折率をnとする。

aは入射光の入射面からはかった入射角度，tbは入射光が屈折した光線の進行方向のプリズム底面に対する傾き，(x1, y1)は全反射する点の座標，tcは反射光の底面に対する傾きの絶対値，tdはプリズムからの出射光の底面からはかった出射角度の勾配，x2は出射点の座標である。

n = 1.6; a = 36.8699;
t1 = Tan[15 Degree]; t2 = Tan[30 Degree]; x3 = Sqrt[3] + 1/Sqrt[3];
tb[n_, a_] := Cos[a Degree]/Sqrt[n^2 - Cos[a Degree]^2]
x1[n_, a_] := Sqrt[3]*t1/(t1 + tb[n, a])
y1[n_, a_] := 1 + Sqrt[3]*t1*tb[n, a]/(t1 + tb[n, a])
tc[n_, a_] := (t2 + tb[n, a])/(1 - t2*tb[n, a])
td[n_, a_] := Sqrt[tc[n, a]^2 - n^2 + 1]/n
x2[n_, a_] := x1[n, a] + y1[n, a]/tc[n, a]

ro = Point[{(1 + 2*t1)/(1 + t1), (1 + 2*t1)/(1 + t1)}];
pol = Polygon[{{0, 0}, {0, 1 + t1}, {1, 1}, {x3, 0}}];

lin[n_, a_] := Line[{{- Sin[a Degree], 1 - Cos[a Degree]}, {0, 1.0}, {x1[n, a], y1[n, a]}, {x2[n, a], 0}, {x3, -td[n, a]*(x3 - x2[n, a])}}]
das[n_, a_] := Line[{{x2[n, a], 0}, {-0.5, -td[n, a] (-0.5 - x2[n, a])}}]

f[n_, a_] := Graphics[{LightGray, pol, Thick, Red, lin[n, a], Dashed, Blue, das[n, a],PointSize[0.03],ro}]

{f[n, a], td[n, a]/Tan[a Degree]}
Manipulate[{f[n, a], n, a, td[n, a]/Tan[a Degree]}, {a, 25, 50}]

図　Pellin-Broca PrismのMathematicaシミュレーション

なお，td[n,a]/Tan[a Degree] = 1ならば入射光と出射光が直交することになる。

［１］定偏角分散プリズムの最適回転中心の位置（小穴純）

球面調和関数（１）

tsujimotterさんが日曜化学：量子力学の基本と球面調和関数の可視化というテーマで python のmatplotlib を使ったコードを書いていたところ，東北大学の堀田さんに褒められていた。

なるほど，これをJuliaでやってみようかと思ったのだが，Plots の三次元散布図のカラーオプションの制御方法がわからなくてとりあえず挫折した。

Mathematicaだと簡単に実例が見つかるんだわこれが。

With[{lmax = 2}, Table[If[Abs[m] <= l,
SphericalPlot3D[
Abs[2 SphericalHarmonicY[l, l-m, t, p]]^2, {t, 0, Pi}, {p, 0,
2 Pi}, PlotRange -> {{-0.5, 0.5}, {-0.5, 0.5}, {-0.5, 0.5}},
BoxRatios -> {1, 1, 1}, PlotPoints -> 30, Axes -> False,
Boxed -> False, Mesh -> False, ColorFunctionScaling -> False,
ColorFunction ->
Function[{x, y, z, t, p, r}, Blend[{Blue, Orange}, (Cos[Arg[SphericalHarmonicY[l, l - m, t, p]]] + 1)/
2]]], Null], {l, 0, lmax}, {m, 0, lmax}]] //
GraphicsGrid[#, AspectRatio -> 1] &

図　球面調和関数の例 by Mathematica（筑波大学　武内修さんより引用）

 

イベント制限

なにがなんでも オリンピックに向かってまっしぐらに進んでいる今日このごろ。

6月20日迄の緊急事態宣言と蔓延防止等重点措置は解除の方向となり，変更されたものの期限は7月11日迄である。

緊急事態宣言10都道府県｛北海道，東京，愛知，京都，大阪，兵庫，岡山，広島，福岡，沖縄｝のうち，沖縄は継続，岡山と広島は解除，その他下線のところは蔓延防止措置に移行。

蔓延防止等重点措置5県｛埼玉，千葉，神奈川，岐阜，三重｝の内，岐阜と三重は解除，その他は継続。したがって，蔓延防止等重点措置は合計10都道府県となる。

イベント制限については，緊急事態宣言と蔓延防止等重点措置が出されている都道府県に対して，現在の定員50%以内かつ上限5000人を解除後は定員50%以内かつ上限10000人に変更し，1ヶ月後程度で問題がなければ，現在のその他地域と同様に，5000人以下か定員50%以内のいずれか大きい方に変更する予定である。

7月23日のオリンピック開始時には，定員の50%がOKになるということだ。関係者が口をそろえて今回のイベント制限変更措置は五輪とは関係ないというのでますます怪しさが露呈している。わかりにくいので，Mathematicaで現在と今後の制限のようすを図示してみる。

現在： g1 = RegionPlot[ {y < 0.5 x && y < 5, y <= 0.5 x || y <= 5}, {x, 0, 40}, {y, 0, 20}, PlotLegends -> "Expressions"]

今後： g2 = RegionPlot[ { y < 0.5 x && y < 10, y <= 0.5 x || y <= 5}, {x, 0, 40}, {y, 0, 20}, PlotLegends -> "Expressions"]

図　イベント制限（左6/20まで，右6/21から，単位は千人）

P. S. 日本経済新聞の朝刊記事の内容をまとめると上のようになったのだけれど，6月17日に開催された，新型コロナウイルス感染症対策本部(第 69 回) によると微妙に解釈が間違っているのかもしれない。その場合，その他地域でも上限の1万人が存在していることになるのだけれど，ちょっと自信がない。どうして，こんなにわかりにくい文書を作れるの？

角運動量代数（２）

角運動量代数（１）からの続き

Mathematicaの関数には，かなり昔からClebshGordan係数が含まれていた。さすがに素粒子物理を専攻していたウルフラムだと思った。Mathematicaには次の関数が実装されている。

ClebschGordan[{j1,m1},{j2,m2},{j,m}]　|jm〉から|j1m1〉|j2m2〉への分解に対する，クレプシュ・ゴルダン係数を与える。
ThreeJSymbol[{j1,m1},{j2,m2},{j3,m3}]　ウィグナーの3-jシンボルの値を与える（m1+m2+m3=0）。
SixJSymbol[{j1,j2,j3},{j4,j5,j6}]　ラカーの6-jシンボルの値を与える。

Juliaではどうかというと，WingerSymbols.jlというパッケージがあり，以下の関数が使える。
wigner3j(j1,j2,j3,m1,m2,m3=-m1-m2)，wigner6j(j1,j2,j3,j4,j5,j6)，clebschgordan(j1,m2,j2,m2,j3,m3=m1+m2)， racahV(j1,m2,j2,m2,j3,m3=-m1-m2)，racahW(j1,j2,J,j3,J12,J23)

いずれにせよ9-j symbolは実装されていないのだった。一方，PythonのSymPyにあるWignerSymbolsには，Wigner 3-j，6-j，9-j，Clebsch-Gordan係数，Racah係数，Gaunt係数が含まれている。Gaunt係数という名前は知らなかったが，換算行列要素の形では頻出する係数だ。それより，JuliaのWignerFamilies.jlというパッケージで，3-j Symbolの配列が計算されてグラフ化されているのが印象的だった。 

using WignerFamilies
# wigner3j for all j fixing j2=100, j3=60, m2=70, m3=-55, m1=-m2-m3
w3j = wigner3j_f(100, 60, 70, -55)
js = collect(eachindex(w3j))
plot(js, w3j.symbols)

図　f(j) = w3j(j,100,60,-15,70,55) （JuliaPackageから引用）

素数の計算（１）

 ソフィー・ジェルマン素数の続き

件の記事の中には，「p=999671で 2p+1=1999343，が約数となる。 2^p−1=2^99671−1 は10進法で300932桁となったが，開発用等ではないスペックのかなり低いパソコンでも2～3分で計算してくれた」とあったので追試しようとして，はまってしまった。

Mathematicaでは，このあたりのソフィー・ジェルマン素数を次のように一瞬で計算した。

In[1]:= Timing[ Do[p = Prime[2^16 + i];  If[Mod[p, 4] == 3 && PrimeQ[2*p + 1],  Print[p, " : ", Mod[2^p - 1, 2*p + 1]]], {i, 12900, 13000}]]
999611 : 0
999623 : 0
999671 : 0
1000151 : 0
1000211 : 0
1000403 : 0
Out[1]= {0.007209, Null}

In[2]:= Timing[Print[p = Prime[10^8], " ", PrimeQ[p]]]
2038074743 True
Out[2]= {0.000185, Null}

そこで，10^8番目の素数をpython3で計算して判定すると，

from sympy import *
import time
start_time=time.time()
p=prime(10**8)
end_time=time.time()
print(p,isprime(p),end_time-start_time)

2038074743 True 2.994921922683716

同様に，10^8番目の素数をjulia 1.5で計算して判定すると， 

using Primes
n=10^8
function sgp(n)
  p=prime(n)
  println(isprime(p)," ",p)
end

@time sgp(n)
true 2038074743
420.535971 seconds (333.45 M allocations: 4.970 GiB, 0.24% gc time)

なお，p=999671は，prime(78476)であり，この計算時間はjuliaでも0.153秒だった。それにしてもこの遅さはいったい何ということでしょう・・・？ 

Mathics

Mahtematicaについてからの続き

Mathematicaを使い始めて30年くらいになるが，アカデミック価格の年間契約からはずれたため，2019年の12.0.0.0でバージョンアップが止まってしまった。最新のMathematicaは12.1だと思われる。Mathematicaは家庭・趣味用デスクトップライセンスでも45,000円なので今となってはかなり高額である。学生用ライセンスこそ21,000円だけれど，大学だとその10倍の210,000円，普通の企業だと，469,000円もするわけで相当厳しいのであった。

で，普通はMathematicaの代替はないかと考えてSegeMathだとかいうことになる。いっときこれをインストールしていたのだけれどかなりサイズが大きいので，SSD領域確保のために消去してしまった。Mathematicaのほうが慣れている分使いやすいことでもあるし。M1マシンとBigSurに移行したらどうなるかはわからない。Jupyter上のJuliaが登場したことによって，かなりの部分はこれで代替できるのだけれど，それでもMathematicaの方がコーディングの負荷が少ない。

さて，そんな中で今日発見したのが，フリーのMathematica文法互換ソフトウェアのMathicsである。円周率は4000桁くらいまでしか計算できなかったので，あくまでも教育用という位置づけになるかもしれない。それでもちょっと安心でWolframAlphaよりはましかもしれない。

インストールは簡単で，(1) pip3 install Mathics3， (2) brew install sqlite (macOSの場合)，(3) \$ mathics でコマンドライン版起動，(4) \$ mathics server でGUI版起動 (localhost:8000) であった。List of computer algebra systems にはリストアップされていないのだけれどどうしてかな？

Mathics 1.1.0 on CPython 3.8.6 (default, Oct 27 2020, 08:57:44) using SymPy 1.6.2, mpmath 1.1.0

Copyright (C) 2011-2020 The Mathics Team.

This program comes with ABSOLUTELY NO WARRANTY. This is free software, and you are welcome to redistribute itunder certain conditions. See the documentation for the full license.

Quit by pressing CONTROL-D

THE MATH(S) FIX

THE MATH(S) FIX はスティーヴン・ウルフラムの弟のコンラッド・ウルフラムによる数学教育革新のための青写真についての本である。

ウェブサイトのはしがきを訳すると，こんな感じだった。

 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

AI時代の教育の青写真

なぜ私たちは皆，生涯のうち何年にもわたって数学を教えられているのか？それは本当に役に立っているのだろうか，あるいは，ほとんどが失敗であり多くの人に学びの力を与えられていないのか？それはAI時代に不可欠なものか，それとも時代遅れの通過儀礼なのだろうか？

"The Math(s) Fix: An Education Blueprint for th AI Age" は，なぜ数学教育が世界的に危機的状況にあるのか，また，どうして根本的に新しい主流科目を創るのことが唯一の解決策なのか，を明らかにする画期的な本だ。

この本は，今日の数学教育が，現代的な計算・データ科学・人工知能（AI）が必要とされる現代社会を発展させるための機能を十分果たしていないということを主張する。その代わりに，学生はコンピュータが得意とするものと競争することを強いられて負けてしまっている。

本書は，「数学が苦手」であることが，学習者のせいというより，科目の失敗が原因であること，すなわち，教育のエコシステムが行き詰まり，学生・親・教師・学校・雇用者・政策立案者が，現実世界の要求に追いつこうと間違った方向に走っていることを説明する唯一の本である。

しかし，この本はさらに先を行くものである。問題を解決して，AI時代の教育の普遍的な改革の種を蒔くために，学校における中核的な計算思考科目としての完全に代替的なビジョンを初めて提示しているものだから。

 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

これがたぶんコンピューテーショナル・シンキング（計算論的思考）の王道なのだと思われる。でも，GIGAスクール構想にみんな出払ってしまっており，プログラミング教育とその周辺の話題はすっかり霞んでしまっていた。

図：THE MATH(S) FIX の表紙（Wolfram Media から引用）

［１］Let's Fix Math Education （computerbasedmath.org）

緊急事態宣言と接触制限モデル（３）

緊急事態宣言と接触制限モデル（２）からの続き

4月7日のtwitterで牧野淳一郎さんがいうには，

一応どれも現実をモデルしているはずのところで、そこまで行動抑制しないといけない割合が、西浦氏:0.2、大橋氏「0.45](R_0>1でいくとして)、佐藤氏 0.02 で3人で20倍違う時点で少なくとも2人は間違っているのは明らかであろう。

西浦氏は，緊急事態宣言と接触制限モデル（１）で取り上げた西浦博さん（北海道大学），大橋氏は，新型コロナウイルス感染症の 流行予測の大橋順さん（東京大学），佐藤氏は，COVID-19情報共有の佐藤彰洋さん（横浜国立大学）である。

上に示されている数字は行動抑制の因子であり，モデルの感染率にかける係数だと思われる。ここでは，どれが正しいかについては考察せずに，ドイツなど欧州の基本再生産数$R_0=2.5$（7日増倍率が7倍）に基づいた西浦さんの計算（我々の評価では$R_0=2.3$だったが）を，東京の現状に合わせて再計算した結果を報告する。

前回述べたように，現在の東京では，7日増倍率が2.5倍（$R_0=1.65$）程度なので，西浦さんが用いている値よりかなり小さいのだ。単純なSIRモデルを用いてこれを再現するパラメタセットを探し，その場合に接触制限の効果がどうなるかを調べよう。

やり方は前回と同じなので，パラメタを示す。β=0.28，初期値はν=0.005とすれば上記の増倍率が再現できた。つまり，さきのMathematicaコードにおいて，
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
f[n_, β_, α_, ν_, p_]:= 
NDSolve[{S'[t]==-(p+(10-p)*Tanh[2*(20-t)])/10*β/n*J[t]*S[t], 
J'[t]==(p+(10-p)*Tanh[2*(20-t)])/10*β/n*J[t]*S[t]-J[t]/α, 
R'[t]==J[t]/α, S[0]==n, J[0]==ν, R[0]==0}, {S,J,R}, {t,0,100}];
n=1400; β=0.28; ν=0.005; p=10;
sol = f[1400, β, 7, 0.005, p];
s[t_] := S[t] /. sol[[1, 1]];
i[t_] := J[t] /. sol[[1, 2]];
cft1[t_]:= (p+(10-p)*Tanh[2*(20-t)])/10*β/n*i[t]*s[t];
・・・
Plot[{cf1[t], cf2[t], cf3[t], cf4[t], cf5[t], cf6[t]},
{t, 10, 30}, PlotRange -> {0, 0.1},
GridLines -> {{19, 20, 21}, {0.005, 0.01, 0.015, 0.02, 0.025}}]
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
などとすると計算結果は次のようになった。

図　SIRモデル＋シグモイド関数制限措置による新規感染者数の変化（東京のデータ）

　（上から順に，青：制限なし，80%，60%，50%，40%，20%に制限）

増加と減少の臨界点は50%あたりにある。すなわち，数字だけでいえば大橋氏の値と近い結果が得られたのかな。大橋氏の資料に対する牧野さんのもろもろの批判にはうなずけるものもあるけれど。そう，このブログも拡散してはいけない情報の仲間である。

まとめ

新型コロナウイルス感染症専門家会議（というか西浦博さん）から出てきた接触8割削減の前提条件は，欧州並の$R_0=2.5$に対応する場合ということだった。それを現在の東京の$R_0=1.65$に対応させると，上記のように接触6割削減で新規感染数は減少させることができる。いずれにせよモデル計算の結果なので，他地域の今後の動向や安全率も勘案して，目標値として接触8割削減を掲げることは意味があるだろう。

P. S. 佐藤彰洋さんのシミュレーションに対する疑義がでてきた。ちょっと安心した。Delay Differential Equationによるモデリングは，それでもまだ自分には理解できていない。

P. P. S. 4/10に牧野淳一郎さんの解説が出てきたので問題点がよく理解できるようになった。ただ，6割削減でも新規感染数が減少に転ずるという牧野さんの解説は，7割削減が必要であるというこちらの結果とは食い違っていた。

［１］人との接触７～８割減、効果は　専門家「感染抑制できる」、全員やれば６割減でも　新型コロナ（朝日新聞）→タイトルと本文があまり整合していない記事

［２］佐藤彰洋氏(横浜市大)の新型コロナ感染予測シミュレーションに関する疑義について（本堂毅・佐野雅己・松下貢）

［３］いろいろなモデル計算（牧野淳一郎，2020.4.10）

［４］交通整理（牧野淳一郎，2020.4.12）

［５］公開質問状に対するコメント（佐藤彰洋，2020.4.13）

［６］交通整理の続き（牧野淳一郎，2020.4.12）

［７］交通整理の続きその２（牧野淳一郎，2020.4.16）

［８］新型コロナウイルス感染症対策専門家会議の見解等（牧野淳一郎，2020.4.18）

緊急事態宣言と接触制限モデル（１）

北海道大学の西浦博さんが4月3日にマスコミで示していた図がある。対策がない状態の新規感染者数が指数関数的な増大をしているときに，人の接触を20%減らした場合と80%減らした場合の新規感染者数の変化を示した図だ。これにより彼は欧米に近い外出制限の必要性を訴えていた。残念ながら，これは専門家会議や政府の共通のコンセンサスとはならず，やがて4月7日の7都府県に対する5月6日までの緊急事態宣言につながっていく。

ここではその図がどのようにして得られたものかを，素人が持っている簡単な道具と知識で理解することを目指す。報道されたニュースや伝わってくる情報は鵜呑みにはせず，自分で考えることが必要だから。集団としては人口1400万人の東京都を想定する。日本経済新聞におけるこの内容に関する報道の文脈でも東京を対象としていることがうかがえる。東京の4月3日における新規感染者数は14人（新規感染数累計は753人）であり，人口の1ppm条件がちょうど達成された頃なので，これ以降をSIRモデルなどの単純な微分方程式系で扱うことは可能だろう。

今，探してもみあたらないのだが，西浦さんはこの報道の後に解説のための動画を公開していた。それによればしばらく前のドイツなみの再生産数を仮定し，現時点から15日目以降に制限措置を行った場合をシミュレートしているとした。なお，彼のモデルでは，計算開始から15日目が現在であり，30日目が制限措置の開始時点となる。

今から3週間前のドイツでは新規感染数累計の7日増倍率が7倍になっていた。仮に6.3倍（1日に1.3倍）を仮定すると，基本再生産数は$R_0 = 1 + 5/7 * \log 6.3 \approx 2.3$となる*。一方現在の東京では新規感染数累計の7日増倍率が2.5倍（1日に1.14倍）程度なので，基本再生産数は$R_0 = 1 + 5/7 * \log 2.5 \approx 1.65$（倍加時間と基本再生産率を参照）である。だからこのドイツ並の仮定はちょっとどうかと思うが，いちおうこれを採用する。
（*西浦さんはドイツなどの欧州の平均基本再生産数が$R_0=2.5$であると話していた）

つまり，現在100人の新規感染者数はこのまま放置すると15日後のt=30（1.3^15≒50）に5000人以上に達するという状況を西浦さんは想定していると思われる。今の東京の水準をそのままあてはめた場合は，15日後のt=30（1.14^15≒7）に東京の新規感染者数は700人程度であることに注意しよう。

SIRモデルに接触制限措置の効果を含めた微分方程式系をMathematicaで書くと次のようになる。nは集団人口（1400，以下人数の単位は1万人当りとなる），βは感染率（0.44），αは感染期間（5），νは感染数の初期値（0.0002），pは制限措置の程度（p=10→制限なし，p=9→80%に制限，p=6→20%に制限）である。なお，接触制限措置は，βにかけるシグモイド関数（幅が1日のオーダー）によってモデル化して階段関数は用いない。t=15が4月3日現在（新規感染者数=100）に対応する。

- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
f[n_,β_,α_,ν_,p_]:=
NDSolve[{S'[t]==-(p+(10-p)*Tanh[2*(30-t)])/10*β/n*J[t]*S[t], 
J'[t]==(p+(10-p)*Tanh[2*(30-t)])/10*β/n*J[t]*S[t]-J[t]/α, R'[t]==J[t]/α, S[0]==n,J[0]==ν,R[0]==0}, {S,J,R}, {t,0,100}]

β=0.44 sol = f[1400, β, 7, 0.0002, 10];
s[t_] := S[t] /. sol[[1, 1]];
i[t_] := J[t] /. sol[[1, 2]];
r[t_] := R[t] /. sol[[1, 3]];
cf1[t_] = (10+0*Tanh[2*(30-t)])/10*β*i[t]*s[t]/1400;

sol = f[1400, β, 7, 0.0002, 9];
s[t_] := S[t] /. sol[[1, 1]];
i[t_] := J[t] /. sol[[1, 2]];
r[t_] := R[t] /. sol[[1, 3]];
cf2[t_] = (9+1*Tanh[2*(30-t)])/10*β*i[t]*s[t]/1400;

sol = f[1400, β, 7, 0.0002, 6];
s[t_] := S[t] /. sol[[1, 1]];
i[t_] := J[t] /. sol[[1, 2]];
r[t_] := R[t] /. sol[[1, 3]];
cf3[t_] = (6+4*Tanh[2*(30-t)])/10*β*i[t]*s[t]/1400;

Plot[{cf1[t], cf2[t], cf3[t]}, {t, 15, 45}, 
PlotRange -> {0, 2}, GridLines -> 
{{29, 30, 31}, {0.1, 0.4, 0.5, 0.6, 0.7, 0.8}}]
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

計算の結果，次のグラフが得られた。西浦さんのものと概ね一致する。

図１　SIRモデル＋シグモイド関数制限措置による新規感染者数の変化

　（縦軸の単位は万人，横軸の単位は日＝3/20を基準とした経過日数,
青：制限なし，オレンジ：80%に制限，緑：20%に制限）

図２　上記に，60%，40%，30%に制限などの場合を加えたもの

なるほど，中途半端な制限ではだめだということか。30%に制限するあたりが増加かどうかの臨界点かもしれない。この度の緊急事態宣言はどの程度の効果があるのだろうか。

感染症の数理シミュレーション（５）

感染症の数理シミュレーション（４）からの続き

ここ１，２週間が瀬戸際といわれ続けて，２週間を迎えようとしている。そこで，対策の効果や時期を感染率$\beta(t)$という１つの時間依存のパラメタに代表させ，これが感染者の増加にどのような影響を及ぼすかを見えるようにしたい。

感染症の数理シミュレーション（３）では，SIIDR2モデルとして，感染対策によって感染率を下げる効果を，感染シミュレーション開始時を基点とする指数関数的な因子$\beta(t) = \beta \exp(-t/\lambda) $によって導入した。これは，感染対策時間$\lambda$の数倍で感染率をほとんどゼロに近づけてしまうが，それは難しいのではないだろうか。

そこで，感染率を1桁程度減少させるとともに，対策導入日$\tau$を新たなパラメータとして加えた新しい対策時間因子を含む感染率を定義する。
\begin{equation}
\beta(t) =\dfrac{\beta}{15} \{ 8 + 7 \tanh \dfrac{-(t - τ)}{\lambda} \}
\end{equation}
感染対策の前後にはそれぞれ$\beta$及び$\beta/15$という一定の感染率を与える（15にはとくに根拠はない，パラメータの定量的な精度は求めておらず，オーダやファクターの評価ができればよいという立場である）。もし，対策開始日が$\tau=0$であれば，シミュレーション開始日の$t=0$の$\beta(0)=8\beta/15$となるので，これまでの$\beta$の定義の半分の値になっていることに注意する。

$\beta(t)$を$\beta$で割って，感染対策時間因子だけを取り出したものをグラフにすると次のようになる。横軸はシミュレーション開始からの日数を表す。

図１　感染対策時間因子β(t) λ=5 左からτ=5,10,15

図２　感染対策時間因子β(t)  τ=20 急峻な方からλ=5,10,15

どうやら政府は情報公開に力を入れるのではなく，マスコミやインターネットにおける情報抑制に全力を投球し始めた。ますます，正確な情報が見えなくなってしまう。

感染症の数理シミュレーション（６）に続く