2014年6月30日月曜日

微分解析機

東京理科大学近代科学資料館で企画展「計算する器械たち」を開催中だ.

そこで毎日微分解析機を実演している. デモは例のサークルテストだが, やはり一周すると2パーセントくらいのエラーが出る.

このデモの様子を, 昔微分解析機を使っていらした, 渡辺勝先生にお見せしたら, 「サークルテストでの反転時のバックラッシュが小さいようにお見受けしております」というメイルを頂いた.

ところで手回し計算機や電動計算機が出回ると, Runge Kuttaという計算法で連立m元常微分方程式を解くようになった. 我々がプログラミングを勉強したケンブリッジ大学のEDSACの本の例題にRunge Kutta Gillがあったので, 当時のコンピュータ屋にはお馴染の方法である. この方法については「伊理正夫,松谷泰行,Runge-Kutta-Gill法について,情報処理,Vol.8,No.2,pp.103-107」に解説がある. 今回はそれを使ってサークルテストを計算してみた. Schemeのプログラムは次のようだ.
(define (runge-kutta-gill m h n)
 ;m 変数の個数, h 分点間隔, n 分点数
(let ((y (make-vector m)) (f (make-vector m))
  (q (make-vector m)) (k 0) (qi 0) (s 0) (r 0)
  (c0 (- 1 (/ 1 (sqrt 2)))) (c1 (+ 1 (/ 1 (sqrt 2)))))
 (define (inity)    ;yの初期値
   (vector-set! y 0 0)
   (vector-set! y 1 0)
   (vector-set! y 2 1))
 (define (calcf)    ;fを計算
   (vector-set! f 0 1)
   (vector-set! f 1 (vector-ref y 2))
   (vector-set! f 2 (- (vector-ref y 1))))
  (define (printy)  ;yの出力
   (display (list 
     (vector-ref y 0)
     (vector-ref y 1)
     (vector-ref y 2))) (newline))
  (newline) (inity)
  (do ((i 0 (+ i 1))) ((= i m)) (vector-set! q i 0))
   ;qiの初期化
  (calcf) (printy)
  (do ((j 0 (+ j 1))) ((= j n))
  (for-each (lambda (e c) 
   (do ((i 0 (+ i 1))) ((= i m))
    (set! k (* h (vector-ref f i)))
    (set! qi (vector-ref q i))
    (set! r (e))
    (set! s (vector-ref y i))
    (vector-set! y i (+ s r))
    (set! r (- (vector-ref y i) s))
    (vector-set! q i (+ qi (* 3.0 r) (* c k))))
   (calcf))
   (list (lambda () (- (* 0.5 k) qi))
     (lambda () (* c0 (- k qi)))
     (lambda () (* c1 (- k qi)))
     (lambda () (/ (- k (* 2.0 qi)) 6.0)))
   (list -0.5 (- c0) (- c1) -0.5))
  (printy))))

(runge-kutta-gill 3 (/ (atan 1) 2.5) 20)
普通とちょっと違うのはy[0]がxであることだ. 従ってf(というのはdy/dxのことだが)のf[0]=1にしてある.

yの初期値では, y[0](xのこと)=0, y[1](yのこと)=0, y[2](dy/dxのこと)=1. fの計算では f[0]=1, f[1]=y[2], f[2]=-y[1]とする.

Runge Kuttaでは分点を1個進めるのに, x, x+h/2, x+h/2, x+hについて4回fの値を計算するが, そのループはfor-eachで回している. ループごとに変る値については, (lambda (e c) ..)の引数で渡す. 特にeの値は, 呼ばれた時に計算する必要があるので, 引数の方は(lambda () ..)の形にし, 使う時に呼出す.

initiy, calcf, printfなどのサブルーチンが, 関数runge-kutta-gillの中で定義してあるのは, 問題ありだが, yやf読み書きするので, これで我慢している.

最後の行で呼び出しているが, 0から2πまでを20分割で計算する. つまりh=2π/20; πを4*(atan 1)で計算したいからh=(atan 1)/2.5としてある.

たったの20分割なのだが, 一周したときのy[2]が.999868...だから1万分の2の程度の精度だ.

その出力を貰ってサークルを描くPostScriptのプログラムは次の通り.

/ps[[0 0 1]
[.3141592653589793 .30899155257892935 .9510578492071948]
[.6283185307179586 .5877376828378169 .8090352529734782]
[.9424777960769379 .8089575954451165 .5878333484965556]
[1.2566370614359172 .9510010098334781 .30910245672629366]
[1.5707963267948966 .9999670230159169 1.2303914619280843e-4]
[1.8849555921538759 .9510645042444495 -.30886434562367165]
[2.199114857512855 .8090808881734999 -.5876187580148391]
[2.5132741228718345 .5879134969774285 -.8088585919500648]
[2.827433388230814 .3092092738117766 -.9509316169859521]
[3.141592653589793 2.4607017746574233e-4 -.9999340319806838]
[3.4557519189487724 -.30873714204448566 -.951071143410806]
[3.7699111843077517 -.5874998314987273 -.8091265072362017]
[4.084070449666731 -.8087595818584492 -.5879936306339407]
[4.39822971502571 -.9508622132841098 -.30931607883671747]
[4.71238898038469 -.9999010268957623 -3.6909309235748844e-4]
[5.026548245743669 -.951077766707203 .30860994184321267]
[5.340707511102648 -.8091721101619074 .5873809032915227]
[5.654866776461628 -.5880737494657694 .8086605651723112]
[5.969026041820607 -.3094228718001786 .9507927987297933]
[6.283185307179586 -4.921078894070952e-4 .9998680077626148]
] def
280 360 translate
/mv {200 mul exch 200 exch mul moveto} def
/ln {200 mul exch 200 exch mul lineto} def
ps 0 get dup 1 get exch 2 get mv
1 1 20 {ps exch get dup 1 get exch 2 get ln} for stroke
図はこのようだ.



2014年6月8日日曜日

鶴亀鴉算

小学生のころ鶴亀算を習った. 例えば鶴亀合せて10匹. 足は合せて30本. 鶴と亀はそれぞれ何匹か. 全部鶴とすると足は20本. 残りの10本は鶴より2本足が多い亀によるから, 亀は5匹. 従って鶴も5匹. 足の本数を検算すると2*5+4*5=30.

矢野健太郎先生がアメリカの子供にこのような話をしたら, 鶴と亀とは見ただけで分り, 足だってまったく違うという反応だったらしい. そこで矢野先生は2ドルのケーキと4ドルのケーキを合せて10個買い, 合計を30ドルにするには2ドルと4ドルのケーキを何個ずつ買うかと話したそうだ.

小学6年の私は近くに住む中学1年の先輩と話をしていた. 先輩は「鶴をx, 亀をyとすると, x+y=10, 2x+4y=30. これを解くとx=5, y=5が得られる」と得意そうに話した. へぇーすごいなぁと思ったが, 1年後には私もその先輩と同じ中学(詳しくいえば7年制高等学校尋常科)に入学し, 同じように解けるようになった.

さて東京理科大学近代科学資料館で, 6月19日から8月8日まで「計算する器械たち —アナログコンピュータ展—」が開催される.

それに門脇廉君(現 九州大学)が作った機械式アナログ計算機の「三元連立方程式求解機」も展示される. 私は展示準備中のその機械を見たが, 非常に美しく作られていて感動した.

企画展開催中に説明を担当する理科大の学生さんたちと三元連立方程式の例題を作って解いてみたりしていたが, 見学に来た小学生に連立方程式の意味を教えるにはどうするかということになった.

そこで私が考えたのが, このブログの始めにあった鶴亀算であった. でもそれは二元であるが, もう1種類の変わった動物を登場させようと思った.

それがサッカーの3本足の鴉である. ボールを掴んでいるのが3本目の足だ.



これが三元連立方程式求解機用問題.

t, 亀k, 鴉cとする. ただし鴉はサッカーJFAのシンボル 3本足のものである. 頭の数が6, 羽の数が10, 足の数が16のとき, 鶴, 亀, 鴉それぞれ何匹か.

つまり連立方程式
t+k+c=6 (0)
2t+0k+2c=10 (1)
2t+4k+3c=16 (2)
を解くわけだ.

早速三元連立方程式求解機のシミュレータでやってみる.



4枚のプレートがあり, 左からt, k, c, 定数に対応している. また制約するテープは内側から順に上の式(0), (1), (2)に対応している.

プレートの傾斜は右端の定数のを10°にしたので, 左からの傾斜は31.395°, 10.0°, 20.322°,10°である. つまり(t, k, c)=(3, 1, 2).

この図で分かるのは定数項の値が係数に較べて大きいので, 右端の定数のプーリーが分散しているのに対し, 係数の方のプーリーがごちゃごちゃと固まっていることである.

これを解決するには上の連立方程式の変数のそれぞれから1, または2を引くのである.

1を引くと
t+k+c=3 (0)
2t+0k+2c=6 (1)
2t+4k+3c=7 (2)

2を引くと
t+k+c=0 (0)
2t+0k+2c=2 (1)
2t+4k+3c=-2 (2)

2を引たので解いたのが下の図である. 各プーリーが分散している.


傾斜が左から10°, -10°, 0°, 10°になり, (t, k, c)=(1, -1, 0)であるから, 鶴が3, 亀が1, 鴉が2と判明した.

却って難しくなったかしら.