python入門: 16. 数値微分
微積分を計算機で行ってみよう。
数値微分
関数f(x)の微分は以下の(1)式または(2)式で定義される(両式は等価):
いずれの式も, 右辺のlimの内側は, 分子がf の変化したぶん(左辺のdfに相当)であり, 分母がxの変化したぶん(左辺のdxに相当)である。
計算機では, 微分を「ほぼ」この定義通りに実行する。「ほぼ」と言ったのはこういうことである: 数学的にはhやx1 - xを限りなくゼロに近づけるという概念操作を行うが, 計算機ではそれを行わないのだ。すなわち, 「限りなくゼロに近づける」かわりに, hやx1 - xをなるべく狭くとるのである。それによって誤差は生じるのだが, とりあえずそのことは気にしない(笑)。
例として, f(x)=sin xを数値微分してみよう(もちろん答はcos xになるはずだが, それを知らぬ体でやるのである)。
まず, xをnumpyの配列(numpy.ndarrayのインスタンス)で作ろう。ここでは0以上7未満とし, 間隔を0.1とする。
import numpy as np x=np.arange(0, 7, 0.1)
そして, y=f(x)=sin xもnumpyの配列(numpy.ndarrayのインスタンス)で作ろう:
y=np.sin(x)
すると, yには, sin(0), sin(0.1), ..., sin(6.8), sin(6.9)という値が並んでいるはずだ。
ここで, 微分係数はyの隣同士の差(右ひく左)をxの隣同士の差(右ひく左)で割ることで求まる。たとえばx=0での微分係数は,
{sin(0.1)-sin(0)}/(0.1-0.0)
である。これは,
(y[1]-y[0])/(x[1]-x[0])
である。実際にPythonで↑これを打ってみると,
0.99833...
となる。厳密には1(=cos 0)になるはずだが, ちょっとズレているのは誤差であり仕方ない。
このような計算を全てのxについて行うには, こうすればよい:
for k in range(x.size-1): print((y[k+1]-y[k])/(x[k+1]-x[k]))
ここでずらっと出てきたのが, xのそれぞれの値(0, 0.1, 0.2, ..., 6.9)における微分係数の値である。
課題16-1: 上の命令で, x.size-1の「-1」を取り去って実行してみよ。エラーが出るだろう。なぜだろう?
上の命令では, 計算結果が逐次表示されるが, それをグラフに描いたりするには, 配列にまとめたい。そのときに便利なのが「リスト内包表記」というテクニックである。こうするのだ:
bibun=[(y[k+1]-y[k])/(x[k+1]-x[k]) for k in range(x.size-1)]
すると, bibunというリストに, f'(x)がまとめて格納される。確認してみよう:
print(bibun)
課題16-2: リスト内包表記とはどのようなものか, 調べよ。
では, これをグラフに描いてみよう:
import matplotlib.pyplot as plt plt.plot(x, y) plt.plot(x[:-1], bibun) plt.plot(x, np.cos(x)) plt.show()
もとの関数(sin x)と数値微分の結果(bibun), そして解析的な結果(cos x)が重ねて表示されただろう。ここで, x[:-1]というのは, xの配列の最初から最後のひとつ前まで(xの最後の項を無視する)という意味である。
課題16-3: 以上の説明を再現せよ。すなわち, sin xの数値微分を行って, もとの関数(sin x)と数値微分の結果(bibun), そして解析的な結果(cos x)を重ねてグラフに描け。
ところで, 上のbibun=で「リスト内包表記」を使ったが, 実はそれを使わなくてももっと簡単にできる。こうやるのだ:
bibun=(y[1:]-y[:-1])/(x[1:]-x[:-1])
このほうがPythonらしくてスマートである。ここで, x[1:]というのは, xの配列の初項の次から最後まで(xの初項を無視する)という意味である。初項の次ならからならx[2:]ではないか?と思うかもしれないが, Pythonも含めて, 計算機では順番は1からではなく0から数え始める。xという配列の初項はx[0]であり, その次はx[1]なのである。
今度は, f(x)=2xを数値微分してみよう。numpyを使うと, 累乗はnp.powerという関数で計算できる。たとえば10の0.3乗は,
np.power(10, 0.3)
である。これを使うと, xが-3以上3未満で, 0.1刻みとするとき, 2xの微分はこうやれば求まって, グラフまで描ける:
x=np.arange(-3, 3, 0.1) y=np.power(2, x) bibun=(y[1:]-y[:-1])/(x[1:]-x[:-1]) plt.plot(x, y) plt.plot(x[:-1], bibun) plt.show()
課題16-4: これを参考にして, f(x)=3xを数値微分してグラフに描け。
Keyword(s):
References:[python入門]