<jupyter_start><jupyter_text>插值<jupyter_code>from sympy import *<jupyter_output><empty_output><jupyter_text>拉格朗日插值<jupyter_code>def lagrange_interpolate(points: list, simplify_result=True, verbose=False):
    """拉格朗日插值

    Args:
        points: list, [(x1, y1), (x2, y2), ..., (xn, yn)]
        simplify_result: bool, 化简最终结果, default True
        verbose: bool, 输出每一步的结果, default False

    Returns: 
        L: sympy object of Symbol('x'), 插值多项式 $L(x)$
    """
    x = Symbol('x')
    L = 0  # 插值多项式
    for i, point in enumerate(points):
        xi, yi = point
        li = 1
        for j in range(len(points)):
            if j == i:
                continue
            xj, yj = points[j]
            li *= (x - xj) / (xi - xj)
        L += yi * li
        if verbose:
            print(f"l_{i}(x) = ", simplify(expand(yi * li)))
        
    if simplify_result:
        L = simplify(L)
    return L
points = [(11, 0.190809), (12, 0.207912), (13, 0.224951)]

L = lagrange_interpolate(points)
print('插值多项式:', L)
print('Hypothesis:  ', L.subs(Symbol('x'), 11.5))
print('Actual vaule:', sin(rad(11.5)).evalf())<jupyter_output>插值多项式: -3.20000000000042e-5*x**2 + 0.0178390000000004*x - 0.00154799999999966
Hypothesis:   0.199368500000004
Actual vaule: 0.199367934417197<jupyter_text>牛顿插值<jupyter_code># 差商缓存
__difference_quotient_cache = {}

def difference_quotient(f, xs: list, verbose=False):
    """求差商 $f[xs...]$
    
    Args:
        f: function, 函数
        xs: list, 要计算的差商的 $f[x0, x1, ...]$ 的参数值：[x0, x1, ...]
        verbose: bool, 打印出每一步计算差商的值, default False

    Returns:
        dq: sympy object, 差商值
    """

    # 尝试从缓存读取
    __key = str([(x, f(x)) for x in xs])
    if __key in __difference_quotient_cache:
        if verbose:
            print(f"cached: f{xs}: {__difference_quotient_cache[__key]}")
        return __difference_quotient_cache[__key]

    if len(xs) == 1: # 0阶
        dq = sympify(f(xs[0]))
    else: # n 阶
        dq_h = difference_quotient(f, xs[:-1], verbose)
        dq_l = difference_quotient(f, xs[1: ], verbose)
        dq = (dq_l - dq_h) / (xs[-1] - xs[0])
    
    if verbose:
        print(f"f{xs}: {dq}")

    # 写入缓存
    __difference_quotient_cache[__key] = dq

    return dq
def newton_interpolate(points: list, N_start=0, points_start=[], simplify_result=True, verbose=False):
    """牛顿插值

    Args:
        points: [(x1, y1), (x2, y2), ..., (xn, yn)]： 插值点

        N_start: a sympy object of Symbol('x'), 起始插值多项式。
        points_start: [(x1, y1), (x2, y2), ..., (xn, yn)]: 计算 N_start 用的插值点
                 该函数中做的插值会「承袭」这个用 points_start 计算的 N_start，在其基础上用新增插值点 points 去改进。
                 default 0 (从头开始构建).

        simplify_result: bool, 化简最终结果, default True
        verbose: 打印出每一步计算差商的值, default False

    Returns: 
        N: a sympy object of Symbol('x'), 插值多项式 $N(x)$
    """
    x = Symbol('x')
    
    # sympify points
    # for i in range(len(points)):
    #     points[i] = sympify(points[i][0]), sympify(points[i][1])
    
    f = lambda x: dict(points + points_start)[x]
    
    xs = [p[0] for p in points_start] # 承袭的插值点
    N = N_start   # 承袭的插值多项式
    for point in points:
        xs.append(point[0])
        N += difference_quotient(f, xs, verbose) * prod([x - xi for xi in xs[:-1]])
    
    if simplify_result:
        N = simplify(N)
    return N
points = [(11, 0.190809), (12, 0.207912), (13, 0.224951)]

N = newton_interpolate(points)
print('插值多项式:', N)
print('Hypothesis:  ', N.subs(Symbol('x'), 11.5))
print('Actual vaule:', sin(rad(11.5)).evalf())

# 承袭之前的 N:
print('承袭:')
new_points = [(11.2, sin(rad(11.2)).evalf()), (11.7, sin(rad(11.7)).evalf())]
N1 = newton_interpolate(new_points, N_start=N, points_start=points, verbose=True)
print('插值多项式:', N1)
print('Hypothesis:  ', N1.subs(Symbol('x'), 11.5).evalf())
print('Actual vaule:', sin(rad(11.5)).evalf())
__difference_quotient_cache<jupyter_output><empty_output><jupyter_text>三弯矩法<jupyter_code>import numpy as np

def spline3_interpolate(points: list, simplify_result=True):
    """三弯矩法插值

    三弯矩法进行三次样条插值。使用自由边界条件。

    Args:
        points: list, [(x1, y1), (x2, y2), ..., (xn, yn)]
        simplify_result: bool, 化简最终结果, default True

    Returns: 
        L: sympy Piecewise object of Symbol('x'), 插值多项式 $L(x)$
    """
    # 排序给定的点
    ps = sorted(points, key=lambda p: p[0])
    n = len(points)

    # points to a function
    _f_dict = dict(ps)
    def f(x):
        return _f_dict[x]

    # $h_k = x_{k+1} - x_k$
    def h(k):
        return ps[k+1][0] - ps[k][0]

    hks = [h(0)]

    # 用方程 D * M = d 解出 M

    D = np.zeros((n, n))
    d = np.zeros(n)

    for k in range(1, n-1):
        # $h_k$, $h_{k-1}$
        hks.append(h(k))
        hk, hks1 = hks[k], hks[k-1]

        # $\mu_k$ -> mu, $\lambda_k$ -> ld
        _fra = hks1 + hk
        mu = hks1 / _fra
        ld = hk / _fra

        # $\mu_kM_{k-1}+2M_k+\lambda_kM_{k+1} = d_k$
        D[k, k-1] = mu
        D[k, k]   = 2
        D[k, k+1] = ld
        d[k] = 6 * difference_quotient(f, [ps[k-1][0], ps[k][0], ps[k+1][0]])

    # 边界条件
    # Natural Boundary
    D[0, 0] = 1
    D[n-1, n-1] = 1

    d[0] = 0
    d[n-1] = 0

    # 解出 M
    M = np.linalg.solve(D, d)

    # 插值函数
    piecewises = []
    for k in range(n-1):
        s  = M[k] * (ps[k+1][0] - _x) ** 3 / (6 * hks[k])
        s += M[k+1] * (_x - ps[k][0]) ** 3 / (6 * hks[k])
        s += (ps[k][1] - M[k] * hks[k]**2 / 6) * (ps[k+1][0] - _x) / hks[k]
        s += (ps[k+1][1] - M[k+1] * hks[k]**2 / 6) * (_x - ps[k][0]) / hks[k]
        if simplify_result:
            s = simplify(s)
        piecewises.append((s, And(_x >= ps[k][0], _x <= ps[k+1][0])))

    return Piecewise(*piecewises)<jupyter_output><empty_output><jupyter_text>题目2．区间 $[-5, 5]$ 作等距划分：...<jupyter_code>def f(x):
    return 1 / (1 + x ** 2)

def get_points(s, n, f):
    """在闭区间 s 作 n 等距划分, 输出 n+1 个插值点 [(x0, f(x0)), (x1, f(x1)), ..., (xn, f(xn))]
    e.g. 
        get_points([-5, 5], 2, lambda x: 1 / (1 + x ** 2))
    output: [(-5.0, 0.038461538461538464), (0.0, 1.0), (5.0, 0.038461538461538464)]
    """
    h = (s[1] - s[0]) / n
    xs = [s[0] + k * h for k in range(n+1)]
    points = [(x, f(x)) for x in xs]
    return points
hps = {}  # 多项式插值
hss = {}  # 三次样条插值
for n in [5, 10, 20]:
    ps = get_points([-5, 5], n, f)
    # 多项式插值
    # h = lagrange_interpolate(ps)
    h = newton_interpolate(ps)
    hps[n] = h

    # 三次样条插值
    s = spline3_interpolate(ps)
    hss[n] = s

print('多项式插值：')
for k in hps:
    print(f'n = {k}:\t{hps[k]}')
print('\n三次样条插值：')
for k in hss:
    print(f'n = {k}:\t{hss[k]}')
# 画图
def show_plot_range(r, hs, title):
    """画图比较插值结果
    r：区间
    hs: 插值结果组 {n: h(x)}
    title: 图标题
    """
    x = Symbol('x')
    colors = (c for c in [0, 'b', 'g', 'r'])

    plt = plot(f(x), (x, *r), show=False, line_color=next(colors), label='f(x)', legend=True, title=title)
    for n in hs:
        p = plot(hs[n], (x, *r), show=False, line_color=next(colors), label=f'n={n}', legend=True)
        plt.extend(p)
    
    plt.show()
# 画图：多项式插值
show_plot_range((-5, 5), hps, 'Polynomial Interpolation')
show_plot_range((-5, -3), hps, 'Polynomial Interpolation')
show_plot_range((-1, 1), hps, 'Polynomial Interpolation')
# 画图：三次样条插值
show_plot_range((-5, 5), hss, 'Cubic Spline Interpolation')
show_plot_range((-5, -3), hss, 'Cubic Spline Interpolation')
show_plot_range((-1, 1), hss, 'Cubic Spline Interpolation')<jupyter_output><empty_output><jupyter_text>拟合<jupyter_code>def normalEquationFit(X, y):
    """最小二乘拟合
    Compute the normal equations:

    $\theta = (X^TX)^{-1}X^Ty$

    Args:
        X: np.array, $$
            X = \left[\begin{array}{c}
            - (x^{(1)})^T - \\
            - (x^{(2)})^T - \\
            \vdots \\
            - (x^{(m)})^T - \\
            \end{array}\right]
        $$
        y: np.array, $y = [y^{(1)},y^{(2)},\dots,y^{(m)}]^T$

    return theta.
    """
    return np.dot(np.dot(np.linalg.pinv(np.dot(X.T, X)), X.T), y)<jupyter_output><empty_output><jupyter_text>题目：3.4.2 数据拟合 1.<jupyter_code>x = np.array([1, 2, 3, 4, 5])
y = np.array([4, 4.5, 6, 8, 8.5])

X = np.c_[np.ones((len(x), 1)), x]

theta = normalEquationFit(X, y)
print(theta)

pred_f = lambda t: theta[0] + theta[1] * t
# 用梯度下降检查

a, b = 0, 0

num_epoch = 10000
learning_rate = 5e-4

for e in range(num_epoch):
    y_pred = a * x + b
    
    grad_a = 2 * (y_pred - y).dot(x)
    grad_b = 2 * (y_pred - y).sum()
    
    a = a - learning_rate * grad_a
    b = b - learning_rate * grad_b

print(a, b)
# 画图
import matplotlib.pyplot as plt

pred_x = [k/10 for k in range(50)]
pred_y = [pred_f(x) for x in pred_x]

plt.plot(x, y, 'bo', label='date')
plt.plot(pred_x, pred_y, 'r-', label='pred')
plt.legend()
plt.show()<jupyter_output><empty_output><jupyter_text>题目：3.4.2 数据拟合 2.<jupyter_code>x = np.array([2, 3, 4, 7, 8, 10, 11, 14, 16, 18, 19])
y = np.array([106.42, 108.2, 109.5, 110, 109.93, 110.49, 110.59, 110.6, 110.76, 111, 111.2])<jupyter_output><empty_output><jupyter_text>(1) $y(x) = c_0 + c_1x + c_2x^2$数据预处理：<jupyter_code># X: [1, x, x**2]
X = np.c_[np.ones((len(x), 1)), x, x ** 2]
print(X)<jupyter_output>[[  1.   2.   4.]
 [  1.   3.   9.]
 [  1.   4.  16.]
 [  1.   7.  49.]
 [  1.   8.  64.]
 [  1.  10. 100.]
 [  1.  11. 121.]
 [  1.  14. 196.]
 [  1.  16. 256.]
 [  1.  18. 324.]
 [  1.  19. 361.]]<jupyter_text>定义预测函数、代价函数：<jupyter_code>def hypothesis(x, theta):
    return np.dot(theta, np.array([1, x, x**2]))

def cost(x, y, hypothesis, theta):
    return np.sum( (hypothesis(x, theta) - y) ** 2 )<jupyter_output><empty_output><jupyter_text>最小二乘正规方程拟合：<jupyter_code># theta: [c0, c1, c2]
theta = normalEquationFit(X, y)
theta<jupyter_output><empty_output><jupyter_text>计算代价：<jupyter_code>cost(x, y, hypothesis, theta)<jupyter_output><empty_output><jupyter_text>可视化拟合：<jupyter_code>x_test = np.linspace(0, 25, num=100, endpoint=True)
y_pred = hypothesis(x_test, theta)

plt.plot(x, y, 'bo', label='date')
plt.plot(x_test, y_pred, 'r-', label='pred')
plt.legend()
plt.show()<jupyter_output><empty_output><jupyter_text>(2) $y(x) = ae^{\frac{b}{x}}$取对数：$\ln(y) = \ln(a) + \frac{b}{x}$数据预处理：<jupyter_code># X: [1, 1/x]
X = np.c_[np.ones((len(x), 1)), 1 / x]
print(X)<jupyter_output>[[1.         0.5       ]
 [1.         0.33333333]
 [1.         0.25      ]
 [1.         0.14285714]
 [1.         0.125     ]
 [1.         0.1       ]
 [1.         0.09090909]
 [1.         0.07142857]
 [1.         0.0625    ]
 [1.         0.05555556]
 [1.         0.05263158]]<jupyter_text>新的预测函数：<jupyter_code>def hypothesis(x, theta):
    a, b = theta
    h = np.exp(a) * np.exp(b / x)
    return h<jupyter_output><empty_output><jupyter_text>拟合：<jupyter_code># theta: [ln(a), b]
theta = normalEquationFit(X, np.log(y))
theta<jupyter_output><empty_output><jupyter_text>代价计算：<jupyter_code>cost(x, y, hypothesis, theta)
x_test = np.linspace(0, 25, num=100, endpoint=True)
y_pred = hypothesis(x_test, theta)

plt.plot(x, y, 'bo', label='date')
plt.plot(x_test, y_pred, 'r-', label='pred')
plt.ylim(90, 120)
plt.legend()
plt.show()<jupyter_output><empty_output>