Maxima-asdfの面白いところは、じぶんで何かMaximaのプログラムを作ったら、簡単にGithubで公開して、世の中に貢献できることです。
と言うわけで一つ作ったのが、オイラー・マクローリン総和パッケージです。
GitHub - YasuakiHonda/euler-maclaurin-sum: Support for Euler Maclaurin Summation Formula
オイラー・マクローリン総和公式をご存知でしょうか。
Wikipedia英語版によればオイラーとマクローリンが独立に発見したそうです。積分の台形公式の精密化したものと考えられます。オイラーは級数の和を求めるのに使い、マクローリンは積分を求めるのに使ったそうです。周期的ベルヌーイ多項式\( \overline{B}_{K}\left(x\right)\)が登場するのですが、あまり恐れずに、まず公式をご覧ください。
$$ \notag \sum_{n=N}^{M}{f\left(n\right)}=\sum_{k=1}^{K-1}{\frac{\left(-1\right)^{k+1}\,B_{k+1}\,\left(\left.\frac{d^{k}}{d\,x^{k}}\,f\left(x\right)\right|_{x=M}-\left.\frac{d^{k}}{d\,x^{k}}\,f\left(x\right)\right|_{x=N}\right)}{\left(k+1\right)!}}+\frac{\left(-1\right)^{K+1}\,\int_{N}^{M}{\overline{B}_{K}\left(x\right)\,\left(\frac{d^{K}}{d\,x^{K}}\,f\left(x\right)\right)\;dx}}{K!}+\int_{N}^{M}{f\left(x\right)\;dx}+\frac{f\left(N\right)}{2}+\frac{f\left(M\right)}{2} $$
左辺が適当な関数の和になっています。左辺は積分などが絡む式です。左辺の総和が、右辺がうまく計算できればわかる、という形になっています。
今回は、まずAPIの使い方を説明します。
(%i1) install_github("YasuakiHonda","euler-maclaurin-sum","master")$
(%i2) asdf_load_source("euler-maclaurin-sum")$
これでインストール終了です。maxima-asdfを導入していない場合には、Githubからeuler-maclaurin-sum.macをダウンロードして読み込めばOKです。
APIと言ってもその核心は以下、変数emsに格納された公式そのものです。Maximaには総和も微分も積分も必要な道具は全て揃っているので、この公式をそのまま打ち込むことができます。それにemsという名前をつけただけです。
(%i3) ems;
$$ \tag{%o3} \sum_{n=N}^{M}{f\left(n\right)}=\sum_{k=1}^{K-1}{\frac{\left(-1\right)^{k+1}\,B_{k+1}\,\left(\left.\frac{d^{k}}{d\,x^{k}}\,f\left(x\right)\right|_{x=M}-\left.\frac{d^{k}}{d\,x^{k}}\,f\left(x\right)\right|_{x=N}\right)}{\left(k+1\right)!}}+\frac{\left(-1\right)^{K+1}\,\int_{N}^{M}{\overline{B}_{K}\left(x\right)\,\left(\frac{d^{K}}{d\,x^{K}}\,f\left(x\right)\right)\;dx}}{K!}+\int_{N}^{M}{f\left(x\right)\;dx}+\frac{f\left(N\right)}{2}+\frac{f\left(M\right)}{2} $$
見ればわかることですが、f(x), N, M, Kなどは利用者が自由に決めていくことができます。例えばf(x)=log(x)として、emsを評価してみます。M,Nについて聞かれるので、M>N-1>0と仮定して答えることにします。
(%i4) ems,f(x):=log(x),nouns;
$$ \tag{*} \verb|Is |N-1\verb| positive, negative or zero?| $$
pos;
$$ \tag{*} \verb|Is |N-M+1\verb| positive, negative or zero?| $$
neg;
$$ \tag{%o4} \sum_{n=N}^{M}{\log n}=\sum_{k=1}^{K-1}{\frac{\left(-1\right)^{k+1}\,B_{k+1}\,\left(\left.\frac{d^{k}}{d\,x^{k}}\,\log x\right|_{x=M}-\left.\frac{d^{k}}{d\,x^{k}}\,\log x\right|_{x=N}\right)}{\left(k+1\right)!}}+\frac{\left(-1\right)^{K+1}\,\int_{N}^{M}{\overline{B}_{K}\left(x\right)\,\left(\frac{d^{K}}{d\,x^{K}}\,\log x\right)\;dx}}{K!}-N\,\log N+\frac{\log N}{2}+N+M\,\log M+\frac{\log M}{2}-M $$
元のemsの中のf(x)にlog(x)を代入して、出来る式変形を全て行ったのが上記の式です。
詳細化すれば、計算も進みます。今度はf(x)=log(x), N=1, K=5, M-1>Nとしてemsを式変形させてみましょう。
(%i5) LOGSUM:ems,f(x):=log(x),N:1,K:5,nouns;
$$ \tag{*} \verb|Is |M-1\verb| positive, negative or zero?| $$
pos;
$$ \tag{%o5} \sum_{n=1}^{M}{\log n}=\frac{\int_{1}^{M}{\frac{\overline{B}_{5}\left(x\right)}{x^5}\;dx}}{5}+M\,\log M+\frac{\log M}{2}-M+\frac{\frac{1}{M}-1}{12}-\frac{\frac{2}{M^3}-2}{720}+1 $$
今度は随分と簡単になりました。ここまでの式変形で近次計算はありませんから、上記式は正確に成立します。
n次の周期的ベルヌーイ多項式periodic_bernpoly(x,n)やその積分の項を処理するために便利な道具を2つ定義してあります。まずは、数値計算などを行うために、単純にperiodic_bernpoly(x,n)に計算可能な定義を与えることと、定義を削除することができます。
(%i6) def_pbp();
$$ \tag{%o6} \left[ \overline{B}_{n}\left(x\right):=\mathrm{bernpoly}\left(x-\left \lfloor x \right \rfloor , n\right) \right] $$
これでperiodic_bernpoly(n,x)にfloor(x)関数を使った定義を与えられます。数値計算したりグラフを書いたりすることができます。
この計算可能な定義があると評価の際にperiodic_bernpoly(n,x)は定義式に置き換えられます。それが不便なこともあるので、その定義を削除できるようにしてあります。
(%i7) undef_pbp();
$$ \tag{%o7} \mathbf{done} $$
これで定義が削除されました。
もう一つの仕組みはemsの評価の過程でよく現れる、
(%i8) integrate(periodic_bernpoly(x,n)/t(x),x,a,b);
$$ \tag{%o8} \int_{a}^{b}{\frac{\overline{B}_{n}\left(x\right)}{t\left(x\right)}\;dx} $$
という形の式を、floor(x)関数も使わない、式変形に便利な積分の総和の形に式変形するパターンマッチ関数emsRemIntです。
(%i9) apply1(%,emsRemInt);
$$ \tag{%o9} \sum_{k=a}^{b-1}{\int_{k}^{k+1}{\frac{\mathrm{bernpoly}\left(x-k , n\right)}{t\left(x\right)}\;dx}} $$
Maximaでパターンによる式変形は、apply1(式、パターンマッチ関数)という形で可能です。式全体を見て、パターンにマッチすれば、この書き換えをやってくれます。周期的ベルヌーイ多項式の周期が1であること、その基本周期である区間[0,1]での値は普通のベルヌーイ多項式で定義されることから、この書き換えが可能になります。
パッケージには今回のセッションのノートブックも含まれていますから、maxima-jupyterとmaxima-asdfをインストールしている人は、上記を実行したりちょっと変えて実験したり、といったことがノートブック上で簡単にできます。今回のセッションをノートブック形式でちょっと見たい、という形はこちらのURLからどうぞ。
