基于python的三次样条插值原理及代码

1 三次样条插值

1.1 三次样条插值的基本概念

        三次样条插值是通过求解三弯矩方程组(即三次样条方程组的特殊形式)来得出曲线函数组的过程。在实际计算中,还需要引入边界条件来完成计算。样条插值的名称来源于早期工程师制图时使用的细长木条(样条),这些木条被固定在样点上,然后自由弯曲以绘制出曲线。

1.2 三次样条插值的条件

        假设在区间[a, b]上有n+1个样本点,即[a, b]区间被划分成n个小区间[(x_{0},x_{1}),(x_{1}, x_{2}), ...,(x_{n-1}, x_{n})],其中x_{0}=a, x_{n}=b。三次样条插值需要满足以下条件:

  1. 在每个分段区间[x_{i},x_{i+1}]上,插值函数S(x) = S_{i}(x)都是一个三次多项式
  2. 满足插值条件,即S(x_i) = y_i \quad (i = 0, 1, ..., n)
  3. 曲线光滑,即S(x)、S'(x)、S''(x)连续

1.3 三次样条插值的公式推导

1. 3.1 三次样条函数的形式

        在每个小区间[x_{i},x_{i+1}]上,三次样条函数S_i(x)可以表示为:

S_i(x) = a_i + b_i(x - x_i) + c_i(x - x_i)2 + d_i(x - x_i)3

        其中,a_i ,b_i,c_i,d_i是待求的系数。

1.3.2 插值条件

        由于S(x_{i}) = y_{i},我们可以得到n+1个方程:

 S_i(x_i) = a_i = y_i \quad (i = 0, 1, ..., n) 

1.3.3 连续性条件

        除了插值条件外,还需要保证曲线在样本点处的一阶导数和二阶导数连续。即:

  • 一阶导数连续:对于任意i(0 < i < n),有

    S_i'(x_{i+1}) = S_{i+1}'(x_{i+1})

  • 二阶导数连续:对于任意i(0 < i < n),有

    S_i''(x_{i+1}) = S_{i+1}''(x_{i+1})

1.3.4 求解未知数

        每个小区间内有4个未知数(a_i ,b_i,c_i,d_i),共有n个小区间,所以一共有4n个未知数。根据插值条件、一阶导数连续和二阶导数连续条件,我们可以列出以下方程组:

  • 插值条件:提供n+1个方程(其中ai=yi已知)。
  • 一阶导数连续:提供n-1个方程(两个端点不参与此条件)。
  • 二阶导数连续:提供n-1个方程(同样,两个端点不参与此条件)。

        这样,我们总共有4n-2个方程,但还需要两个额外的方程来求解所有未知数。这两个额外的方程由边界条件给出。

1.3.5 边界条件

        常用的边界条件有三种:

  • 自然边界(Natural Spline):指定端点二阶导数为0,即

    S''(x_0) = 0 = S''(x_n)

  • 固定边界(Clamped Spline):指定端点一阶导数,即

    S'(x_0) = A, \quad S'(x_n) = B 

  • 非扭结边界(Not-a-Knot Spline):使两端点的三阶导与这两端点的邻近点的三阶导相等,即

    S'''(x_0) = S'''(x_1), \quad S'''(x_{n-1}) = S'''(x_n) 

1.3.6 求解过程

        以自然边界条件为例,即 S''(x_0) = 0, \quad S''(x_n) = 0,我们需要通过这些条件来求解每个小区间上的二阶导数 m_{i}=S''(x_i)(通常称为样条曲线的“弯矩”或“斜率”)。

        步骤一:设置方程

        在每个小区间[x_{i},x_{i+1}]上,三次样条S_i(x)的二阶导数 S''(x_n) 是一个常数 m_{i}。因此,三次样条S_i(x) 可以表示为:

S_i(x) = \frac{m_i}{6}(x_{i+1} - x)3 + \frac{m_{i+1}}{6}(x - x_i)3 + a_i(x_{i+1} - x) + b_i(x - x_i) + c_i

其中 a_i ,b_i,c_i​ 是需要通过插值条件和其他连续性条件求解的系数。

        步骤二:利用插值条件

        由于 S_i(x_{i}) = y_{i},S_i(x_{i+1}) = y_{i+1},我们可以列出 n+1 个插值方程。但由于 a_{i}=b_{i}由 S_i(x_{i}) = y_{i} 直接得出),我们实际上只需要解出a_{i},b_{i}​(或者更直接地,解出 m_{i})。

        步骤三:利用一阶导数连续性

        在 x=x_{i}+1处,有 S_i{}'(x_{i}) = S_{i+1}{}'(x_{i+1})。对 S_i{}(x_{i}) = S_{i+1}{}(x) 求导并设置相等,我们可以得到一个关于 m_{i}和 m_{i+1}​ 的方程。由于有 n−1 个这样的内部节点,因此我们可以得到 n−1 个这样的方程。

        步骤四:利用二阶导数连续性(自然边界条件)

        在 x=x_{0} 和 x=x_{n} 处,有S''(x_0) = 0, \quad S''(x_n) = 0,即 m_{0}=0 和 m_{n}=0这两个条件提供了额外的两个方程。

        步骤五:解线性方程组

        现在,我们有了 n−1 个一阶导数连续性方程,加上两个自然边界条件方程,共 n+1 个方程,用于求解 n 个未知数 m_{1},m_{2},...m_{n-1}​(注意 m_{0}和 m_{n}​ 已知为0)。这是一个线性方程组,可以通过高斯消元法、LU分解等方法求解。

        步骤六:计算系数 a_i ,b_i,c_i

        一旦我们得到了所有的 mi​,就可以通过插值条件和其他连续性条件回代求解出每个小区间上的 a_i ,b_i,c_i​。具体来说,可以使用以下公式:

  • a_{i}=b_{i}(已知)
  • b_{i},c_{i}可以通过 S_i(x_{i+1}) = y_{i+1}和 S_i{}'(x_{i}) = S_{i+1}{}'(x_{i+1}) 的联立方程求解。

        步骤七:构建完整的三次样条函数

        最后,我们得到了每个小区间上的三次多项式 S_i(x),它们共同构成了完整的三次样条插值函数 S(x)。在任意点 x∈[x_{i},x_{i+1}] 上,S(x) = S_{i}(x)

1.4 总结

        三次样条插值通过求解一系列线性方程组来确定每个小区间上的三次多项式,这些多项式在样本点处满足插值条件,并且具有连续的一阶导数和二阶导数。自然边界条件、固定边界条件或非扭结边界条件等不同的边界条件选择会影响求解过程和最终的插值曲线形状。三次样条插值因其光滑性和易于实现的特性,在数值分析、数据可视化、计算机辅助设计等领域有着广泛的应用。

2 代码

import numpy as np
import matplotlib.pyplot as plt
from scipy.interpolate import CubicSpline# 定义数据点
x = np.array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15])
y1 = np.array([0, 4, 6, 8, 9, 11, 13, 16, 19, 21, 23, 25, 28, 30, 33, 34])
y2 = np.array([0, 1, 2, 4, 5, 6, 7, 8, 10, 12, 14, 16, 20, 25, 29, 34])
y3 = np.array([0, 8, 12, 18, 20, 26, 30, 32, 33, 34, 34, 36, 36, 35, 35, 34])# 创建三次样条插值函数
cs1 = CubicSpline(x, y1)
cs2 = CubicSpline(x, y2)
cs3 = CubicSpline(x, y3)# 生成插值后的x坐标
x_new = np.linspace(0, 15, 100)# 计算插值后的y坐标
y_new1 = cs1(x_new)
y_new2 = cs2(x_new)
y_new3 = cs3(x_new)# 设置颜色选择器和颜色映射
select2 = (0.1, 0.4, 0.5, 0.3, 0.8, 1.0) # 非连续型色组在0-(色组长度-1)之间选择
colors = plt.get_cmap('rainbow')(select2) # 从色组里选择颜色# 绘制原始数据点和插值曲线
# 调整图表大小和位置
plt.figure(figsize=(10, 5))
plt.plot(x, y1, 'o', color=colors[0], label='ori_y1')
plt.plot(x, y2, 'o', color=colors[2], label='ori_y2')
plt.plot(x, y3, 'o', color=colors[4], label='ori_y3')
plt.plot(x_new, y_new1, '-', color=colors[3], label='y1_interp_curve')
plt.plot(x_new, y_new2, '-', color=colors[1], label='y2_interp_curve')
plt.plot(x_new, y_new3, '-', color=colors[5], label='y3_interp_curve')# 设置图例和标题
plt.legend()
plt.xticks(np.arange(0, 16, 1))
plt.title('Cubic Spline Interpolation Example')
plt.xlabel('x')
plt.ylabel('y')# 显示图表
plt.show()

3 插值结果

图3-1 插值结果

本文来自互联网用户投稿,该文观点仅代表作者本人,不代表本站立场。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如若转载,请注明出处:http://www.mzph.cn/pingmian/47472.shtml

如若内容造成侵权/违法违规/事实不符,请联系多彩编程网进行投诉反馈email:809451989@qq.com,一经查实,立即删除!

相关文章

oracle 经营范围 设计

在Oracle数据库中设计经营范围通常涉及创建相关的数据库表来记录和管理经营范围内的数据。以下是一个简单的例子&#xff0c;展示了如何设计一个经营范围表&#xff1a; CREATE TABLE business_units (bu_id NUMBER PRIMARY KEY,bu_name VARCHAR2(100),parent_bu_id NUMBER,CO…

探索Node.js中的node-xlsx:将Excel文件解析为JSON

在Node.js开发中&#xff0c;处理Excel文件是一个常见需求&#xff0c;特别是在需要导入大量数据或生成报表的场景中。node-xlsx 是一个强大的库&#xff0c;它提供了Excel文件的解析和生成功能。本文将深入探讨 node-xlsx 的使用&#xff0c;并通过一个案例演示如何将Excel文件…

算法——双指针(day2)

目录 202.快乐数 题目解析&#xff1a; 算法解析&#xff1a; 代码&#xff1a; 11.盛最多水的容器 题目解析&#xff1a; 算法解析&#xff1a; 代码&#xff1a; 202.快乐数 力扣链接&#xff1a;202.快乐数 题目解析&#xff1a; 本文中最重要的一句话就是重复平方和…

AI自动生成PPT哪个软件好?高效制作PPT优选这4个

7.15初伏的到来&#xff0c;也宣告三伏天的酷热正式拉开序幕~在这个传统的节气里&#xff0c;人们以各种方式避暑纳凉&#xff0c;享受夏日的悠闲时光。 而除了传统的避暑活动&#xff0c;我们还可以用一种新颖的方式记录和分享这份夏日的清凉——那就是通过PPT的方式将这一传…

班迪录屏Bandicam使用详解

Bandicam是一款功能强大的视频录制工具&#xff0c;录制出来的视频体积较小且内容清晰度较高&#xff0c;平时录屏、录游戏都非常合适。可以全屏幕录制&#xff0c;也可以自定义录制区域&#xff0c;还可以在录制时添加自定义的logo&#xff0c;并且有个绘制模式&#xff0c;适…

学习008-01-02 Define the Data Model and Set the Initial Data(定义数据模型并设置初始数据 )

Define the Data Model and Set the Initial Data&#xff08;定义数据模型并设置初始数据 &#xff09; This topic explains how to implement entity classes for your application. It also describes the basics of automatic user interface construction based on a da…

cookies和session的区别【面试】

一、 共同点&#xff1a; 目的&#xff1a;Cookie和Session都是用来跟踪浏览器用户身份的会话方式。 二、 工作原理&#xff1a; 1. Cookie的工作原理 浏览器端第一次发送请求到服务器端。服务器端创建Cookie&#xff0c;包含用户信息&#xff0c;然后将Cookie发送到浏览器…

基于AT89C51单片机的多功能自行车测速计程器(含文档、源码与proteus仿真,以及系统详细介绍)

本篇文章论述的是基于AT89C51单片机的多功能自行车测速计程器的详情介绍&#xff0c;如果对您有帮助的话&#xff0c;还请关注一下哦&#xff0c;如果有资源方面的需要可以联系我。 目录 选题背景 原理图 PCB图 仿真图 代码 系统论文 资源下载 选题背景 美丽的夜晚&…

力扣刷题(自用)

哈希 128.最长连续序列 128. 最长连续序列 - 力扣&#xff08;LeetCode&#xff09; 这个题要求O(n)的时间复杂度&#xff0c;我一开始想的是双指针算法&#xff08;因为我并不是很熟悉set容器的使用&#xff09;&#xff0c;但是双指针算法有小部分数据过不了。 题解给的哈…

JavaScript:移除元素

这是原题&#xff1a;给你一个数组 nums 和一个值 val&#xff0c;你需要 原地 移除所有数值等于 val 的元素。元素的顺序可能发生改变。然后返回 nums 中与 val 不同的元素的数量。 假设 nums 中不等于 val 的元素数量为 k&#xff0c;要通过此题&#xff0c;您需要执行以下操…

【PyTorch][chapter 26][李宏毅深度学习][attention-2]

前言&#xff1a; Multi-Head Attention 主要作用&#xff1a;将Q,K,V向量分成多个头&#xff0c;形成多个子语义空间&#xff0c;可以让模型去关注不同维度语义空间的信息 目录&#xff1a; attention 机制 Multi-Head Attention 一 attention 注意力 Self-Attention&#x…

三分钟速通银行家算法

银行家算法&#xff08;Bankers Algorithm&#xff09;是一种用于避免死锁的经典算法&#xff0c;广泛应用于操作系统、数据库管理系统及分布式系统中。下面将结合实战场景&#xff0c;详细介绍银行家算法的实现过程和应用。 一、算法背景 银行家算法由荷兰计算机科学家Edsge…

什么是im即时通讯?WorkPlus im即时通讯私有化部署安全可控

IM即时通讯是Instant Messaging的缩写&#xff0c;指的是一种实时的、即时的电子信息交流方式&#xff0c;也被称为即时通讯。它通过互联网和移动通信网络&#xff0c;使用户能够及时交换文本消息、语音通话、视频通话、文件共享等信息。而WorkPlus im即时通讯私有化部署则提供…

防火墙--双机热备

目录 双击热备作用 防火墙和路由器备份不同之处 如何连线 双机 热备 冷备 VRRP VGMP&#xff08;华为私有协议&#xff09; 场景解释 VGMP作用过程 主备的形成场景 接口故障的切换场景 整机故障 原主设备故障恢复的场景 如果没有开启抢占 如果开启了抢占 负载分…

对红酒品质进行数据分析(python)

http://t.csdnimg.cn/UWg2S 数据来源于这篇博客&#xff0c;直接下载好csv文件。 这篇内容均在VScode的jupyter notebook上完成&#xff0c;操作可以看我的另一篇博客&#xff1a;http://t.csdnimg.cn/69sDJ 一、准备工作 1. 导入数据库 #功能是可以内嵌绘图&#xff0c;并…

如何查看Linux中某个项目是否在Docker中运行

方法一&#xff1a;检查进程的 cgroup Docker 容器的进程运行在特定的 cgroup 中。你可以通过检查进程的 cgroup 信息来判断它是否在 Docker 容器中运行。 找到项目的进程 ID (PID)&#xff1a; 假设你知道项目的进程名称&#xff0c;例如 my_project&#xff0c;你可以使用 p…

纯硬件一键开关机电路的工作原理

这是一个一键开关机电路: 当按一下按键然后松开&#xff0c;MOS管导通&#xff0c;VOUT等于电源电压; 当再次按一下按键然后松开&#xff0c;MOS管关闭&#xff0c;VOUT等于0; 下面来分析一下这个电路的工作原理。上电后&#xff0c;输入电压通过R1和R2给电容充电&#xff0c;最…

继承和多态常见的面试问题

文章目录 概念问答 概念 下面哪种面向对象的方法可以让你变得富有( A) A: 继承 B: 封装 C: 多态 D: 抽象 (D )是面向对象程序设计语言中的一种机制。这种机制实现了方法的定义与具体的对象无关&#xff0c; 而对方法的调用则可以关联于具体的对象。 A: 继承 B: 模板 C: 对象的…

Java的Cpp本地库调用

叠甲&#xff1a;以下文章主要是依靠我的实际编码学习中总结出来的经验之谈&#xff0c;求逻辑自洽&#xff0c;不能百分百保证正确&#xff0c;有错误、未定义、不合适的内容请尽情指出&#xff01; 文章目录 1.设置项目文件结构2.生成 Cpp 头文件3.编写 Cpp 程序实现4.生成本…

TS相较于JS有什么优缺点

TypeScript&#xff08;TS&#xff09;是JavaScript的一个超集&#xff0c;它添加了静态类型检查和编译时的强大功能&#xff0c;目的是提高代码质量和维护性。相较于JavaScript&#xff0c;TS的主要优点和缺点如下&#xff1a; 优点&#xff1a; 类型安全性&#xff1a;通过…