偏微分を数値計算で求めたい

今回は数値微分の話です。

偏微分を数値で計算することが必要な場面は日常ではあまりないのですが、常微分を数値計算したことがあったので、偏微分はどうなんだろ？

ということで計算してみました。

簡単な関数で考えます。

f(x,y) = x^2+y^2

いわゆる円の方程式です。

√(f)

で半径が出ます。

これをxで偏微分すると

fx(x,y) = 2x

yで偏微分すると

fy(x,y) = 2y

点(2,3)においては

fx(2,3) = 4

fy(2,3) = 6

です。

これを数値計算で求めてみようという話になります。

参考にしたのは以下。

堀之内總一 酒井幸吉　榎園茂（2015）『Cによる　数値計算法入門（第2版）新装版』森北出版株式会社.

偏微分の計算方法の数式は書くと面倒なので書籍などを参考にしてください。

作成したプログラムは以下。

def f1(x):
    return x[0]**2 + x[1]**2
def bibun1(step,x,k):
    high = []
    low = []
    for i in range(len(x)):
        if(i != k):
            high.append(x[i])
            low.append(x[i])
        else:
            high.append(x[i] + step)
            low.append(x[i] - step)
    return (f1(high) - f1(low))/(2*step)

step = 0.01
x = [2,3]
print(f1(x))
print("1:",bibun1(step,x,0))
print("2:",bibun1(step,x,1))

そして出力は以下。

13

1: 3.9999999999999147

2: 5.999999999999872

ちゃんと4と6が出ました。

数値微分ってコンピュータでどうやるんだろうと疑問だったのですがスッキリしました。

なんで数値計算する必要があるのかっていうと、元の関数が数式で微分できるならいいんですけど、そうじゃないシチュエーションがあるんです。

そういうときに必要なんでしょうね。

これで偏微分方程式も解けるようになるんですけど、波動方程式とか熱拡散方程式とか日常で使わないのでやりません。

まあ使う人には入り口になるかなという感覚で書いてみました。