三次样条插值是要解决,在模拟函数的时候,插值节点处不平滑的问题
S(x)是f(x) (x0 .... xn)的三次样条插值函数,则S(x)要在区间(xj,xj+1)上满足S(xj-0) = S(xj+1+0) , S'(xj-0) = S'(xj+1 + 0), s''(xj - 0) = S''(xj+1 +0),
这样在每一个这样的区间建立一个等式条件,同时在两个端点出处得到两个确定值
(S'(x0) = f'(x0),S'(xn) = f'(xn) ,S''(x0) = f''(x0) S''(xn) = f''(xn))
通过上述的等式可以得出S(x)的最终形式
因为S(x)是三次的,导数两次之后就是线性的,在两点之间可以求出两次导数后的S''(x)在往回积分得到下面的式子
S(x) = Mj*(xj+1 - x)^3/(6*hj) + Mj+1(x - xj)^3/(6 * hj) + (yj - Mj*hj^2/6)*(xj+1 - x)/hj+(yj+1 - Mj+1 * hj^2)*(x - xj)/hj j = 0,1,......n-1
hj 就是 xj+1 - xj,现在就是需要求 Mj (Mj是S''(xj)的值)
由求得的S(x)导数一次 + S'(xj-0) = S'(xj+1+0) 可以得到
uj * Mj-1 + 2Mj +lj * Mj+1 = dj
uj = hj-1/(hj-1 + hj) lj = 1 - uj
S(x)是f(x) (x0 .... xn)的三次样条插值函数,则S(x)要在区间(xj,xj+1)上满足S(xj-0) = S(xj+1+0) , S'(xj-0) = S'(xj+1 + 0), s''(xj - 0) = S''(xj+1 +0),
这样在每一个这样的区间建立一个等式条件,同时在两个端点出处得到两个确定值
(S'(x0) = f'(x0),S'(xn) = f'(xn) ,S''(x0) = f''(x0) S''(xn) = f''(xn))
通过上述的等式可以得出S(x)的最终形式
因为S(x)是三次的,导数两次之后就是线性的,在两点之间可以求出两次导数后的S''(x)在往回积分得到下面的式子
S(x) = Mj*(xj+1 - x)^3/(6*hj) + Mj+1(x - xj)^3/(6 * hj) + (yj - Mj*hj^2/6)*(xj+1 - x)/hj+(yj+1 - Mj+1 * hj^2)*(x - xj)/hj j = 0,1,......n-1
hj 就是 xj+1 - xj,现在就是需要求 Mj (Mj是S''(xj)的值)
由求得的S(x)导数一次 + S'(xj-0) = S'(xj+1+0) 可以得到
uj * Mj-1 + 2Mj +lj * Mj+1 = dj
uj = hj-1/(hj-1 + hj) lj = 1 - uj
dj = 6f[xj-1,xj,xj+1]均差
其中 d0 = 6/h0 * (f[x0,x1] - f'(x0)) dn = 6/hn-1 * (f'(xn) - f[xn-1,xn])
可以根据这些写出矩阵
2 l0
u1 2 l1 *M = d[i =0 ....n]
u2 2 l2
.....
un 2
这样可以得到M了,其他需要的也可以计算了.
画图的时候用到了matplotlib
import scipy.interpolate as itp
import matplotlib.pyplot as plt
import numpy as np
''''
线性插值
'''
def line_insert(xn,yn):
def line(x):
for value in range(len(xn)-1):
if xn[value] <= x <= xn[value+1]:
result1 = ((x-xn[value+1])/(xn[value] - xn[value+1]))*yn[value]
result2 = ((x-xn[value])/(xn[value+1] - xn[value]))*yn[value+1]
return result1+result2
return line
这个是线性插值,返回一个直接计算的函数
接下来是实现的三次样条的茬插值
'''
三次样条
'''
last_trax = []
#均差 求d的时候
def sub_spl(x,y):
if len(x) == 2:
return (y[1] -y[0])/(x[1] - x[0])
return (sub_spl(x[1:],y[1:]) - sub_spl(x[:len(x)-1],y[:len(y) - 1]))/(x[len(x) - 1] - x[0])
def three_spline_trax(tup_x,tup_y,s0,sn):
trax = []
h = []
for i in range(len(tup_x) - 1):
h.append(tup_x[i+1] - tup_x[i])
d = [0 for i in range(len(tup_x))]
d[0] = 6/h[0] * (sub_spl(tup_x[:2],tup_y[:2]) - s0)
d [-1] = 6/h[-1] * (sn - sub_spl(tup_x[-2:],tup_y[-2:]))
u = [0 for i in range(len(tup_x)-1)]
for j in range(1,len(tup_x) - 1):
d[j] = 6 * sub_spl(tup_x[j-1:j+2],tup_y[j-1:j+2])
u[j] = h[j-1]/(h[j-1] + h[j])
l = [1-i for i in u]
u.append(1) #un = 1
for i in range(len(tup_x)):
trax.append([0 for j in range(len(tup_x))])
trax[i][i] = 2
for i in range(len(tup_x) - 1):
trax[i][i+1] = l[i]
trax[i+1][i] = u[i+1]
trax[i].append(d[i])
trax[-1].append(d[-1])
return trax,h
#对得到的矩阵简化上三角
def return_trax(trax,m,n,num):
if num is 1:
last_trax.append(trax[n-1][n-1:])
return
else:
#寻找主元
Max = abs(trax[0][0])
t1 = 0
for i in range(1,m):
if abs(trax[i][0]) > Max:
Max = trax[i][0]
t1 = i
trax[0],trax[t1] = trax[t1],trax[0]
tmp2 = []
last_trax.append(trax[0])
for i in range(1,n):
#倍数
tmp = -trax[i][0]/trax[0][0]
def cal(tup):
return tup[0] + tup[1] * tmp
tmp2.append(map(cal,zip(trax[i][1:],trax[0][1:])))
trax = tmp2
return return_trax(trax,m-1,n-1,num-1)
def s_um(tup):
return reduce(lambda x,y:x*y,tup)
#这个函数得到M0 - Mn
def calculate(trax,m):
#倒着循环
l = [0 for i in range(m)]
for i in range(m-1,-1,-1):
if i == m-1:
#直接求出最后一个参数的数值
l[i] = trax[i][m-i]/trax[i][0]
else:
tmp3 = trax[i][m-i]
tmp3 -= sum(map(s_um,zip(l[i:],trax[i][:m-i])))
l[i] = tmp3/trax[i][0]
return l
def res(l,tup_x1,tup_y1,h):
def outcome(x):
for j in range(len(tup_x1) - 1):
if tup_x1[j] <= x <= tup_x1[j+1]:
sum_func = l[j]*(tup_x1[j+1] - x)**3/(6*h[j]) + l[j+1]*(x - tup_x1[j])**3/(6*h[j])+(tup_y1[j] - l[j]/6*h[j]**2)*(tup_x1[j+1] - x)/h[j] + (tup_y1[j+1]-l[j+1]/6*h[j]**2)*(x - tup_x1[j])/h[j]
return sum_func
return outcome
在计算矩阵的乘法的时候,直接用上次的列主元消去....
最后的res函数就是在构造S(x)
three_line_trax 就是返回上面写的那个带2的矩阵,h就是相邻x的插值
最后再用return_trax返回上三角矩阵,calculate计算就可以
可能中间有些地方说的不准确,如果发现请指正