三次样条插值详解(附代码实现)
文章目录
- 前言
- 引入
- 二次样条的原理
- 二次样条代码实现
- 三次样条的原理
- 三次样条代码实现
前言
当已知某些点而不知道具体方程时候,最经常遇到的场景就是做实验,采集到数据的时候,我们通常有两种做法:拟合或者插值。拟合不要求方程通过所有的已知点,讲究神似,就是整体趋势一致。插值则是形似,每个已知点都必会穿过,但是高阶会出现龙格库塔现象,所以一般采用分段插值。今天我们就来说说这个分段三次样条插值。
引入
首先我们先抛开众多的回归算法不谈, 我们对于给出如下的离散的数据点,现在想根据如下的数据点来推测 x=6 时的值,我们应该采用什么方法呢?
用于拟合样条函数的数据
x | f(x) |
---|---|
3.0 | 2.5 |
4.5 | 1.0 |
7.0 | 2.5 |
9.0 | 0.5 |
我们知道在平面上两个点确定一条直线,三个点确定一条抛物线(假设曲线的类型是抛物线),那么现在有四个点,我们很自然的会想到,既然两个点确定一条直线,那么最简单的方法就是,两个点之间连一条线,两个点之间连一条线,最后得到的一种折线图如下:这样我们只要确定 x=6 时的直线,把自变量的值带进去,就显然会得到预测的函数值。但是,这种方法显然具有很明显的缺陷:曲线不够光滑,连接点处的斜率变化过大。可能会导致函数的一阶导数不连续。那么我们应该如何解决这个问题呢?
二次样条的原理
我们会想到既然直线不行,那么我们就用曲线来近似的代替和描述。最简单的曲线是二次函数,如果我们用二次函数:ax2+bx+cax^2 + bx + cax2+bx+c 来描述曲线,最后的结果可能会好一点,下图中一共有4个点,可以分成3个区间。每一个区间都需要一个二次函数来描述,一共需要9个未知数。下面的任务就是找出9个方程。
如下所示:一共有x0,x1,x2,x3x_0,x_1,x_2,x_3x0,x1,x2,x3四个点,三个区间,每个区间上都有一个方程。
曲线方程在节点处的值必须相等,即函数在x1,x2两个点处的值必须符合两个方程,这里一共是4个方程:
a1x12+b1x1+c1=f(x1)a_1x_1^2 + b_1x_1 + c_1 = f(x_1)a1x12+b1x1+c1=f(x1)
a2x12+b2x1+c2=f(x1)a_2x_1^2 + b_2x_1 + c_2 = f(x_1)a2x12+b2x1+c2=f(x1)
a2x22+b2x2+c2=f(x2)a_2x_2^2 + b_2x_2 + c_2 = f(x_2)a2x22+b2x2+c2=f(x2)
a3x22+b3x2+c3=f(x2)a_3x_2^2 + b_3x_2 + c3_ =f (x_2)a3x22+b3x2+c3=f(x2)
第一个端点和最后一个端点必须过第一个和最后一个方程:这里一共是2个方程
节点处的一阶导数的值必须相等。这里为两个方程。
2a1x1+b1=2a2x1+b22a_1x_1 + b_1 = 2a_2x_1 + b_22a1x1+b1=2a2x1+b2
2a2x2+b2=2a3x2+b32a_2x_2 + b_2 = 2a_3x_2 + b_32a2x2+b2=2a3x2+b3
在这里假设第一个方程的二阶导数为0:这里为一个方程:
a1=0a_1 = 0a1=0
上面是对应的9个方程,现在只要把九个方程联立求解,最后就可以实现预测 x=6 处节点的值。
下面是写成矩阵的形式,由于a1=0,所以未知数的个数少了一个
二次样条代码实现
下面是二次样条的python实现,最后将结果绘制在图上。
import numpy as np
import matplotlib.pyplot as plt
from pylab import mpl
"""
二次样条实现:
函数x的自变量为:3, 4.5, 7, 9因变量为:2.5, 1 2.5, 0.5
"""
x = [3, 4.5, 7, 9]
y = [2.5, 1, 2.5, 0.5]"""一共有三个区间,用二次样条求解,需要有9个方程""""""
功能:完后对二次样条函数求解方程参数的输入
参数:要进行二次样条曲线计算的自变量
返回值:方程的参数
"""
def calculateEquationParameters(x):#parameter为二维数组,用来存放参数,sizeOfInterval是用来存放区间的个数parameter = []sizeOfInterval=len(x)-1;i = 1#首先输入方程两边相邻节点处函数值相等的方程为2n-2个方程while i < len(x)-1:data = init(sizeOfInterval*3)data[(i-1)*3]=x[i]*x[i]data[(i-1)*3+1]=x[i]data[(i-1)*3+2]=1data1 =init(sizeOfInterval*3)data1[i * 3] = x[i] * x[i]data1[i * 3 + 1] = x[i]data1[i * 3 + 2] = 1temp=data[1:]parameter.append(temp)temp=data1[1:]parameter.append(temp)i += 1#输入端点处的函数值。为两个方程,加上前面的2n-2个方程,一共2n个方程data = init(sizeOfInterval*3-1)data[0] = x[0]data[1] = 1parameter.append(data)data = init(sizeOfInterval *3)data[(sizeOfInterval-1)*3+0] = x[-1] * x[-1]data[(sizeOfInterval-1)*3+1] = x[-1]data[(sizeOfInterval-1)*3+2] = 1temp=data[1:]parameter.append(temp)#端点函数值相等为n-1个方程。加上前面的方程为3n-1个方程,最后一个方程为a1=0总共为3n个方程i=1while i < len(x) - 1:data = init(sizeOfInterval * 3)data[(i - 1) * 3] =2*x[i]data[(i - 1) * 3 + 1] =1data[i*3]=-2*x[i]data[i*3+1]=-1temp=data[1:]parameter.append(temp)i += 1return parameter"""
对一个size大小的元组初始化为0
"""
def init(size):j = 0;data = []while j < size:data.append(0)j += 1return data"""
功能:计算样条函数的系数。
参数:parametes为方程的系数,y为要插值函数的因变量。
返回值:二次插值函数的系数。
"""def solutionOfEquation(parametes,y):sizeOfInterval = len(x) - 1;result = init(sizeOfInterval*3-1)i=1while i<sizeOfInterval:result[(i-1)*2]=y[i]result[(i-1)*2+1]=y[i]i+=1result[(sizeOfInterval-1)*2]=y[0]result[(sizeOfInterval-1)*2+1]=y[-1]a = np.array(calculateEquationParameters(x))b = np.array(result)return np.linalg.solve(a,b)"""
功能:根据所给参数,计算二次函数的函数值:
参数:parameters为二次函数的系数,x为自变量
返回值:为函数的因变量
"""
def calculate(paremeters,x):result=[]for data_x in x:result.append(paremeters[0]*data_x*data_x+paremeters[1]*data_x+paremeters[2])return result"""
功能:将函数绘制成图像
参数:data_x,data_y为离散的点.new_data_x,new_data_y为由拉格朗日插值函数计算的值。x为函数的预测值。
返回值:空
"""
def Draw(data_x,data_y,new_data_x,new_data_y):plt.plot(new_data_x, new_data_y, label="拟合曲线", color="black")plt.scatter(data_x,data_y, label="离散数据",color="red")mpl.rcParams['font.sans-serif'] = ['SimHei']mpl.rcParams['axes.unicode_minus'] = Falseplt.title("二次样条函数")plt.legend(loc="upper left")plt.show()result=solutionOfEquation(calculateEquationParameters(x),y)
new_data_x1=np.arange(3, 4.5, 0.1)
new_data_y1=calculate([0,result[0],result[1]],new_data_x1)
new_data_x2=np.arange(4.5, 7, 0.1)
new_data_y2=calculate([result[2],result[3],result[4]],new_data_x2)
new_data_x3=np.arange(7, 9.5, 0.1)
new_data_y3=calculate([result[5],result[6],result[7]],new_data_x3)
new_data_x=[]
new_data_y=[]
new_data_x.extend(new_data_x1)
new_data_x.extend(new_data_x2)
new_data_x.extend(new_data_x3)
new_data_y.extend(new_data_y1)
new_data_y.extend(new_data_y2)
new_data_y.extend(new_data_y3)
Draw(x,y,new_data_x,new_data_y)
二次样条函数运行之后的结果如下,从图像中,我们可以看出,二次样条在函数的连接处的曲线是光滑的。这时候,我们将x=5输入到函对应的函数端中,就可以预测相应的函数值。但是,这里还有一个问题,就是二次样条函数的前两个点是直线,而且函数的最后一个区间内看起来函数凸出很高。我们还想解决这些问题,这时候,我们想是否可以用三次样条函数来进行函数的模拟呢?
三次样条的原理
三次样条的原理和二次样条的原理相同,我们用函数 ax3+bx2+cx+dax^3 + bx^2 + cx + dax3+bx2+cx+d 这个函数来进行操作,这里一共是4个点,分为3个区间,每个区间一个三次样条函数的话,一共是12个方程,只要我们找出这12个方程,这个问题就算解决了。
内部节点处的函数值应该相等,这里一共是4个方程。
函数的第一个端点和最后一个端点,应该分别在第一个方程和最后一个方程中。这里是2个方程。
两个函数在节点处的一阶导数应该相等。这里是两个方程。
两个函数在节点处的二阶导数应该相等,这里是两个方程。
端点处的二阶导数为零,这里是两个方程。
a1=0a_1 = 0a1=0
b1=0b_1 = 0b1=0
三次样条代码实现
import numpy as np
import matplotlib.pyplot as plt
from pylab import mpl
"""
三次样条实现:
函数x的自变量为:3, 4.5, 7, 9因变量为:2.5, 1 2.5, 0.5
"""
x = [3, 4.5, 7, 9]
y = [2.5, 1, 2.5, 0.5]"""
功能:完后对三次样条函数求解方程参数的输入
参数:要进行三次样条曲线计算的自变量
返回值:方程的参数
"""
def calculateEquationParameters(x):#parameter为二维数组,用来存放参数,sizeOfInterval是用来存放区间的个数parameter = []sizeOfInterval=len(x)-1;i = 1#首先输入方程两边相邻节点处函数值相等的方程为2n-2个方程while i < len(x)-1:data = init(sizeOfInterval*4)data[(i-1)*4] = x[i]*x[i]*x[i]data[(i-1)*4+1] = x[i]*x[i]data[(i-1)*4+2] = x[i]data[(i-1)*4+3] = 1data1 =init(sizeOfInterval*4)data1[i*4] =x[i]*x[i]*x[i]data1[i*4+1] =x[i]*x[i]data1[i*4+2] =x[i]data1[i*4+3] = 1temp = data[2:]parameter.append(temp)temp = data1[2:]parameter.append(temp)i += 1# 输入端点处的函数值。为两个方程, 加上前面的2n - 2个方程,一共2n个方程data = init(sizeOfInterval * 4 - 2)data[0] = x[0]data[1] = 1parameter.append(data)data = init(sizeOfInterval * 4)data[(sizeOfInterval - 1) * 4 ] = x[-1] * x[-1] * x[-1]data[(sizeOfInterval - 1) * 4 + 1] = x[-1] * x[-1]data[(sizeOfInterval - 1) * 4 + 2] = x[-1]data[(sizeOfInterval - 1) * 4 + 3] = 1temp = data[2:]parameter.append(temp)# 端点函数一阶导数值相等为n-1个方程。加上前面的方程为3n-1个方程。i=1while i < sizeOfInterval:data = init(sizeOfInterval * 4)data[(i - 1) * 4] = 3 * x[i] * x[i]data[(i - 1) * 4 + 1] = 2 * x[i]data[(i - 1) * 4 + 2] = 1data[i * 4] = -3 * x[i] * x[i]data[i * 4 + 1] = -2 * x[i]data[i * 4 + 2] = -1temp = data[2:]parameter.append(temp)i += 1# 端点函数二阶导数值相等为n-1个方程。加上前面的方程为4n-2个方程。且端点处的函数值的二阶导数为零,为两个方程。总共为4n个方程。i = 1while i < len(x) - 1:data = init(sizeOfInterval * 4)data[(i - 1) * 4] = 6 * x[i]data[(i - 1) * 4 + 1] = 2data[i * 4] = -6 * x[i]data[i * 4 + 1] = -2temp = data[2:]parameter.append(temp)i += 1return parameter"""
对一个size大小的元组初始化为0
"""
def init(size):j = 0;data = []while j < size:data.append(0)j += 1return data"""
功能:计算样条函数的系数。
参数:parametes为方程的系数,y为要插值函数的因变量。
返回值:三次插值函数的系数。
"""def solutionOfEquation(parametes,y):sizeOfInterval = len(x) - 1;result = init(sizeOfInterval*4-2)i=1while i<sizeOfInterval:result[(i-1)*2]=y[i]result[(i-1)*2+1]=y[i]i+=1result[(sizeOfInterval-1)*2]=y[0]result[(sizeOfInterval-1)*2+1]=y[-1]a = np.array(calculateEquationParameters(x))b = np.array(result)for data_x in b:print(data_x)return np.linalg.solve(a,b)"""
功能:根据所给参数,计算三次函数的函数值:
参数:parameters为二次函数的系数,x为自变量
返回值:为函数的因变量
"""
def calculate(paremeters,x):result=[]for data_x in x:result.append(paremeters[0]*data_x*data_x*data_x+paremeters[1]*data_x*data_x+paremeters[2]*data_x+paremeters[3])return result"""
功能:将函数绘制成图像
参数:data_x,data_y为离散的点.new_data_x,new_data_y为由拉格朗日插值函数计算的值。x为函数的预测值。
返回值:空
"""
def Draw(data_x,data_y,new_data_x,new_data_y):plt.plot(new_data_x, new_data_y, label="拟合曲线", color="black")plt.scatter(data_x,data_y, label="离散数据",color="red")mpl.rcParams['font.sans-serif'] = ['SimHei']mpl.rcParams['axes.unicode_minus'] = Falseplt.title("三次样条函数")plt.legend(loc="upper left")plt.show()result=solutionOfEquation(calculateEquationParameters(x),y)
new_data_x1=np.arange(3, 4.5, 0.1)
new_data_y1=calculate([0,0,result[0],result[1]],new_data_x1)
new_data_x2=np.arange(4.5, 7, 0.1)
new_data_y2=calculate([result[2],result[3],result[4],result[5]],new_data_x2)
new_data_x3=np.arange(7, 9.5, 0.1)
new_data_y3=calculate([result[6],result[7],result[8],result[9]],new_data_x3)
new_data_x=[]
new_data_y=[]
new_data_x.extend(new_data_x1)
new_data_x.extend(new_data_x2)
new_data_x.extend(new_data_x3)
new_data_y.extend(new_data_y1)
new_data_y.extend(new_data_y2)
new_data_y.extend(new_data_y3)
Draw(x,y,new_data_x,new_data_y)
参考自:https://zhuanlan.zhihu.com/p/62860859
三次样条插值详解(附代码实现)相关推荐
- 动态规划---01背包问题--Dp(详解附代码)
一.动态规划 代表一类问题(最优子结构或子问题最优性)的一般解法,是设计方法或者策略,不是具体算法 本质:递推,核心是找到状态转移的方式,写出dp方程. 解决问题:交叉,重叠子问题(最优子问题) 形式 ...
- 【排序】堆排序详解 附代码
按照国际惯例,开篇前先简单介绍(吹一波)堆排序(Heapsort).Heapsort是一种优秀的排序算法(个人感觉基本排序算法中仅次于快速排序),时间复杂度为O(nlgn),同时,Heapsort具有 ...
- 各种进制转换(二,八,十,十六进制间转换)详解附代码
进制转换 原理 进制转换是人们利用符号来计数的方法.进制转换由一组数码符号和两个基本因素"基数"与"位权"构成. 基数是指,进位计数制中所采用的数码(数制中用来 ...
- 前序遍历、中序遍历、后序遍历层序遍历详解附代码(数据结构C语言)
目录 (1)前序遍历 (DLR) 递归算法 (2)中序遍历 (LDR) 递归算法 (3)后序遍历 (LRD) 递归算法 (4)层序遍历 队列实现方法 层序遍历的定义: 实现方法: 代码实现 结果截图 ...
- Numpy学习笔记(二):argmax参数中axis=0,axis=1,axis=-1详解附代码
文章目录 1.argmax和max函数区别 2.axis=0/axis=1/axis=-1的区别 3.具体代码分析 ---3.1一维数组 ---3.2二维数组 ---3.3三维数组 1.argmax和 ...
- 随机分布嵌入(RDE)框架详解附代码
介绍 研究了好一阵子马欢飞老师在PNAS上发的文章,下面附上个人的研究心得与代码与大家讨论. 在基于非线性系统的理论基础上,延迟嵌入理论以及广义嵌入理论等相空间重构的理论基础上,观察者便有可能从一个观 ...
- c++实现贪吃蛇详解(附代码)
文章目录 前言 一.运行界面 二.类的大致抽象 三.关于一些问题的思考 四.最后一些想说的 五.代码 前言 经过一个多月的学习,又加深了对c++的理解,所以接下来,就和大家分享一下,一个月学习c++的 ...
- 【大道至简】机器学习算法之EM算法(Expectation Maximization Algorithm)详解(附代码)---通俗理解EM算法。
☕️ 本文来自专栏:大道至简之机器学习系列专栏
- 目标检测模型的评估指标mAP详解(附代码)
https://zhuanlan.zhihu.com/p/37910324 对于使用机器学习解决的大多数常见问题,通常有多种可用的模型.每个模型都有自己的独特之处,并随因素变化而表现不同. 每个模型在 ...
- 数学建模——一维、二维插值模型详解Python代码
数学建模--一维.二维插值模型详解Python代码 一.一维插值 # -*-coding:utf-8 -*- import numpy as np from scipy import interpol ...
最新文章
- FreeSWITCH IVR中lua调用并执行nodejs代码
- 共享可写节包含重定位_PE结构学习01-DOS头-NT头-节表头
- Python中的sorted函数以及operator.itemgetter函数
- 一文读懂spring boot 和微服务的关系
- VS2012 中 c++项目中的各个选项介绍
- 《通过C#学Proto.Actor模型》之 HelloWorld
- 编程语言的“别样”编年史
- 2018-2019-1 《信息安全系统设计基础》教学进程
- 人脸解锁除了要穿衣服,还有什么秘密?
- 【Julia】Julia v1.5.1 更改Pkg存放位置
- Git、Github、Gitlab、Gitee、Git-ce的区别
- 龙芯CPU芯片介绍说明
- postgres关键字、常量和数据类型
- 彼岸花的传说——彼岸繁花,开一千年,落一千年,花叶不相见。情不为因果,缘注定生死。...
- APIS(BOM)——Window对象、本地存储
- CSS控制,彩色图片变灰色
- 从沟通的一般模型想到互联网,再想到数字媒体,最后想到信息世界
- Rsa前后端加密交互 带demo
- 本科毕业设计(论文)开题报告模板1
- 我发的文章变成了0和1,那0和1是怎么发送给你的?计算机网络(二)物理层
热门文章
- 为什么手机浏览器打不开html文件,win7浏览器打不开本地html文件的原因及解决方法...
- asp网站如何设置默认页_IIS 7.5 在 Windows Server(R) 2008 R2
- XP下grub4dos硬盘安装和启动FreeBSD-8.0-i386+GNOME桌面
- Asp.Net MVC4.0 官方教程 入门指南之六--查看Edit方法和Edit视图
- 计算机用户全部删除,电脑用户怎么删除:批量删除计算机用户方法
- 文曲星猜数字用c语言编程,文曲星中的猜数字游戏,要猜一个四位数,有什么通用公式?...
- 一个屌丝程序猿的人生(四十三)
- debezium系列之:理解database.server.name和database.history.kafka.topic
- 解决The kernel appears to have died. It will restart automatically问题
- 工业蒸汽_到底什么是蒸汽机,我想要一个吗?