とらりもん - python入門: 16. 数値微分と数値積分 Diff
- Added parts are displayed like this.
- Deleted parts are displayed
like this.
微積分を計算機で行ってみよう。
!数値微分
関数~~f~~(~~x~~)の微分は以下の(1)式または(2)式で定義される(両式は等価)：
http://pen.envr.tsukuba.ac.jp/~torarimon/etc/deriv.png
いずれの式も, 右辺のlimの内側は, 分子が~~f~~ の変化したぶん（左辺の~~df~~に相当）であり, 分母が~~x~~の変化したぶん（左辺の~~dx~~に相当）である。
計算機では, 微分を「ほぼ」この定義通りに実行する。「ほぼ」と言ったのはこういうことである: 数学的には~~h~~や~~x~~__1__ - ~~x~~を限りなくゼロに近づけるという概念操作を行うが, 計算機ではそれを行わないのだ。すなわち, 「限りなくゼロに近づける」かわりに, ~~h~~や~~x~~__1__ - ~~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~~)=2""~~x~~""を数値微分してみよう。numpyを使うと, 累乗はnp.powerという関数で計算できる。たとえば10の0.3乗は,
np.power(10, 0.3)
である。これを使うと, xが-3以上3未満で, 0.1刻みとするとき, 2""~~x~~""の微分はこうやれば求まって, グラフまで描ける:
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~~)=3""~~x~~""を数値微分してグラフに描け。
!数値微分
関数~~f~~(~~x~~)の微分は以下の(1)式または(2)式で定義される(両式は等価)：
http://pen.envr.tsukuba.ac.jp/~torarimon/etc/deriv.png
いずれの式も, 右辺のlimの内側は, 分子が~~f~~ の変化したぶん（左辺の~~df~~に相当）であり, 分母が~~x~~の変化したぶん（左辺の~~dx~~に相当）である。
計算機では, 微分を「ほぼ」この定義通りに実行する。「ほぼ」と言ったのはこういうことである: 数学的には~~h~~や~~x~~__1__ - ~~x~~を限りなくゼロに近づけるという概念操作を行うが, 計算機ではそれを行わないのだ。すなわち, 「限りなくゼロに近づける」かわりに, ~~h~~や~~x~~__1__ - ~~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]なのである。
np.power(10, 0.3)
である。これを使うと, xが-3以上3未満で, 0.1刻みとするとき, 2""~~x~~""の微分はこうやれば求まって, グラフまで描ける:
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~~)=3""~~x~~""を数値微分してグラフに描け。