sympyで定義した2変数関数をマクローリン展開してプロットする

テイラー展開やマクローリン展開は1変数関数では有名ですが、2変数関数の場合でもよく使います。ここではsympyを使って記号ベースで2変数関数をマクローリン展開し、matplotlibによるプロットをする所までをやってみます。

こんにちは。wat(@watlablog)です。ここでは2変数関数のテイラー展開/マクローリン展開を簡単に説明し、sympy/matplotlibで展開とプロットを行います！

目次（項目クリックでジャンプできます）

2変数関数のテイラー展開について

1変数関数の場合のおさらい

1変数の場合のテイラー変換は「Python/sympyでテイラー展開した結果をグラフ化する方法」に記載していますが、以下に式(1)と、式(1)をまとめた式(2)をn次までのテイラー展開として再掲します。

\[ f(x)=f(a)+\frac{1}{1!}{f}'(a)(x-a)+\frac{1}{2!}{f}''(a)(x-a)^{2}+\frac{1}{3!}{f}'''(a)(x-a)^{3}+\cdots+\frac{1}{n!}{f}^{(n)}(a)(x-a)^{n}+\cdots\tag{1} \]

\[ f(x)=\sum_{n=0}^{\infty }\frac{f^{(n)}(a)}{n!}(x-a)^{n}\tag{2} \]

2変数関数の場合

2変数関数の場合のn次までのテイラー変換は少し微分が複雑になり、式(3)となります。

\[ \begin{align} f(a+h, b+k)=f(a, b)&+\frac{1}{1!}(h\frac{\partial }{\partial x}+k\frac{\partial }{\partial x})f(a, b)\\&+\frac{1}{2!}(h\frac{\partial }{\partial x}+k\frac{\partial }{\partial x})^{2}f(a, b)\\\cdots&+\frac{1}{n!}(h\frac{\partial }{\partial x}+k\frac{\partial }{\partial x})^{n}f(a, b)+\cdots\tag{3} \end{align} \]

僕は数学があまり得意ではなく、この辺になると下手な説明をして誤解を生む可能性があるので、以下の参考書P242にお任せします。

微分積分(改訂版)

微分積分学の教科書にも使われる。矢野健太郎、石原繁著

Amazon

楽天

ただ、本ページでは以下の微分の関係式を使ってプログラムを書いていきますので、備忘録程度にメモしておきます（上記参考書のP234あたりの話です。どの微分の参考書でも載っていると思いますので、詳しくは書籍をご覧下さい）。

その内容ですが、\(z=f(x, y), x=a+ht, y=b+kt\)として、合成関数\(z=f(a+ht, b+kt)\)は、

\[ \frac{\mathrm{d} z}{\mathrm{d} t}=\frac{\partial z}{\partial x}\frac{\partial x}{\partial t}+\frac{\partial z}{\partial y}\frac{\partial y}{\partial t}\tag{4} \]

という式(4)の合成関数の微分公式を利用すれば、

\[ \frac{\mathrm{d} z}{\mathrm{d} t}=h\frac{\partial z}{\partial x}+k\frac{\partial z}{\partial y}\tag{5} \]

と、式(5)となるのは前提から自明なので、2階微分についても同様に、

\[ \begin{align} \frac{\mathrm{d} ^2z}{\mathrm{d} t^2}=\frac{\mathrm{d} }{\mathrm{d} t}(\frac{\mathrm{d} z}{\mathrm{d} t})&=\frac{\mathrm{d} }{\mathrm{d} t}(h\frac{\partial z}{\partial x}+k\frac{\partial z}{\partial y})\\ &=h\frac{\partial }{\partial x}(h\frac{\partial z}{\partial x}+k\frac{\partial z}{\partial y})+k\frac{\partial }{\partial y}(h\frac{\partial z}{\partial x}+k\frac{\partial z}{\partial y})\\ &=h^{2}\frac{\partial^2 z}{\partial x^2}+2hk\frac{\partial^2 z}{\partial x\partial y}+k^{2}\frac{\partial^2 z}{\partial y^2}\tag{6} \end{align} \]

と、式(6)を導くことが出来ます。これは先ほどのn次まで2変数関数のテイラー展開結果、2次の項で使用します。

2変数関数の任意点\((a, b)\)周りのテイラー展開は式(3)に示しましたが、任意点を\((0, 0)\)，\(h, k\)を\(x, y\)とすることで2変数関数のマクローリン展開（式(7)）を得ることが出来ます。

\[ \begin{align} f(x, y)=f(0,0)&+\frac{1}{1!}(x\frac{\partial }{\partial x}+y\frac{\partial }{\partial y})f(0, 0)\\&+\frac{1}{2!}(x\frac{\partial }{\partial x}+y\frac{\partial }{\partial y})^{2}f(0, 0)\\\cdots&+\frac{1}{n!}(x\frac{\partial }{\partial x}+y\frac{\partial }{\partial y})^{n}f(0, 0)+\cdots\tag{7} \end{align} \]

今回はこの2変数関数のマクローリン展開をPythonでプログラミングし、グラフに近似曲面をプロットしてみます。

sympyを使った2変数関数のマクローリン展開コード

\(f(x, y)=e^{x+2y}\)

以下にsympyを使った2変数関数のマクローリン展開のコードを示します。
マクローリン展開部分はmclaughlin_2varというdef関数にまとめており、この関数に展開する元の式と、グラフプロットするためのgridmesh形式のXとYを引数として渡しています。

def関数の中身としては、
①f_taylor：sympyの微分(diff)と階乗(factorial)を使ってテイラー展開を2次精度までする。
②f_taylorに(0, 0)を代入(subsを2回使って1つずつ代入)する。
③hとkにgridmesh形式のXとYを渡してプロット用の曲面Zを得る。
という流れで計算をしています。

from sympy import*
import numpy as np
from matplotlib import pyplot as plt


def mclaughlin_2var(f, X, Y):
    ''' 2変数/2次精度のマクローリン展開をする関数:(f:sympyの式, p:任意点(a,b), XY:np.gridmeshの軸) '''

    # 2変数の場合のテイラー展開(2次精度まで)
    f_taylor = f + (1/factorial(1)) * (h * diff(f, x) + k * diff(f, y))\
                 + (1/factorial(2)) * (h ** 2 * diff(f, x, 2) + 2 * h * k * diff(diff(f, x), y) + k ** 2 * diff(f, y, 2))

    # 点(0,0)を代入(原点周りの近似)しマクローリン展開にする
    subs_zero = f_taylor.subs(x, 0).subs(y, 0)
    print('Original equation=', f)
    print('Series equation=', subs_zero)

    # 1次元に変換
    X_1d = X.ravel()
    Y_1d = Y.ravel()
    Z = np.zeros(len(X_1d))

    # hをx, kをyとしマクローリン展開から曲面を得る
    for i in range(len(X_1d)):
        Z[i] = subs_zero.subs(h, X_1d[i]).subs(k, Y_1d[i])

    return Z

def plot(X, Y, Z):
    ''' グラフプロット '''

    # フォントの種類とサイズを設定する。
    plt.rcParams['font.size'] = 14
    plt.rcParams['font.family'] = 'Times New Roman'

    # グラフの入れ物を用意する。
    fig = plt.figure()
    ax1 = fig.add_subplot(111, projection='3d')

    # 軸のラベルを設定する。
    ax1.set_xlabel('x')
    ax1.set_ylabel('y')
    ax1.set_zlabel('z')

    # データプロットする。
    ax1.plot_wireframe(X, Y, Z, label='Theory', color='black')
    ax1.plot_wireframe(X, Y, np.reshape(result, X.shape), label='Series:O=2')
    plt.legend()

    # グラフを表示する。
    plt.show()
    plt.close()
    # -------------------------------------------------------


if __name__ == '__main__':
    ''' Main '''

    # 理想波形を生成
    np_x = np.arange(-1, 1.2, 0.2)
    np_y = np.arange(-1, 1.2, 0.2)
    X, Y = np.meshgrid(np_x, np_y)
    Z = np.exp(X + 2 * Y)
    
    # sympy変数、展開する式の定義と関数実行
    x, y, h, k = symbols('x y h k')
    f = exp(x + 2 * y)
    result = mclaughlin_2var(f, X, Y)

    # プロット
    plot(X, Y, Z)

from sympy import*

import numpy as np

from matplotlib import pyplot as plt

def mclaughlin_2var(f, X, Y):

''' 2変数/2次精度のマクローリン展開をする関数:(f:sympyの式, p:任意点(a,b), XY:np.gridmeshの軸) '''

# 2変数の場合のテイラー展開(2次精度まで)

f_taylor = f + (1/factorial(1)) * (h * diff(f, x) + k * diff(f, y))\

+ (1/factorial(2)) * (h ** 2 * diff(f, x, 2) + 2 * h * k * diff(diff(f, x), y) + k ** 2 * diff(f, y, 2))

# 点(0,0)を代入(原点周りの近似)しマクローリン展開にする

subs_zero = f_taylor.subs(x, 0).subs(y, 0)

print('Original equation=', f)

print('Series equation=', subs_zero)

# 1次元に変換

X_1d = X.ravel()

Y_1d = Y.ravel()

Z = np.zeros(len(X_1d))

# hをx, kをyとしマクローリン展開から曲面を得る

for i in range(len(X_1d)):

Z[i] = subs_zero.subs(h, X_1d[i]).subs(k, Y_1d[i])

return Z

def plot(X, Y, Z):

''' グラフプロット '''

# フォントの種類とサイズを設定する。

plt.rcParams['font.size'] = 14

plt.rcParams['font.family'] = 'Times New Roman'

# グラフの入れ物を用意する。

fig = plt.figure()

ax1 = fig.add_subplot(111, projection='3d')

# 軸のラベルを設定する。

ax1.set_xlabel('x')

ax1.set_ylabel('y')

ax1.set_zlabel('z')

# データプロットする。

ax1.plot_wireframe(X, Y, Z, label='Theory', color='black')

ax1.plot_wireframe(X, Y, np.reshape(result, X.shape), label='Series:O=2')

plt.legend()

# グラフを表示する。

plt.show()

plt.close()

# -------------------------------------------------------

if __name__ == '__main__':

''' Main '''

# 理想波形を生成

np_x = np.arange(-1, 1.2, 0.2)

np_y = np.arange(-1, 1.2, 0.2)

X, Y = np.meshgrid(np_x, np_y)

Z = np.exp(X + 2 * Y)

# sympy変数、展開する式の定義と関数実行

x, y, h, k = symbols('x y h k')

f = exp(x + 2 * y)

result = mclaughlin_2var(f, X, Y)

# プロット

plot(X, Y, Z)

以下に実行結果（printの結果)を示します。Original equationが元の式、Series equationがsympyで計算されたマクローリン展開結果です。

Original equation= exp(x + 2*y)
Series equation= h**2/2 + 2*h*k + h + 2*k**2 + 2*k + 1

1 2	Original equation= exp(x + 2y) Series equation= h2/2 + 2hk + h + 2k*2 + 2k + 1

\(f(x, y)=e^{x+2y}\)は各微分形が

\[ f_{x}=f_{xx}=e^{x+2y}, f_{y}=f_{xy}=2e^{x+2y}, f_{yy}=4e^{x+2y} \]

と表せるので、2次までのマクローリン展開は式(8)となります。

\[ 1+x+2y+\frac{x^2}{2}+xy+2y^{2} \tag{8}\]

先ほどのSeries equationと同じ形をしていることが確認できました。

上記コードを実行すると以下の図も出力されます。この図は元の式Theoryと展開した式を図示したSeriesを重ね描きプロットしているもので、原点(0, 0)あたりがよく一致していることが確認できると思います。

\(f(x, y)=\sin(x+y^{2})\)

関数を変えてやってみましょう。以下部分のコードで、式のみを\(\sin(x+y^{2})\)と変更してみます。

# 理想波形を生成
np_x = np.arange(-1, 1.2, 0.2)
np_y = np.arange(-1, 1.2, 0.2)
X, Y = np.meshgrid(np_x, np_y)
Z = np.sin(X + Y ** 2)

# sympy変数、展開する式の定義と関数実行
x, y, h, k = symbols('x y h k')
f = sin(x + y ** 2)
result = mclaughlin_2var(f, X, Y)

# 理想波形を生成

np_x = np.arange(-1, 1.2, 0.2)

np_y = np.arange(-1, 1.2, 0.2)

X, Y = np.meshgrid(np_x, np_y)

Z = np.sin(X + Y ** 2)

# sympy変数、展開する式の定義と関数実行

x, y, h, k = symbols('x y h k')

f = sin(x + y ** 2)

result = mclaughlin_2var(f, X, Y)

以下が実行結果としてコンソールに表示される式です。

Original equation= sin(x + y**2)
Series equation= h + k**2

1 2	Original equation= sin(x + y2) Series equation= h + k2

グラフは下図になります。コード上で元の式を変更しただけですがやはり原点周りはよく一致しているようです。

まとめ

本記事では2変数関数にテイラーの定理を適用させ、原点周りの近似を2次まで行うマクローリン展開の式とPythonコードを紹介しました。

実際に指数関数と三角関数の2種類に適用させ、グラフプロットにて原点周りの近似が出来ていることを確認しました。

（本当は任意個数の変数、任意点周り、n次のテイラー展開までやりたかったのですが、sympyがうまく扱えず、まずはマクローリンで妥協しました…）

前回の1変数関数のテイラー展開ではsympy.seriesでテイラー展開自体を自動で計算していましたが、今回のような多変数関数に対応していないようでした。もっと簡単に多変数関数、任意点周りの近似が得られれば様々な開発案件で役に立つものとなると思いますが、先は少し長そうです。

まずは2変数関数ですが、曲面を原点周りで近似する方法を1つ習得しました！
Twitterでも関連情報をつぶやいているので、wat(@watlablog)のフォローお待ちしています！

2変数関数のテイラー展開について

1変数関数の場合のおさらい

2変数関数の場合

sympyを使った2変数関数のマクローリン展開コード

\(f(x, y)=e^{x+2y}\)

\(f(x, y)=\sin(x+y^{2})\)

まとめ

コメントを残すコメントをキャンセル

最近の投稿

アーカイブ

カテゴリー

2変数関数のテイラー展開について

1変数関数の場合のおさらい

2変数関数の場合

sympyを使った2変数関数のマクローリン展開コード

\(f(x, y)=e^{x+2y}\)

\(f(x, y)=\sin(x+y^{2})\)

まとめ

SNSでもご購読できます。

コメントを残す コメントをキャンセル

最近の投稿

アーカイブ

カテゴリー

コメントを残すコメントをキャンセル