とらりもんHOME  Index  Search  Changes  Login

python入門: 16. 数値微分

微積分を計算機で行ってみよう。

数値微分

関数f(x)の微分は以下の(1)式または(2)式で定義される(両式は等価):

deriv.png

いずれの式も, 右辺のlimの内側は, 分子がf の変化したぶん(左辺のdfに相当)であり, 分母がxの変化したぶん(左辺のdxに相当)である。

計算機では, 微分を「ほぼ」この定義通りに実行する。「ほぼ」と言ったのはこういうことである: 数学的にはhx1 - xを限りなくゼロに近づけるという概念操作を行うが, 計算機ではそれを行わないのだ。すなわち, 「限りなくゼロに近づける」かわりに, hx1 - 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を数値微分してグラフに描け。

Last modified:2023/01/03 10:05:02
Keyword(s):
References:[python入門]