Simpson 則で数値積分。

世間並みに，いや，それ以上に GW な日々を過ごしている。ところによっては 5 月 1 日と 2 日は平日扱いであったろうが，私はその 2 日間が休みになっただけでなく，週末の 6 日も併せてまるっと一週間，否，期間の両端の日曜日を含めれば実に 8 日間もお休みという，名実ともにゴールデンなウィークを享受している。

私は GW が始まる前から寝正月ならぬ寝ゴルウィクを過ごすと心に決めており，今のところその計画は順調に実施できている。

とはいえさすがに起きている時間がないわけではないので，森口繁一，伊理正夫，武市正人編『Full BASIC による算法通論』（東京大学出版会）に掲げられたお題で，私にもできそうなものに手を付けて見ることにした。

この本は 1991 年に出版されたもので，「まえがき」によると東京大学教養学部の第 4 学期に，工学部へ進学予定の学生向けに開講された科目「算法通論」のための教科書として作成されたとある。

主要コンテンツは第 I 部と第 II 部で，合わせるとちょうど 15 週分の授業内容になっている。今でも東大は 1 学期の授業週が 15 週なのかどうか知らないのだが，かつての 90 分授業ではなく，1 コマ 100 分で 14 週という流行りの授業形態だとそのまま使用するわけにはいかず，アレンジが必須となる。

それはさておき，かねてよりとある連立微分方程式の数値趣味レーションをするのに，Euler 法を卒業して Runge-Kutta 法に鞍替えしようと願っていたが，ウォーミングアップに定積分の数値計算法の一つである Simpson 則で腕試しをしようと，この書の § 6 関数，6.0 数値積分法に目を付けた次第である。

ところが，この本はどうやら森口先生の『JIS FORTRAN 入門 [上]』第 3 版が手元にないとどうにも使いづらい仕様になっている。そちらのテキストの例題の Full BASIC 版のプログラムを『算法通論』の方で例題として示しているからである。そのため，例題でどのような関数の数値積分を行うかはプログラムを読まないと分からないようになっており，仕方がないので自分で勝手に関数をでっち上げてその近似計算を行うことにした。

使用する処理系はもちろん白石和夫先生の十進 BASIC である。久々なのでバージョンアップしていないか公式サイトを見に行ったら，案の定手持ちのもの (7.8.5.4) より新しいバージョン　 (7.8.6.4) がリリースされていた。

さっそくインストーラなしの ZIP ファイルをダウンロードして展開し，BASIC を立ち上げてヘルプで関数定義の仕方を確認しようとしたところ，目次は表示されるものの肝心の本文が表示されない。

公式サイトに記されたトラブルシューティングを参考に，解凍したフォルダ内の BASIC.chm をダブルクリックして「常に確認する」のチェックを外し，一旦ヘルプと BASIC 本体を閉じて再度 BASIC を立ち上げ直したらヘルプの本文も無事に表示されるようになった。わーい！

今回のプロジェクトの主眼は，『算法通論』に記された Simpson 則の公式の捉え方があっているかどうかの検証である。そこには，和の最初の 3 項と，最後の 4 項が記されているが，中間の一般項がどのような規則で定められているのか，自分で類推するしかない。どうやら区間の両端の f(a) と f(b) は重み 1 で，区間の分割数 n は偶数とし，区間の両端をそれぞれ第 0 項と第 n 項とすると，中間の項のうち，奇数項は重み 4 で，偶数項は重み 2 で足し合わせるらしい。

今現在，私の身の回りには数値積分の公式が記載されている書籍が何冊かあるし，ネットで検索しても山ほどヒットするだろうが，この類推があっているかどうかを，結果が分かっている定積分で検証してしまおうというワケである。

せっかくなので円周率でも求めようと思い，1/(1+x^2) の不定積分が arctan(x) であることを利用することにした。
今この瞬間にようやく気付いたことであるが，0 から 1 までの積分が円周率の 4 分の 1 であるから，すべてが有理数の範囲内で済む話であるが，休みボケのせいか，1/√3 から √3 までの積分を使おうとしており，そうすると無理数の √3 があちこちに出てきてヤダなと日和り，スケール変換をして √3 を一回だけ乗じる形に書き換えた形で計算することにした。この記事を書き終えたら区間 [0,1] での積分にソッコーで直そう・・・。

拙いソースコードは以下の通りである。
初めは分点の個数を 2^10=1024 個で様子を見たのだが，調子が良さそうだったのでお気楽に 2^20 個に増やして実行したら結果が出るまでしばらく待たされた。
単純に考えて，1000 倍の手間が掛かることになったわけだから，計算時間は 1000 倍になったはずである。しかも分点を増やした方が精度が下がるという現象まで起きる始末であった。

原因は私にはさっぱりわからないが，きっとデフォルトの桁数では 2^20 では細かすぎて却って精度が落ちたのだと考え，十進 BASIC の便利機能の一つである 1000 桁モードに変更して実行し直してみた。

ついでに，ヘルプを参照して TIME 関数を応用した実行時間の計測も行うこととした。

参考までに，1000 桁モードで n=2^10 としたときの実行時間は 0.3 秒だったのに対し，n=2^20 だと 39.5 秒かかった。1000 倍になるという理屈は当たらずとも遠からずである。
十進 BASIC の組み込み定数としての PI の値と比較すると，小数点以下 24 桁までしか合っていない。

実はこっそり区間 [0,1] 版に変えてみたのだが，n=2^20 で実行時間が 60.46 秒で，小数点以下 35 桁くらいまでは合っているようだった。

誤差が大きいのか小さいのか私にはまるで見当が付かない。『算法通論』には理論的に予想される誤差の表式も与えられているが，使用した関数の 3 次導関数もしくは 4 次関数を必要とするため，ちょっとやる気が起きない。ただ，今回使用した f(x)=1/(3+x^2) ないしは f(x)=1/(1+x^2) は導関数の漸化式を作れば要領よく計算できそうなので，そういった方面からの検証もよい演習になりそうである。

ちなみに，答え合わせ（？）として『通論』のプログラムを見てみると，私のものよりうんと洗練されたプロの香りが強く漂う（その道の大家である森口先生やそのお弟子さんの手になるものであろうから，こんな感想を述べる事すら恐れ多いことであるが）。

ちらっと頭をよぎったことではあるが，和の中間項はやはり奇数項と偶数項をセットにして足し合わせるのね。そして第 n-1 項と第 n 項は最後に別で処理するわけね。
お手本を参考に手直しすれば実行時間が短くなる可能性を感じる。伝説の大先生方の胸をお借りして改良しようと思う。

というわけで，プログラム名に "draft" という釈明を忍ばせておいた。

REM Simpson's Rule [draft]
LET t0=TIME
DEF f(x)=1/(3+x^2)
LET n=2^20
LET a=1
LET b=3
LET h=(b-a)/n
LET S=f(a)
FOR k=1 TO n-1
   LET p=MOD(k,2)
   LET a=a+h
   LET S=S+(2+2*p)*f(a)
NEXT k 
LET S=S+f(b)
PRINT S*h/3
PRINT 2*SQR(3)*S*h
PRINT PI
PRINT TIME-t0;"秒かかりますた。"
END

[付記：2023 年 5 月 5 日]
戸川隼人『数値計算』（情報処理入門コース7，岩波書店，1991 年）8-2 では Simpson 則（シンプソンの公式と呼んでいる）をより一般的な Newton-Cotes の公式の特別な例として取り上げている。そこでは計算例として f(x)=4/(1+x²) の区間 [0,1] における定積分が挙げられており，あらかじめ関数 f(x) の分子に 4 を掛けておくという実にもっともな取り扱いがなされていた。

上記の私の拙いプログラムをそのように修正して末尾の PRINT 文のうち不要になったものを省いて 1000 桁計算したところ，59.92 秒で小数第 37 位まで正確な数値を得た。分点はもちろん 2²⁰ のままである。

戸川氏の本ではプログラム例は FOTRAN で示されているので，一応それも学んだ上で，拙プログラムを改良するのに参考になる点があれば取り入れたいと思う。