Playing with Polynomials

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)])

$\displaystyle \frac{7 x^{3}}{12} - \frac{28 x^{2}}{3} + \frac{499 x}{12} - \frac{99}{2}$

interpolate([(random.randint(1, 100), random.randint(1, 10)) for _ in range(5)])

$\displaystyle - \frac{57187 x^{4}}{9355254048} + \frac{220445 x^{3}}{228176928} - \frac{462110449 x^{2}}{9355254048} + \frac{136997773 x}{148496096} + \frac{8757387}{5303432}$

对函数的近似

$y = \sin(x)$

xs = np.arange(-1, 1, 0.1)
ys = np.sin(xs * 10)
interpolate(list(zip(xs, ys)), draw_points=False)

$\displaystyle - 32.4687383202836 x^{19} + 2.79396974735164 \cdot 10^{-9} x^{18} + 213.314283646643 x^{17} - 1.49011669608375 \cdot 10^{-8} x^{16} - 710.678128842264 x^{15} + 1.49011678604879 \cdot 10^{-8} x^{14} + 1579.66070281342 x^{13} - 1.86264919496656 \cdot 10^{-9} x^{12} - 2497.36639823671 x^{11} + 1.86264654430386 \cdot 10^{-9} x^{10} + 2754.31533635175 x^{9} - 1.16415598560152 \cdot 10^{-10} x^{8} - 1983.97996999821 x^{7} + 1.45519457601099 \cdot 10^{-11} x^{6} + 833.325373093016 x^{5} - 1.71725907896593 \cdot 10^{-18} x^{4} - 166.666480422499 x^{3} + 4.10284968316589 \cdot 10^{-20} x^{2} + 9.99999879997535 x - 2.66458998985356 \cdot 10^{-22}$

Lagrangian Interpolation

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)

$\displaystyle \frac{x \left(y_{0} - y_{1}\right)}{x_{0} - x_{1}} + \frac{x_{0} y_{1} - x_{1} y_{0}}{x_{0} - x_{1}}$

lagrangian_interpolate(3)

$\displaystyle \frac{x^{2} \left(- x_{0} y_{1} + x_{0} y_{2} + x_{1} y_{0} - x_{1} y_{2} - x_{2} y_{0} + x_{2} y_{1}\right)}{x_{0}^{2} x_{1} - x_{0}^{2} x_{2} - x_{0} x_{1}^{2} + x_{0} x_{2}^{2} + x_{1}^{2} x_{2} - x_{1} x_{2}^{2}} + \frac{x \left(x_{0}^{2} y_{1} - x_{0}^{2} y_{2} - x_{1}^{2} y_{0} + x_{1}^{2} y_{2} + x_{2}^{2} y_{0} - x_{2}^{2} y_{1}\right)}{x_{0}^{2} x_{1} - x_{0}^{2} x_{2} - x_{0} x_{1}^{2} + x_{0} x_{2}^{2} + x_{1}^{2} x_{2} - x_{1} x_{2}^{2}} + \frac{x_{0}^{2} x_{1} y_{2} - x_{0}^{2} x_{2} y_{1} - x_{0} x_{1}^{2} y_{2} + x_{0} x_{2}^{2} y_{1} + x_{1}^{2} x_{2} y_{0} - x_{1} x_{2}^{2} y_{0}}{x_{0}^{2} x_{1} - x_{0}^{2} x_{2} - x_{0} x_{1}^{2} + x_{0} x_{2}^{2} + x_{1}^{2} x_{2} - x_{1} x_{2}^{2}}$

Fast Fourier Transformation

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)])

$\displaystyle \frac{x^{4} \cdot \left(12 y_{0} - 48 y_{1} + 72 y_{2} - 48 y_{3} + 12 y_{4}\right)}{288} + \frac{x^{3} \left(- 168 y_{0} + 624 y_{1} - 864 y_{2} + 528 y_{3} - 120 y_{4}\right)}{288} + \frac{x^{2} \cdot \left(852 y_{0} - 2832 y_{1} + 3528 y_{2} - 1968 y_{3} + 420 y_{4}\right)}{288} + \frac{x \left(- 1848 y_{0} + 5136 y_{1} - 5616 y_{2} + 2928 y_{3} - 600 y_{4}\right)}{288} + 5 y_{0} - 10 y_{1} + 10 y_{2} - 5 y_{3} + y_{4}$

尝试 2：$1, 1/2, 1/3, \ldots, 1/n$

xs = [sympy.Rational(1, i) for i in range(1, n + 1)]
try_basis(xs)

$\displaystyle \frac{9765625 x^{4} \cdot \left(\frac{12 y_{0}}{15625} - \frac{48 y_{1}}{15625} + \frac{72 y_{2}}{15625} - \frac{48 y_{3}}{15625} + \frac{12 y_{4}}{15625}\right)}{288} + \frac{9765625 x^{3} \left(- \frac{216 y_{0}}{78125} + \frac{816 y_{1}}{78125} - \frac{1152 y_{2}}{78125} + \frac{144 y_{3}}{15625} - \frac{168 y_{4}}{78125}\right)}{288} + \frac{9765625 x^{2} \cdot \left(\frac{1428 y_{0}}{390625} - \frac{4992 y_{1}}{390625} + \frac{6552 y_{2}}{390625} - \frac{768 y_{3}}{78125} + \frac{852 y_{4}}{390625}\right)}{288} + \frac{9765625 x \left(- \frac{4104 y_{0}}{1953125} + \frac{12864 y_{1}}{1953125} - \frac{15552 y_{2}}{1953125} + \frac{1728 y_{3}}{390625} - \frac{1848 y_{4}}{1953125}\right)}{288} + 15 y_{0} - 40 y_{1} + 45 y_{2} - 24 y_{3} + 5 y_{4}$

尝试 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)

$\displaystyle \frac{x^{4} \cdot \left(4 y_{0} + 2 y_{0} e^{- \frac{3 i \pi}{5}} - 2 y_{0} e^{\frac{i \pi}{5}} - y_{0} e^{\frac{2 i \pi}{5}} - y_{0} e^{\frac{4 i \pi}{5}} - y_{0} e^{- \frac{4 i \pi}{5}} - 2 y_{0} e^{- \frac{i \pi}{5}} + 3 y_{0} e^{\frac{3 i \pi}{5}} - 3 y_{1} + 3 y_{1} e^{- \frac{2 i \pi}{5}} + 3 y_{1} e^{- \frac{i \pi}{5}} - y_{1} e^{\frac{2 i \pi}{5}} + 2 y_{1} e^{\frac{4 i \pi}{5}} + 2 y_{1} e^{\frac{3 i \pi}{5}} - 4 y_{1} e^{- \frac{3 i \pi}{5}} - y_{2} - 3 y_{2} e^{\frac{3 i \pi}{5}} - 3 y_{2} e^{\frac{i \pi}{5}} + 3 y_{2} e^{- \frac{i \pi}{5}} + y_{2} e^{- \frac{3 i \pi}{5}} + y_{2} e^{\frac{2 i \pi}{5}} + 2 y_{2} e^{\frac{4 i \pi}{5}} - 2 y_{2} e^{- \frac{4 i \pi}{5}} - 2 y_{2} e^{- \frac{2 i \pi}{5}} - y_{3} - 4 y_{3} e^{\frac{4 i \pi}{5}} + 3 y_{3} e^{- \frac{4 i \pi}{5}} + y_{3} e^{- \frac{2 i \pi}{5}} - y_{3} e^{- \frac{i \pi}{5}} + y_{3} e^{\frac{3 i \pi}{5}} + 2 y_{3} e^{\frac{i \pi}{5}} - 3 y_{3} e^{- \frac{3 i \pi}{5}} + y_{4} + 4 y_{4} e^{- \frac{3 i \pi}{5}} - 3 y_{4} e^{\frac{3 i \pi}{5}} + y_{4} e^{\frac{4 i \pi}{5}} + y_{4} e^{\frac{2 i \pi}{5}} - 3 y_{4} e^{- \frac{i \pi}{5}} + 3 y_{4} e^{\frac{i \pi}{5}} - 2 y_{4} e^{- \frac{2 i \pi}{5}}\right)}{8 + 10 e^{- \frac{3 i \pi}{5}} - 10 e^{\frac{i \pi}{5}} + 6 e^{- \frac{2 i \pi}{5}} - 8 e^{\frac{4 i \pi}{5}} + 3 e^{\frac{2 i \pi}{5}} - 7 e^{- \frac{4 i \pi}{5}} - 11 e^{- \frac{i \pi}{5}} + 13 e^{\frac{3 i \pi}{5}}} + \frac{x^{3} \left(- 4 y_{0} e^{\frac{4 i \pi}{5}} + 2 y_{0} e^{- \frac{2 i \pi}{5}} - 2 y_{0} e^{\frac{i \pi}{5}} + y_{0} e^{- \frac{3 i \pi}{5}} + y_{0} e^{- \frac{4 i \pi}{5}} - y_{0} e^{- \frac{i \pi}{5}} + 2 y_{0} e^{\frac{2 i \pi}{5}} + 3 y_{0} e^{\frac{3 i \pi}{5}} + 5 y_{1} - 2 y_{1} e^{\frac{i \pi}{5}} + 2 y_{1} e^{- \frac{i \pi}{5}} + y_{1} e^{- \frac{3 i \pi}{5}} - y_{1} e^{\frac{2 i \pi}{5}} - y_{1} e^{\frac{3 i \pi}{5}} + 3 y_{1} e^{\frac{4 i \pi}{5}} - 3 y_{1} e^{- \frac{4 i \pi}{5}} - 4 y_{1} e^{- \frac{2 i \pi}{5}} - y_{2} + 3 y_{2} e^{- \frac{3 i \pi}{5}} - 2 y_{2} e^{\frac{3 i \pi}{5}} + y_{2} e^{- \frac{4 i \pi}{5}} - 2 y_{2} e^{- \frac{i \pi}{5}} + 2 y_{2} e^{\frac{4 i \pi}{5}} + 2 y_{2} e^{\frac{i \pi}{5}} + 2 y_{2} e^{\frac{2 i \pi}{5}} - 3 y_{2} e^{- \frac{2 i \pi}{5}} - 3 y_{3} + 5 y_{3} e^{- \frac{2 i \pi}{5}} - 3 y_{3} e^{\frac{2 i \pi}{5}} + 2 y_{3} e^{- \frac{i \pi}{5}} - 2 y_{3} e^{- \frac{4 i \pi}{5}} + 3 y_{3} e^{\frac{4 i \pi}{5}} - 2 y_{3} e^{- \frac{3 i \pi}{5}} - y_{4} - 4 y_{4} e^{\frac{4 i \pi}{5}} + 3 y_{4} e^{- \frac{4 i \pi}{5}} - y_{4} e^{- \frac{i \pi}{5}} + 2 y_{4} e^{\frac{i \pi}{5}} - 3 y_{4} e^{- \frac{3 i \pi}{5}}\right)}{8 + 10 e^{- \frac{3 i \pi}{5}} - 10 e^{\frac{i \pi}{5}} + 6 e^{- \frac{2 i \pi}{5}} - 8 e^{\frac{4 i \pi}{5}} + 3 e^{\frac{2 i \pi}{5}} - 7 e^{- \frac{4 i \pi}{5}} - 11 e^{- \frac{i \pi}{5}} + 13 e^{\frac{3 i \pi}{5}}} + \frac{x^{2} \cdot \left(2 y_{0} e^{- \frac{3 i \pi}{5}} + 2 y_{0} e^{- \frac{2 i \pi}{5}} - 2 y_{0} e^{\frac{i \pi}{5}} - y_{0} e^{\frac{4 i \pi}{5}} - y_{0} e^{- \frac{4 i \pi}{5}} + y_{0} e^{\frac{2 i \pi}{5}} - 2 y_{0} e^{- \frac{i \pi}{5}} + 3 y_{0} e^{\frac{3 i \pi}{5}} - y_{1} - 3 y_{1} e^{\frac{2 i \pi}{5}} - 2 y_{1} e^{\frac{3 i \pi}{5}} + 2 y_{1} e^{- \frac{2 i \pi}{5}} - 3 y_{1} e^{\frac{4 i \pi}{5}} + 2 y_{1} e^{- \frac{4 i \pi}{5}} - 2 y_{1} e^{- \frac{i \pi}{5}} + 3 y_{1} e^{\frac{i \pi}{5}} - 2 y_{1} e^{- \frac{3 i \pi}{5}} + 3 y_{2} + 3 y_{2} e^{- \frac{i \pi}{5}} - y_{2} e^{\frac{2 i \pi}{5}} + y_{2} e^{- \frac{2 i \pi}{5}} + 4 y_{2} e^{\frac{3 i \pi}{5}} - 4 y_{2} e^{- \frac{3 i \pi}{5}} - y_{3} + 3 y_{3} e^{- \frac{3 i \pi}{5}} - 2 y_{3} e^{\frac{3 i \pi}{5}} + y_{3} e^{- \frac{4 i \pi}{5}} - 2 y_{3} e^{- \frac{i \pi}{5}} + 2 y_{3} e^{\frac{4 i \pi}{5}} + 2 y_{3} e^{\frac{i \pi}{5}} + 2 y_{3} e^{\frac{2 i \pi}{5}} - 3 y_{3} e^{- \frac{2 i \pi}{5}} - y_{4} - 3 y_{4} e^{\frac{3 i \pi}{5}} - 3 y_{4} e^{\frac{i \pi}{5}} + 3 y_{4} e^{- \frac{i \pi}{5}} + y_{4} e^{- \frac{3 i \pi}{5}} + y_{4} e^{\frac{2 i \pi}{5}} + 2 y_{4} e^{\frac{4 i \pi}{5}} - 2 y_{4} e^{- \frac{4 i \pi}{5}} - 2 y_{4} e^{- \frac{2 i \pi}{5}}\right)}{8 + 10 e^{- \frac{3 i \pi}{5}} - 10 e^{\frac{i \pi}{5}} + 6 e^{- \frac{2 i \pi}{5}} - 8 e^{\frac{4 i \pi}{5}} + 3 e^{\frac{2 i \pi}{5}} - 7 e^{- \frac{4 i \pi}{5}} - 11 e^{- \frac{i \pi}{5}} + 13 e^{\frac{3 i \pi}{5}}} + \frac{x \left(3 y_{0} e^{- \frac{3 i \pi}{5}} + 2 y_{0} e^{- \frac{2 i \pi}{5}} - y_{0} e^{\frac{i \pi}{5}} - y_{0} e^{\frac{4 i \pi}{5}} + y_{0} e^{\frac{3 i \pi}{5}} + 2 y_{0} e^{\frac{2 i \pi}{5}} - 4 y_{0} e^{- \frac{i \pi}{5}} - 4 y_{0} e^{- \frac{4 i \pi}{5}} - y_{1} + 5 y_{1} e^{- \frac{4 i \pi}{5}} + 2 y_{1} e^{- \frac{3 i \pi}{5}} - y_{1} e^{\frac{3 i \pi}{5}} + y_{1} e^{- \frac{i \pi}{5}} - y_{1} e^{\frac{4 i \pi}{5}} + 3 y_{1} e^{\frac{2 i \pi}{5}} - 4 y_{1} e^{- \frac{2 i \pi}{5}} - y_{2} - 3 y_{2} e^{\frac{2 i \pi}{5}} + 3 y_{2} e^{- \frac{2 i \pi}{5}} - 3 y_{2} e^{\frac{4 i \pi}{5}} + 2 y_{2} e^{- \frac{4 i \pi}{5}} - y_{2} e^{\frac{3 i \pi}{5}} - 2 y_{2} e^{- \frac{i \pi}{5}} + 3 y_{2} e^{\frac{i \pi}{5}} - 2 y_{2} e^{- \frac{3 i \pi}{5}} + 5 y_{3} - 2 y_{3} e^{\frac{i \pi}{5}} + 2 y_{3} e^{- \frac{i \pi}{5}} + y_{3} e^{- \frac{3 i \pi}{5}} - y_{3} e^{\frac{2 i \pi}{5}} - y_{3} e^{\frac{3 i \pi}{5}} + 3 y_{3} e^{\frac{4 i \pi}{5}} - 3 y_{3} e^{- \frac{4 i \pi}{5}} - 4 y_{3} e^{- \frac{2 i \pi}{5}} - 3 y_{4} + 3 y_{4} e^{- \frac{2 i \pi}{5}} + 3 y_{4} e^{- \frac{i \pi}{5}} - y_{4} e^{\frac{2 i \pi}{5}} + 2 y_{4} e^{\frac{4 i \pi}{5}} + 2 y_{4} e^{\frac{3 i \pi}{5}} - 4 y_{4} e^{- \frac{3 i \pi}{5}}\right)}{8 + 10 e^{- \frac{3 i \pi}{5}} - 10 e^{\frac{i \pi}{5}} + 6 e^{- \frac{2 i \pi}{5}} - 8 e^{\frac{4 i \pi}{5}} + 3 e^{\frac{2 i \pi}{5}} - 7 e^{- \frac{4 i \pi}{5}} - 11 e^{- \frac{i \pi}{5}} + 13 e^{\frac{3 i \pi}{5}}} + \frac{4 y_{0} + 2 y_{0} e^{- \frac{3 i \pi}{5}} - 3 y_{0} e^{\frac{i \pi}{5}} - y_{0} e^{\frac{2 i \pi}{5}} - y_{0} e^{\frac{4 i \pi}{5}} - 2 y_{0} e^{- \frac{i \pi}{5}} - 2 y_{0} e^{- \frac{4 i \pi}{5}} + 3 y_{0} e^{\frac{3 i \pi}{5}} + 3 y_{1} e^{- \frac{3 i \pi}{5}} + 3 y_{1} e^{- \frac{2 i \pi}{5}} - y_{1} e^{\frac{i \pi}{5}} - y_{1} e^{\frac{4 i \pi}{5}} + 2 y_{1} e^{\frac{3 i \pi}{5}} + 2 y_{1} e^{\frac{2 i \pi}{5}} - 4 y_{1} e^{- \frac{i \pi}{5}} - 4 y_{1} e^{- \frac{4 i \pi}{5}} + 2 y_{2} e^{- \frac{3 i \pi}{5}} - 2 y_{2} e^{\frac{i \pi}{5}} + y_{2} e^{- \frac{2 i \pi}{5}} - y_{2} e^{\frac{4 i \pi}{5}} - y_{2} e^{- \frac{4 i \pi}{5}} + y_{2} e^{\frac{2 i \pi}{5}} - 2 y_{2} e^{- \frac{i \pi}{5}} + 2 y_{2} e^{\frac{3 i \pi}{5}} - 4 y_{3} e^{\frac{4 i \pi}{5}} - 2 y_{3} e^{\frac{i \pi}{5}} + y_{3} e^{- \frac{3 i \pi}{5}} + y_{3} e^{- \frac{2 i \pi}{5}} + y_{3} e^{- \frac{4 i \pi}{5}} - y_{3} e^{- \frac{i \pi}{5}} + 2 y_{3} e^{\frac{3 i \pi}{5}} + 2 y_{3} e^{\frac{2 i \pi}{5}} + 4 y_{4} + 2 y_{4} e^{- \frac{3 i \pi}{5}} - 2 y_{4} e^{\frac{i \pi}{5}} - y_{4} e^{\frac{2 i \pi}{5}} + y_{4} e^{- \frac{2 i \pi}{5}} - y_{4} e^{\frac{4 i \pi}{5}} - y_{4} e^{- \frac{4 i \pi}{5}} - 2 y_{4} e^{- \frac{i \pi}{5}} + 4 y_{4} e^{\frac{3 i \pi}{5}}}{8 + 10 e^{- \frac{3 i \pi}{5}} - 10 e^{\frac{i \pi}{5}} + 6 e^{- \frac{2 i \pi}{5}} - 8 e^{\frac{4 i \pi}{5}} + 3 e^{\frac{2 i \pi}{5}} - 7 e^{- \frac{4 i \pi}{5}} - 11 e^{- \frac{i \pi}{5}} + 13 e^{\frac{3 i \pi}{5}}}$

同时求多个值

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)
]

$\displaystyle \left[ A + B + C + D + E, \ A + B e^{\frac{i \pi}{5}} + C e^{\frac{2 i \pi}{5}} + D e^{\frac{3 i \pi}{5}} + E e^{\frac{4 i \pi}{5}}, \ A + B e^{\frac{2 i \pi}{5}} + C e^{\frac{4 i \pi}{5}} + D e^{- \frac{4 i \pi}{5}} + E e^{- \frac{2 i \pi}{5}}, \ A + B e^{\frac{3 i \pi}{5}} + C e^{- \frac{4 i \pi}{5}} + D e^{- \frac{i \pi}{5}} + E e^{\frac{2 i \pi}{5}}, \ A + B e^{\frac{4 i \pi}{5}} + C e^{- \frac{2 i \pi}{5}} + D e^{\frac{2 i \pi}{5}} + E e^{- \frac{4 i \pi}{5}}\right]$