Python 三次样条插值法

代码

'''

本函数通过三次样条插值法进行函数值计算

'''

# 三次样条插值
import numpy as np

# 用于存放x，y，m的值
x = np.array([1,2,4,5])
y = np.array([1,3,4,2])
m = np.array([17/8,None,None,-19/8])
lens = len(x)
x_f = 3.0     # 待插值点
# 用于存放几种数据
h = np.zeros(lens-1)
a = np.zeros(lens)
b = np.zeros(lens)
n = np.zeros([lens,lens])

for i in range(len(h)):
    h[i] = x[i+1] - x[i]

for i in range(1,len(a)-1):
    a[i] = h[i-1]/(h[i]+h[i-1])
    b[i] = 3*(((1-a[i])*(y[i]-y[i-1]))/h[i-1]+(a[i]*(y[i+1]-y[i]))/h[i])
a[0] = 1
b[0] = (3*(y[1]-y[0]))/h[0]
b[-1] = (3*(y[-1]-y[-2]))/h[-1]

for i in range(len(n)):
    n[i][i] = 2
    if i != len(n) - 1:
        n[i][i + 1] = a[i]
    if i != 0:
        n[i][i-1] = 1 - a[i]

for i in range(len(n)-1):
    if i == 0:
        m[i+1] = (b[i] - n[i][i]*m[i])/a[i]
    else:
        m[i+1] = (b[i] - n[i][i]*m[i] - n[i][i-1]*m[i-1])/a[i]

p = 0
for i in range(len(x)-1):
    if x[i] < x_f <= x[i + 1]:
        p = i
        break

if p == 0:
    print('待求值点超出函数范围！')
    exit()

fun_1 = (1+2*(x_f-x[p])/(x[p+1]-x[p]))*(((x_f-x[p+1])/(x[p]-x[p+1]))**2)*y[p]
fun_2 = (1+2*(x_f-x[p+1])/(x[p]-x[p+1]))*(((x_f-x[p])/(x[p+1]-x[p]))**2)*y[p+1]
fun_3 = (x_f-x[p])*(((x_f-x[p+1])/(x[p]-x[p+1]))**2)*m[p]
fun_4 = (x_f-x[p+1])*(((x_f-x[p])/(x[p+1]-x[p]))**2)*m[p+1]

fun = fun_1 + fun_2 + fun_3 + fun_4

print('函数在点 {0} 处的插值函数值为: {1} '.format(x_f,fun))