插值-樣條插值

本文轉載自查看原文 2019-07-16 12:04 5080 數學與算法

百度百科定義

插值：在離散數據的基礎上插補連續函數，使得這條連續曲線經過全部離散點，同時也可以估計出函數在其他點的近似值。

樣條插值：一種以可變樣條來作出一條經過一系列點的光滑曲線的數學方法。插值樣條是由一些多項式組成的，每一個多項式都是由相鄰的兩個數據點決定的，這樣，任意的兩個相鄰的多項式以及它們的導數在連接點處都是連續的。

樣條插值法

簡單理解，就是每兩個點之間確定一個函數，這個函數就是一個樣條，函數不同，樣條就不同，所以定義中說可變樣條，然后把所有樣條分段結合成一個函數，就是最終的插值函數。

思路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)