はじめての数式処理ソフト


はじめての数式処理ソフト CD-ROM付 (ブルーバックス) (新書)竹内薫『はじめての数式処理ソフト』(講談社ブルーバックス)

数式処理ソフトを使ったことのない人に、その強力さを垣間見させてくれる入門書。

小説の登場人物にも使わせるほどMathematica好きの竹内さんが、Maximaを題材にした本を出すなんて意外だなあ、と思って読んでみたら、何のことはない。竹内さんはMathematicaで書いていて、それをMaximaに翻訳してもらったとのこと。

ペンローズのねじれた四次元―時空をつくるツイスターの不思議

まあ、Mathematicaは高いからねえ、ブルーバックス読んで興味を持って、ちょっと試すってわけにはいかないし、しょうがないか。「数式処理ソフトを使ったことのない人」向けには十分だと思うよ。実際、『ペンローズのねじれた四次元』の表紙の絵なんかも、Maximaでがんばって描いているし(ちなみに、あの絵はPlotPointsを調整しないで描いているからちょっとおかしくて、ちゃんとやればああいうギザギザにはならない。「調整していない」ということは、前に竹内さん本人から聞いた)。

MathematicaとMaxima、両方使った経験から言うと(と言ってもMaximaはちょっとだけ)、Maximaがフリーソフトウェアだということを考慮しても、Mathematicaに軍配が上がると思う(数式処理ソフトとしてはもちろん、何でもありのプログラミング言語としても)。高いと言っても学生版なら3万円くらいだし(3万円で得られる生産性という意味では驚異的)。

とはいえ、こういうところからMaximaの人気が高まって、びびったWolframがMathematicaの一部でもオープンソース化してくれたらいいなあ、などと思う。ていうか、そういう話があったじゃないですか、どうなったんですか、その後。

Maximaは簡単にダウンロードできるのだから、CDなんて付けないで、その分安くすればよかったのに、と思う。

追記:ついに無料になったMathematica

数値積分、それが目的ではないとき


学生向け

数値積分の簡単な方法を教えてください

と訊かれました(もう自分でできるはずなんだけど)。いつもどおり、低級言語で実装すれば数値計算の勉強になります。答えが欲しいだけなら、Maximaのような高級言語を使いましょう(Mathematicaは使える人が少なそうなので、フリー・ソフトウェアのMaximaでまず実装することにしました)。

例として、f(x)=1+2 Pi sin(2 Pi x)を0<=x<=1の範囲で数値積分します。

/* romberg積分 */
(%i1) romberg(sin(x),x,0,1);
(%o1) 0.45969769038987

/* 別の方法(この問題についてはこっちのほうがいい) */
/* Newton-Cotesの8次多項式による求積法(ルーチンは適応型) */

(%i2) load("qq");
(%o2) C:/PROGRA~1/MAXIMA~1.0/share/maxima/5.10.0/share/numeric/qq.lisp
(%i3) quanc8(sin(x),x,0,1);
(%o3) 0.45969769413186

/* もちろんこの問題は、記号的に解けるわけで */
(%i4) integrate(sin(x),x,0,1);
(%o4) 1-cos(1)

/* 数値にすると */
(%i5) bfloat(%);
(%o5) 4.596976941318603b-1

ブラックボックスで価格も高いけど高機能なMathematicaでは次のようになります(最後の「%」以外はUMMでも動きます)。

In[1]:= NIntegrate[Sin[x], {x, 0, 1}]

Out[1]= 0.459698

In[2]:= NIntegrate[Sin[x], {x, 0, 1},Method->Trapezoidal]

NIntegrate::slwcon:
   Numerical integration converging too slowly; suspect one of the following:
    singularity, value of the integration being 0, oscillatory integrand, or
    insufficient WorkingPrecision. If your integrand is oscillatory try using
    the option Method->Oscillatory in NIntegrate.

Out[2]= 0.459698

In[3]:= Integrate[Sin[x], {x, 0, 1}]

Out[3]= 1 - Cos[1]

In[4]:= %//N

Out[4]= 0.459698

指定できる数値積分法はいろいろあります。(参照:数値積分の実装

高次方程式、それが目的ではないとき


学生向け(WolframAlphaがなかった頃の話です。)

高次方程式を数値的に解く簡単な方法を教えてください

と訊かれました。低級言語で実装すれば数値計算の勉強になります。答えが欲しいだけなら、Mathematicaのような高級言語を使いましょう。次のように、簡単に解けます(UMMでも動きます)。

例:

In[3]:= With[{f = Sum[Random[Integer, {0, 9}] x^i, {i, 1, 8}] - 10},
          {f, NSolve[f == 0, x] // TableForm}]

                        3      4      5      6
Out[3]= {-10 + 2 x + 5 x  + 4 x  + 4 x  + 4 x , x -> -1.43626             }

                                                x -> -0.586233 - 1.04486 I

                                                x -> -0.586233 + 1.04486 I

                                                x -> 0.380996 - 1.13446 I

                                                x -> 0.380996 + 1.13446 I

                                                x -> 0.846732

計算機プログラムの構造と解釈Mathematicaは高いので、Maximaのようなフリー・ソフトウェアで試すとよいかもしれません。高次方程式は次のように説きます。

allroots(-10+5*x+9*x^2+4*x^3+4*x^4+6*x^5+9*x^6+3*x^8);

もう少しアルゴリズムを楽しみたいのなら、Schemeを使うのはどうでしょうか。Mathematicaのような高級言語とJavaのような低級言語の間に位置するLispの方言です(この文脈では)。Structure and Interpretation of Computer Programsに解説があります(オンライン版日本語訳

行列の計算、それが目的ではないとき(Maxima版)


学生向け

Maximaの使い方がわからないのですが

と言われました。

別に私も詳しくないのですが。

MaximaでCayley‐Hamiltonの定理を確認する例:

(%i1) n:5;
(%o1)                                        5
(%i2) h[i,j]:=1+n-max(i,j);
(%o2)                     h     := 1 + n - max(i, j)
                           i, j
(%i3) m:genmatrix(h,n,n);
                               [ 5  4  3  2  1 ]
                               [               ]
                               [ 4  4  3  2  1 ]
                               [               ]
(%o3)                          [ 3  3  3  2  1 ]
                               [               ]
                               [ 2  2  2  2  1 ]
                               [               ]
                               [ 1  1  1  1  1 ]
(%i4) f:expand(determinant(m-diagmatrix(length(m),x)))
/* or f:expand(charpoly(m,x)); */;
                       5       4       3       2
(%o4)               - x  + 15 x  - 35 x  + 28 x  - 9 x + 1
(%i5) sum(coeff(f,x,i)*m^^i,i,0,nterms(f)-1);
                               [ 0  0  0  0  0 ]
                               [               ]
                               [ 0  0  0  0  0 ]
                               [               ]
(%o5)                          [ 0  0  0  0  0 ]
                               [               ]
                               [ 0  0  0  0  0 ]
                               [               ]
                               [ 0  0  0  0  0 ]

はじめてのMaximaはじめてのMaximaが参考になるかもしれません。Knoppix/Mathの紹介ということで、Maxima以外にも無料で(一部は自由に)使える数学ソフトウェアがたくさん紹介されています。

Maximaに限定するならMaxima簡易マニュアルも有用です。

参考:行列の計算、それが目的ではないとき(Mathematica版)