（十三）二次样条插值

 

#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()