feature image

2017年4月9日 | ブログ記事

Pythonで数学の課題を解こう【新歓ブログリレー2017 6日目】

こんにちは。hukuda222です。情報通信系の学部2年生です。

バスゼミに行った2,5類の方は初めての休日、他の類の方は連休2日目ですね。

一年前の僕は2類のバスゼミで徹夜で遊んで疲れたので昼まで寝ていました。

・前置き

近年、ハードウェアの進歩によって機械学習の分野に注目が集まっており、その影響でその辺に強いプログラミング言語であるところのPython、Rあたりの需要が高まっています。Pythonは言語別年収で1位らしいですし、これはマスターしておけばきっと将来安泰でしょう[要出典]

Pythonにはバージョン2系と3系があり、2系にしか対応してないライブラリも多いので併用して使われています。いずれは、全てのライブラリが最新バージョンに対応するようなのでここでは最新の安定バージョンである3.6.0(4/9現在)を使います。

それから、未習範囲の数学の話も出てきますが、予習がてら目を通すなり、授業で課題が出て「この回答であってるかな?」となった時にでも見てください。

・準備

Python最新版をインストール

……するのは面倒かもしれないので、数学の検算だけできればいいやって場合はこのサイトを使ってください。
http://live.sympy.org/

Macに入れたい場合はこのサイトを参考にすると良いと思います。
http://www.task-notes.com/entry/20141215/1418612400

windowsに入れたい場合はこのサイトを参考にすると良いと思います。
windowsでのPythonのインストール方法

Pythonのインストール以降は、基本的にコマンドプロントでの操作になります。

SymPyのインストール

SymPyは、PythonでSymbol、つまり変数を扱うことのできるライブラリです。
これを用いることで、大学数学の計算機を用いることを想定していない計算問題を一瞬で解くことができます。
ブラウザで動くsympyを使う場合は、以下の操作は必要ありません。

pipが使える場合は、

pip install sympy

でインストールできます。

SymPyをダウンロード(GitHubに飛びます)

Macの方なら、sympy-1.0.tar.gz、Windowsの方ならsympy-1.0-win32.exeを上のリンクからダウンロードしてください。

ダウンロードしたものを解凍したsympy-1.0フォルダ内のsetup.pyファイルがある場所で

python setup.py install

を実行することで、インストールできます。


Sympyで遊ぼう

今回は対話モードで実験したいと思います。

pythonとコンソールにタイプすると、入力待ち画面になるはずなのでそこで、

>>import sympy
>>from sympy import *

とタイプします。(ブラウザのSympyを使う場合はここまでのコマンドは必要ありません)

最初のコマンドは"Sympy"っていうのを使うよ、2つ目はSympyの全ての機能を使うよという意味です。

※以下ではpprintという関数を使っていますが、ブラウザ版のSympyだとpprintを使わなければ、latexで綺麗に表示してくれるので使わない方がいいです。

・変数の宣言

>>x = Symbol('x')
>>y = Symbol('y')

これで、xとyという変数を使うことができます。

・四則演算とべき乗

 >>pprint(x+y)
  x+y
 >>pprint(x-y)
  x-y
 >>pprint(x*y)
  xy
 >>pprint(x/y)
  x
  -
  y

yという感じで表示されます。x/yの結果がちょっとわかりにくいかもしれませんが、ちゃんとyが分母でxが分子の分数として表示されています。

 >>pprint(x**2)
  ㅤ 2
  x

見にくいですが、x**nでxのn乗を表します。
縦長ですが、普段使う表記で表示してくれます。

・有名な定数と関数

>>init_session()

この関数を実行することでよく使う定数や関数を使えるようになります。

 >>pprint(exp(cos(x**I)+sin(E*pi)))

       ⎛   i⎞
   cos⎝x ⎠ + sin(e⋅π)
e

かなり見辛いですが、"E"は自然対数e、"I"は虚数単位i、"pi"で円周率πを示せる他、expやlog、三角関数はそのまま使えます。(exp(n)で自然対数eのn乗を表します)

・方程式を解く

 >>solve(x**2-2*x+1,x)
  [1]
 >>solve([x+y-4,x-y-2],[x,y])
 {y:1,x:3}

こんな感じで、=0となる式と求めたい変数を入れると計算してくれます。
あとの具体的に使いそうな機能は、手元にある実際に授業で使った教科書の計算問題で紹介します。


1.の導関数を求めよ

>>diff((x**2+1/((x-1)**2))**(1/3),x,1)

 
2/3を丸めて計算しているので0.6666666667という表記になってしまっていますが、ちゃんと計算できています。
ただ、問題によっては計算ができているものの約分ができてなかったりすることがあります。

2.を漸化展開を用いて書き表せ

数学的な説明は長くなるので省きますが、式が極限まで0に近づくとどうなるのかを考えると、分子分母がそれぞれ0になってしまいます。
そこで、漸化展開をすることによって分子と分母を約分できるような形に近似させます。
こうすることで、xが0に限りなく近い場合にとる値も求めることができます。

 >>f1 = sin(x)-x*E**x+x**2
 >>f2 = x*(cos(x)-1)
 >>limit(f1/f2,x,0)
  3
  -
  4

心配しなくてもSympyが勝手に漸化展開して計算してくれます。漸化展開した結果は、series(f1,x)で見ることができます。

3.をガンマ関数を用いて求めよ

 >>integrate(x**3*(E**(-x**(0.5))),(x,0,oo))
 >>10080.0000000000

無限は、ooで表記します。ガンマ関数も必要に応じて勝手に使ってくれます。

4. この行列Aを簡約化せよ

 >>A = Symbol('A')
 >>A = Matrix([[1,2,3,2],[1,2,1,1],[1,2,-1,0]])
 >>A.rref()
 (Matrix([
 [1, 2, 0, 1/2],
 [0, 0, 1, 1/2],
 [0, 0, 0,   0]]), [0, 2])

一つ目の要素のMatrix([[1, 2, 0, 1/2],[0, 0, 1, 1/2],[0, 0, 0,   0]])がAを簡約化した行列です。
二つ目の要素は、Aをどのように簡約化したかを示すピポットというものです。
詳しくはググるか数学の先生に聞いてください。
慣れるとこんな簡単な問題ならすぐにできるのですが、簡約化は計算ミスの温床になりやすいのでチェックは入念にしたいところです。

それから、省きましたが行列同士の演算は商以外は実数と同様に計算できます。

5. この行列Bの行列式を求めよ

 >>var("a:z") 
 (a, b, c, d, e, f, g, h, i, j, k, l, m, n, o, p, q, r, s, t, u, v, w, x, y, z)
 >> B=Matrix([[a,0,0,b],[c,d,0,0],[e,f,g,0],[0,0,h,i]])
 >> B.det()
 a*d*g*i - b*c*f*h + b*d*e*h

var("a:z")で、a~zの変数を一気に宣言することができます。
変数を含む行列の行列式を出すのも余裕です。

6. この行列Cの余因子行列を用いて逆行列を求めよ

 >>C=Symbol("C")
 >> C=Matrix([[x-2,1,1][0,2*x-1,x-1][-2,1,1]])
 >>C.adjugate()/C.det()
 >>Matrix([
 [            1/x,    0,                     -1/x],
 [(-2*x + 2)/x**21/x,  -(-x + 1)*(-x + 2)/x**2],
 [ (4*x - 2)/x**2, -1/x, (-2*x + 1)*(-x + 2)/x**2]])

C.adjugate()で余因子行列が出せます。逆行列というのは、余因子行列を行列式で割ったものなのでそのままかくと答えが出ます。

 >> C.inv()
  Matrix([
  [1/(x - 2) - 2*(-(x - 1)/((x - 2)*(2*x - 1)) + 1/(x - 2))*(2*x**2 - 5*x + 2)/(x**2*(x - 2)), -1/((x - 2)*(2*x - 1)) + (1 + 2/(x - 2))*(-(x - 1)/((x - 2)*(2*x - 1)) + 1/(x - 2))*(2*x**2 - 5*x + 2)/(x**2*(2*x - 1)), -(-(x - 1)/((x - 2)*(2*x - 1)) + 1/(x - 2))*(2*x**2 - 5*x + 2)/x**2],
  [ -2*(x - 1)*(2*x**2 - 5*x + 2)/(x**2*(x - 2)*(2*x - 1)), 1/(2*x - 1) + (1 + 2/(x - 2))*(x - 1)*(2*x**2 - 5*x + 2)/(x**2*(2*x - 1)**2), -(x - 1)*(2*x**2 - 5*x + 2)/(x**2*(2*x - 1))],
  [ 2*(2*x**2 - 5*x + 2)/(x**2*(x - 2)), -(1 + 2/(x - 2))*(2*x**2 - 5*x + 2)/(x**2*(2*x - 1)), (2*x**2 - 5*x + 2)/x**2]])
  >>simplify(C.inv())
 Matrix([
 [             1/x,    0,              -1/x],
 [ (-2*x + 2)/x**21/x, -1 + 3/x - 2/x**2],
 [2*(2*x - 1)/x**2, -1/x,  2 - 5/x + 2/x**2]])

C.inv()でもそのまま逆行列を出すことができます。
この際、そのまま書くとぐちゃぐちゃなことになるのでsimplifyの中に入れることで見やすい形に変形させることができます。simplifyは様々な式に対して使えます。

・終わりに

実際に数学の課題を見るとわかりますが、単純な計算よりも難しい問題はたくさんあります。
証明問題とかε-δ問題とか、どちらかというと、計算するだけの段階にまで落とし込めれば正解したも同然という問題が多くあります。
そこで、立式まではあってるのに最後の最後の計算ミスによって、数学の成績を下げてしまうと大変もったいないのです。

数学の単位は1,2Qの必修で4単位、3,4Qの選択で6単位なのですごーく重要です。

試験の場合は何回もチェックするくらいしかできないですが、課題の場合は何を使ってもいいので、計算機や友人の回答と照らし合わせて是非是非高得点を目指して、希望の系に行けるように頑張ってください。

明日のブログは、「ドット絵体験会案内」とyuuさんの「AndroidアプリのGoogle Drive連携」の二本立てです。乞うご期待。

・参考文献

・公式リファレンス(英語)
http://docs.sympy.org/latest/index.html

・Scipy Lecture Notesのおまけ記事(日本語)
http://www.turbare.net/transl/scipy-lecture-notes/packages/sympy.html

・プログラミング速報(日本語)
http://programming.blogo.jp/python/sympy/introduction

hukuda222 icon
この記事を書いた人
hukuda222

主にプログラムを書いてます。

この記事をシェア

このエントリーをはてなブックマークに追加
共有

関連する記事

2017年11月14日
IBIS2017参加報告
Keijan icon Keijan
2022年10月11日
アルゴリズム班にKaggle部を設立し、初心者向けデータ分析体験会を開催しました!
abap34 icon abap34
2021年4月18日
ベズー係数とN項の拡張ユークリッドの互除法
0214sh7 icon 0214sh7
2017年11月4日
文章をよしなに分散表現しよう
David icon David
2022年11月5日
ピックの定理の多次元拡張 -エルハート多項式-
0214sh7 icon 0214sh7
2022年3月27日
予想外のバグから生まれた、最高に楽しい課題の話
logica icon logica
記事一覧 タグ一覧 Google アナリティクスについて