（十三）二次樣條插值

 

#encoding=utf-8
from numpy import *
import  numpy as  np

import matplotlib.pyplot as plt
plt.close()
fig=plt.figure()
plt.grid(True)
plt.xlabel("X")
plt.ylabel("Y")
plt.title(u"二次樣條", fontproperties='SimHei')
#二次樣條
#設已知n+1個點 設方程式為yi = ai*x**2+bi*x+ci,則需要求解3n個參數，則需要3n個條件。
# 由曲線過點得2n 個條件 ，由內點處一次求導相等得n-1個條件 ，令第一個節點處二階導為0 得一個條件 a1 = 0


# #n+1個點
# xi = [3,4.5,7,9]#存儲x的值
# n=len(xi)-1
# yi = [2.5,1,2.5,0.5]#存儲y 的值

# 輸入相關數值
n=input("請輸入取到樣點數目：")-1
xi = []
yi = []
for t in range(0,n+1):
    inputX = "請輸入第"+str(t+1)+"個點 X："
    inputY = "請輸入第"+str(t+1)+"個點 Y："
    xi.append(input(inputX))
    yi.append(input(inputY))



xi_2 = []# 存儲 x**2 的值
#向 x**2 中存入數據
for i in range(0,len(xi)):
    xi_2.append(xi[i]**2)
#將方程組轉化為矩陣 使 mat1 * mat2 = mat3
# 可以得到行列為3n 的矩陣
m1 = [[0 for a in range(3*n)] for b in range(3*n)]
m3 = [[0 for a in range(1)] for b in range(3*n)]
#若 mat2 內的值定為 a1,b1 ,c1,a2,b2,c2,a3,b3,c2... 便可以確定mat1的值
#向 m1 中存入數據
#定義變量 p 用於記錄向矩陣插入的行數，及后續插入
p=0
for j in range(0,n):

#從0 到n-1 的 n 個點代入
    m1[p][3*j]=xi_2[j]
    m1[p][3*j+1]=xi[j]
    m1[p][3*j+2]=1
    m3[p][0]=yi[j]
    p=p+1
#從 1 到 n 的 n 個點代入
    m1[p][3*j]=xi_2[j+1]
    m1[p][3*j+1]=xi[j+1]
    m1[p][3*j+2]=1
    m3[p][0]=yi[j+1]
    p=p+1
# 中間節點的一階導 相等 2ai + bi - 2a(i+1)-b(i+1) = 0
#從 1 到n-1 的 n-1 個 點處代入
for k in range(1,n):
    m1[p][3*(k-1)]= 2*xi[k]
    m1[p][3*(k-1)+1]= 1
    m1[p][3*k]= -2*xi[k]
    m1[p][3*k+1]= -1
    p=p+1
#  代入條件 二階導為零a1 = 0
m1[p][0] = 2
p=p+1
#將list轉化為 矩陣
mat1 = np.matrix(m1)
mat3 = np.matrix(m3)
#求mat2
_mat1 = mat1.getI()
mat2=_mat1*mat3
# mat2=mat3/mat1
#將矩陣轉化為list提取數據
m2=mat2.tolist()
#整理求得的曲線
line=[]
#用於收集區間
interval=[]
for q in range(0,n):
    interval.append(np.linspace(xi[q],xi[q+1]))
    a=m2[q*3]
    b=m2[q*3+1]
    c=m2[q*3+2]
    line.append(a*interval[q]**2+b*interval[q]+c)
    plt.plot(interval[q],line[q])
plt.show()