技术宅的结界

 找回密码
 立即注册→加入我们

QQ登录

只需一步,快速开始

搜索
热搜: 下载 VB C 实现 编写
查看: 1426|回复: 2
收起左侧

【Julia1.5.2】三次样条插值

[复制链接]

4

主题

7

帖子

392

积分

用户组: 中·技术宅

UID
6285
精华
4
威望
43 点
宅币
208 个
贡献
71 次
宅之契约
0 份
在线时间
4 小时
注册时间
2020-10-2
发表于 2020-12-1 17:44:33 | 显示全部楼层 |阅读模式

欢迎访问技术宅的结界,请注册或者登录吧。

您需要 登录 才可以下载或查看,没有帐号?立即注册→加入我们

x
本帖最后由 流星雨 于 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,不包含首尾项)。另为方便调试,未封装“数据读入部分(读取表格数据)”和“计算设置部分(设置插值点并获取插值函数值)”。
QQ图片20201201173805.png

算法原理:
1.png 2.png 3.png 4.png

使用说明:
安装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)
5.png

评分

参与人数 3威望 +13 宅币 +46 贡献 +17 收起 理由
0xAA55 + 2 + 10 + 2
Ink_Hin_fifteen + 1 + 6 + 5
watermelon + 10 + 30 + 10 赞!

查看全部评分

回复

使用道具 举报

1088

主题

2606

帖子

7万

积分

用户组: 管理员

一只技术宅

UID
1
精华
236
威望
474 点
宅币
21358 个
贡献
45937 次
宅之契约
0 份
在线时间
2058 小时
注册时间
2014-1-26
发表于 2020-12-11 01:57:15 | 显示全部楼层
我表示这些数学公式我看不下去,但样条线我还是懂的

29

主题

331

帖子

1999

积分

用户组: 上·技术宅

UID
3808
精华
11
威望
105 点
宅币
1238 个
贡献
165 次
宅之契约
0 份
在线时间
378 小时
注册时间
2018-5-6
发表于 2020-12-28 09:07:35 | 显示全部楼层
我来个python的!
[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[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)])
    
    def Interpolation(self):
        h = np.array([self.x[i+1]-self.x[i] for i in range(len(self.x)-1)])
        u = np.array([h[i] / (h[i]+h[i+1]) 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[0, 1] = my_lambda[0]
        Matrix[np.size(d)-1, np.size(d)-2] = u[np.size(d)-1]
        for i in range(1, np.size(d)-1):
            Matrix[i, i-1] = u[i]
            Matrix[i, i+1] = my_lambda[i]
        # solve the function group
        M = np.linspace(0, 1, len(self.x))
        M[0] = 0
        M[1:len(self.x)-1] = np.linalg.solve(Matrix, d)
        M[len(self.x)-1] = 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[i][0] = M[i] / (6*h[i])
            coeff[i][1] = M[i+1] / (6*h[i])
            coeff[i][2] = (self.y[i] - h[i]**2*M[i]/6)/h[i]
            coeff[i][3] = (self.y[i+1] - h[i]**2*M[i+1]/6)/h[i]
        
        # get the x in which section
        location = 0
        for i in range(len(self.x)-1):
            if self.x[i] <= x <= self.x[i+1]:
                break
            location += 1
        
        result = (self.x[location+1]-x)**3*coeff[location][0]
        result += (x-self.x[location])**3*coeff[location][1]
        result += (self.x[location+1]-x)*coeff[location][2]
        result += (x-self.x[location]) * coeff[location][3]
        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[i], x[i+1], 5)
        test_y = []
        for j in range(0, 5):
            test_y.append(C.GenerateFunc(test_x[j]))
        plt.scatter(test_x[0], test_y[0])
        plt.scatter(test_x[1:5], test_y[1:5], 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的点为插值结果,每两个圆点代表了插值的分段区间。
Figure 2020-12-28 090709.png
Passion Coding!

本版积分规则

QQ|申请友链||Archiver|手机版|小黑屋|技术宅的结界 ( 滇ICP备16008837号 )|网站地图

GMT+8, 2021-10-22 05:35 , Processed in 0.039459 second(s), 37 queries , Gzip On.

Powered by Discuz! X3.4

Copyright © 2001-2021, Tencent Cloud.

快速回复 返回顶部 返回列表