第 14 回 (2005.08.01) -- Mathematica 入門 2

さて,今回は Mathematica を普通に使うために必要な 「プログラミング」と 「グラフ描画」 についてごく簡単に学ぶ.

プログラミングというと難しく聞こえるが, 「具体的にはどうやって計算するか」をコンピュータに 「具体的に教える」だけの話である. 計算したいものが同じでも具体的な計算方法は異なっても良いことから, プログラムも個人によって異なるのは当然であり, そこに「上手,下手」が現れてくる.
プログラミングに習熟することは「思考過程を言語化」することに慣れるという側面もあるので, 早くから慣れておくのが良い.

グラフについては,gnuplot に既に学習したので不要に思うかも知れないが, 対話計算を主とする Mathematica などでは結果を「その場で目でみる」ことが重要であり, 習熟不足のためにグラフ化が面倒になるようでは Mathematica を使えるとは言えない.

簡単なプログラミングを用いた関数定義

Mathematica ではちょっとややこしい計算は名前をつけて関数にしてしまい, それを単位として扱うのが全体を見通し良くするコツの一つである. それには自分の思考,計算過程をいかにコンピュータの言語で表すかという技術が必要になる. ここでは,ごく簡単に計算の方法を Mathematica で表す方法について述べる.

なお,ごく簡単な関数をどうやって定義するかについては前回の授業の 簡単な関数の定義 を参照すること.

条件判断

まずは条件分岐について触れておこう. 一般に条件分岐は「もし A ならば処理 1 を,B ならば 処理 2 を」というように条件によって動作を変えるのに使われる. これによって,プログラムの動作を制御することができるという意味で初歩の初歩の技術であるので,しっかり学習するように.

条件判断 … If [ 条件, 条件が真の時の式, 条件が偽の時の式 ]

まずは一つだけの条件で,その条件も yes か no という答えしか返さない 一番簡単な例について学習しよう.

例えば,BMI(Body Mass Index): 体重(kg) / (身長(m)^2) を見て太り過ぎかどうかの判定をする関数を考えてみよう. (日本肥満学会では BMI が 25以上で肥満とする)

    In[1]:= BMI[w_, h_] := w/(h^2) ← w:体重(kg), h:身長(m)

    In[2]:= BMI[64, 1.68]
    Out[2]= 22.6757

    In[3]:= himan[x_] := If[ 25 <= x, 
              Print["BMI=",x, " is HIMAN"], 
              Print["BMI=",x, " is not himan"]
              ]

    In[4]:= himan[BMI[64, 1.68]]
    Out[4]= BMI=22.6757 is not himan

    In[5]:= himan[BMI[75, 1.68]]
    Out[5]= BMI=26.5731 is HIMAN
    

ただし,途中で

表示 … Print[表示したいもの(カンマで区切って並べる)]

というコマンドを用いている.

注意 Print コマンドは「画面に出力はする」が,「関数の値にはならない」. つまり,Print コマンドを使った出力は他の関数の入力に使うための再利用が難しい.
よって,Print コマンドはあくまで人間の「確認用」に使うだけにして, 関数の結果出力に直接使わないようにしよう.
そういう意味で,この himan 関数はあまり良くない書き方になっている. 改善方法については,後述の Himando 関数のところで述べているのでそこを見よ.

条件判断(多値) … Which [ 評価式1, 式1, 評価式2, 式2, 評価式3, 式3,... ]

評価式n を順に比較して,最初に True (真)になる最初の n に対して 式n を実行する.
(次の Switch と似たような動作をするので,そちらもみておくこと)

このコマンドを用いる例を考えてみよう. 例えば,日本肥満学会では BMI に応じて

BMI < 18.5 18.5≦ 25≦ 30≦ 35≦ 40≦
分類 低体重 普通 肥満(1度) 肥満(2度) 肥満(3度) 肥満(4度)

と肥満症を分類しているので,この分類を行う関数を作ってみよう.

    In[1]:= Himando[x_] := Which[
              x < 18.5, Print["yase"],
              x < 25, Print["futsuu"],
              x < 30, Print["himan 1"],
              x < 35, Print["himan 2"],
              x < 40, Print["himan 3"],
              True, Print["himan 4"]    ← ここまで来たら必ず成り立つよう書かれている.True というのは常に「真」となる特別な式であると思えば良い.
              ]

    In[2]:= Himando[BMI[65,1.68]]
    Out[2]= futsuu

    In[3]:= Himando[BMI[120,1.68]]
    Out[3]= himan 4
    

注意 さて,上の Himando 関数を他の関数から再利用したい,ということを考えると, 先にも述べた理由によって Print コマンドは使わない方がよい. つまり,Print 命令は「画面に出力はする」が,「関数 Himando の値にはならない」からである.
そう考えて改良すると,例えば次のように「関数 Himando の値として数字や文字列を返す」などとするのがよい.

    In[4]:= Himando[x_] := Which[
              x < 18.5, -1,
              x < 25,   0,
              x < 30,   1,
              x < 35,   2,
              x < 40,   3,
              True,     4
              ]

    In[5]:= Himando[BMI[65,1.68]]
    Out[5]= 0                        ← かえって分かりにくくなったような気が最初はするが…
    

こうすると,Himando を他の関数から簡単に利用できる. Module の説明の部分で Advice という関数を作ってその例を示すのでみておくと良い.

条件判断(多値) … Switch [ 評価式, 評価1, 式1, 評価2, 式2, 評価3, 式3,... ]

評価式と評価n を順に比較して,評価式が評価n にマッチする最初の n に対して 式n を実行.
先の Which 命令と似ているが,より問題がシンプルなケースに使える.

複雑な関数の定義 … Module

さて,関数の内部でやや複雑な計算が必要となるような場合にどうやって関数を定義して用いるのか, について学習しよう.

なお,複雑な関数をどうやって定義するかについては前回の授業の 複雑な関数の定義 … Module を参照すること.

さてここで,先の修正した Himando 関数がどう「再利用されるか」の例を作成してみよう. これをみれば「関数の再利用」がどういうのものか,そして,その便利性がよく分かるはずである.

    In[6]:= Advice[x_] := Module[{h, alist, r},
              alist = {
                "yasesugi desu. motto tabemasyo",
                "futsuu desu. konomama ikimasyou",
                "himando 1 desune. sukosi yasemasyou",
                "himando 2 desu. kanari yabaidesu. yasero.",
                "himando 3 da. yabai. yaseroyo.",
                "himando 4. futorisugi. iikagen yasero."
                };
              h = Himando[x];
              r = alist[[h+2]];  ← リストは常に「1番」から始まるので +2 として合わせておく.
              Return[r]
              ]    

    In[7]:= Advice[BMI[75, 1.68]]
    Out[7]= himando 1 desune. sukosi yasemasyou
    

注意 こうすることで,「肥満度の判定」と「アドバイス」を分離できるということが重要なのである. つまり,肥満度の判定の基準が多少変わっても Himando 関数を変更するだけでよく, Advice 関数は変更しないで済む.

反復 I … 基本

コンピュータに計算をさせる基本的なテクニックに「反復」がある.
この反復が理解できないと,手で計算するのと機械に計算させるのとで手間 が変わらないことになりかねないので,まずは理解できるように努力しよう.

まずは次に示す基本的な,通常のコンピュータ言語にも必ず実装されている と言って良いコマンドから習熟すると良い.
ただし,後でも示すように,Mathematica にはこれら基本コマンドよりも便利なコマンドが多くあり, かつ,それらを使った方が計算が早かったりするので, 基本コマンドだけで済ませようと思わない方が良い.

反復 … Do[式, {i, min_i, max_i, step_i}]

step_i おきに i を min_i から max_i まで変えながら式を繰り返す.

まずはこの簡単な Do や次の For, While などを理解してから次へ進むのが良いだろう.

例えば,月の利息が 1% の借金をして,12ヶ月たつといくらになるか? を計算する関数をつくると, 次のようになる.

    In[1]:= Debt1[x_] := Module[{y},
              y = x;
              Do[y = y*1.01, {i, 1, 12}];  ← y に新しい y を代入して更新している.
              Return[y]
              ]
  

注意 上の例にもある「代入による変数の値の更新」は, コンピュータを用いた計算の基本テクニックなのでよく慣れておくこと.

反復 … For[初期条件, ループ条件, 式2, 式1]

初期条件を評価した後に,ループに入る. 各ループでは,ループ条件が満たされれば式1,式2 の順序で評価され,もう一度ループに入る. ループ条件が満たされない場合,ループ終り.

Continue コマンドとの整合性のためにこのような仕様になってい るが分かりにくい(^-^).
ちなみに,Continue コマンドは使わない方が良い代表的なコマンドである. よって紹介しない.

例えば,月の利息が 1% の借金をして,12ヶ月たつといくらになるか? を計算する関数をつくると, 次のようになる.

    In[2]:= Debt2[x_]  := Module[{y},
              y = x;
              For[i = 1, i <= 12, i++, y = y*1.01];
              Return[y]
              ]
  

注意 i++ は(値を 1増やして元の値を返す)の省略形である. 他にも++i, i--, --i, i+=a, i-=a, i*=a, i/=a などの省略した形式がある. これらの意味を Help で調べてみよ.

反復 … While[条件, 式]

条件が真ならば式を評価し,また条件をチェック,を繰り返す. 条件が偽になったら終り.

例えば,月の利息が 1% の借金をして,12ヶ月たつといくらになるか? を計算する関数をつくると, 次のようになる.

    In[3]:= Debt3[x_] := Module[{y, i},
            y = x;
            i = 1;
            While[i <= 12, y = y*1.01; i++];   ← y=..;i++ は大きな一つの式と見なされている
            Return[y]
            ]
  

反復 II … 高階関数

Mathematica には上に示した基本的な反復コマンドよりも便利なコマンドがある.
それが
高階関数 … (Mathematica では)関数を引数にとる関数のこと
である. これは数式処理ソフト特有のコマンド群であり, 使いこなせると上の反復の基本コマンドを使うよりも便利なことが多い. Mathematica を使いこなそうと思うレベルになったらまずはこれらのコマンドに習熟すべきである.

リストの要素全てに関数を適用する … Map[関数,リスト]

Map で明示しなくても良いケースも多い.

    In[1]:= a = {1,2,3,4}

    In[2]:= Map[Sin, a]
    Out[2]= {Sin[1],Sin[2],Sin[3],Sin[4]}

    In[3]:= a // Sin
    Out[3]= {Sin[1],Sin[2],Sin[3],Sin[4]}

    In[4]:= Sin[a]
    Out[4]= {Sin[1],Sin[2],Sin[3],Sin[4]}

    In[5]:= % // N
    Out[5]= {0.841471, 0.909297, 0.14112, -0.756802}
  

式の(レベル指定された)「頭部」を関数に置き換える … Apply[関数,式, レベル]

レベル指定を {1} として,リストの中身を書き換えるのに便利.

    In[1]:= a = Table[Sin[n], {n,4}]
    Out[1]:= {Sin[1], Sin[2], Sin[3] , Sin[4]}

    In[2]:= Apply[Cos, a, {1}]
    Out[2]= {Cos[1], Cos[2], Cos[3], Cos[4]}

    In[3]:= Apply[Plus, a]               
    Out[3]= Sin[1]+Sin[2]+Sin[3]+Sin[4]  
         ← レベル指定しないと,{0},すなわち最初の頭部と見なされる.
            a は正確には List[Sin[1],Sin[2],Sin[3],Sin[4]] のことであるので,
            これは Plus[Sin[1],Sin[2],Sin[3],Sin[4]] = ... と置き換えられる.
  

条件にあう要素をリストから取り出す … Select[リスト,条件関数]

条件関数をリストの要素に適用して,True となるものだけを取り出す.

    In[1]:= a = {1.0, -3, 5, 4, 3.5, -1.2}

    In[2]:= Select[a, Positive]    
    Out[2]= {1.,5,4,3.5}
  

関数を引数に繰り返して適用 … Nest[関数,引数, 適用回数]

例えば,Nest[f, x, 5] は f[f[f[f[f[x]]]]] を意味する.
とても便利なコマンドなので是非とも使えるようになっておきたい.

月の利息が 3% の預金をするとして,100万円を預けると 2年後には複利で幾らになるか…

    In[1]:= Debt[x_] := x * 1.03

    In[2]:= Nest[Debt, 100, 24]
    Out[2]= 203.279
    

関数を引数に繰り返して適用 … Fold[関数,引数, リスト]

Nest に似ているが,関数の第二引数がリストの要素になる.
例えば,Fold[f, x, {2,1,5}] は f[f[f[x, 2], 1], 5] を意味する.

例えば,((1.32)0.7)1.2 を計算するには…

    In[1]:= Fold[Power, 1.3, {2, 0.7, 1.2}]
    Out[1]= 1.55391
    

グラフ描画

こうしたツールでグラフを描くというと, 「データをプロットする」という場合と, 「数式を描く」という場合がある.
データをプロットする方法は前回 リストをグラフで見る で学習したので,そちらを参考にすると良い. より詳しくは Help をみるのが良いだろう.

関数の 2次元グラフを描く … Plot[関数, 変数の範囲]

複数のグラフを同時にに描くには, Plot[{関数1, 関数2, ... }, 変数の範囲] とする.

10日借りると 1割の利息が複利でつく借金をしたとして, x 日後に借金はいくらになっているか,関数を作ってグラフ化してみよう.

    In[1]:= Pay[x_, gankin_] := gankin * (1.1^(x/10))   ← 元金を gankin としての x 日後の借金総額

    In[2]:= Table[{n, Pay[n, 100]}, {n, 0, 100, 10}] // TableForm
    Out[2]//TableForm=            ← 100(万円)借りたとして借金がどうなるか10日毎に100日までの様子の確認. 
        0    100
        10   110.
        20   121.
        30   133.1
        40   146.41
        50   161.051
        60   177.156
        70   194.872
        80   214.359
        90   235.795
        100  259.374

    In[3]:= Plot[Pay[x,100], {x,0,400}]  ← 100(万円)借りたとして400日までの様子をグラフ化. 
  

Graph of Debt
↑ かなりすごいことになっている

注意 これはグラフの両軸が「そのまま」であるが,片方,もしくは両方の軸を対数表示する 「対数グラフ」も Mathematica では簡単に書ける. それには, LogLinearPlotLogLogPlot などのコマンドを使えば良い. 詳しくは Help で調べるか,先週示した参照文献を調べると良いだろう.

グラフをファイルとして出力する

Mathematica で得られたグラフをファイルに出力したい,というときは,

そのグラフを右クリック → 編集 → 選択範囲の形式保存 → EPS

としていったん eps ファイルにグラフを保存する. その後, Postscript ファイルから他の画像ファイルへの変換 を参照しながら png ファイルなどに変換してから web に掲載すればよい.

レポート課題

以下に示す課題について 能う限り賢明な調査と考察を行い,

InformationLiteracy-04

という題名をつけて e-mail にて教官宛にレポートとして提出せよ(教官のメールアドレスは授業中に口頭で伝える).
なお,レポートを e-mail の代わりに TeX で作成した書面にて提出してもよい.

注意 なお,本レポートは単位認定のための必須要件「ではない」.
ただし,提出することにより成績に加点されるので,単位があやしい, 成績が気になるという者は積極的に提出する方がよいだろう.

  1. ■ 3x+1 問題 ■ (別名: Collatz 問題, Syracuse 問題,角谷問題)
    まず,(1以上の)自然数 n に対して次のような関数 f(n)
      f(n) = n/2   : n が偶数の場合
      f(n) = 3n+1  : n が奇数の場合  
    を考える.
    これに対し,任意の自然数 n に対して,関数 f を繰り返し適用するといつか必ず 1 になるという「予想」がある. この予想は真か偽か… というのがこの問題である.

    さて,この予想の真偽を体感すべく,上の f(n) を Mathematica で 定義し,n=1 から n=100 程度まで予想が成り立つかどうか試してみよ.

    注意 こうした課題に対して「試してみました」と一行だけ書いて提出しても無意味だ… わかるな?
    また,こうした間抜けなレポートは晒す(公開する)こともあるので,覚悟しておくこと(^-^).

  2. 上の 3x+1 問題の予想では自然数 n に対して関数 f を何回か適用すると 1 になる. そこで,fm(n) = 1 なる最小の m を K(n) と書いたとき, K(n) がいくつになるか,n=1 から n=100 程度まで調べ,グラフにせよ.

  3. sin(x)+sin(1.1 x) を x ∈ [0, 200] の範囲でグラフ化し, そのグラフを自分の web に掲載せよ. その具体的な方法については,すぐ上の グラフをファイルとして出力する を参照せよ.
    さらに,このグラフがなぜそんな形になっているか解説せよ.

  4. Plot の例で示した Pay という関数を片対数グラフ両対数グラフで示せ. また,どちらのグラフが直線状になっているかを答え,なぜ直線になるのか, その理由を数式で解説せよ.




最終更新日 … $Date: 2005-05-04 06:25:03+09 $
Valid CSS! Valid XHTML 1.1!