数学わかんねーからプログラムでズルしてみた 微分編（Python：Sympy）

微分ってなんぞ？

知るわきゃねーだろ！
コンビニの会計で使うくらい一般的になって出なおせや！

Wikipediaを読む

数学、とくに解析学における微分法（びぶんほう、differentiation, derivation）は、
空間やその上に定義される関数・写像を各点の近傍で考え、
その局所的な振舞いを調べることによって、それらの特徴を記述する方法である。

微分法 - Wikipedia

ふむ、局所的な振る舞いが空間で写像するわけだな。
振る舞い関数の近傍を記述する方法であると。

はい、お手上げ！

とりあえず微分する

In [4]: diff(3*x, x)
Out[4]: 3

変数が一つなので、何も指定しなくても良い。

In [3]: diff(x**2)
Out[3]: 2⋅x

定数は消える。

In [3]: diff(x**2 + 10000)
Out[3]: 2⋅x

こっからは偏微分（変数が2個以上あると、一部だけ微分するので偏微分というらしい）

In [7]: a = x**2 * sin(x + y)

In [8]: a
Out[8]: 
 2           
x ⋅sin(x + y)

In [9]: diff(a, x)
Out[9]: 
 2                            
x ⋅cos(x + y) + 2⋅x⋅sin(x + y)

In [10]: diff(a, y)
Out[10]: 
 2           
x ⋅cos(x + y)

変数が2つあるので、明示的に何で微分するか渡す必要がある。
第2引数を渡さないとエラー。

なにやら正解っぽいのが出たが、これっぽっちも微分がなんたるかを理解していない。
ググったところ、ぐにゃぐにゃしてるグラフのある一点での傾きを求めるものっぽい。

さっき定数が消えたのは、傾きと関係ないからか？

newton_methodなるものが、微分使ってがんばるらしいのでやってみた。

newton_methodとは

りんごが落ちそうな雰囲気が漂ってる。
Newton-Raphson method（ニュートン・ラフソン法）とも呼ばれるらしい。
名称はアイザック・ニュートンとジョセフ・ラフソンに由来するらしいが、ラフソンさん省略されることがあるとかかわいそう。

Wikipediaを読む

方程式系を数値計算によって解くための反復法による求根アルゴリズムの1つである。

http://ja.wikipedia.org/wiki/%E3%83%8B%E3%83%A5%E3%83%BC%E3%83%88%E3%83%B3%E6%B3%95

方程式系の系ってなんだよ系って。
細かいことは気にせず先に進もう。

やってみる

Wikipediaにの近似値の求め方が書いてあったので、こいつを題材にしよう。

1.41421356（ヒトヨヒトヨニヒトミゴロ）を求めるのが目標。
9桁までそれっぽくでりゃいいっしょ。

In [1]: sqrt(2)
Out[1]: 
   ___
╲╱ 2 

In [2]: sqrt(2).evalf(9)
Out[2]: 1.41421356

よっしゃぁ！！！
クリアぁああああ！！

じゃなくて、newton_method使わないとね、せっかくだし。

まず、平方根をグラフで表現してみる。
方程式にするとこんな感じ。

In [1]: plot(x**2)
Out[1]: <sympy.plotting.plot.Plot at 0x108989dd0>

x=1の時、y=1
x=1.4142... y=2
x=2の時、y=4
x=3の時、y=9 ...
となる。

で、こんなグラフ書いてどうするわけ？ってとこなんだが、どうにもならない・・・
この方程式をちょっといじれば、newton_methodで解ける形になる。

がその方程式。

Cは切片。C=-2としてプロットしてみる。

In [1]: plot(x**2 - 2, (x, -5, 5))
Out[1]: <sympy.plotting.plot.Plot at 0x106f5af10>

この関数のyに0を代入すると



となり、yが0の時に、求めたい値が出るって仕組み。

こいつをsympyで見てみると

In [1]: C = symbols('C')

In [2]: solve(x**2 - C, x)
Out[2]: 
⎡   ___    ___⎤
⎣-╲╱ C , ╲╱ C ⎦

になる。
ずれてるってレベルじゃねーぞ。

ここからが本題。
Wikipediaの導入の項目を見てみる。

この方法の考え方は以下のようである：まず初めに、予想される真の解に近いと思われる値をひとつとる。
次に、そこでグラフの接線を考え、その x 切片を計算する。

http://ja.wikipedia.org/wiki/%E3%83%8B%E3%83%A5%E3%83%BC%E3%83%88%E3%83%B3%E6%B3%95#.E5.B0.8E.E5.85.A5

予想される真の解・・・1.4142...ってわかってんだけど、今回は2にしてみる。

のグラフと、左記方程式のx座標が2の時の接線を同時にプロットしてみる。

まず、x座標が2の時の方程式を求める。

In [5]: Eq(y, x ** 2 - 2)
Out[5]: 
     2    
y = x  - 2

を微分すると

In [6]: diff(x ** 2 - 2)
Out[6]: 2⋅x

x座標が2の時の傾きは

In [8]: (2 * x).subs(x, 2)
Out[8]: 4

これぞ微分パワー！

点(2, 2)を通る傾き4の1次方程式のグラフは
bを切片とすると

In [10]: b = symbols('b')

In [11]: Eq(2, 4 * 2 + b)
Out[11]: 2 = b + 8

これをbについて解くと

In [12]: solve(Eq(2, 4 * 2 + b))
Out[12]: [-6]

※戻り値がlistなのは、解が複数ある方程式があるからね。とか

これで傾きと切片が出せたので、1次方程式が書ける。

と一緒にプロットすんぜ！

In [13]: plot(x**2-2, 4*x-6, (x, -4, 4))
Out[13]: <sympy.plotting.plot.Plot at 0x1078a2890>

ちなみに赤文字はmac付属のソフトでちまちま加筆してます・・・

√2に近い何かが何かというと・・・
のyが0の時なので

In [31]: Eq(y, 4 * x -6)
Out[31]: y = 4⋅x - 6

In [32]: Eq(y, 4 * x -6).subs(y, 0)
Out[32]: 0 = 4⋅x - 6

In [33]: solve(Eq(y, 4 * x -6).subs(y, 0), x)
Out[33]: [3/2]

In [34]: solve(Eq(y, 4 * x -6).subs(y, 0), x)[0].evalf()
Out[34]: 1.50000000000000

ようやく出たよに近い何か！
1.5だった！！確かに1.41421356...に近い！！

Wikipediaの導入の項目の続きをを読むと

このx切片の値は、予想される真の解により近いものとなるのが一般である。
以後、この値に対してそこでグラフの接線を考え、同じ操作を繰り返していく。

http://ja.wikipedia.org/wiki/%E3%83%8B%E3%83%A5%E3%83%BC%E3%83%88%E3%83%B3%E6%B3%95#.E5.B0.8E.E5.85.A5

って書いてあるけど、最初に使ったx切片が2だね。
同じ操作を繰り返していくってのは、次はx切片を1.5にしてやってみようって事。
これを繰り返すとに近づきまっせってことだ。

定式化しないと使い物にならんので、定式化する。

だが、ここでまさかの気力切れ。

あるじゃない、同じ事解説してる動画。
めっちゃわかりやすい。
しかもイケメンボイス。

さて、手を抜いたところで最後にプログラム貼り付けて終わり。

def sqrt_newton_method(value, cnt):
	 # x切片
	preval = value

	for i in range(cnt):
		preval = 1/2 * (preval + value/preval)

	return preval

if "__main__" == __name__:
	# √2を1回計算
	print(sqrt_newton_method(2, 1)) #=>1.5
	# √2を2回計算
	print(sqrt_newton_method(2, 2)) #=>1.4166666666666665
	# √2を4回計算
	print(sqrt_newton_method(2, 4)) #=>1.4142135623746899

	print()
	# √3を4回計算
	print(sqrt_newton_method(3, 4)) #=>1.7320508100147274