Mathics（３）

Mathics（２）・wolframscript からの続き

図：jupyter上のmathics3実行画面

アレコレたいへんのとき，こんなもんいるかと勢いで消したうちの一つがWolframの Mathematica 的なフリーソフトのmathicsだった。そのすぐあとで，Mathics3が登場したというニュースをみかけて，再度インストールする気になったが，前回消したつもりの残骸が残っていてうまくいかない。

中国製の生成AIは，自然科学系が強いという噂を思い出して，Qwenに相談してみた。途中まではよかったのだけれど，jupyterでも使えるようにするところでつまづいてしまい，その度にエラーコードをいれるのだが，にっちもさっちもいかない。相談相手をChatGPTにスイッチするとすいすいとインストールできた。

% brew update

% brew install python@3.12 sqlite3 llvm@18

% python3.12 -m venv venv_mathics3

% source ~/venv_mathics3/bin/activate

> python -m pip install -U pip setuptools wheel

> pip install Mathics3

> mathics3

OK Quit[]

> pip install mathicsscript

> mathicsscript

Mathicscript: 10.0.0, Mathics3 10.0.1

on CPython 3.12.13 (main, Mar  3 2026, 12:39:30) [Clang 21.0.0 (clang-2100.0.123.102)]

Using:

SymPy 1.14.0, mpmath 1.3.0, numpy 2.4.4

cython Not installed, 

No matplotlib installed,

Asymptote version 3.09

Copyright (C) 2011-2026 The Mathics3 Team.

This program comes with ABSOLUTELY NO WARRANTY.

This is free software, and you are welcome to redistribute it

under certain conditions.

See the documentation for the full license.

Quit by evaluating Quit[] or by pressing CONTROL-D.

In[1]:=

> python -m ipykernel install --user --name mathics3-py --display-name "Python + Mathics3"

> jupyter lab

(1) カーネルとして Python + Mathics3 を選択する。

(2) ノートブックの最初のセルで以下を実行する。

%load_ext mathics3_kernel.frontend.jupyter

さらに，旧Mathicsと現Mathics3の違いを尋ねるとおおむね以下のようなことだった。

観点　　　　　　｜ 旧 Mathics　　　　　　　 ｜Mathics3

基本性格　　　　｜ Mathematica風CAS　　　  ｜その後継・再設計版

Python　　　　　｜ Python 2時代的な設計が残 ｜Python 3前提

構造　　　　　　｜ 一体型に近い　　　　　　｜カーネル、CLI、Web、Jupyter等を分離

インストール　　｜ 旧環境依存が出やすい　　｜現代Python環境向け

Jupyter連携　　　｜ 周辺的　　　　　　　　　｜Mathics3-notebook-frontends等で対応

Mathematica互換性 ｜ 部分互換 　　　　　　｜部分互換。完全互換ではない

現在使うべきもの ｜ 基本的には非推奨 　　｜こちらを使うべき

単にコマンドラインでMathematicaを使うのが目的ならば，wolframscriptでよい。もしMathematicaを持っていない人に奨める場合はMathics3もよい。Mathics3では，超越数は3500桁まで，代数方程式の解である冪乗根は7000桁まで可能に見える。例えば，N[Pi+E,3500]やN[Sqrt[3],7000]は可能でそれを超える桁は無理だ。

P. S. なぜか以前は存在していたwikipediaのmathics項目が全言語から消されていた。なぜ？waybackmachineをたどっても，List of computer algebra systems に転送されるだけで，そこにもmathics項目がない。完全に抹殺されている。なぜ？

積分（３）

中嶋慧さんがある積分について語っていた。

$\displaystyle \int_0^1 \int_0^1  \int_0^1  \sqrt{x^2+y^2+z^2} \ dx dy dz $ を求めるというものだ。

ChatGPT-5 Pro で 15分20秒かかったらしい。Hypercube Point Picking とよばれるものだ。

(1) Mathematica 14.3 の場合

Mathematicaが14.3にバージョンアップしていたので，試しにやらせてみた。だめだ。$\displaystyle \int_0^1 \int_0^1  \sqrt{x^2+y^2} \ dx dy $は解けるが，3次元にすると5分かかって undefined ときた。

仕方がないので，1変数ずつ順次やらせてみた。

Integrate[Sqrt[x^2 + y^2 + z^2], {x, 0, 1}]

1/2 (Sqrt[1 + y^2 + z^2] + (y^2 + z^2) ArcTanh[1/Sqrt[1 + y^2 + z^2]])

Integrate[1/2 (Sqrt[1 + y^2 + z^2] + (y^2 + z^2) ArcTanh[1/Sqrt[1 + y^2 + z^2]]), {x,0,1}]

1/24 (8 Sqrt[2 + z^2] + 2 (5 + 9 z^2) ArcCoth[Sqrt[2 + z^2]] +  Log[1 - I z] + Log[1 + I z] 

- 2 Log[1 + Sqrt[2 + z^2]] + z^2 (4 z (ArcCot[z Sqrt[2 + z^2]] - 3 ArcTan[z] + 

 3 ArcTan[(1 + z^2 - Sqrt[2 + z^2])/z]) + 9 Log[1 - I z] + 9 Log[1 + I z] 

- 6 (2 Log[-1 + Sqrt[2 + z^2]] + Log[1 + Sqrt[2 + z^2]])))

Integrate[ 同上, {z,0,1}]

1/72 (18 Sqrt[3] - 3 \[Pi] + 24 ArcCoth[Sqrt[3]] -  4 Log[1351 - 780 Sqrt[3]]) 

NIntegrate[Sqrt[x^2 + y^2 + z^2], {x, 0, 1}, {y, 0, 1}, {z, 0, 1}]

0.960592

数値はあっているが，中嶋さんの答えと表式が異なっている。なぜ？

(2) ChatGPT-5 Thinking の場合

9分24秒かけて

$\displaystyle \int_0^1 \int_0^1  \int_0^1  \sqrt{x^2+y^2+z^2} \ dx dy dz =\dfrac{\sqrt{3}}{4} -\dfrac{\pi}{24} - \frac{1}{2} \log 2 + \log(1+\sqrt{3})$ 

この食い違いについてChatGPT-5に尋ねてみると，

まず，${\rm ArcCoth}[\sqrt{3}] = \frac{1}{2}\log \dfrac{\sqrt{3}+1}{\sqrt{3}-1} = \log(\sqrt{3} + 1)-\frac{1}{2} \log{2} $

つぎに，${\rm Log}[1351 - 780 {\rm Sqrt} [3]] = \log (2-\sqrt{3})^6 = 6 \log (2-\sqrt{3}) $

$= 6 \log \dfrac{1}{2+\sqrt{3}} = -6 (2 \log(\sqrt{3}+1) - \log 2) $

$\therefore 1/72 (24 {\rm ArcCoth}[{\rm Sqrt} [3]] -  4 {\rm Log}[1351 - 780 {\rm Sqrt} [3]])  =  \log(\sqrt{3}+1) - \frac{1}{2} \log 2 $

これで無事に一致することが確かめられた。

パデ近似（３）

パデ近似（２）からの続き

では，$\log x$ の近似式はどうして導かれたのか。ChatGPT o3に教えてもらうことにした。

(1) $x=t^2$ とおいて，$f(x) = \log x = 2 \log t = g(t)$

(2) $h=t-1$ とおいて，$t=1 \ (x=1)$のまわりでパデ近似を行う。

(3) テイラー展開すると， $g(t) = 2 \log(1+h) = 2 (h - \frac{1}{2} h^2 + \frac{1}{3} h^3 -\frac{1}{4} h^4  + \cdots)$

(4) これを $h$について パデ近似して$P_{2,2}(h)$ を求めると$P_{2,2}(h) = \dfrac{2h+h^2}{1+h+\frac{1}{6} h^2}$

(5) これを $h=t-1$ で $t$の関数に戻すと，$g(t) \approx \dfrac{6(t^2-1)}{1+4t+t^2}$ となる。

(6) さらに $t=\sqrt{x}$ で $x$の関数に戻すと，$f(x) \approx \dfrac{6(x-1)}{1+4 \sqrt{x} + x}$となる。

こうして$\log x \approx  \dfrac{6(x-1)}{1+4 \sqrt{x} + x}$ が得られた。なるほどね。

じゃあ，次数を上げて計算し，自分が本当に理解したかどうかチェックしてみよう。

(7) $g(t) = 2 \log(1+h) = 2 (h - \frac{1}{2} h^2 + \frac{1}{3} h^3 -\frac{1}{4} h^4  + \frac{1}{5} h^5 -\frac{1}{6} h^6 + \cdots)$

(8) $P_{3,3}(h) = \dfrac{2h+2h^2+\frac{11}{30}h^3}{1+\frac{3}{2}h+\frac{3}{5}h^2+\frac{1}{20}h^3}$

(9) $g(t) \approx \dfrac{2 (-11 - 27 t + 27 t^2 + 11 t^3)}{3 (1 + 9 t + 9 t^2 + t^3)}$

(10) $f(x) \approx \dfrac{2 (-11 - 27 \sqrt{x} + 27 x + 11 x \sqrt{x}}{3 (1 + 9 \sqrt{x} + 9 x + x \sqrt{x})}$

グラフをみるといちおうもっともらしい。

図：$\sqrt{x}$を含んだパデ近似

実はMathematicaにはパデ近似を求める関数があるので，原点以外でも簡単に求まる。

PadeApproximant[2(t-1)-(t-1)^2+(2(t-1)^3)/3-(t-1)^4/2, {t, 1, {2, 2}}] 

まあ，さきほどのチェックでも連立方程式を解いてもとの変数に戻すところはMathematicaの力を借りている。3元連立方程式を間違えずに解くのは筋トレ以上にたいへんだ。そういえば，中学生のときは2元連立方程式がちょっと苦手だったかもしれない。

PLUコード

バナナについてるバーコードで，有機農法か無農薬かを見分けられるというので，そんなバカなことはないでしょうと，鼻であしらった。どうせ，環境保護派のフェイクニュースに違いないと思ったからだ。

それでも念のために調べてみると，なんとそれはいくぶん正しい説だった。ちなみに，スーパーマーケットで自分でスキャンするバーコード（広義のPLUコード = Price Look-Up Code）は，JANコード（日本の場合）と呼ばれるコードだ。これは今回の鍵となる狭義のPLUコードとは直接関係ない。以下に[2]から狭義のPLUコードの定義を引用する。

① PLU コード（番号）

通常バラで流通する野菜や果物などには、流通段階での商品管理を行うため、4

桁または 5 桁の PLU コード（Price Look-Up code）が付されることがある。この

番号は、野菜や果物の種類を示すもので、種類によっては大きさ（重さ）などでさらに細分される場合もある（基本的にはどの生産者が出荷する場合でも同じ種類であれば同じ番号）。

この番号を管理しているのは、IFPC（International Federation For Produce

Coding）という団体で、米国、カナダのほか、南米等でも用いられている。番号の使用は任意であるが、卸売業者やスーパーマーケットなどは商品管理のために導入しているところが多く、IFPC によれば、全米に流通する 8 割の野菜や果物にはこの番号が添付されているのではないかということであった（2004 年 10 月 31 日現在）。

PLU コードは、通常、出荷する生産者（団体）または卸売業者がシールを作製し、商品を流通させる際に個々の個体に貼り付けるか、商品が梱包された箱に商品の個数分のシールを同封し、小売店等で個別に貼り付けることが多いが、小売店でシールを作製する場合もある。

番号は、通常４桁で、オーガニック（有機）の場合にはこの 4 桁の数字の頭に

“９”を、遺伝子組換えの場合は“８”を加えることとされている。

IFPC は、PLU コードのシールを作製する場合の注意事項として以下の事項を推

奨している。

○ できるなら枠で囲うようにし、数字の頭には“＃”のマークをつける。

○ 読みやすさの観点から、個々の数字は 14 ポイント以上の大きさで、大きけれ

ば大きいほどよい。文字の幅は、高さの 3 分の１以下であってはならない。数

字は背景と対照的な色合いとすること（理想は、白地に黒の数字）。  

図：PLUコードの例（Monique Shouten から引用）

さらにびっくりしたのは，Mathematica（Wolfram言語）では，PLUコードを調べる関数が備わっていることだ。まだまだ知らないことが多すぎる世界だ。

図：Mathematicaの例（こんな使い方ができるのか！）

［１］International Federation for Produce Standards

［２］平成 16 年度 農林水産物貿易円滑化推進事業 貿易情報海外調査報告書

［３］PLU(Price Look Up)コードを調べる（Wolfram言語）

火星19

弾道ミサイルの軌道（５）からの続き

ロイターが伝えた朝鮮中央通信のニュースによれば，10月31日に発射したのは大陸間弾道ミサイル「火星19」だとのこと。ICBMの最終完結形態で世界最強のミサイルをうたっている。知らんけど。飛行距離は1001.2km，飛行時間は5156秒，最高高度は7687.5km だ。

無意味なアラートを出すよりも，こういった基本的な情報がほしいところだ。しかし，防衛省が韓国軍や米軍からこれらの情報を入手してもことさらにマスコミに流す必要はないと考えているのだろう。日本政府にしてみれば，アラートを出すだけで危機意識の醸成という目的は達成される。

2年前の「火星17」については，同じロフテッド軌道で，飛行距離 999.2 km ，飛行時間 4135 秒，最高高度 6040.9 km だった。これに比べると，今回は飛行時間が1000秒伸びて，最高高度も1600kmほど高くなっている。

そこで，忘れかけていたMathematicaのコードを引っ張り出して再計算してみた。目の子でパラメータサーチをした。単純なモデルなので，そんなにピッタリ再現できるわけではないが，ほぼほぼ雰囲気は得られる。飛行時間を5156秒=86分あたりに固定すると，最高高度は170kmほど高くなって7850km（+2.1%）になってしまう。最高高度を7600km台にするためには，飛行時間を2分ほど減らして84分にする（-2.3%）必要がある。燃焼時加速度a，燃焼時間τ，燃料重量比pなどを変えても，飛行時間で制約がかかる範囲では，飛行高度を再現するには至らなかった。

これらのパラメタのままで，投射角度だけを45度にすれば（飛行時間Tは，地表に到達する点で調整する），大円到達距離は16500kmを越えてくる。

g = 0.0098; R = 6350; τ = 90.8; p = 0.75; a = 0.045; s = 86.7 Degree; T = 5160;

g = 0.0098; R = 6350; τ = 90.8; p = 0.75; a = 0.045; s = 45 Degree; T = 7050; 

(* g=0.0098; R=6350; τ=90.4; p=0.75; a=0.045;s=86.7Degree;T=5160-120; *)

(* g=0.0098; R=6350; τ=90.4; p=0.75; a=0.045;s=45 Degree;T=6850; *)

fr[t_, τ_] := a*Sin[s]*HeavisideTheta[τ - t]

ft[t_, τ_] := a*Cos[s]*r[t]*HeavisideTheta[τ - t]

fm[t_, τ_] := -p/(τ - p*t)*HeavisideTheta[τ - t]

sol = NDSolve[{r''[t] == -fm[t, τ]*r'[t] + h[t]^2/r[t]^3 - g R^2/r[t]^2 + fr[t, τ], 

r[0] == R, r'[0] == 0,  h'[t] == -fm[t, τ]*h[t] + ft[t, τ], h[0] == 0}, {r,  h}, {t, 0, T}]

f[t_] := r[t] /. sol[[1, 1]]

d[t_] := h[t] /. sol[[1, 2]]

Plot[{6350, f[t]}, {t, 0, T}]

Plot[{f[t + 1] - f[t], d[t]*R/f[t]^2, d[t]/f[t]}, {t, 0, T}, 

 PlotRange -> {-8, 8}]

{f[T/2] - 6350, tyx[T], v[T]/0.340}

→ {7851.42, {5160, -46.3658, 1007.58}, 24.5585}

g0 = ParametricPlot[{NIntegrate[R d[t]/f[t]^2 , {t, 0, tt}], f[tt] - R}, {tt, 0, T}]

→　{8231.43, {7050, 121.626, 16657.9}, 25.6423}

g9 = ParametricPlot[{NIntegrate[R d[t]/f[t]^2 , {t, 0, tt}], f[tt] - R}, {tt, 0, T}]

図：北朝鮮のミサイルの軌道（ロイター記事から引用）

OpenAI o1（４）

OpenAI o1（３）からの続き

引き続き，o1 のプログラミング能力について検証してみよう。数学体感教室2024の問題3を一般化したものについて，o1-mini に Julia プログラムを作成させた。1回目は関数定義が重複するエラーが出ていたので，そのエラーコードを読ませると，2回目には正しい答えが得られた。

Ｑ（koshix）：

$a_i \ (i=1 ... n) $が正の整数であり，$n$個の$a_i$の組 $(a_1, a_2,  ..., a_n)$ 最大公約数は1であるとする。

$a_1/a_2=k_1,\  2*a_2/a_3=k_2, \ 3*a_3/a_4=k_3,\  ... (n-1)*a_{n-1}/a_n \\ =k_{n-1},  n*a_n/a_1 = k_n $において，$(k_1, k_2, ... k_n)$ は整数となって，$k_1*k_2*k_3*...*k_{n-1}*k_n = n!$ が成り立つ。

$n$を与えたときに，上記の条件を満足する$(k_1, k_2, ..., k_n)$ の組の数 $b_n$ を求める Julia プログラムを作成せよ。

このプログラムを用いて，$ b_1, b_2, ..., b_8$ の値を計算せよ。

Ａ（OpenAI o1-mini 思考時間29秒）：

以下に、与えられた条件を満たす $(k_1, k_2, ..., k_n)$の組の数 $b_n$ を計算する Julia プログラムを示します。このプログラムは、各 $n$に対して $n!$ の素因数分解を行い、各素数の指数を用いて組み合わせの数を計算します。

アプローチ

１．素数の列挙：各 $n$ に対して，$n!$ の素因数分解に必要な素数を列挙します。

２．指数の計算：各素数 $p$ について，$n!$ 中の $p$ の指数 $e_p$ を計算します。

これは，$\lfloor \dfrac{n}{p} \rfloor + \lfloor \dfrac{n}{p^2} \rfloor +\lfloor \dfrac{n}{p^3} \rfloor +\cdots $で求められます・

３．組み合わせの計算：各素数 $p$ に対して，$\binom{n+e_p-1}{e_p}$を計算します。これにいおり，$e_p$ 個の指数を $k_i$ に分配する方法の数が求まります。

４．総数の計算：すべての素数についての組み合わせの数を掛け合わせることで、総数 $b_n$ を得ます。

Julia プログラム（修正版）

以下のJuliaプログラムは、上記のアプローチに基づいて $b_n$ を計算します。

# Primes.jl パッケージのインポート

import Primes: primes

# n! における素数 p の指数を計算する関数

function exponent_in_factorial(n::Int, p::Int)::Int

    exponent = 0

    power = p

    while power <= n

        exponent += div(n, power)

        power *= p

    end

    return exponent

end

# b_n を計算する関数

function calculate_b(n::Int)::Int

    if n == 0

        return 1

    end

    

    # n >= 2 の場合のみ素数を取得

    if n >= 2

        primes_list = Primes.primes(2, n)

    else

        primes_list = Int[]  # 空の配列

    end

    

    total = 1

    for p in primes_list

        e_p = exponent_in_factorial(n, p)

        # 組み合わせの計算

        total *= binomial(n + e_p - 1, e_p)

    end

    return total

end

# b_1 から b_8 までを計算して表示

for n in 1:8

    bn = calculate_b(n)

    println("b_$n = $bn")

end

プログラムの説明

・exponent_in_factorial 関数：与えられた $n$ と素数 $p$ に対して，$n!$ の中の $p$ の指数を計算します。

・calculate_b 関数：与えられた $n$ に対して，$b_n$ を計算します。まず，$n$ 以下の素数を取得し、各素数についての組み合わせの数を計算して掛け合わせます。

・メインループ： $n=1$ から $n=8$ までの $b_n$ を計算し，結果を表示します。

実行結果

上記のプログラムを実行すると、以下のように $b_1$ から $b_8$ までの値が得られます。

$b_1 = 1$

$b_2 = 2$

$b_3 = 9$

$b_4 = 80$

$b_5 = 875$

$b_6 = 15876$

$b_7 = 288120$

$b_8 = 7907328$

これにより、各 $n$ に対して条件を満たす $(k_1, k_2, ..., k_n)$ の組の数 $b_n$ を効率的に計算することができます。

このプログラムは，Mathematicaだとより簡潔に表現できる。これは，もとの関数を２つ与えた上で，ChatGPT-4oに作ってもらった。エラーが出ていたが，3回目の修正で正しい答えに到達した。

f[n_] := FactorInteger[n!]

g[n_, k_] := Binomial[n + k - 1, k]

h[n_] := 

 Module[{factors = f[n], product = 1}, 

  Do[product *= g[n, pair[[2]]], {pair, factors}];

  product]

a = Table[h[k], {k, 1, 8}]

ListPlot[Log[a]] 

{1, 2, 9, 80, 875, 15876, 288120, 7907328} 

図：n=1-8 に対するh[n] の対数のListPlot（Mathematica）

２６

$\displaystyle 26 = \Sigma_{n=1}^\infty \dfrac{n^3}{2^n}$ という式が急に現われてちょっとビックリした。自分で考える力が衰えている老人はさっそくPerplexity PROに聞いてみる。

Ｑ（koshix）：

Mathematicaで f[k_]:=Sum[n^k/2^n, {n, 1, Infinity}]  と定義したとき f[k]の一般式を与えてください。 

 Ａ（Perplexity PRO）： 

Mathematicaで定義された関数  f[k]:=Sum[n^k/2^n, {n, 1, Infinity}] の一般式は以下のように与えられます：f(k)=Li_−k (1/2)。ここで、Li_−k (1/2) はポリログ関数（polylogarithm function）を表しています。この結果は、無限級数  ∑_{n=1}^∞ n^3/2^n の閉形式解を示しています。ポリログ関数は、このような形の無限級数を表現するのに適した特殊関数です。ポリログ関数 Li_−k (1/2)は一般に次のように定義されます：　Li_s (z)=∑_{n=1}^∞ n^s/z^n この場合、s=−k および z=1/2 となっています。この一般式は、任意の実数 k に対して成り立ちます。特定の k の値に対しては、より簡単な形に簡略化できる場合もありますが、一般的にはこのポリログ関数の形で表現されます。この結果は、複雑な無限級数を単一の特殊関数で表現できることを示しており、数学的に非常に興味深い結果です。

正しいのだけれど，単に名前がついただけなので理解が深まらない。ChatGPT-4oに別の聞き方をしてみよう。

Ｑ（koshix）：

高校数学までの知識で Σn=1,∞ n^3/2^n の値を導出する方法を教えてください

Ａ（ChatGPT-4o）：

（もろもろ略）

Ｑ（koshix）：

いい線いっているようなのですが，ちょっと違うのではないですか。

Ａ（ChatGPT-4o）：

（もろもろ略）

Ｑ（koshix）：

惜しいですね，もう一息ではないでしょうか。

Ａ（ChatGPT-4o）：

（もろもろ略）

しかし、ここでの正確な計算をすると、最終的には $14$ ではなく、 ${\sum }_{n=1}^{\mathrm{\infty }}\frac{n^3}{2^n}$ の値は正しくは $\boxed{{\displaystyle 26}}$ です。

導出過程の考え方は正しかったのだけれど，肝腎の代数計算が間違っている。何度かやり直してもうまく計算できず，最終的には，答えだけ無理矢理あわせてきた。夏休みの宿題に四苦八苦しているスーちゃん状態である。

ChatGPT-4oに提案された方法は，$|x| < 1 $として，$\displaystyle \Sigma_{n=1}^{\infty} x^n = \dfrac{1}{1-x}$という恒等式から出発し，この両辺を何度か微分して変形した式に $x=\dfrac{1}{2}$を代入するというものだ。

一回微分して両辺に$x$をかけると，$\displaystyle \Sigma_{n=1}^{\infty} n x^n = \dfrac{x}{(1-x)^2}$

この手続きをもう一度繰り返すと，$\displaystyle \Sigma_{n=1}^{\infty} n^2 x^n = \dfrac{x(1+x)}{(1-x)^3}$

さらにこれをもう一度繰り返すと，$\displaystyle \Sigma_{n=1}^{\infty} n^3 x^n = \dfrac{x(1+4x+x^2)}{(1-x)^4}$

これに，$x=\frac{1}{2}$を代入すると，$\displaystyle \Sigma_{n=1}^{\infty} \dfrac{n^3}{2^n} = 26$

都道府県の長さ

ACTIVE GALATTICさんが，次のように書いていた。

中核都市の金沢と被害が激しい輪島は東京と水戸くらい離れていると石川県の広さを説明する投稿を先日見かけたが、石川県の広さは34~35位とランキング下位なので、むしろ東京や大阪が都道府県としては異様に狭いと認識すべきなのだろう

面積と距離が微妙に絡まっていて，ややモニョる。ある図形が与えられたときに，その面積ではなくて，長さを図形の大きさの指標として与えることはできるだろうか。地図上の領域だと東西南北の端点を考えて，南北間距離，東西間距離を抽出できるが，これでは2つの数字の組になる。あるいは図形の周長も考えられるが，フラクタル的な図形だと面積が等しくても周長はいくらでも長くできる。

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

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

まず，方法１のプログラムを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+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を除く。

［１］地域内距離（腰塚武志）

円軌道はむずかしい

万里鏡１号・弾道ミサイルの軌道（２）からの続き

北朝鮮の弾道ミサイルの簡単なシミュレーションコードをMathematicaで作っていた。これを少しアレンジすれば人工衛星を軌道に投入するところまでできそうな気がする。

早速，以前のコードを修正してみた。まずは通常の加速直後に角度方向だけに加速度を加えるようにしたがうまくいなかい。打ち上げ加速は投射角の方向になっているので，動径速度成分が大きく残っているうえ，加速すれば軌道は膨らむ。このため，離心率の大きな長円軌道になって地表にぶつかってしまうか，地球の重力圏から脱出してしまうのだ。

次に，打ち上げ加速の直後に空白時間をおいて，動径速度成分が小さくなったところで角度方向に加速できるようにした。それでもうまくいかない。簡単な試行錯誤では周回軌道にのせるのが難しい。そもそも角度方向に加速するということは面積速度すなわち角運動量をふやし，動径方向の微分方程式で軌道半径を膨らませる方向に作用してしまう。

そこで，後期加速では衛星をその速度ベクトルの方向に加速することにした。$t=0$で速度ベクトルをとりだすところに発散があったので，これを回避するため，地球の自転による2倍面積速度$h(t)$の初期値として，$h(0) = R^2\omega=R^2 \frac{2\pi}{24*3600}=2930$ km$^2/$sを与えた。初期加速度はこれまでの$\ a=0.0446$として$30$秒加速する。その後，800秒程度休止した後に，後期加速度$\ b=0.1445$（ここを微調整した）で$250$秒加速すると，なんとか軌道に投入することができた。投射角は$s=45$度，初期加速における燃料比は$p=0.85$であり，加速方向の角度には$0.3$をかけて動径成分を抑えた。

なかなか難易度の高いゲームである。衛星の軌道高度が1200km程度の準円軌道となっている。これを500kmにしなさいといわれても，こんな単純な2段階制御ではちょっと難しい。なお，プログラムの検証のため，$r=6850$kmの宇宙空間で第一宇宙速度に相当する$v=\sqrt{gr}=8.2$ km/sを角度方向の初速度として与えると，正確に円軌道を描くことが確かめられた。

g = 0.0098; R = 6350; τ = 30; τs = τ*27; τt = 250; p = 0.85;

 a = 0.0446; b = 0.1445 a; s = 45 Degree; T = 15400; 

fs[t_] := 0.3*ArcTan[r[t]*r'[t]/ h[t]]

fr[t_, τ_] :=  a*Sin[s]*HeavisideTheta[τ - t] + 

   b*Sin[fs[t]]*HeavisideTheta[t-τs-τ]

   *HeavisideTheta[τ+τs+τt-t]

ft[t_, τ_] :=  a*Cos[s]*r[t]*HeavisideTheta[τ - t] + 

   b*Cos[fs[t]]*r[t]*HeavisideTheta[t-τs-τ]

   * HeavisideTheta[τ+τs+τt-t]

fm[t_, τ_] := -p/(τ - p*t)*HeavisideTheta[τ - t]

sol = NDSolve[{r''[t] == -fm[t, τ]*r'[t] 

   +h[t]^2/r[t]^3 -g R^2/r[t]^2 +fr[t, τ],

   r[0] == R, r'[0] == 0, 

   h'[t] == -fm[t, τ]*h[t] + ft[t, τ], 

   h[0] == 2930 + 0*Sqrt[g] R^(3/2)}, {r, h}, {t, 0, T}]

f[t_] := r[t] /. sol[[1, 1]]

d[t_] := h[t] /. sol[[1, 2]]

Plot[{6350, f[t]}, {t, 0, T}]

Plot[{f[t + 1] - f[t], d[t]*R/f[t]^2, d[t]/f[t]},

 {t, 0, T}, PlotRange -> {-5, 15}]

 

図：苦労すると有難みがわかる衛星の準円軌道のグラフ

P. S. もう少しがんばると，軌道高度650km（r=6980km）の準円軌道まで達成できた。

g = 0.0098; R = 6350; τ = 25; τs = τ*15.3; τt = 350;  p = 0.85;

a = 0.0446; b = 0.1275 a; s = 45 Degree; T = 15400; 

出生時の体重分布

赤ちゃんの出生時の体重分布がないのか調べてみた。500gきざみはすぐ見つかったのだが，もう少し精度がないかと探したところ，統計でみる日本 e-Stat にあった。人口動態調査 人口動態統計 確定数 保管統計表 出生 ／ 出生数，性・出生時の体重（100g階級）・出生時の身長別 というものだ。

これをExcelでグラフにするが，普通の正規分布と比較しやすくするため，縦軸の100g区間ごとの出生数を全体の出生数に対する割合（%）に直した。なお，これは2021年（最新）のデータであり，男子の（出生数，平均値，標準偏差）は（415,836人，3.047 kg，0.390 kg），女子のそれは（395,687人，2.962 kg，0.358 kg）となっている。

図１：女子（青）と男子（オレンジ）の出生数割合（横軸100gビン数）

正規分布の方は，Mathematicaの組み込み関数を使うのだけれど，100g刻みにしているデータと比較するため，μを10倍，σを10倍にしたものを使っている。

In[1]:= f[\[Mu]_, \[Sigma]_] := 

 Integrate[ Exp[-(x - \[Mu])^2/(2*\[Sigma]^2)]/

  Sqrt[2*Pi*\[Sigma]^2], {x, -Infinity, Infinity}]

In[2]:= f[29.62, 3.58]

Out[2]= 1. - 8.59978*10^-15 I

In[3]:= g[\[Mu]_, \[Sigma]_] := 

 Plot[Exp[-(x - \[Mu])^2/(2*\[Sigma]^2)]/Sqrt[2*Pi*\[Sigma]^2],

 {x, 0, 20 \[Sigma]}, PlotRange -> {0, .12}]

In[4]:= g1 = g[29.62, 3.58];

In[5]:= g2 = g[30.47, 3.90];

In[6]:= Show[g1, g2]

Out[6]= -Graphics-

In[7]:= Plot[

 {PDF[NormalDistribution[29.62, 3.58], x] // Evaluate, 

  PDF[NormalDistribution[30.47, 3.90], x] // Evaluate},

 {x, 0, 50}, Filling -> Axis, PlotRange -> {0, 0.13}]

Out[7]= -Graphics-

ほぼそれらしいような気もするが，微妙に正規分布からずれているのかもしれない。例えば，青の4 kg の位置では累積割合で 99.5% になる。この位置は正規分布では，2.86σ相当に対応する。それで積分すると99.8%になってしまうので，少しずれているようだ。

図２：平均値と標準偏差をそろえた正規分布

折り紙で３次方程式（２）

折り紙で３次方程式（１）からの続き

Mathematicaで折り紙の3次方程式問題をplotする関数を書いてみた。

最初に，Plot[{-x - 3, x + 1, 2 x, (x + 2)^2/4,   Sqrt[8 (x)] - 1, -Sqrt[8 (x)] - 1}, {x, -8, 8},  AspectRatio -> Automatic, PlotRange -> {-8, 8}]　こんな風に問題ごとにプロットしてみた。係数を混同していてなかなか正解にたどりつかなかったが，どうにか安定して答えがでてきた。ただし，接線の方程式は目の子であてはめている。

係数，a,b,c,dを与えて，焦点，準線，放物線，接線をすべて描画するように改善したものが次のコードである。最終ステップで，三次方程式を解いた解が3つリストとして出てくるので，これを実数の個数によって場合分けして描画するのはさぞ面倒だろうと心配していたのだが，リストの成分表示からリストに直しただけで求めていた結果が得られてしまった。虚数解が入ってきても何の問題もなくスキップしてくれた。やはりMathematicaは賢い。

そのコードは次のようなものである。

fold3[a_, b_, c_, d_] := 

{g0 = Graphics[{Point[{b, a}], Point[{d, c}], 

      Line[{{-10, -a}, {10, -a}}], Line[{{-d, -10}, {-d, 10}}]}, 

     Axes -> True]; 

 g1 = Plot[{(x - b)^2/(4*a), Sqrt[4 d x] + c, -Sqrt[4 d x] + c},

    {x, -10, 10}, PlotRange -> {-10, 10}, AspectRatio -> Automatic]; 

 sol = Solve[a t^3 + b t^2 + c t + d == 0, t]; 

 t1 = t /. sol; 

 s1 = -( a t1*t1 + b t1);

 g2 = Plot[t1 x + s1, {x, -10, 10}, PlotRange -> {-10, 10}, 

     AspectRatio -> Automatic, PlotStyle -> Red];

 Show[g0, g1, g2]}[[1]]

この関数で，fold3[1, 2, -2, -1]，fold3[1, 1, 1, 1]，fold3[1, 0, -1, 0.0001] を実行した結果を以下に示す。傾きが0の場合は，極限値として得られることになる。

図：fold3の実行結果

ブランコの物理

Physical Reveiw E に，Initial phase and frequency modulations of pumping a playground swing という日本のグループの論文が掲載された。ブランコを漕いで揺らすコツを解明したというものだ，たぶん。アブストラクトをDeepLにかけてGPT-4で整えると次のようになる。

遊具のブランコは、ブランコ本体とそれを揺らす人間からなる動的な結合振動子系である。本研究では、上半身の自然な動きの初期段階がブランコの連続的なポンピングに与える影響を分析するモデルを提案し、3種類の長さのブランコで10人の参加者が実際にポンピングした際の運動データを用いて検証する。

研究の結果から、スイングが垂直位置（中点）にあり、振幅が小さい時に前方へ移動する際、最大リーンバックの初期段階が発生し、最大振幅のポンピングが予測された。振幅が大きくなるにつれて、最適な初期段階はサイクルの初期段階、すなわちスイングの軌道の後端に向かって徐々に移動する。モデルの予測に従って、参加者全員が、スイングの振幅が大きくなるにつれて上半身の動作の初期位相が早くなることが確認された。

この結果から、スイングする人は、上半身の動作の頻度と初期位相を適切に調整することで、ブランコを効果的にポンピングすることができると考えられる。

人間の体を3つの線質量要素で表現し，上半身と下肢の角度をそれぞれ変数として導入したラグランジアンを考えている。このとき，これらの身体の角度については，何通りかのあらかじめ定められたモデルの周期運動だとしている。

現実に近いのだけれど面倒なモデルなので今一つピンと来ない。そこで，人体を振り子につながって距離が固定された２質点で表し，その間で周期的な質量移動をすることでスウィングするというモデルを考えてみた。質量移動は，位相と角振動数の２つの変数で制御できる。速度に比例する減衰振動のある単振子をこの質量移動で励振するというものだ。

Do[

 s0 = Pi/6*0; w0 = Pi; m1 = 1; m2 = 1;

 λ = Pi/12; ν = 0.05/Pi; ε = 0.95 + i*0.01; δ = Pi/8*0;

 m[t_] := (m1 + m2)/20 * Cos[ε w0 t + δ];

 sol = NDSolve[

       {(m1 + m2) s''[t] == -w0^2 * 

        ( (m1 + m[t]) * (Sin[s[t] + λ + ν * s'[t]) 

        + (m2 - m[t]) * (Sin[s[t] - λ + ν * s'[t])), 

        s[0] == s0, s'[0] == 0}, s, {t, 0, 60}];

  s /. sol[[1]];

  f[t_] := s[t] /. sol[[1]];

  g[t_] := s'[t] /. sol[[1]];

  h[i] = Plot[f[t], {t, 0, 30}, AspectRatio -> 1/2, 

   PlotRange -> {-0.6, 0.6}, PlotStyle -> Hue[0.08*i]],

 {i, 0, 7}]

もとの振動子の初期位相は s0 であり，固有角振動数は w0 (=√g/ℓ)，２質点の距離は角度 λ で表現している。速度に比例した減衰係数は ν，質量移動の角振動数と固有角振動数の比が ε，質量移動の位相のズレが δ である。以下の図は，ε = 0.95 + i*0.01 としたもの（上）と，δ = Pi/8*i; としたもの（下）である。固有角振動数比が 0.97のときが最も励振が大きく，2%ふえるだけで減衰に移行してしまう。

図：振幅の質量移動振動への依存性，振動数（上）と位相（下）

凹面鏡

軸対称な凹面鏡の断面の曲線が放物線ならば，鏡の対称軸＝主軸に平行に入射する光線は反射して焦点に集まることが良く知られている。この曲線が一般の形の場合は主軸近傍では2次関数で近似できるので同様に焦点に光が集まる。近傍から外れた場合に生ずる焦点からのズレが収差となる。ところで，焦点を結ぶのは2次関数に限られるのだろうか？

$x$軸上の原点Oに凹面鏡の底を接地させ，$y$軸上にある焦点Fと原点Oの距離を$f$とする。原点近傍の凹面鏡の断面の曲線上の点P$(x, y(x))$について，直線FPの傾きは$-\frac{f}{x}$と近似できる。したがって，点Pにおける鏡面の傾きは，$\frac{x}{2f}$となる。つまり，$\frac{dy}{dx}= \frac{x}{2f}$という微分方程式が成り立つ。つまり，少なくとも原点近傍では $y(x) = \frac{x^2}{4f}$という2次関数でなければならない。

Mathematicaでのシミュレーションコードを書いてみた。

f2[a_] := a^2

t2[a_] := ArcTan[D[f2[x], x]] /. x -> a 

s2[a_] := Tan[Pi/2 + 2*t2[a]]

g2[y_, a_] := s2[a]*(y - a) + f2[a]

TrigExpand[f2[a] - a*s2[a]] // Simplify

Out[-]= 1/4

p0 = Plot[f2[x], {x, 0, 1}, PlotStyle -> {Red}]; 

p2 = Table[Plot[g2[y, 0.1*n], {y, 0, 0.1*n}], {n, 1, 8}];

q2 = Table[b = 0.1*n; 

   Graphics[ Line[{{b, 1}, {b, f2[b]}}, 

   VertexColors -> {Green, Blue}]],{n, 1, 8}];

Show[{p0, p2, q2}, PlotRange -> All, AspectRatio -> Automatic]

図：凹面鏡（放物面鏡）の断面図と焦点への結像

先ほどの議論を一般化する。凹面鏡の断面の曲線を$y(x)$とする。点P$(x,y(x))$における接線の傾きは，$y'(x)$であり，点Pにおける入射光線と法線のなす角度は，$\theta = \arctan{ y'(x)}$。そこで，点Pにおける反射光線の傾きは，$m = \tan(\frac{\pi}{2} +2 \theta) = -\frac{\cos 2 \theta}{\sin 2 \theta}  = \frac{y'(x)^2-1}{2 y'(x)}$となる。

さて，反射光線の方程式は，$Y-y(x)=m(X-x)$である。したがって，焦点の位置 を$(X,Y)=(0,f)$とすると，$f = y(x) - x \cdot  \frac{y'(x)^2-1}{2 y'(x)}$という微分方程式で定まる。先ほどの，$y(x)=\frac{x^2}{4f}$はこの方程式を満足している。

wolframscript

Wolfram Engine からの続き

ChatGPTは，因数分解ができないので，Wolfram Alphaと連携すればよいという話について。

ChatGPTは数学が苦手である。Stephen Wolframもそう言っている。BingChatに聞いてみると，WolframAlphaのAPIを使う必要があるようだ…というかそんなことは先刻承知。それにはライセンスが必要なので少々ためらう。何か他の方法はないものかと眺めていると，wolframscriptというWolfram言語のコマンドライン版が見つかった。

とりあえずインストールしてみようと（すでに当初の目的から話がずれつつある），確かめてみれば，自分のMacBook Air M1には既にインストール済みだった。以前，JupyterでWolfram言語を使うためにWolfram Engineを導入した時にインストールしていたらしい。うかつなことだ。

最近，Mathematica本体を立ち上げるのが億劫な場合，コマンドラインでフリーのMathicsを使うことが多かった。ただ機能はかなり制限されている。今回，はじめて wolframscript を実行してみると，Mathicsよりは大分制限が少ないというか，インストール済みのMathematica 12.0 の機能がフルに使えるのかもしれない。そこで次のような単独で実行可能なスクリプトも組めるのだった。拡張子は .wls である。

#!/usr/local/bin/wolframscript

(* ::Package:: *)

n=ToExpression[$ScriptCommandLine[[2]]]

Print[RandomInteger[{1,6},n^2]]

コマンドライン引数（文字列変数）が，$ScriptCommandLine[[2]]以下で取得できる。これを整数型に変換するために，ToExpression[] 関数を使った。これで，プログラム外から与えた変数に対して，Wolfram言語での処理が可能になる。

［１］LangChainのWolfram Alpha連携を試す（npaka）

［２］WolframScript（Wolfram言語＆システム ドキュメントセンター）

弾道ミサイルの軌道（６）

弾道ミサイルの軌道（４）からの続き

H3ロケットの発射がフェイルセーフで中止された次の日，なんだか間の悪いことに北朝鮮のICBM級ミサイルが平壌近傍から発射され，北海道渡島大島西200kmに落下した。

飛行時間は66分，飛行距離は900km，最高高度は5700km，落下速度はマッハ20程度と推測される。ペイロードを調整すれば，射程が14,000kmでアメリカ本土が射程に入ると報道された。

それでは，以前作ったのプログラムで確認してみよう。一生懸命Juliaのフォルダを探していたが，Mathematicaのコードだった。昨年11月のコードのパラメータで，飛行時間と飛行距離を微調整した結果が次の通り。

g = 0.0098; R = 6350; τ= 86; p = 0.75; 

a = 0.0446; s = 86.5 Degree; T = 3960; 

(* g = 0.0098; R = 6350; τ = 87; p = 0.75; a

 = 0.0446; s = 45 Degree; T = 5100; *)

fr[t_, τ_] := a*Sin[s]*HeavisideTheta[τ - t]

ft[t_, τ_] := a*Cos[s]*r[t]*HeavisideTheta[τ - t]

fm[t_, τ_] := -p/(τ - p*t)*HeavisideTheta[τ - t]

sol = NDSolve[{r''[t] == -fm[t, τ]*r'[t] + h[t]^2/r[t]^3 - 

 g R^2/r[t]^2 + fr[t, τ], r[0] == R, r'[0] == 0, 

 h'[t] == -fm[t, τ]*h[t] + ft[t, τ], h[0] == 0}, {r, h},

 {t, 0, T}]

f[t_] := r[t] /. sol[[1, 1]]

d[t_] := h[t] /. sol[[1, 2]]

Plot[{6350, f[t]}, {t, 0, T}]

Plot[{f[t + 1] - f[t], d[t]*R/f[t]^2, d[t]/f[t]}, 

{t, 0, T}, PlotRange -> {-4, 8}]

tyx[T_] := {T, f[T] - R, NIntegrate[R d[t]/f[t]^2 , {t, 0, T}]}

v[T_] := Sqrt[(f[T] - f[T - 1])^2 + (R d[T]/f[T]^2)^2]

{tyx[T], v[T]/0.340}

g0 = ParametricPlot[{NIntegrate[R d[t]/f[t]^2 , {t, 0, tt}], 

  f[tt] - R}, {tt, 0, T}]

燃焼時間τが86秒，燃料重量比が0.75，燃焼加速度が 0.0446km/s^2，投射角が86.5度である。これで飛行時間を与えると，飛行距離と最高高度と落下時速度が概ね再現できる。万有引力は距離の2乗に反比例し，コリオリ力や空気抵抗は無視，地球の形を考慮して2次元極座標で質量が変化する1段ロケットの微分方程式を解いている。

燃焼時間τを87秒にして，投射角を45度にすれば，飛行時間が5100秒で飛行距離は14,600km，落下時速度はマッハ24になる。このときペイロードはほとんど変化させていない。燃料重量比は同じで燃焼時間を1秒のばしただけだ。

図：投射角86度のロフテッド軌道と投射角45度の弾道軌道

軸対称電荷分布（５）

軸対称電荷分布（４）からの続き　　

クーロンポテンシャルを場合分けによって書くことができたが，実際に計算する場合には，いちいち分けるのは面倒だ。といっても高々六類型なのでどうってことはないといえばない。

軸対称電荷分布（２）でやったように，ヘヴィサイドの階段関数を使えば，$r_<,\ r_> \ $が表現できるので，クーロンポテンシャルのルジャンドル成分$\ \widetilde{V}_n(r)\ $を一つの関数を使って表わせる。そこで，$R(t)=\dfrac{1}{\sqrt{1+t^2}}, \ a=1.0,\ c= 0.707107$の場合について具体的に計算してみた。

rs[s_, r_] := s* HeavisideTheta[r - s] + r*HeavisideTheta[s - r];

rl[s_, r_] := r* HeavisideTheta[r - s] + s*HeavisideTheta[s - r];

R[t_] := 1/Sqrt[1 + t^2];

v[n_, r_] := 3*NIntegrate[ 1 + 0 * LegendreP[n, t] 

    *rs[s, r]^n * s^2/rl[s, r]^(n + 1), {t, 0, 1}, {s, 0, R[t]}]

v0 = v[0, 0]

2.64412

Plot[{v0, v[0, r]}, {r, 0, 3}]

Plot[{v0, v[2, r]}, {r, 0, 3}, PlotRange -> {2.55, 2.65}]

Plot[{v0, v[4, r]}, {r, 0, 3}, PlotRange -> {2.64, 2.665}]

Plot[{v0, v[6, r]}, {r, 0, 3}, PlotRange -> {2.638, 2.648}]

ここで，被積分関数の中に 1 を加えたのは，積分値が 0 になって計算が不安定になるのを防ぐためである。その部分の定数$\ v_0\ $を取り除くために，ルジャンドル関数の前に 0 をかけて積分することで$\ v_0\ $を求め，実際の計算ではこれを取り除く。

その結果は次のようになった。$\widetilde{V}_4(r),\ \widetilde{V}_6(r)$が，$r<c\ $でゼロになるのは，$\int ds\ $からでてくる$\ R(t)^{n-2}\ $がルジャンドル関数の$n$次未満の線形結合になるためである。

図：クーロンポテンシャルのルジャンドル成分 $\widetilde{V}_n(r)\  (n=0,2,4,6)$

軸対称電荷分布（３）

軸対称電荷分布（２）からの続き　

今から25年前，1997年度の出世・松井・渡会の卒業論文では，変形核の軸対称電荷分布の形として原点からの距離を$R(\cos\theta)$として，2種類の関数形を採用した。なお，$\theta$は対称軸である$z$軸正方向からの方位角を表わし，$x-y$平面に対して対称な分布を考えるので，$R(\cos\theta)$は$\cos\theta$の偶関数に限定される。$t=\cos\theta$とおいて，電荷密度$\rho_0$は，$Ze = 4 \pi \rho_0 \int_0^1 \int_0^{R(t)} s^2 ds\ dt $で与えられる。

(1) 回転楕円体：$R(\cos\theta) = \dfrac{ac}{\sqrt{c^2+(a^2-c^2)\cos^2\theta}}$である。

この場合の電荷密度は，$\rho_0=\dfrac{3Ze}{4\pi a c^2}$となる。

(2) 2次ルジャンドル関数：$R(\cos\theta) = a+(c-a) \cos^2\theta = \alpha + \beta P_2(\cos\theta)$である。

ただし，$\alpha = \dfrac{2a+c}{3c},\quad \beta = \dfrac{2(c-a)}{3}$。

この場合の電荷密度は，$\rho_0=\dfrac{105 Ze}{4 \pi (16 a^3+8a^2c + 6 a c^2 + 5 c^3)}$となる。

ところで，ルジャンドル関数は直交多項式であり，回転楕円体の断面の楕円はルジャンドル関数による展開で近似できる。この近似がどのくらい妥当するかを，長軸と短軸の比が1:2や2:1の場合について確かめてみる。楕円の離心率は，$e=\sqrt{\dfrac{a^2-c^2}{a^2}}$で与えられる，$\varepsilon=\dfrac{e^2 a^2}{c^2}$とおくと，$R(t) = \dfrac{a}{\sqrt{1+ \varepsilon t^2}}$となる。まず，$n=4$までの展開係数を求めておく。

$\dfrac{1}{\sqrt{1+ \varepsilon t^2}} \approx 1 - \dfrac{\varepsilon}{2}t^2 + \dfrac{3}{8}\varepsilon^2 t^4 = \alpha P_0 (t)+ \beta P_2(t) + \gamma P_4(t)$

$= \alpha + \beta \dfrac{3 t^2-1}{2} + \gamma \dfrac{35t^4-30t^2+3}{8}=1-\dfrac{\beta}{2}+\dfrac{3\gamma}{8} +\Bigl( \dfrac{3\beta}{2}-\dfrac{15\gamma}{4}\Bigr) t^2+\dfrac{35\gamma}{8}t^4$

これを解いて，$\alpha = 1 -\dfrac{\varepsilon}{6} + \dfrac{3\varepsilon^2}{40}$，$\beta = -\dfrac{\varepsilon}{3} + \dfrac{3 \varepsilon^2}{14}$，$\gamma = \dfrac{3\varepsilon^2}{35}$

a = 0.5; c = 1.0;

(* a = 1.0; c = 0.5; *)

d[a_,c_,t_] := a/Sqrt[1 + (a^2-c^2)/c^2 * t^2]

4 Pi * Integrate[Integrate[r^2, {r,0,d[t,a,c]}], {t,0,1}]

g0 = Plot[d[t,a,c], {t,0,1}, PlotStyle->Hue[0]]

f[n_,a_,c_] := Integrate[LegendreP[n,t]*d[a,c,t], {t,0,1}]

b = Table[f[n,a,c], {n,0,16,2}]

g[m_] := Plot[Sum[b[[n/2+1]]*LegendreP[n,t]*(2*n+1), {n,0,m,2}],

         {t,0,1}, PlotStyle-> Hue[(m+1)/10]]

Show[Table[{g0, g[k]}, {k, 4, 6, 2}], PlotRange -> {0.4,1}]

図：軸対称電荷分布（左 oblate，右 prolate，横軸 $t=\cos\theta$，縦軸 $R(t)$）

扁平率 $f=\dfrac{a-c}{a}\ $が-1〜1/2の範囲では，4次のルジャンドル展開で回転楕円体をほぼ表現できる。6次ならば十分だ。

軸対称電荷分布（２）

軸対称電荷分布（１）からの続き

前回の$z$軸対称で，$x-y$平面にも対称な電荷分布のつくるクーロンポテンシャルは，

$\displaystyle V(r, s) = \dfrac{-e}{4 \pi \varepsilon_0} \sum_{n=0}^\infty  P_n(s)  \int_{0}^{1}  \int_0^{R(t)} \rho(r', t) P_n (t) \Bigl\{   \dfrac{1}{r_>} \bigl (\dfrac{r_<}{r_>}\bigr)^n \Bigr\}  4\pi r'^2 dr' d t\ $であった（ただし，$s=\cos\theta$，$t=\cos\theta'$）。

最初に球対称の場合で，上記の式を確認する。このとき$R(t) =R$， $\rho(r', t) = \rho(r') = H(R-r') \rho_0$である。$H(x)\ $はヘヴィサイドの階段関数，$\dfrac{4\pi}{3} R^3 \rho_0 = Ze $である。

このとき$\ n=0\ $の項だけが残り，$\displaystyle V(r) =  \dfrac{-e}{4 \pi \varepsilon_0}  \int_0^{R} \rho_0  \dfrac{1}{r_>}  4\pi r'^2 dr' $

(1) $r>R\ $の場合：

$\displaystyle V(r) =  \dfrac{-e}{4 \pi \varepsilon_0}  \int_0^{R} \rho_0  \dfrac{1}{r}  4\pi r'^2 dr'  =  \dfrac{-e}{4 \pi \varepsilon_0}  4\pi  \rho_0  \dfrac{1}{r} \int_0^{R}  r'^2 dr'  =   \dfrac{-Ze^2}{4 \pi \varepsilon_0} \dfrac{1}{r} $

(2) $r<R\ $の場合：

$\displaystyle V(r) =  \dfrac{-e}{4 \pi \varepsilon_0}  4 \pi \rho_0 \Bigl\{ \int_0^r   \dfrac{1}{r}  r'^2 dr'  + \int_r^R  \dfrac{1}{r'} r'^2 dr'  \Bigr\} $

$\displaystyle  =  \dfrac{-e}{4 \pi \varepsilon_0}  4 \pi \rho_0 \Bigl\{ \dfrac{1}{r} \dfrac{r^3}{3}  +   \left [ \dfrac{r'^2}{2} \right ]_r^R  \Bigr\} =  \dfrac{-Ze^2}{4 \pi \varepsilon_0} \Bigl\{\dfrac{r^2}{R^3} + \dfrac{3(R^2-r^2)}{2R^3}\Bigl\}$

$\displaystyle  =  \dfrac{-Ze^2}{4 \pi \varepsilon_0 }\dfrac{1}{R} \Bigl\{\dfrac{3}{2} - \dfrac{r^2}{2R^2}\Bigl\}$

これで，球対称一様電荷分布によるクーロンポテンシャルが再現できた。

Mathematicaでは次のように計算する。

v[r_] := -3＊NIntegrate[s^2/

       (r*HeavisideTheta[r-s] + s*HeavisideTheta[s-r]), {s, 0, 1}]

Plot[v[r], {r, 0, 5}, PlotRange -> {0, -1.6}]

図：球対称一様電荷分布のクーロンポテンシャル(R=1.0, V(R)=-1.0)

Mathics（２）

Mathics（１）からの続き

ちょうど2年前の今ごろ，Mathicsについての記事をかいていたが，完全に記憶の彼方に跳んでいた。別の記事を復習していてたまたま目に付いたのである。

Mathicsは，Pythonベースのフリーでオープンソースの汎用計算機代数システムであり，Mathematica®と互換性のある構文と関数を特徴としている。早速，MacBook Air M1で動かしてみると動かない。pip3 install Mathics で簡単に復活した。

Mathics 5.0.2

on CPython 3.10.8 (main, Oct 13 2022, 09:48:40) [Clang 14.0.0 (clang-1400.0.29.102)]

using SymPy 1.10.1, mpmath 1.2.1, numpy 1.23.4

Copyright (C) 2011-2022 The Mathics Team.

This program comes with ABSOLUTELY NO WARRANTY.

This is free software, and you are welcome to redistribute it

under certain conditions.

See the documentation for the full license.

前回より，かなりバージョンがあがっている。 GUIの画面も実装しているようなので試してみたが，これがなかなか難渋する。どうやら，Three.jsというJavascriptによる三次元グラフックシステムを利用しているらしい。これも7-8年前に流行ったようだが全く知らなかった。MathicsのGUI版の方は，Mathics-Djangoという名前で，PythonのウェブフレームワークのDjangoの上で動いていた。

pip3 install Mathics-Django 

cd mathics-django-master/mathics_django

*You have 21 unapplied migration(s). 

python3 manage.py migrate

python3 manage.py runserver 

*error AttributeError: module 'scipy.linalg' has no attribute 'inv'

pip3 install scipy-stack

python3 manage.py runserver

*error sqlite3.OperationalError: unable to open database file

*error django.db.utils.OperationalError: unable to open database file

mkdir ~/.local/var/mathics 

djangoは，作業ディレクトリ（~/mathics-django-master/mathics_django）で，manage.py runserver を実行すると，ローカルのウェブサーバー（ http://127.0.0.1:8000/）が立ち上がる仕組みになっている（エラーが出たら指示に従ってディレクトリに必要なファイルを持ってくる）。最初のエラーは，scipy-stackのインストールで，次のエラーは，データベースの場所の問題だったが，~/.local/var/mathicsというデフォルトのディレクトリを作成することでようやく解決することができた。

早速実行してみたが，相変わらず，計算量の制限は大きく，1500!（10進 4114桁）くらいまでしか計算できない。本物のMathematicaならば，10000000!（10進 65 657 059桁）でも平気で計算できるのだけれど。

なお，List of computer algebra systems には無事に掲載されていた。

図：MacBook Air M1にインストールしたMathics GUI版

弾道ミサイルの軌道（５）

弾道ミサイルの軌道（４）からの続き

NHKの報道によれば，11月18日（金）午前10時14分，北朝鮮のICBMが発射され，午前11時23分に北海道渡島大島の西200kmの海上に落下した。飛行時間は69分，飛行距離は1,000km，最高高度は6,100km，速度はマッハ22とされている。

また，弾頭質量によっては15,000kmを越える射程となり，米国本土を狙えるとあった。しかし，今回の実験もロフテッド軌道だったので射出角は90度にかなり近いはずだ。これをより低い角度にすれば質量を変化させなくても十分到達できると思うけれどどうなのか。とりあえず，自分のモデルで計算して確かめてみよう。

g = 0.0098; R = 6350; τ = 87; p = 0.75; a = 0.0446; s =  45 Degree; T = 5100; 

g = 0.0098; R = 6350; τ = 87; p = 0.75; a = 0.0446; s =  86.2 Degree; T = 4140; 

fr[t_, τ_] := a*Sin[s]*HeavisideTheta[τ - t]

ft[t_, τ_] := a*Cos[s]*r[t]*HeavisideTheta[τ - t]

fm[t_, τ_] := -p/(τ - p*t)*HeavisideTheta[τ - t]

sol = NDSolve[{r''[t] == -fm[t, τ]*r'[t] + h[t]^2/r[t]^3 - g R^2/r[t]^2 + fr[t, τ], 

r[0] == R, r'[0] == 0, h'[t] == -fm[t, τ]*h[t] + ft[t, τ], h[0] == 0}, {r, h}, {t, 0, T}]

f[t_] := r[t] /. sol[[1, 1]]

d[t_] := h[t] /. sol[[1, 2]]

Plot[{6350, f[t]}, {t, 0, T}]

Plot[{f[t + 1] - f[t], d[t]*R/f[t]^2, d[t]/f[t]}, {t, 0, T},  PlotRange -> {-4, 8}]

tyx[T_] := {T, f[T] - R, NIntegrate[R d[t]/f[t]^2 , {t, 0, T}]}

v[T_] := Sqrt[(f[T] - f[T - 1])^2 + (R d[T]/f[T]^2)^2]

{tyx[T], v[T]/0.340}

{{4140, 41.9866, 1036.42}, 22.9628}

g0 = ParametricPlot[{NIntegrate[R d[t]/f[t]^2 , {t, 0, tt}],  f[tt] - R}, {tt, 0, T}]

このモデルの仮定は次のようなものである。(1) 距離の二乗に反比例する万有引力が働く極座標での運動方程式を解く。(2) コリオリの力や空気抵抗は考えない。(3) 推進剤の燃焼噴射によって，投射体の速度ベクトルの方向に等加速度（a=0.0446 km/s^2）を発射時から一定時間加える。これにより投射体の質量は燃料が0になるまで一定の割合で減少する。(4) 多段式ロケットは考えず，一段式として全質量に対する燃料比（p=0.75）を与える。このとき，運動方程式は全質量に依存しない。(5) 飛行時間の観測値を与え，投射角度と推進加速時間を主パラメータとして，飛行時間と最高高度の観測値が再現されるように調整する。

結論としては，弾頭質量（というか全質量に対する燃料比）を変化させなくても，射出角を86.2度から45度に変更すれば，飛行時間5100秒（85分）で到達距離が14,600kmになる。なお，ロフテッド軌道の場合の弾道ミサイルの最高速度はマッハ23となり，報告されたマッハ22をほぼ再現した。

図：同一パラメタで射出角と飛行時間だけを変えた場合の弾道軌道

P. S. 11月19日の北朝鮮の朝鮮中央通信によれば，このロフテッド軌道の弾道ミサイルは，新型大陸間弾道弾「火星17」の試験発射であり，最高高度 6,040.9 km ，飛行距離 999.2 km ，飛行時間 4135 秒とのことだ。なお，同じ火星17は，今年の3月24日には，最高高度 6,248.5 km ，飛行距離 1090 km ，飛行時間 4052 秒で飛行している。

［１］【詳報】北朝鮮ICBM級ミサイル発射 北海道の西 EEZ内落下か（NHK）

［２］北朝鮮のミサイル等関連情報（続報第２報，防衛省・自衛隊）

　　（防衛省からは，わかった部分から順に発表されている。発表時刻を入れるべきでは）