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にお任せします。

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

ただ、本ページでは以下の微分の関係式を使ってプログラムを書いていきますので、備忘録程度にメモしておきます(上記参考書の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変数関数のマクローリン展開コード

Advertisements

\(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を得る。
という流れで計算をしています。

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

\(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)あたりがよく一致していることが確認できると思います。

exp(x+2y)のマクローリン展開結果

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

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

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

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

sin(x+y^2)のマクローリン展開結果

まとめ

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

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

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

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

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

SNSでもご購読できます。

コメントを残す

*