微分解析機

機械式の微分解析機では, 変数同士の乗算は, 積分機を2台と加算でするのだと思っていた. ところが図面に乗算卓というのがあり, 積分をせずとも乗算しているらしかった.

幾何学での乗算法はよく知られている.



点Pから半直線を引き, PA=a, PB=bとなる点A, Bをその上に取る. またPから別の半直線を引き, PC=1の点Cを取る. A,B,Cを通る円Oを描き, PCの延長上で円との交点をDとするとPDの長さxがa*bになっている.

これは方冪の定理による. その説明が次の図である.



PO²-AO²を考える. 正三角形POMでPythagorasの定理を使うと
PO²=PM²+OM²
同様に
AO²=AM²+OM²
辺々相引くと
PO²-AO²=PM²-AM²=(PM-AM)*(PM+AM)=PA*PB
つまりまえの図のPA*PB=a*b=PO²-AO²
=PO²-半径²=一定

である.

しかし微分解析機でこれをやっているか疑問であった. ところである文献に微分解析機の使い方の説明があり, 乗算卓の記号が下の図のようになっていた.



元々の記号に赤や青の色はない.

微分解析機からs, t, xの3本の軸で接続されている, s*t=xと出力が出るらしい.

x軸の右にはハンドルのような絵がある. そこで考えてみて, こうなっていると推測した. t軸が廻転すると, 軸に接続されている赤色のT型の足が円の中心を軸に振れる. 一方s軸の廻転で青い棒が左右に移動する. オペレータはハンドルを廻し, 青軸上のカーソルを赤線の上に保つのである. そのハンドルの廻転がx軸で出力される.

その説明が次の図だ.



円の中心がAの点である. tもsもAから左右に計測する. tが増えるとB点は左へ移動し, ABと直角についている腕, ADも一緒に廻転する. EDの長さをxとすると, 三角形ABCと三角形ADEは相似であり, ACの長さを1とすると, 1/t=s/x 従ってx=s*t である.

そういう乗算卓を使っていたとは面白い.

微分解析機

その昔, 微分解析機という低速のアナログ計算機があった. 1階の常微分方程式を解く機械式の積分機を何台かと, 入出力テーブルを持ち, その間の情報は多くの廻転軸で伝達した.

Wikipediaによれば, 原理を発明したのは1876年, James Thomsonで, 実用機の製作は1927年, MITのVannever Bushによる.

我が国にも2, 3台はあったらしい. 神楽坂にある東京理科大学の近代科学資料館にいくと, 展示されている. 私は1950年代に, 東大生産技術研究所にあった大型機を見学したことがある.

ないよりは数段ましの計算機だが, 大阪大学にいらした, 清水辰次郎先生(1887-1992)が微分解析機の使用経験を書かれたものを読むと, これを使うのは, 結構大変であったらしい.

清水先生の話は

d²x/dt²-μ(1-x²)dx/dt+x=0

を積分機2台の微分解析機でどう解いたかというものである. この微分方程式はVan der Polの式といって, 非線形屋さんの好きなものである.

2階の式は2つの変数に分け, 1階の連立にするのが常道である.
dx/dt=y      積分機0番
dy/dt=g(x,y)    積分機1番
where g(x,y)=μ(1-x²)y-x

xの初期値を1, yのを0, μを0.1にして, この結果をxとyとの図にしたのを下に示す.


今ではPCで簡単に図が描けるが, 昔は手間がかかったので, こういうきれいな図が得られると, 額にいれた飾っていた人もいた.


さて, 阪大の微分解析機は積分機が2台あったから, これで解けそうに見えるが, 問題はμ(1-x²)y-xの計算である. 微分解析機の変数は前述のように廻転軸なので, 最後の引き算は差分ギアを使えばよいし, 最初のμの掛け算はμが定数なので, ギア比で実行できる. xの自乗やそれにyを掛けるように, 変数同士の乗算があると, 微分解析機は積の微分の式を使い, 2台の積分機を使うことになるのだが, 阪大には余っている積分機はない.

そこで聞くも涙の話になるのだが, 清水先生は, 固定したμと, ある範囲のxとyとに対して, g(x,y)をあらかじめ計算し, 表か図にした.

適当な初期値から, (例えばx=1, y=0から,) 起動するとき, それに対応するg(x,y)を積分機1番に入力する. するとx, yが計算され, 出力テーブルにプロットされる; その値を読んで1人が叫ぶと, もう1人がそれに見合うg(x,y) を入力するという, 恐ろしいことをやった. この遅れが誤差になるのは承知の上だが, 傾向をみるだけにはそうするしかなかったのであろう.

私の興味は, g(x,y)はどういうものかにあった. 微分解析機では, f(x)のような1変数の入力なら, 入力テーブルにこの関数形を描いておき, 人手でトレースするなり, 生研の機械なら, 光電管で曲線を自動追尾するなりして, 入力出来るが, gは2変数関数なので, 3次元の入力装置を作れば自動化できたかも知れない. それで関数の形が知りたかったのである.

-4≤x,y≤4について, 0.2おきにg(x,y)を描いたのが下である. もちろん3次元には表示できないので, 高さを向こうへ倒して表現してある. つまりy=4は一番上の線で, 基準が赤線であり, x=-4のとき, g(-4,4)=-2と読む. 他の3つの隅の例では, g(-4,-4)=10, g(4,4)=-10, g(4,-4)=2である.


これならなんとか立体模型をつくり, 入力装置も出来たのではないかと想像する.

この後, この図を等高線で描いてみた. 等高線のプログラムは, 天気図を描きたいという人などが, 昔から挑戦しているが, 私は始めてプログラムしてみた.

それが次の図である.


さらに, 横倒しの図と等高線図を重ねたものも示す.



積分だけでなく, こういう描画も楽になったとつくづく思う.

彩色立方体

東急ハンズで1辺が2センチの木の立方体を売っていたので, 早速買い求めた. 50個で500円程度であった.

これに白, 赤, 緑, 橙, 青, 黒の紙を張って30個の彩色立方体をこしらえた.



作り方は次の通り.

まず30個の立方体のそれぞれの各面に, 色の名を鉛筆でつける.
白で1辺が3センチ程度のほぼ正方形の紙を30枚用意する.
糊で立方体の白の面にその紙を張る.
糊が乾いたら, カッターで余白を切り取る.
白の作業が終わったら, 他の色でも同様に行う.

と書けば簡単だが, 180枚についてこの作業をするのは, 一仕事であった.

今回の反省は, やはり緑, 青など色が濃すぎた; 糊が隣りの面にはみ出て汚れたところがあった; カッターの刃が傾き, 稜を削った個所があった; など.

プログラムの結果と合わせるための, 立方体の番号を探すカタログを用意する必要も感じている.

彩色立方体

Martin Gardnerの本に30 color cubesという話題があった.

立方体には面が6つある. それに相異なる6色を塗る. 廻転して同じになるようなものは作らないとすると, 30通りの塗り方がある.

6色を0,1,2,...,5とする. 立方体の面にも名前をつける. 島内剛一先生の「ルービック・キューブと数学パズル」にあったように, 上をtopのt, 下をbottomのb, 北をn, 東をe, 南をs, 西をwとするのは, 常識だろう. これを(t n e s w b)の順に並べることにする.

まずtを基準として色0を割当てる. するとbには, 他の5色のいずれかが来る. そして残りの4色で側面を塗るが, 4色の最小の色をnに割当てると, 後は3色の置換の6通りになり, 都合30通りというわけだ.

こんなプログラムを書いた.

(define (gencubes)
 (let ((t 0) (neswb '(1 2 3 4 5)))
  (apply append (map (lambda (b)
   (let* ((nesw (delete b neswb)) (n (apply min nesw))
    (esw (delete n nesw)))
    (map (lambda (esw) (append (list t n) esw (list b)))
     (permutation esw)))) neswb))))

2行目でtを0, neswbを(1 2 3 4 5)とし, 3行目のmapのlambda変数でそのいずれかをbとする. neswbからbを除いたものがneswで, その最小をn, 残りをeswとした. それの辞書式順の置換を作り, 前に(t n), 後に(b)をappendして30通りを作っている.

((0 2 3 4 5 1) (0 2 3 5 4 1) (0 2 4 3 5 1) (0 2 4 5 3 1)
 (0 2 5 3 4 1) (0 2 5 4 3 1) (0 1 3 4 5 2) (0 1 3 5 4 2)
 (0 1 4 3 5 2) (0 1 4 5 3 2) (0 1 5 3 4 2) (0 1 5 4 3 2)
 (0 1 2 4 5 3) (0 1 2 5 4 3) (0 1 4 2 5 3) (0 1 4 5 2 3)
 (0 1 5 2 4 3) (0 1 5 4 2 3) (0 1 2 3 5 4) (0 1 2 5 3 4)
 (0 1 3 2 5 4) (0 1 3 5 2 4) (0 1 5 2 3 4) (0 1 5 3 2 4)
 (0 1 2 3 4 5) (0 1 2 4 3 5) (0 1 3 2 4 5) (0 1 3 4 2 5)
 (0 1 4 2 3 5) (0 1 4 3 2 5))

これを辞書式順にソートしたものを, 0は白, 1は赤, 2は緑, 3は青, 4は橙, 5は黒として表示したのが, 以下の図である.

左上の0番で説明すると, それぞれは2つの6角形が連接した形になっていて, 左の6角形は立方体を北東の上方から眺めたもので, t,n,eの3つの面がそれぞれ菱形で表示してある. 上にt=0, 右下にn=1, 左下にe=2の面が見える. 一方, 右の6角形は,南西の下方から眺めたもので, 右上にs=3, 左上にw=4, 下にb=5の面がある. wの4とnの1は, 図でつながっているだけでなく, 実物でもつながっている. 立方体の側面は, 1,2,3,4と北から右まわりにつながっている様子が描いてある.

そこでクイズなのだが, 仮にいま話題にした0番を見本とすると, 他の立方体から8個を選び, 2x2x2のちょうど2倍の立方体を作り, その立方体の外方の面を見本と同じ色にし, 立方体の内側で相接している面同士は, 同じ色にするというのである.

通常なら, 積み木を用意し, とっかえひっかえやってみることになるわけだが, そこはやはりプログラムで全解探索せざるべからずである.

そしてやってみた結果が下の図で, 解は2通りあった.

この図の見方は, それぞれの解で, 大型の菱形に配置した4個の立方体が上下2段に描いてある. 上の4個がそのまま上段に, 下の4個が下段になる. それぞれの菱形の下(0と4)が北東, 左(1と5)が南東, 上(2と6)が南西, 右(3と7)が北西を占める. それぞれの解で, 上段の4個(0,1,2,3)のtは0, 下段の4個(4,5,6,7)のbは5, 北の4個(0,3,7,4)のnは1,...のように, 外の6面は4個ずつ見本と同じである. また相接する面も, 例えば上の解の1のb=4は5のt=4で一致しているし, 1のn=5は0のs=5で一致している.

プログラムは通常の深さ優先の探索木を辿るものである. 解の図のように, 8個の立方体を0,1,2,...の順で候補を探す. そのため, ある立方体をあらゆる方向から見た24個を作るプログラムと, (0 1 2 () () ())のようなパターンと, ある立方体の色の順が, 合っているかを見るプログラムを用意した.

探し方は以下の通り.

    探すパターン 決った色
c0 (t n e - - -) s0 w0 b0
c1 (t s0 e s - -) w1 b1
c2 (t  - w1 s w -) n2 b2
c3 (t n w0 n2 w -) b3
c4 (b0 n e - - b) s4 w4
c5 (b1 s4 e s - b) w5
c6 (b2 - w5 s w b) n6
c7 (b3 n w4 n6 w b)

30通りの立体の組を全候補とする. その中から見本を選ぶ. 上の例では0番. 見本からt, n, e, s, w, bを設定する. 全候補から見本を除いたものをc0用候補とし, そこからc0のパターンに合ったものをc0とする. c0が決るとs0, w0, b0も決る. c0用候補からからc0を除いたものをc1用候補とし, そこからc1のパターンに合ったものをc1とする. c1が決るとw1, b1も決る. 以下同様にしてc7まで決ればc0,c1,...c7を出力する.

プログラムはこれだけである.

30種の立方体は, 全て対称的に出来ているから, 0番を見本とする解が得られれば, 1番,...,29番を見本とする解もまったく同様に作れる. 実際プログラムがあるのだから, 全部やってみた. 全く同じような解が2つずつ得られた.

追記 上と下の解で使われている立方体の番号(上の図参照)は, それぞれ

        6               4
    22     7         1    16
       19              12
       12              19
    16     1         7    22
        4               6

であり, 2つの解で同じ組み合せであった.

二進木の置換表現

図のような二進木があるとする. ノードの番号は中順序で付けてある.

二進木だから左リンクと右リンクがある. つまり

   k      1  2  3  4  5  6  7  8  9 10 11 12 13 14 15
LLINK[k]  0  0  1  0  0  4  0  0  8  0  7  6  0 13  0
RLINK[k]  2  0 12  5  0 11  9  0 10  0  0 14  0 15  0

この2組のリンクを1組に圧縮する方法がTAOCP V4F4にあった. permutation representation of binary treeという.

ノードがn個あれば, 左と右でリンクは2nあるわけだが, 意味有るリンクは1を子にするもの, 2を子にするもの, (3はルートだから, 子にするものはない,)4を子にすもの,... という n-1個しかない. あとは子なしのnilのリンクだ. n=15の上の表でも, 0 は16ある.

この0は無駄なので, それを使おうと考えるのが, この方法の味噌である.

ノード番号は中順序なので, LLINK[k]<k かつ LLINK[k]!=0ならLLINK[k]は重要な情報である.

LLINK=0の場所がうまくRLINKに転用出来るかが, 問題である. つまり図で,
1, 3, 4, 6, 7, 9, 12, 14 (1)
から出ている右リンクの情報を,
1, 2, 4, 5, 7, 8, 10, 13, 15 (2)
の左リンクの場所を使って表現出来るかということだ.

ありがたいことに, 右部分木には, 必ずもっとも左のノードがあり, もっとも左だから, 左リンクは空いていて, しかも中順序だと, ノード番号が1だけ多い. つまり上の(2)の列は, 先頭の1を除いて(2)の列より1ずつ多い.

例えはノード3の右部分木は12をルートとするものだが, そのもっとも左のノードは4. ノード4の右部分木のもっとも左のノードは5. ノード6の右部分木は11をルートとするものだが, そのもっとも左のノードは7である.

これを活用することにすると, 3の右リンク12を4に置かせて貰うが, 部分木のルート自身がもっとも左でなければ, いまの場合がそうだが, 4<12になり, 4の右リンク5を5に置かせて貰う場合, 右部分木のルートがもっとも左なので, 5≤5となる.

要するに, こういう情報だけを保存すると,

   k      1  2  3  4  5  6  7  8  9 10 11 12 13 14 15
LINK[k]   0  2  1 12  5  4 11  9  8 10  7  6 14 13 15

で済むのである. 1のLINKが空いているのは, もったいないようにも思うが, 1はこの二進木全体のもっとも左のノードなので, 二進木がスーパーノード0の右部分木であるとすれば, 0の右リンク3を置かせてあげればよく, 従って1のLINKは3とするのが正解である.

   k      1  2  3  4  5  6  7  8  9 10 11 12 13 14 15
LINK[k]   3  2  1 12  5  4 11  9  8 10  7  6 14 13 15

上の二進木を中順序で辿る通常のアルゴリズムは,

(define (traverse k)
 (if (not (eq? (llink k) '())) (traverse (llink k)))
 (display k) (display " ")
 (if (not (eq? (rlink k) '())) (traverse (rlink k))))

(define (llink k)
 (list-ref '(() () 1 () () 4 () () 8 () 7 6 () 13 ())
  (- k 1)))
(define (rlink k)
 (list-ref '(2 () 12 5 () 11 9 () 10 () () 14 () 15 ())
  (- k 1)))

(traverse 3) => 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15

リンクを

(define (link k)
 (list-ref '(3 2 1 12 5 4 11 9 8 10 7 6 14 13 15) (- k 1)))

とした場合のllinkとrlinkは

(define (llink k) (let ((l (link k)))
  (if (< l k) l '())))
(define (rlink k) (let ((l (link (+ (modulo k 15) 1))))
  (if (> l k) l '())))

(traverse (link 1)) => 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15

なるほど.

ドーナッツを截る

前のblogでドーナッツの断面の平面図を描いたが, ドーナッツがどう見えるかは, あれからでは, よく分からない人が多いのではないか.

そこでもう少し違う角度の絵を描いてみることにした.



この図は, 前の最後の絵に似ている. ただ楕円が2つあり, それぞれ直かに楕円として描いた. 青の楕円は, 長軸の上が(0,0.5), 下が(0,-1.5), 短軸の右が(√3/2,0), 左が(-√3/2,0)なので, 座標の原点を0.-0.5に移動し, スケールを(√3/2,1)に設定して, 半径1の円を描いている.

さて, 任意の鉛直面P,P'でドーナッツを切ったとき, その面と交わる楕円の位置は(y+0.5)²+4x²/3=1とy=kxの交点だから, EとFは簡単にえられる.

P,P'がドーナッツを切る2つの円を母線ということにすると, EやFがその母線と交わる点も分かり, 母線の円の中心からの角度も分かる.

例えばAでは青楕円も緑楕円も120度になり, その点でのドーナッツは切断面と接するので, 母線全体は残る. y'の方向では, 青とは0度, 緑とは180度で, ドーナッツの下半分が残り, 上半分は切り取られる.

これを細かい間隔の母線で決定し, それを描いたのが, 次の図である.



最初の図は, 切り取る部分がまだ置いてあって, 切り取る部分の母線には, 色をつけた. 色はxからy'が赤, y'からx'が青, x'からyが緑, yからxが橙になっている.



もう一つの図が上半分をどけた図で, 切断面は破線で描いてある. 隠面消去していないので, 描いてはみたが, やはりちらついて, 納得出来ていない. 母線の一部の円弧は, みな下向きに垂れ下がっていると思って見る必要がある.

ドーナッツを截る

The New Margin Gardner Mathematical LibraryにSliced Doughnutという話題があった. ドーナッツを平面で截ると切口はどうなるかというのだ. 中心を通る水平面で截ると, 同心円になる. 中心を通る鉛直面で截ると, あなの幅だけ離れた2つの円になる. さらにある面で截ると, 切口は交差する2つの真円になる. これはどういう場合かという問題である. ドーナッツの直径は3インチ, 穴の直径は1インチとする.

立体の切口を想像するのは, 結構難しい. 円錐を截って楕円が現れることなど, 証明がなければ直ぐには信じられない.

まぁこんなことだろうと, 解答を見たら, 一応合ってはいたが, これもあてずっぽだったで, もう少し確かめたいと思った.

立体の断面の作図は, 駒場の学生のころ, 得意中の得意だったので, それをやってみたら, 存外簡単であった.

上の図はドーナッツの平面図と側面図である. 側面図に鉛直面で截ったときの, 断面の2つの円が見えるが, それに接する, 赤線の平面で截ると, 断面が真円になるらしい.



中の図が, 断面の描き方である. 側面図の中央に座標原点があるとする. 高さtの水平面で截ったときの, ドーナッツとの境界が, C₀, C₃と, C₁, C₂で, それを半径をする同心円を, 平面図に示した. 高さtの水平面と斜めの断面との交点をxとすると, ドーナッツの切口は平面図のy₀, y₁, y₂やy₃を通る.



従って, tをドーナッツの最下面から最上面まで, 少しずつ動かしながら, yの点を求めてつなげればよい. こういう作業にはPostScriptは適している.

下の図がそうやって描いた断面で, 楕円をしている. もう1本は, 上下対称のはずなので, 省略した. A点は, 截る面がドーナッツに接するところで, 平面図においてAから上を通りBまでは, y>0の点をつなぐ. Bを過ぎ, C,D,Aはy<0の点をとる. Cから下を通ってCまでは, 外側の同心円との交点である.




この楕円の上下の長さは2インチである. 赤線の勾配は30度, 赤線を斜辺をする直角三角形の高さは1なので, 斜辺の長さはちょうどは2インチになり, 真円であることが分かる.