LagrangePolynomial_ErrorAnalysis.png

import numpy as np 
import matplotlib.pyplot as plt 
import sympy as sb

# 森正武, "数値解析", 共立出版株式会社, 第2版, 2002.
# 4.12 ラグランジュ補間公式の標本点の分布
# 等間隔な標本点の分布(P170 - P173)

def LagrangePolynomial(X,Y):
    def l(j,N):
        f = 1
        for i in range(N+1):
            if i is not j:
                f *= (x - xN[i]) / (xN[j] - xN[i]) 
        return f

    def L(N):
        f = 0
        for j in range(N+1):
            f += yN[j] * l(j,N)
        return f

    N = len(X)-1

    p = L(N)
    for i in range(N+1):
        p = p.subs(xN[i],X[i])
        p = p.subs(yN[i],Y[i])
    pp = sb.lambdify(x,p)

    return pp
   
def Error(X):
    def Fn(x,N):
        y = sb.prod((x-xN[i]) for i in range(N))
        return y

    # 誤差の特性関数
    def Φ(N,x,z):
        y = Fn(x,N) / ( (z-x) * Fn(z,N) )
        return y
        
    # 点数
    N = len(X)

    # 特異点
    ζj = list(sb.singularities(F,z))

    # 留数
    Rj = []
    for ζ in ζj:
        R = sb.residue(F,z,ζ)
        Rj.append(R)
    
    # 誤差
    ε = 0
    for i,R in enumerate(Rj):
        ε += - Φ(N,x,ζj[i]) * R
    
    # ノードの代入
    for i in range(N):
        ε = ε.subs(xN[i],X[i])

    εε = sb.lambdify(x,ε)

    return εε

# 変数
xN = sb.IndexedBase("xN")
yN = sb.IndexedBase("y")
z = sb.Symbol('z')
x = sb.Symbol('x')

# 関数
F = 1/(1+25*z**2)
f = sb.lambdify(z,F)

# 区間[a,b]
a = -1
b =  1

# 点数
N = 4

# 等間隔ノード
k = np.arange(-N,N+1)
x_ES = k/N
y_ES = f(x_ES)

xx = np.linspace(a,b,1000)
yy = f(xx)

ε = Error(x_ES)
p = LagrangePolynomial(x_ES,y_ES)

yy_LP = p(xx)

# 誤差
E1 = yy - yy_LP
E2 = ε(xx).real

plt.close('all')
plt.figure(figsize=(6,6))
plt.plot(xx,yy,'k',lw=3,label='$f(x)=\\frac{1}{1+25x^2}$')
plt.plot(xx,yy_LP,lw=1,label='$p_{%d}(x)$' % (2*N+1))
plt.plot(x_ES,y_ES,'C3o')

plt.plot(xx,E1,'k--',lw=2,label='$f(x)-p_{%d}(x)$' % (2*N+1))
plt.plot(xx,E2,'--',lw=1,label='$\\varepsilon(x)$')
plt.plot(x_ES,np.zeros_like(x_ES),'ko')


plt.xlabel('$x$')
plt.ylabel('$y$')
plt.grid()
plt.legend(loc=1,fontsize=12)

plt.ylim([-1.5,1.5])
plt.tight_layout()

plt.savefig('LagrangePolynomial_ErrorAnalysis.png')