Python/sympyでテイラー展開した結果をグラフ化する方法

  • このエントリーをはてなブックマークに追加

テイラー展開は任意点周りの近似計算をすることが出来るため、様々な工学的・物理的場面で活躍します。ここではsympyを使って面倒なテイラー展開をPythonにやらせ、グラフプロットまで自動で行うことを目標とします。

こんにちは。wat(@watlablog)です。ここではPython/sympyを使ってテイラー展開を自動で行ってプロットする所までをやってみます

テイラー展開の式

Pythonのsympyでテイラー展開を計算する前に、まずは式の確認をしていきましょう。

このページでは\(f(x)\)と関数の変数が\(x\)とただ1つだけの場合を考えます。
任意点\(x=a\)周りのテイラー展開を式(1)で示します。

\[ 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 (1) \]

式(1)は無限項の足し算になっており、式(2)とまとめて書くことができます。

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

この式は一見すると複雑で難しい印象を受けます。

しかし、1次近似までを見ると高校で習う接線の公式になっていることに気付くと思います。

テイラー展開の目的は任意点の情報から近傍の状態を知る(近似する)ことにあります。

テイラー展開は接線の公式を延長し、2次、3次と無限次まで精度を高めると元の関数と一致するということを意味しています。複雑な式の意味も元の式を多項式で近似しているということですね。

sympyの使い方:テイラー展開を体験する

sympyとは?

sympyとは、記号(Symbol)を記号のまま数学的な処理をすることが可能なPythonライブラリです。

例えば、\(y=x^{2}\)を微分して値を求める場合、通常のプログラミングであればプログラマが\(y'=2x\)を記述し、そこに値を入れるか、近似的な計算手法を使って求めるかの2通りが考えられます。

しかしsympyを使えば\(x\)と\(y\)をSymbolとして定義し、プログラム的に\(y'=2x\)を求めることが可能です。近似的手法に頼ることなく解が得られる等の様々なメリットが考えられますね。

インストール

sympyはpipでインストールすることが可能です。Windowsであればコマンドプロンプトに以下のコードを打ち込むことでインストールが完了します。

import文

import文の書き方は人により様々ですが、import*とアスタリスクを付けた形で書く人が多いかも知れません。 sympyはsinやcos, expといった関数はそれぞれimportすることも出来ますが、import*と書いておけば個別に書く必要がありません。

記号(シンボル)を定義する

sympyではこれから扱いたい\(x\)や\(y\)といった変数を以下のコードで示すようにSymbolを使って定義します。

テイラー展開を計算する

以下がテイラー展開を計算するコードです。指数関数\(f(x)=e^{x}\)を\(x=0\)点で第3項\(n=3\)まで展開しています(原点周りなので、正確にはマクローリン展開と呼びます)。

\(e^{x}\)をマクローリン展開すると、以下の式(3)となります(上の式(2)に当てはめて計算してみて下さい)。

\[
f(x)=1+x+\frac{1}{2}x^{2}+\cdots (3)
\]

上記コードを実行すると、以下の結果を得ます。式(3)と同じですね。
O()は剰余項を意味していて、nで指定した項以外がまとめられているものだそうです。

なんとなくsympyの使い方がわかってきました!それにしても、記号を処理してくれるのは地味にすごいですね!

ただ、これだけだと「計算結果を表示させてみました」に過ぎず、特に面白くありません。

そこで、以下にsympyで得られた式をnumpy形式で出力させ、matplotlibでグラフへプロットするということをやってみます。

sympyで自動テイラー展開してプロットをするコード

\(f(x)=e^{x}\)

以下のコードは\(e^{x}\)をseriesでテイラー展開し、計算された多項式をlambdifyを使ってnumpyの関数化させています。solは関数化させたオブジェクトtaylor_yに横軸cal_xを代入してnumpy形式で配列を得ています。

このコードではfor文を使ってn(項数)違い、つまり複数の精度違いで結果をプロットしています。

実際に値を入れて計算するには、removeO()で剰余項を外す必要がある所がキーポイントです。removeO()を付けないと以下のエラーが発生します。
NameError: name 'O' is not defined

以下は任意点\(a\)違いの結果を2つ例として比較しています。左が\(a=0\)点周りのマクローリン展開、右が\(a=1\)周りのテイラー展開の結果です。

expのテイラー展開

コード上で\(a\)の値を変えるだけで、\(a\)点周りの近似曲線を得ることが出来ました!

\(f(x)=\sin x\)

続いて\(\sin x\)もテイラー展開してみましょう。以下は上記expをsinに変更しただけです。

結果は下図です。同じく\(a=0\)と\(a=1\)を比較してみました。同じように任意点周りの近似がされていることがわかりました。

sin(x)のテイラー展開

まとめ

本記事では\(f(x)\)のテイラー展開の概要を説明し、sympyの簡単なテイラー展開の使い方を紹介しました。

ただコンソール画面に展開結果を表示させるだけではなく、numpy形式としてグラフ化可能な曲線を得る方法も調査してみました。

sympyにも簡易プロッティング機能はあるそうですが、Python使いはnumpyとmatplotlibで出来た方が他の処理と連携しやすく有用でしょう。

最後にテイラー展開でよく例題として出される\(e^{x}\)と\(\sin x\)を実際にsympyで展開し、教科書でお馴染みの結果を得ることが出来ました。

テイラー展開を手計算でやらなくても良いというのは地味に画期的!?
Twitterでも関連情報をつぶやいているので、wat(@watlablog)のフォローお待ちしています!

SNSでもご購読できます。

コメント

コメントを残す

*