【Julia1.5.2】三次样条插值
本帖最后由 流星雨 于 2020-12-1 18:29 编辑参考数目：
【1】《数值分析》(李乃成、梅立泉编著，科学出版社，ISBN：978-7-03-032192-3)
背景：
当得到n个x和对应的点y时，如果需要得到一条插值曲线，并且曲线在每个点上二阶导数连续——曲线光滑，早期的做法是用压铁将一根富有弹性的细条(样条，Spline)固定在这些点上，然后沿着细条画出一根光滑的曲线【1】。
结构：
采用Julia语言(1.5.2)编制脚本Spline.jl，主要功能为求解三类边界条件(双侧二阶导数已知、双侧一阶导数已知、双侧周期)约束下的三次样条插值曲线。
脚本包含5个函数：cal_k(计算多项式(x±x_i )^n的系数)、seletbound(选择边界条件)、cal_M(计算每个点的二阶导数M)、Spline_out(插值函数输出部分)、cal_hμλd(计算hμλd，不包含首尾项)。另为方便调试，未封装“数据读入部分(读取表格数据)”和“计算设置部分(设置插值点并获取插值函数值)”。
算法原理：
使用说明：
安装Julia编译环境JuliaPro-1.5.2-1，可在URL：https://juliacomputing.com/products/juliapro/下载最新版的JuliaPro-1.5.3-1。安装完毕后用Juliapro打开文件Spline.jl。在REPL窗口中输入pwd()查看当前工作路径，可以选择在当前工作路径下放置输入文件xy.csv和boundary.csv，分别为xy点集和边界条件。也可在脚本文件Spline.jl中更改：
ID = "C:\\Users\\hp\\xy.csv"
ID_b = "C:\\Users\\hp\\boundary.csv"
其中C:\\Users\\hp为工作路径。在工作路径下放置xy.csv和boundary.csv，格式见附件。特别地，当选定“周期”边界条件时，应确保所给的数据满足式(3-2)，此时在“周期”列随便填入一个内容(不能是none，建议填1)，其余列填none。
单击▶即可运行。UI右侧Workspace窗口可实时显示各变量值。
计算设置部分(153行后)采用两种插值方式：单值和等距取点。其中变量xt表示单值插值点，单击▶计算后右侧窗口St为对应的插值函数值。number表示在插值区间内取等距点的数量，单击▶计算后右侧窗口Stest为对应的插值函数值集合。
运行结束后会在工作目录下生成插值函数表out.csv，第一列为插值点，第二列为对应的插值函数值。
附上参考书【1】给出的使用建议：
三次样条插值具有良好的收敛性、计算稳定性和二阶光滑性，当线性插值不可满足精度时可优先考虑；
慎用高次插值，避免龙格现象。
我表示这些数学公式我看不下去，但样条线我还是懂的 我来个python的！
# -*- coding: utf-8 -*-
"""
Created on Sun Dec 27 19:30:57 2020
Copyright@Xi'an Jiaotong University - NuTHeL
@author: Zhou Xingguang
"""
import numpy as np
import matplotlib.pyplot as plt
class CubicSpline(object):
def __init__(self, x, y):
self.x = x[:]
self.y = y[:]
@staticmethod
def SecondOrderDiff(x, y):
return np.array([((y-y)/(x-x)-(y-y)/(x-x))/(x-x) for i in range(len(x)-2)])
def Interpolation(self):
h = np.array(-self.x for i in range(len(self.x)-1)])
u = np.array( / (h+h) for i in range(len(h)-1)])
my_lambda = 1 - u
d = 6 * CubicSpline.SecondOrderDiff(self.x, self.y)
# generate the M matrix
Matrix = 2*np.identity(np.size(d), dtype=float)
Matrix = my_lambda
Matrix = u
for i in range(1, np.size(d)-1):
Matrix = u
Matrix = my_lambda
# solve the function group
M = np.linspace(0, 1, len(self.x))
M = 0
M = np.linalg.solve(Matrix, d)
M = 0
return M, h
def GenerateFunc(self, x):
M, h = self.Interpolation()
# generate the coefficient matrix
coeff = np.zeros((np.size(self.x)-1, 4), dtype=float)
for i in range(0, np.size(self.x)-1):
coeff = M / (6*h)
coeff = M / (6*h)
coeff = (self.y - h**2*M/6)/h
coeff = (self.y - h**2*M/6)/h
# get the x in which section
location = 0
for i in range(len(self.x)-1):
if self.x <= x <= self.x:
break
location += 1
result = (self.x-x)**3*coeff
result += (x-self.x)**3*coeff
result += (self.x-x)*coeff
result += (x-self.x) * coeff
return result
if __name__ == '__main__':
x = np.linspace(-1, 1, 11)
y = 1/(1+25*x**2)
C = CubicSpline(x, y)
for i in range(10):
test_x = np.linspace(x, x, 5)
test_y = []
for j in range(0, 5):
test_y.append(C.GenerateFunc(test_x))
plt.scatter(test_x, test_y)
plt.scatter(test_x, test_y, marker='x')
# draw the analytic solution
plt.plot(test_x, test_y)
Graphic_x = np.linspace(-1, 1, 100)
Graphic_y = 1/(1+25*Graphic_x**2)
plt.plot(Graphic_x, Graphic_y)
plt.show()
插值结果如下所示，其中橘色的实线为解析解y=1/(25*x^2)的图像，圆点和打x的点为插值结果，每两个圆点代表了插值的分段区间。
页:
[1]