火星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）

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

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

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

弾道ミサイルの軌道（３）・ロフテッド軌道（３）からの続き

11月3日文化の日，午前7時50分に例のけたたましいJ-アラートが新潟県・山形県・宮城県に発出された。7時40分にミサイルが発射されたということだったが，続報では，午前7時48分に日本上空を通過したとなり，その後ミサイルは日本海上空で消失したなどと情報は二転三転する。

防衛省が11月3日の落ち着いた段階で発表したものは次の通り。

北朝鮮は本日7時台から8時台にかけ，ICBM級の可能性があるものも含め，少なくとも3発の弾道ミサイルを，東方向に向けて発射しました。詳細については現在分析中ですが，落下したのはいずれも我が国の排他的経済水域（EEZ）外であり，飛翔距離等については以下の通りと推定しています。

① 7時39分頃，北朝鮮西岸付近から発射し，最高高度約2000 km程度で，約750 km程度飛翔し，朝鮮半島東側の日本海に落下。当該ミサイルはICBM級の可能性があります。

② 8時39分頃，北朝鮮内陸部から発射し，最高高度約50 km程度で，約350 km程度飛翔し，朝鮮半島東岸付近に落下。

③ 8時48分頃，北朝鮮内陸部から発射し，最高高度約50 km程度で，約350 km程度飛翔し，朝鮮半島東岸付近に落下。

なお，日本列島を超えて飛翔する可能性があると探知したものについては，その後，当該情報を確認したところ，探知したものは日本列島を超えず，日本海上空にてレーダーから消失したことが確認されました。

航空自衛隊が設置している警戒管制レーダーのJ/FPS-5（ガメラレーダー），J/FPS-7は弾道ミサイルを探知することができるはずだ。

実際に弾道ミサイルの飛来を探知した場合は，まず航空総隊/BMD統合任務部隊にデータが送信され，防衛省を経て内閣総理大臣および各省庁にデータが送られる。総務省は，全国瞬時警報システム（J-アラート）を通じてマスメディアや地方公共団体に警報を伝達する。国民には，報道や防災無線を通じて屋内退避などの自己防衛が指示される（Wikipedia J/FPS-5から引用）。

ということならば，問題は，Jアラート側というより，航空自衛隊の警戒管制レーダーの精度と分析能力にあることになる（あるいは発出の政治的判断を含めて）。こんな調子で，防衛費を今後5年間で48兆円に倍増しても大丈夫なのかな。

なお，失敗ともいわれているICBM級ミサイルの軌道は高高度のロフテッド軌道に準ずるものだ。前回の4600kmとんだものの高度は1000kmだったので，失敗としてもその2倍に達している。先に作ったコードを用いて，加速パラメタを共通（0.0446 km/s^2）にしておけば，投射角度と燃焼時間を調整することで，成功したロフテッド軌道（最大高度6250 km，到達距離1090 km，到達時間 4000 秒）と，今回の（失敗かもしれない）ロフテッド軌道（最大高度 1980 km，到達距離 760 km，到達時間 不明）をそれぞれ再現することはできそうだった。

このとき，11月3日午前7時39分に発射された弾道ミサイルの飛行時間は27分前後=1650秒となるので，落下時間は午前8時6分ごろである。多段式だとすれば，今のモデルではうまく記述できていないのだが，第1段の落下時刻はこれより早かったと思われる。

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

s = 86 Degree; T = 4000; *)

g = 0.0098; R = 6350; τ = 77; p = 0.65; 

a = 0.0446; s = 84 Degree; T = 1650;

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[TT_] := {TT, f[TT] - R, NIntegrate[R d[t]/f[t]^2 , {t, 0, TT}]}

tyx[T]

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

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

図：11月3日朝の弾道ミサイルのロフテッド軌道

P. S. 1  成功したロフテッド軌道（最大高度6250 km，到達距離1090 km，到達時間 4000 秒）の条件で，投射角度を45度にすると，到達時間 4850秒で到達距離 14300 kmとなる。北米大陸は射程内だ。

P. S. 2  アメリカの軍事偵察衛星からでは発射直後の赤外線を感知できるが，日本のレーダーでは，高度100 kmに達しないと感知できないという説明がサンデーモーニングであった。もしそうならば，上記の11月3日朝のロフテッド軌道が感知できたのは，発射後66秒後となる。

P. S. 3  仮に上記弾道ミサイルの第1段が，同じ条件で燃料噴射後60秒で放物運動に至る物体だとすれば（切り離された第1段とは正確には一致しないが），それは（最大高度 1000 km，到達距離 470 km，到達時間 1100 秒）で飛行することになる（高度100 km に達するのは 発射後 64秒後）。

［１］ミサイル防衛について（防衛省・自衛隊）

［２］弾道ミサイル防衛（平成20年3月，防衛省）

［３］国民保護ポータルサイト（内閣官房）

［４］北朝鮮 弾道ミサイル「火星17型が飛行に失敗か」（NHK）

［５］北朝鮮からの射程距離（地図蔵）

ドーナツ地球（３）

ドーナツ地球（２）からの続き 

12個の球で近似するくらいならば，いっそのことドーナツ型の内部に一様分布する無数の質点モデルをとったほうがよいかもしれない。$\bm{r}$の位置における単位質量に対する万有引力のポテンシャルは，$V(\bm{r}) = \sum_{i=1}^{N} \dfrac{G m_i }{\sqrt{(\bm{r}-\bm{r}_i )^2}}$である。$N$ 個の質点は，質量 $m_i = \frac{M}{N}$，位置ベクトル$\ \bm{r}_i\ $で与えられ，空間に一様分布している。

トーラス（ドーナツ）は，円環の半径が$R$で，断面が半径$a$の円形であるとする。$\ \bm{r}=(x,y,z)\ $とすると，$x=\cos\phi (R + r \cos\theta)$，$ y=\sin\phi (R + r \sin\theta)$，$z = r \sin\theta )$ と表わすことができる。ただし，$\phi$はトーラス中心に対し円環上で$x$軸から測った方位角，$\theta$はトーラス断面円の中心に対し$x-y$平面から測った仰角で，トーラス中心に遠い方から測ったものである。

トーラス断面円の中心からの距離が$r$であり，$(x-R \cos\phi)^2+(y-R \sin\phi)^2 +z^2 = r^2\ $である。そこで，トーラス面上の点は，$(x-R \cos\phi)^2+(y-R \sin\phi)^2 +z^2 = a^2\ $を，トーラス内部の点は不等式$\ (\sqrt{x^2+y^2}-R)^2 + z^2 < a^2\ $を満足する。

m = 1000000; R = 2.0; r = 1.0;

t = Table[{RandomReal[{-3, 3}], RandomReal[{-3, 3}], 

    RandomReal[{-3, 3}]}, m];

(* s = Select[t,(#[[1]]^2+#[[2]]^2+#[[3]]^2)<1 &]; *)

s = Select[t, ((Sqrt[#[[1]]^2 + #[[2]]^2] - R)^2 + #[[3]]^2) < 1 &]; 

smax = Length[s]

n[w_] := 1/Sqrt[w[[1]]^2 + w[[2]]^2 + w[[3]]^2 + 1/smax]

ListPointPlot3D[s, BoxRatios -> Automatic]

図１：トーラス内部に一様分布する質点

p = {4 - z, 0, 0};

d = Table[p - s[[i]], {i, 1, smax}];

w = 1/smax*Sum[n[d[[k]]], {k, 1, smax}];

v[z_] := w

g[1] = Plot[v[z], {z, 0, 8}]

p = {0, 0, z - 4};

d = Table[p - s[[i]], {i, 1, smax}];

w = 1/smax*Sum[n[d[[k]]], {k, 1, smax}];

v[z_] := w

g[2] = Plot[v[z], {z, 0, 8}]

図２：万有引力ポテンシャルの値（左x軸上，右z軸上）

不連続な質点モデルで計算しているため，たまたま観測点が質点と重なると大きなスパイクが現れる。この現象を緩和するため，分母にくる距離の項 n[d[[k]]] において，1/smaxという微小項を加えている。

ドーナツ地球（２）

ドーナツ地球（１）からの続き

前回のモデル式をMathematicaで書くと次のようになる。地球の半径を$R \rightarrow r$としており，重力加速度の大きさは，${\rm [g]}=\dfrac{GM}{R^2}$を単位として得られたものになる。変数$t$が，重力の大きさを検討している大円上の角度になる。

v = Table[{2 r Cos[Pi/6 i], 2 r Sin[Pi/6 i], 0}, {i, 1, 12}];

d = Table[{2 r + r Cos[t], 0, r Sin[t]} - v[[i]], {i, 1, 12}];

n[w_] := w/Sqrt[w[[1]]^2 + w[[2]]^2 + w[[3]]^2]

m[i_] := (1 + Mod[i + 1, 2])/2

f = Sum[n[d[[i]]]*m[i], {i, 1, 12}];

g[t_] := f /. r -> 1.0

Plot[g[t], {t, 0, Pi}]

Plot[{Sqrt[g[t][[1]]^2 + g[t][[3]]^2], 

    180/Pi*ArcTan[g[t][[3]]/g[t][[1]]]}, {t, 0, Pi}];

vp = Table[{{Cos[t], Sin[t]}, {-g[t][[1]], -g[t][[3]]}}, 

    {t, 0., Pi, Pi/30}];

g1 = ListVectorPlot[vp, AspectRatio -> 0.5];

g2 = Graphics[{White, Disk[{0, -0.05}, 0.95, {0, Pi}]}];

Show[g1, g2] 

 

図１：重力加速度の方向依存性（大きさ，x成分，z成分 vs 大円上の角度 t ）

図２：断面における重力加速度ベクトル場の様子（x = -2.0 がドーナツの中心）

ドーナツ地球の外側で8G，内側で2G（ただしドーナツの中心を向く）となった。ということは，断面図が円形の状態では安定な物質分布とならないわけか。

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

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

10月4日の弾道ミサイルは，本来グアムを標的とする射程5,000kmの火星12型相当の中距離弾道弾なので，日本を対象とした射程1,500kmのノドン改良型とは違う。が，面倒なので，前回求めた解の推進加速度を7割程度に落とし加速角度を30°にして，日本を狙ったときの到達時間が適当な値となる解があるかを確認してみた。

g = 0.0098; R = 6350; τ = 60; p = 0.75; a = 0.032; s =  30 Degree;

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, 360}]

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

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

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

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

(f[t] - R) /. FindRoot[D[f[t], t] == 0, {t, 200}]

　86.6527

FindRoot[D[f[t], t] == 0, {t, 300}]

　　{t -> 199.993}

{d[τ]/f[τ], f[τ+1] - f[τ]}

　　{4.14703, 0.971161}

{Sqrt[(d[τ]/f[τ])^2 + (f[τ+1] - f[τ])^2], 

 Sqrt[(d[τ]/f[τ])^2 + (f[τ+1] - f[τ])^2]/.34}

　　{4.25923, 12.5271}

{NIntegrate[R d[t]/f[t]^2 , {t, 0, 200}], 

 NIntegrate[R d[t]/f[t]^2 , {t, 0, 360}]}

　　{651.9, 1305.61}

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

図：弾道ミサイルの軌道（横軸 水平距離 km，縦軸 高度 km）

60秒加速後の速度がマッハ 12.5となり，360秒（これは初期値として与えた）で1300 km離れた日本に到達する。最高高度は90 km弱となる。このあたりは推進加速度とその角度を調整すればなんとでもなる。いずれにせよ，カップ麺ができる3分以内に対応する必要があるのだけれど，Ｊアラートにインプットする情報分析や対応措置は時間内に間に合うのだろうか。