数值计算 利用Python实现三次样条插值(不调用函数库)

三次样条插值介绍

由于之前写的被误删了，而且公式很多很难打，这里就不赘述了，详见《数值分析》李庆扬(第五版)

代码实现

这里仅给出第一类边界条件的代码，其余两种边界情况可以自行修改。

import numpy as np
import matplotlib.pyplot as plt
from pylab import mpl


def SI_dd(x: list, f: list, n: int):  # n-节点数   只要二阶均差
    dd = [0 for _ in range(0, n - 1)]  # divided difference

    for i in range(0, 2):
        for j in range(0, n - 1 - i):  # 0, 1, 2, 3
            if i == 0:
                dd[j] = (f[j + 1] - f[j]) / (x[j + 1] - x[j])
            else:
                dd[j] = (dd[j + 1] - dd[j]) / (x[j + 2] - x[j])
    # print(dd)
    return dd


def solve_m(x: list, f: list, n: int, f_: list):
    h = [x[_ + 1] - x[_] for _ in range(0, n - 1)]
    hh = [(f[_ + 1] - f[_]) / h[_] for _ in range(0, n - 1)]
    # print(hh)
    p = np.array(np.eye(n, n)) * 2
    d = np.array(SI_dd(x, f, n)[0: -1])
    d_0 = (hh[0] - f_[0]) / h[0]
    d_n = (f_[1] - hh[n - 2]) / h[n - 2]
    d = np.insert(d, 0, d_0)
    d = np.append(d, d_n)
    # print(d_0, d_n)

    lmd = [h[_ + 1] / (h[_] + h[_ + 1]) for _ in range(0, n - 2)]
    lmd.insert(0, 1)

    miu = [h[_] / (h[_] + h[_ + 1]) for _ in range(0, n - 2)]
    miu.append(1)
    for i in range(0, n - 1):
        p[i][i + 1] = lmd[i]
    for j in range(1, n):
        p[j][j - 1] = miu[j - 1]
    m = np.dot(np.linalg.inv(p), 6 * d)

    return m


def get_function(data: np.ndarray, start: float, end: float, i: int, x: list, f: list, n: int, f_: list):
    h = [x[_ + 1] - x[_] for _ in range(0, n - 1)]
    m = solve_m(x, f, n, f_)

    res = m[i] / (6 * h[i]) * pow(end - data, 3) + m[i + 1] / (6 * h[i]) * pow(data - start, 3) + \
          (f[i + 1] - m[i + 1] / 6 * pow(h[i], 2)) * (data - start) / h[i] + \
          (f[i] - m[i] / 6 * pow(h[i], 2)) * (end - data) / h[i]
    return res


def draw(x: list, f: list, n: int, f_: list):
    h = [x[_ + 1] - x[_] for _ in range(0, n - 1)]

    for i in range(0, n - 1):
        if i != n - 2:
            plt.plot(np.arange(x[i], x[i + 1] + 0.001, 0.001),
                     get_function(np.arange(x[i], x[i + 1] + 0.001, 0.001), x[i], x[i] + h[i], i, x, f, n, f_),
                     color="blue")
        else:
            plt.plot(np.arange(x[i], x[i + 1] + 0.001, 0.001),
                     get_function(np.arange(x[i], x[i + 1] + 0.001, 0.001), x[i], x[i] + h[i], i, x, f, n, f_),
                     label="三次样条拟合曲线",
                     color="blue")

    plt.scatter(x, f, label="节点数据", color="red")

    plt.title("三次样条插值结果")
    mpl.rcParams['font.sans-serif'] = ['SimHei']
    mpl.rcParams['axes.unicode_minus'] = False
    plt.legend(loc="upper right")

使用实例

考虑龙格函数：$$f(x)=\frac{1}{1+25x^2}$$
在$[-1,1]$上取11个等距节点，对这些节点进行三次样条插值。
并使用第一类边界条件:$f^{\prime }(-1)=50/(1+25{)}^2=-f^{\prime }(1)$

import SI_first
import numpy as np
import matplotlib.pyplot as plt
import newton as nt


def draw(x: list, f: list, n: int, f_: list):
    xx = list(np.arange(-1, 1.01, 0.01))
    ff = [1 / (1 + 25 * pow(num, 2)) for num in xx]
    plt.plot(xx, ff, label="龙格函数", color="black")

    SI_first.draw(x, f, n, f_)
    plt.show()


if __name__ == "__main__":
    n = 11
    x = list(np.arange(-1, 1 + 2 / (n - 1), 2 / (n - 1)))
    f = [1 / (1 + 25 * pow(num, 2)) for num in x]
    f_ = [50 / pow(1 + 25, 2), -50 / pow(1 + 25, 2)]
    draw(x, f, n, f_)

得到图像：