三次(自然)样条插值及其Python代码实现

本文介绍了如何使用三次样条插值方法对给定点进行插值,并通过Python的numpy和matplotlib库实现。通过计算系数矩阵和常数矩阵,求解插值函数的参数,最后绘制了插值曲线。

摘要生成于 C知道 ,由 DeepSeek-R1 满血版支持, 前往体验 >

前言:

利用一个未知理论推出的结果来证明该未知理论是一种荒谬的行为……

理论:

参考《数值分析》这本书(为什么不参考课上PPT,是不想吗QAQ)

代码:

(强烈建议先看一看理论部分的讲解,尝试自己写一写,其实并不难的)

# 三次样条插值
import numpy as np
from matplotlib import pyplot as plt

''' 
定义插值函数为: S_i(x_i) = y_i + b_i(x - x_i) + c_i(x - x_i)^2 + d_i(x - x_i)^3
'''

#  输入需要插值的点
ls = [(-1.00, -14.58), (0.00, 0.00), (1.27, 0.00), (2.55, 0.00), (3.82, 0.00), (4.92, 0.88), (5.02, 11.17)]
x_1 = [tuple_[0] for tuple_ in ls]  # 插值点的横坐标
y_1 = [tuple_[1] for tuple_ in ls]  # 插值点的纵坐标
delta = np.zeros(len(ls) - 1)
Delta = np.zeros(len(ls) - 1)

for i in range(len(ls) - 1):
    delta[i] = ls[i + 1][0] - ls[i][0]
    Delta[i] = ls[i + 1][1] - ls[i][1]

A = np.zeros([len(ls), len(ls)])
B = np.zeros(len(ls))

# n个方程求解n个未知变量c_i

# 列出系数矩阵A

for i in range(len(ls)):
    if i == 0:
        A[i][0] = 1
        continue
    if i == len(ls) - 1:
        A[i][i] = 1
        continue
    A[i][i - 1] = delta[i - 1]
    A[i][i] = 2 * (delta[i - 1] + delta[i])
    A[i][i + 1] = delta[i]

# 常数矩阵B

for i in range(len(ls)):
    if i == 0:
        B[i] = 0
        continue
    if i == len(ls) - 1:
        B[i] = 0
        continue
    B[i] = 3 * (Delta[i] / delta[i] - Delta[i - 1] / delta[i - 1])

# 求解线性方程组,得到系数c_i
c = np.linalg.solve(A, B)

b = np.zeros(len(ls))
d = np.zeros(len(ls))

# 求解b_i,d_i
for i in range(len(ls) - 1):
    b[i] = Delta[i] / delta[i] - delta[i] / 3 * (2 * c[i] + c[i + 1])
    d[i] = (c[i + 1] - c[i]) / (3 * delta[i])

# 绘制散点图
plt.scatter(x_1, y_1, color='blue')
for i in range(len(ls) - 1):
    x = np.linspace(x_1[i], x_1[i + 1], 100)
    y = y_1[i] + b[i] * (x - x_1[i]) + c[i] * (x - x_1[i]) ** 2 + d[i] * (x - x_1[i]) ** 3
    plt.plot(x, y, color='green')
plt.show()

结语:

好好学习,天天向上!

### Python 实现三次样条插值算法 #### 示例代码与解释 为了实现三次样条插值,在Python中通常会利用`scipy.interpolate`库中的`splrep`和`splev`函数来构建并评估样条曲线。下面展示了一个完整的例子,包括数据准备、计算以及绘图部分。 ```python import numpy as np from scipy.interpolate import splrep, splev import matplotlib.pyplot as plt # 定义原始的数据点集 x = np.linspace(0, 10, 11) # 创建自变量序列 y = np.sin(x) # 对应的因变量取正弦值作为示例 # 计算B-spline表示形式所需参数tck=(t,c,k),其中k=3代表三阶多项式即立方样条 tck = splrep(x, y, k=3) # 插入更多点用于平滑绘制图形 xi = np.linspace(min(x), max(x), 100) yi = splev(xi, tck) # 绘制原离散点及其对应的光滑样条拟合曲线 plt.figure(figsize=(8, 6)) plt.plot(x, y, 'o', label='Original Data Points') plt.plot(xi, yi, '-', lw=2, label='Cubic Spline Interpolation') plt.legend() plt.title('Cubic Spline Interpolation Example') plt.show() ``` 这段程序首先定义了一组简单的输入输出对(这里用的是正弦波),接着调用了`splrep()`函数得到描述这条特定路径所需的控制顶点(tension points)、系数(coefficient)及阶次(degree)[^1]的信息;最后通过`splev()`函数基于这些信息生成新的预测值,并将其可视化出来以便直观理解效果如何。 对于更复杂的场景下可能还需要考虑边界条件等问题,则可以根据具体需求调整上述过程中的设置项或是直接编写更加底层的操作如手动建立线性方程组求解节点处的一阶导数值等。
评论 2
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值