【Julia1.5.2】三次样条插值

流星雨

参考数目：
【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】给出的使用建议：
三次样条插值具有良好的收敛性、计算稳定性和二阶光滑性，当线性插值不可满足精度时可优先考虑；
慎用高次插值，避免龙格现象。
三次样条插值.rar (128.51 KB, 下载次数: 0)

2020-12-1 17:44 上传

点击文件名下载附件

0xAA55

我表示这些数学公式我看不下去，但样条线我还是懂的

watermelon

我来个python的！

1. # -*- coding: utf-8 -*-
   
2. """
   
3. Created on Sun Dec 27 19:30:57 2020
   
4. Copyright@Xi'an Jiaotong University - NuTHeL
   
5. @author: Zhou Xingguang
   
6. """
   
7. 
   
8. import numpy as np
   
9. import matplotlib.pyplot as plt
   
10. 
    
11. class CubicSpline(object):
    
12.     def __init__(self, x, y):
    
13.         self.x = x[:]
    
14.         self.y = y[:]
    
15.    
    
16.     @staticmethod
    
17.     def SecondOrderDiff(x, y):
    
18.         return np.array([((y[i+2]-y[i+1])/(x[i+2]-x[i+1])-(y[i+1]-y[i])/(x[i+1]-x[i]))/(x[i+2]-x[i]) for i in range(len(x)-2)])
    
19.    
    
20.     def Interpolation(self):
    
21.         h = np.array([self.x[i+1]-self.x[i] for i in range(len(self.x)-1)])
    
22.         u = np.array([h[i] / (h[i]+h[i+1]) for i in range(len(h)-1)])
    
23.         my_lambda = 1 - u
    
24.         d = 6 * CubicSpline.SecondOrderDiff(self.x, self.y)
    
25.         # generate the M matrix
    
26.         Matrix = 2*np.identity(np.size(d), dtype=float)
    
27.         Matrix[0, 1] = my_lambda[0]
    
28.         Matrix[np.size(d)-1, np.size(d)-2] = u[np.size(d)-1]
    
29.         for i in range(1, np.size(d)-1):
    
30.             Matrix[i, i-1] = u[i]
    
31.             Matrix[i, i+1] = my_lambda[i]
    
32.         # solve the function group
    
33.         M = np.linspace(0, 1, len(self.x))
    
34.         M[0] = 0
    
35.         M[1:len(self.x)-1] = np.linalg.solve(Matrix, d)
    
36.         M[len(self.x)-1] = 0
    
37.         return M, h
    
38.    
    
39.     def GenerateFunc(self, x):
    
40.         M, h = self.Interpolation()
    
41.         # generate the coefficient matrix
    
42.         coeff = np.zeros((np.size(self.x)-1, 4), dtype=float)
    
43.         for i in range(0, np.size(self.x)-1):
    
44.             coeff[i][0] = M[i] / (6*h[i])
    
45.             coeff[i][1] = M[i+1] / (6*h[i])
    
46.             coeff[i][2] = (self.y[i] - h[i]**2*M[i]/6)/h[i]
    
47.             coeff[i][3] = (self.y[i+1] - h[i]**2*M[i+1]/6)/h[i]
    
48.         
    
49.         # get the x in which section
    
50.         location = 0
    
51.         for i in range(len(self.x)-1):
    
52.             if self.x[i] <= x <= self.x[i+1]:
    
53.                 break
    
54.             location += 1
    
55.         
    
56.         result = (self.x[location+1]-x)**3*coeff[location][0]
    
57.         result += (x-self.x[location])**3*coeff[location][1]
    
58.         result += (self.x[location+1]-x)*coeff[location][2]
    
59.         result += (x-self.x[location]) * coeff[location][3]
    
60.         return result
    
61. 
    
62. 
    
63. if __name__ == '__main__':
    
64.     x = np.linspace(-1, 1, 11)
    
65.     y = 1/(1+25*x**2)
    
66.     C = CubicSpline(x, y)
    
67.     for i in range(10):
    
68.         test_x = np.linspace(x[i], x[i+1], 5)
    
69.         test_y = []
    
70.         for j in range(0, 5):
    
71.             test_y.append(C.GenerateFunc(test_x[j]))
    
72.         plt.scatter(test_x[0], test_y[0])
    
73.         plt.scatter(test_x[1:5], test_y[1:5], marker='x')
    
74.    
    
75.     # draw the analytic solution
    
76.     plt.plot(test_x, test_y)
    
77.    
    
78.     Graphic_x = np.linspace(-1, 1, 100)
    
79.     Graphic_y = 1/(1+25*Graphic_x**2)
    
80.     plt.plot(Graphic_x, Graphic_y)
    
81.     plt.show()

复制代码

插值结果如下所示，其中橘色的实线为解析解y=1/(25*x^2)的图像，圆点和打x的点为插值结果，每两个圆点代表了插值的分段区间。