Mathematica 入門

数値計算

まず，Mathematica に数の計算をさせてみましょう．まずは足し算から．

In[1]:= 
123456789+987654321
Out[1]= 
1111111110
掛け算は，2つの数を空白で区切って並べればOKです．

In[2]:= 
142857142857 7
Out[2]= 
999999999999
では次に，100 の階乗を計算してみましょう．

In[3]:= 
100!
Out[3]= 
933262154439441526816992388562667004907159682643816214\
 
  6859296389521759999322991560894146397615651828625369\
 
  7920827223758251185210916864000000000000000000000000
最後の桁まで正しく(? だれか確認してください :-) 求まっています． 1000! でも大丈夫! 上の 100! のところでマウスをクリックして書換えてみてください．

2 の 100 乗なんていうのでも，へっちゃらです．

In[4]:= 
2^100
Out[4]= 
1267650600228229401496703205376
加減乗除以外にも，Mathematica にはさまざまな関数が用意されています． 例えば，10! の平方根は，以下のようにすれば計算できます． ここで，Sqrt[x] という記法で √x を表していることに注意してください．

In[5]:= 
Sqrt[10!]
Out[5]= 
720 Sqrt[7]
√7 が Sqrt[7] となっていますね． Mathematica は結果をできるだけ厳密なままで表現しようとしますので， 結果を勝手に数値に直すことはしません．

例えば，Log[2] (2の自然対数)を計算させてみても，もとのままです．

In[6]:= 
Log[2]
Out[6]= 
Log[2]
次のようにすれば，この結果を数値に直すことができます．

In[7]:= 
N[%]
Out[7]= 
0.693147
ここで，N[x] はxを数値に変換する Mathematica の関数です． また % は直前に実行した計算結果を表します．

関数Nの第2引数に精度を与えることによって， より高い精度の数値を求めることもできます．

In[8]:= 
N[Log[2], 40]
Out[8]= 
0.6931471805599453094172321214581765680755
同様に，円周率(Mathematica では Pi で表す)の最初の100桁を求めたければ， 以下のようにします．

In[9]:= 
N[Pi, 100]
Out[9]= 
3.1415926535897932384626433832795028841971693993751058\
 
  20974944592307816406286208998628034825342117068
NやSqrtのように Mathematica にあらかじめ用意されている関数を 組込関数と呼びます． Mathematica には750近くの関数が組み込まれています．

例えば，n 番目の素数を求める関数 Prime[n] なんていうのも，用意されています．

In[10]:= 
Prime[303]
Out[10]= 
1999

数式処理

まずは，簡単な計算から始めてみましょう．

In[11]:= 
2a + 3a
Out[11]= 
5 a
数式の展開を計算するには， Expand を使います．

In[12]:= 
Expand[(a+b)^5]
Out[12]= 
 5      4         3  2       2  3        4    5
a  + 5 a  b + 10 a  b  + 10 a  b  + 5 a b  + b
因数分解を行うには， Factor を使います．

In[13]:= 
Factor[x^4+64]
Out[13]= 
            2              2
(8 - 4 x + x ) (8 + 4 x + x )
その他，微積分や方程式の求解などさまざまな計算が可能です． たとえば Integrate[expr, x] は式 expr を変数 x の もとで積分し，D[expr, x] は偏微分します． Simplify[expr] は式 expr を簡単化した結果を返します．

In[14]:= 
Integrate[ x^2 Sin[x]^2, x ]
Out[14]= 
   3                                  2
4 x  - 6 x Cos[2 x] + 3 Sin[2 x] - 6 x  Sin[2 x]
------------------------------------------------
                       24
In[15]:= 
D[%, x]			(* 上の結果の偏微分 *)
Out[15]= 
    2       2
12 x  - 12 x  Cos[2 x]
----------------------
          24
In[16]:= 
Simplify[%]		(* 上の結果の簡単化 *)
Out[16]= 
 2       2
x  Sin[x]
方程式を解くには，Solve を使います．

In[17]:= 
Solve[ x^2 - 7x + 3a == 0, x ]
Out[17]= 
       7 - Sqrt[49 - 12 a]
{{x -> -------------------}, 
                2
 
        7 + Sqrt[49 - 12 a]
  {x -> -------------------}}
                 2
以下のようにして，連立方程式もSolveで解くことができます (答の中の I は虚数単位を表しています)．

In[18]:= 
Solve[ {x^3+y^3==1, x+y==2}, {x,y} ]
Out[18]= 
              I          12 + 2 I Sqrt[6]
{{x -> 1 - -------, y -> ----------------}, 
           Sqrt[6]              12
 
               I          12 - 2 I Sqrt[6]
  {x -> 1 + -------, y -> ----------------}}
            Sqrt[6]              12
微分方程式にはDSolveを使います (答の中の C[1], C[2]は積分定数です)．

In[19]:= 
DSolve[ y''[x]-k y[x]==1, y[x], x ]
Out[19]= 
            1       C[1]       Sqrt[k] x
{{y[x] -> -(-) + ---------- + E          C[2]}}
            k     Sqrt[k] x
                 E

2次元グラフィックス

まず簡単なグラフを書かせてみましょう． sin x^2 のグラフを x が -π から π の範囲でプロットします．

In[20]:= 
Plot[ Sin[x^2], {x,-Pi,Pi} ]
Out[20]= 
-Graphics-
Plotにはさまざまなオプションを指定できます． どのようなオプションが指定できるかは，??Plot により表示されます． ためしに色々なオプションを指定してみましょう．

In[21]:= 
Plot[ Sin[x^2], {x,-Pi,Pi},
		Axes->None,                      (* 軸を消す *)
		Frame->True,                     (* 枠を付ける *)
		FrameLabel->{"x","sin(x^2)"},    (* ラベルを付ける *)
		GridLines->Automatic,            (* グリッドを付ける *)
		PlotStyle->{{Thickness[0.03]}}   (* 線の太さの指定 *)
]
Out[21]= 
-Graphics-
実は Mathematica では，任意の2次元図形を表示することが簡単にできます． 詳しい説明は省略しますが，例えば，こんな絵も描けます．

In[22]:= 
Show[Graphics[{
		Circle[{0,0}, 1],
		Disk[{0.2,0.1}, {0.05,0.1}],
		Disk[{-0.2,0.1}, {0.05,0.1}],
		Circle[{0,0}, 0.7, {-0.85Pi,-0.15Pi}] }],
	 AspectRatio->Automatic ]
Out[22]= 
-Graphics-

3次元グラフィックス

3次元グラフィックスも2次元グラフィックスと同様に簡単にプロットできます． まずは z =sin x y のプロットです．

In[23]:= 
Plot3D[ Sin[x y], {x,-Pi,Pi}, {y,-Pi,Pi} ]
Out[23]= 
-SurfaceGraphics-
さらに，PlotPointsオプションで計算する点(デフォールトは15*15)を増やし， ViewPointオプションで視点の座標を指定してみます．

In[24]:= 
Plot3D[ Sin[x y], {x,-Pi,Pi}, {y,-Pi,Pi},
		PlotPoints->30,           (* 30*30点で計算 *)
		Axes->False,              (* 軸を消す *)
		Boxed->False,             (* 枠を消す *)
		ViewPoint->{1.5, -1, 3}   (* 視点の座標を指定 *)
]
Out[24]= 
-SurfaceGraphics-
また，2次元グラフィックスと同様に任意の図形を簡単に表示できます． 例えば，手のあるてるてる坊主(?)を書いてみましょう．

In[25]:= 
Needs["Graphics`Shapes`"];
Show[Graphics3D[{
		TranslateShape[Sphere[0.7], {0,0,1}],
		Cone[],
		RotateShape[Cylinder[0.2,1,8], 0, Pi/2, 0] }],
	Boxed->False, ViewPoint->{2,0,1}
]
Out[25]= 
-Graphics3D-

アニメーション

アニメーションを表示させよう．

図のすぐ右にある線のさらに少し右側をマウスでクリックしてから， Ctrl-Yを押すと絵が動き出します． マウスをクリックすると終わります．

In[26]:= 
Needs["Graphics`Animation`"];
Animate[ Plot[Sin[n x], {x,0,2Pi}], {n,1,6,0.5} ]
In[27]:= 
Needs["Graphics`Animation`"];
Animate[
	Plot3D[Sin[n x] Sin[n y], {x,0,2Pi}, {y,0,2Pi},
			PlotRange->{-1,1} ],
	{n,1,2,0.1} ]

最小2乗法

Mathematica で最小2乗法によりデータを解析する方法を説明します． これが，皆さんの勉強の助けになれば幸いです．

データの入力

まずデータを下のように入れていきます． 新たにデータを入れる場合には，各行の終わりで Return を入力してください． 最後まで入れ終ったところで Shift-Return を入力します．

In[28]:= 
data = {
	{0.05, 0.04},
	{0.10, 0.06},
	{0.15, 0.10},
	{0.20, 0.14},
	{0.25, 0.22},
	{0.30, 0.26},
	{0.35, 0.34},
	{0.40, 0.38},
	{0.45, 0.46}, 
	{0.50, 0.58}, 
	{0.55, 0.62},
	{0.60, 0.74},
	{0.65, 0.82},
	{0.70, 0.86},
	{0.75, 0.94},
	{0.80, 1.06},
	{0.85, 1.18},
	{0.90, 1.22},
	{0.95, 1.34},
	{1.00, 1.42}
	};
それぞれの行は，x座標の値とy座標の値の組を表します． たとえば最初のデータ点の x座標値は 0.05 で y座標値は 0.04 です．

データのプロット

では，読み込んだデータをプロットしてみます． データのプロットするには， ListPlot 関数を使います．

In[29]:= 
ListPlot[data]
Out[29]= 
-Graphics-
x-座標とy-座標のきざみを同一にするには AspectRatio オプションを Automatic と指定します (Mathematica は通常，グラフの外形 を自動的に整形し， 黄金分割比を持つ横長の長方形になるようにします)．

In[30]:= 
ListPlot[ data, AspectRatio->Automatic ]
Out[30]= 
-Graphics-
これだと点が小さくて見えないので，点を大きくしてみます．

In[31]:= 
ListPlot[ data, AspectRatio->Automatic,
		        PlotStyle->PointSize[0.02] ]
Out[31]= 
-Graphics-
このグラフは後から使うので，変数 graph1 に代入しておきます．

In[32]:= 
graph1 = %
Out[32]= 
-Graphics-

最小2乗法による近似

では，このデータに対して，最小2乗法で多項式を当てはめてみましょう． データを最小2乗法で多項式近似するには，関数Fitを用います． 2次多項式で近似するには次のようにします．

In[33]:= 
Fit[ data, {1,x,x^2}, x ]
Out[33]= 
                                    2
-0.0384737 + 0.882201 x + 0.597403 x
3次多項式で近似するには {1,x,x^2,x^3} とすれば OK です．

近似した結果は後から使うので，関数 f として定義しておきます．

In[34]:= 
f[x_] = %
Out[34]= 
                                    2
-0.0384737 + 0.882201 x + 0.597403 x

グラフのプロット

では，近似した関数のグラフを Plot で描いて見ましょう．

In[35]:= 
Plot[ f[x], {x,0,1}, AspectRatio->Automatic ]
Out[35]= 
-Graphics-
このグラフも後から使うので，変数 graph2 に代入しておきます．

In[36]:= 
graph2 = %
Out[36]= 
-Graphics-
2つのグラフ graph1 と graph2 を重ねて表示するには， Show を使います．

In[37]:= 
Show[ graph1, graph2 ]
Out[37]= 
-Graphics-

誤差の計算

2乗誤差を求めるには次のようにします．

In[38]:= 
Sum[ ( f[data[[i,1]]] - data[[i,2]] )^2, {i,1,Length[data]} ]
Out[38]= 
0.00919605
ここで Sum[ xi, {i,m,n} ] は，変数 i を m から n まで変化させた時の xi の総和，つまり 
          n
          Σ xi
         i=m
を求めています． また Length[data] は data の長さ，つまりデータの個数を求めています． data[[i,1]] は i 番目のデータ点の x 座標値， data[[i,2]] は y 座標値を取り出しています．

[ 神戸大学Mathematicaホームページへ ]