插值-样条插值

插值-样条插值百度百科定义插值:在离散数据的基础上插补连续函数,使得这条连续曲线经过全部离散点,同时也可以估计出函数在其他点的近似值。样条插值:一种以可变样条来作出一条经过一系列点的光滑曲线的数学方法。插值样条是由一些多项式组成的,每一个多项式都是由相邻的两个数据点决定的,这样,任意的两个相邻的多项式以及它们的导数在连接点处都是连续的。样条插值法简单理解,就是每两个点之间确定…

大家好,欢迎来到IT知识分享网。

百度百科定义

插值:在离散数据的基础上插补连续函数,使得这条连续曲线经过全部离散点,同时也可以估计出函数在其他点的近似值。

样条插值:一种以 可变样条 来作出一条经过一系列点的光滑曲线的数学方法。插值样条是由一些多项式组成的,每一个多项式都是由相邻的两个数据点决定的,这样,任意的两个相邻的多项式以及它们的导数在连接点处都是连续的。

 

样条插值法

简单理解,就是每两个点之间确定一个函数,这个函数就是一个样条,函数不同,样条就不同,所以定义中说 可变样条,然后把所有样条分段结合成一个函数,就是最终的插值函数。

 

思路1 – 线性样条

两点确定一条直线,我们可以在每两点间画一条直线,就可以把所有点连起来。

插值-样条插值

显然曲线不够光滑,究其原因是因为连接点处导数不相同。

 

思路2 – 二次样条

直线不行,用曲线代替,二次函数是最简单的曲线。

假设4个点,x0,x1,x2,x3,有3个区间,需要3个二次样条,每个二次样条为  ax^2+bx+c,故总计9个未知数。

插值-样条插值

1. x0,x3两个端点都有一个二次函数经过,可确定2个方程

2. x1,x2两个中间点都有两个二次函数经过,可确定4个方程

3. 中间点处必须连续,需要保证左右二次函数一阶导相等

    2*a1*x1+b1=2*a2*x1+b2

            2*a2*x2+b2=2*a3*x2+b3

可确定2个方程,此时有了8个方程。

4. 这里假设第一方程的二阶导为0,即 a1=0,又是一个方程,共计9个方程。   【见补充】 

联立即可求解。

 

python 实现

# encoding:utf-8
import numpy as np
import matplotlib.pyplot as plt
from pylab import mpl
"""
二次样条实现
"""
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*3)
        data[(i-1)*3]=x[i]*x[i]
        data[(i-1)*3+1]=x[i]
        data[(i-1)*3+2]=1
        data1 =init(sizeOfInterval*3)
        data1[i * 3] = x[i] * x[i]
        data1[i * 3 + 1] = x[i]
        data1[i * 3 + 2] = 1
        temp=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] = 1
    parameter.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] = 1
    temp=data[1:]
    parameter.append(temp)
    #端点函数值相等为n-1个方程。加上前面的方程为3n-1个方程,最后一个方程为a1=0总共为3n个方程
    i=1
    while i < len(x) - 1:
        data = init(sizeOfInterval * 3)
        data[(i - 1) * 3] =2*x[i]
        data[(i - 1) * 3 + 1] =1
        data[i*3]=-2*x[i]
        data[i*3+1]=-1
        temp=data[1:]
        parameter.append(temp)
        i += 1
    return parameter

"""
对一个size大小的元组初始化为0
"""
def init(size):
    j = 0
    data = []
    while j < size:
        data.append(0)
        j += 1
    return data


"""
功能:计算样条函数的系数。
参数:parametes为方程的系数,y为要插值函数的因变量。
返回值:二次插值函数的系数。
"""

def solutionOfEquation(parametes,y):
    sizeOfInterval = len(x) - 1
    result = init(sizeOfInterval*3-1)
    i=1
    while i<sizeOfInterval:
        result[(i-1)*2]=y[i]
        result[(i-1)*2+1]=y[i]
        i+=1
    result[(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=u"拟合曲线", color="black")
        plt.scatter(data_x,data_y, label=u"离散数据",color="red")
        mpl.rcParams['font.sans-serif'] = ['SimHei']
        mpl.rcParams['axes.unicode_minus'] = False
        plt.title(u"二次样条函数")
        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)

可以看到 y 是多段二次函数拼接而成。

 

输出

插值-样条插值

二次样条插值连续光滑,看起来效果还行。

只是前两个点之间是条直线,因为假设a1=0,二次函数变成b1x+c1,显然是直线;

而且最后两个点之间过于陡峭 。

 

思路3 – 三次样条

二次函数最高项系数为0,导致变成直线,那三次函数最高项系数为0,还是曲线,插值效果应该更好。

三次样条思路与二次样条基本相同,

同样假设4个点,x0,x1,x2,x3,有3个区间,需要3个三次样条,每个三次样条为  ax^3+bx^2+cx+d,故总计12个未知数。

1. 内部节点处的函数值应该相等,这里一共是4个方程。

2. 函数的第一个端点和最后一个端点,应该分别在第一个方程和最后一个方程中。这里是2个方程。

3. 两个函数在节点处的一阶导数应该相等。这里是两个方程。

4. 两个函数在节点处的二阶导数应该相等,这里是两个方程。       【见补充】

5. 假设端点处的二阶导数为零,这里是两个方程。          【见补充】

           a1=0 

           b1=0

python 实现

# encoding:utf-8
import numpy as np
import matplotlib.pyplot as plt
from pylab import mpl
"""
三次样条实现
"""
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] = 1
        data1 =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] = 1
        temp = 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] = 1
    parameter.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] = 1
    temp = data[2:]
    parameter.append(temp)
    # 端点函数一阶导数值相等为n-1个方程。加上前面的方程为3n-1个方程。
    i=1
    while 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] = 1
        data[i * 4] = -3 * x[i] * x[i]
        data[i * 4 + 1] = -2 * x[i]
        data[i * 4 + 2] = -1
        temp = data[2:]
        parameter.append(temp)
        i += 1
    # 端点函数二阶导数值相等为n-1个方程。加上前面的方程为4n-2个方程。且端点处的函数值的二阶导数为零,为两个方程。总共为4n个方程。
    i = 1
    while i < len(x) - 1:
        data = init(sizeOfInterval * 4)
        data[(i - 1) * 4] = 6 * x[i]
        data[(i - 1) * 4 + 1] = 2
        data[i * 4] = -6 * x[i]
        data[i * 4 + 1] = -2
        temp = data[2:]
        parameter.append(temp)
        i += 1
    return parameter



"""
对一个size大小的元组初始化为0
"""
def init(size):
    j = 0
    data = []
    while j < size:
        data.append(0)
        j += 1
    return data

"""
功能:计算样条函数的系数。
参数:parametes为方程的系数,y为要插值函数的因变量。
返回值:三次插值函数的系数。
"""
def solutionOfEquation(parametes,y):
    sizeOfInterval = len(x) - 1
    result = init(sizeOfInterval*4-2)
    i=1
    while i<sizeOfInterval:
        result[(i-1)*2]=y[i]
        result[(i-1)*2+1]=y[i]
        i+=1
    result[(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=u"拟合曲线", color="black")
        plt.scatter(data_x,data_y, label=u"离散数据",color="red")
        mpl.rcParams['font.sans-serif'] = ['SimHei']
        mpl.rcParams['axes.unicode_minus'] = False
        plt.title(u"三次样条函数")
        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)

 

输出

插值-样条插值

 

补充

微分连续性

插值-样条插值

s 代表三次样条,s’是一阶导,s”是二阶导

 

端点条件

上面我们对端点处的样条进行了假设,为什么呢?其实端点可以有多种不同的限制,常见有3种。

自由边界 Natural

首尾两端没有受到任何使他们弯曲的力,二次样条就是 s’=0,三次样条就是 s”=0

固定边界 Clamped

首尾两端点的微分值被指定

非节点边界 Not-A-Knot

把端点当做中间点处理,三次函数不做假设,即

插值-样条插值

插值-样条插值

 

不同边界的比较

插值-样条插值

 

 

 

参考资料:

https://blog.csdn.net/flyingleo1981/article/details/53008931  三次样条插值(Cubic Spline Interpolation)及代码实现(C语言)

https://blog.csdn.net/deramer1/article/details/79034201  三次样条插值法

转载于:https://www.cnblogs.com/yanshw/p/11194058.html

免责声明:本站所有文章内容,图片,视频等均是来源于用户投稿和互联网及文摘转载整编而成,不代表本站观点,不承担相关法律责任。其著作权各归其原作者或其出版社所有。如发现本站有涉嫌抄袭侵权/违法违规的内容,侵犯到您的权益,请在线联系站长,一经查实,本站将立刻删除。 本文来自网络,若有侵权,请联系删除,如若转载,请注明出处:https://yundeesoft.com/22192.html

(0)

相关推荐

发表回复

您的邮箱地址不会被公开。 必填项已用 * 标注

关注微信