流星雨 发表于 2020-12-1 17:44:33

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

0xAA55 发表于 2020-12-11 01:57:15

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

watermelon 发表于 2020-12-28 09:07:35

我来个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]
查看完整版本: 【Julia1.5.2】三次样条插值