Copyright (C) 2020 Andreas Kloeckner
import numpy as np
import matplotlib.pyplot as plt
x = np.linspace(-1, 1, 5)
xx = np.linspace(-1, 1, 100)
def f(x):
return np.sin(5/(1.6-x))
plt.plot(xx, f(xx))
plt.plot(x, f(x), "o")
dmat = x.reshape(-1, 1) - x
np.fill_diagonal(dmat, 1)
dmat[0, 2], x[0]-x[2]
w = 1/dmat.prod(axis=1)
w
def bary_first_form(z, f):
z = z[:, np.newaxis]
ell = (z-x).prod(axis=-1)
return ell * (w / (z - x) * f(x)).sum(axis=-1)
plt.plot(xx, bary_first_form(xx, f))
plt.plot(xx, f(xx))
def bary_second_form(z, f, weights=None):
if weights is None:
weights = w
z = z[:, np.newaxis]
num = (weights / (z - x) * f(x)).sum(axis=-1)
denom = (weights / (z - x)).sum(axis=-1)
return num/denom
wr = w.copy()
wr[:] = np.random.randn(len(x))
plt.plot(xx, bary_second_form(xx, f))
plt.plot(xx, bary_second_form(xx, f, weights=wr))
plt.plot(xx, f(xx))
plt.plot(x, f(x), "o")
plt.ylim([-2, 2])