ラグランジュ補間を通して学ぶSympyライブラリの使い方 ~計算物理学I (朝倉書店)を参考にpythonを使って~

前置き

　今回はラグランジュ補完多項式を作り出すプログラムでSympyライブラリを使用してみたので今後の参考のためにまとめます。

ラグランジュ補間とは

　ラグランジュ補間とは、データ点に式をフィットする方法の一つで、与えられた$i = 1, \cdots ,n$のデータ点($x_i$,$g_i$)を使用して、全てのデータ点を通る$(n-1)$次の補間多項式を作る方法のことです。
　この方法で求められる$(n-1)$次多項式の$\bar{g}(x)$は、以下のように表されます。

\bar{g}(x) = \sum_{i=1}^n g_i \lambda_i
\\
\lambda_i = \prod_{j(\ne i)=1}^{n} \frac{x - x_i}{x_i - x_j}

　各データ点がほとんど誤差を含まず正確なものであるならば、与えられたデータ点の区間内では信頼して使用することができます。ただ、データ点の区間外で補間多項式を使用しようとするとかなりおかしな結果がでることがあるので注意が必要です。

Sympyライブラリを使用してラグランジュ補間の補間多項式を作ってみよう

　今回はSympyライブラリを使用してラグランジュ補間多項式を求めるプログラムを作製しますが、まずは今回使用するSympyライブラリの簡単な説明をします。

Sympyライブラリとは

　SympyライブラリとはPythonの代数計算のライブラリで、代数の$x$や$y$そのものを計算の中に組み込み使用することができます。具体的には以下のような計算が実行できるのです。

>>>import sympy

#xとyを代数として宣言して使用
>>> x = sympy.Symbol('x')
>>> y = sympy.Symbol('y')
>>> x+y+y
x + 2*y

#x=2として計算
>>>(x+y+y).subs(x,2)
2*y + 2

#代数表現を展開
>>>sympy.expand((x + y)**4)
x**4 + 4*x**3*y + 6*x**2*y**2 + 4*x*y**3 + y**4

より詳細な内容は https://www.sympy.org/en/index.html 等をご参考ください。
　今回は補間多項式そのものの形を知りたかったので、代数計算に強いSympyライブラリを使用するプログラムを作製してみました。

##ラグランジュ補間多項式を作り出すプログラム
　さてさてお待ちかねの(?)、与えられたデータ点からラグランジュ補間多項式を作り出すプログラムは以下のようになりました。あまりひねりが無く数式のまんまですね。ライブラリは先ほどのSympyと、Numpyを使用していますね。

import sympy
import numpy as np

#データ点list_x及びlist_gを与えることで、ラグランジュ補間多項式gbarを代数式の形で求める。
def Lagrange_interpolation(list_x, list_g):
    n = list_x.shape[0]
    x = sympy.Symbol('x')   #xを代数として使用することを宣言
    list_lambda = np.array([])
    
    #補間多項式のλを計算
    for i in range(n):
        lambda1 = 1
        for j in range(n):
            if i==j:
                lambda1 *= 1
            else:
                lambda1 *= (x - list_x[j])/(list_x[i]-list_x[j])
        list_lambda = np.append(list_lambda,lambda1)
    #各λとlist_gから補間多項式gbarを作製
    gbar = sympy.expand(np.sum(list_lambda*list_g))
    return gbar

プログラムの実践例

　さて、せっかく作ったプログラムなので実践していきます。
　今回は、$ g(x) = \ln(x^2 +1) \times \sin(x) $という関数の$x = 1,2, \cdots 10$の10個のデータ点を補間多項式によって補ってみます。
　まずはデータ点を作ります。

>>>import numpy as np
>>>list_x = np.arange(1,11,1)
>>>list_g = np.log(list_x**2+1)*np.sin(list_x)
>>>print(list_x, '\n',list_g)
[ 1  2  3  4  5  6  7  8  9 10] 
 [ 0.58326324  1.46345775  0.32494083 -2.14418293 -3.12426786 -1.00894643
  2.57014669  4.12996447  1.81609046 -2.51072299]

　作り出したデータ点を図にすると以下のようになります。

　上記のデータ点を使用して、ラグランジュ補間多項式を作ってみます。

>>>gbar = Lagrange_interpolation(list_x, list_g)
>>>gbar
3.41706115246062e-6*x**9 - 0.000238318344027363*x**8 + 0.00643158890937556*x**7 - 0.0882402687361221*x**6 + 0.666958583912063*x**5 - 2.80385499635639*x**4 + 6.48008209579984*x**3 - 8.78330812338027*x**2 + 8.05300325374583*x - 2.94757399196828

　どうやら目論見通り9次の多項式が作り出されたようですが、この多項式はちゃんとデータ点を補間できているのでしょうか。
　補間多項式とデータ点を合わせてプロットした図が以下になります。

結果としては問題なく補間できました。もともとの式が単純なためによく再現できたようです。(データ点の配置によっては、うまく補間できない場合があります。実際はもう少しデータ点を狭い範囲で区切って適宜補間していく必要があるかもしれません。)
　ところで、多項式を作り出すために用いたデータ点の外では、補間多項式はどのように振る舞うでしょう。データ点の外を含む領域、$x = 0$から$x = 11$の範囲で$ g(x) = \ln(x^2 +1) \times \sin(x) $と補間多項式を比較したのが以下の図となります。(青い丸が$g(x)$、赤いラインが補間多項式です。)

用いたデータ点の範囲の外では$g(x)$と補間多項式が大きくずれていることが分かります。このように、補間多項式を使ってデータを外挿するのは難しいようです。

まとめ

　Sympyライブラリを使用してラグランジュ補間多項式を作り出すプログラムを作製してみました。関数形が露わになるのでプログラムの確認にも使用できますね。また、Sympyライブラリは代数の微分にも対応しているので、関数の極大極小を求めることもできます。求めた多項式の詳細を知りたい場合にはとても相性がよいライブラリなのではないでしょうか。

参考書籍
R.H.Landau・他 (2018)『実践Pythonライブラリー 計算物理学I -数値計算の基礎/HPC/フーリエウェーブレット解析』 (小柳義夫・他訳) 朝倉書店