import sympy
import random
import numpy as np
import matplotlib.pyplot as plt
import sympy.polys
plt.rc('text', usetex=True)
x = sympy.symbols('x')
def plot(f, points=[], draw_points=True):
"""Plot a polynomial f."""
xmin = min([x_ for x_, _ in points]) - 0.1
xmax = max([x_ for x_, _ in points]) + 0.1
xs = np.arange(xmin, xmax, (xmax - xmin) / 100)
ys = [f.subs(x, x_) for x_ in xs]
plt.grid(True)
plt.plot(xs, ys)
if draw_points:
plt.scatter(
[x_ for x_, y_ in points],
[y_ for x_, y_ in points],
)
for x_, y_ in points:
plt.text(x_, y_, f'$({x_},{y_})$', va='bottom', ha='center')
plt.title(f'$y = {sympy.latex(f)}$')
def interpolate(points, draw_points=True):
n = len(points) - 1
vars = sympy.symbols([f'a{i}' for i in range(n + 1)])
cons = []
for x_, y in points:
cons.append(
sum(
vars[k] * (x_ ** (n-k)) for k, _ in enumerate(points)
) - y
)
sol = list(sympy.linsolve(cons, vars[:n+1]))[0]
f = sum(
coef * (x ** (n-k))
for k, coef in enumerate(sol)
)
plot(f, points, draw_points)
return f
interpolate([(2, 1), (3, 7), (5, -2), (6, -10)])
interpolate([(random.randint(1, 100), random.randint(1, 10)) for _ in range(5)])
$y = \sin(x)$
xs = np.arange(-1, 1, 0.1)
ys = np.sin(xs * 10)
interpolate(list(zip(xs, ys)), draw_points=False)
def lagrangian_interpolate(n):
xs = [sympy.symbols(f'x{i}') for i in range(n)]
ys = [sympy.symbols(f'y{i}') for i in range(n)]
vars = [sympy.symbols(v) for v in 'abcdefghijklmnopqrst'[:n]]
power = list(range(n))
cons = []
for x_, y in zip(xs, ys):
cons.append(
sum(
v * (x_ ** k) for v, k in zip(vars, power)
) - y
)
sol = list(sympy.linsolve(cons, vars))
assert len(sol) == 1
sol = sol[0]
f = (sum(
v * (x ** k) for v, k in zip(sol, power)
))
for x_, y_ in zip(xs, ys):
assert sympy.simplify(f.subs(x, x_)) == y_
return f
lagrangian_interpolate(2)
lagrangian_interpolate(3)
def try_basis(xs):
n = len(xs)
f = lagrangian_interpolate(n)
for i, x_ in enumerate(xs):
f = f.subs(f'x{i}', x_)
return f
尝试 1:$1, 2, 3, \ldots, n$
n = 5
try_basis([i + 1 for i in range(n)])
尝试 2:$1, 1/2, 1/3, \ldots, 1/n$
xs = [sympy.Rational(1, i) for i in range(1, n + 1)]
try_basis(xs)
尝试 3: $\omega^0, \omega^1, \ldots, \omega^n$
w = sympy.E ** (sympy.I * sympy.pi / n)
xs = [w**i for i in range(n)]
try_basis(xs)
同时求多个值
coefs = [sympy.symbols(a) for a in list('ABCDEFGHIJKLMNOPQRSTUVWXYZ')]
f = 0
for i in range(n):
f += coefs[i] * (x**i)
[
f.subs(x, w**i) for i in range(n)
]