三次样条曲线插值的基本原理及其C#实现
声明：本人空间的所有文章，若无特别声明，皆为本人原创，可自由转载，但要注明原作者和原始出处，不可作为商业用途。
下面的内容是直接从Word文档复制粘贴出来的，有很多内容丢失，完整的PDF版本可到百度网盘下载：
链接：https://pan.baidu.com/s/1aEc0htnERV3I75hLq4fFsA 密码：wf6t
三次样条曲线插值的基本原理及其C#实现
作 者: simxpert
2018.5.14
作者简介：
本人长期从事结构强度仿真计算方面的学习和工作，可以为您解决各种结构强度仿真分析以及专业计算软件定制方面的工作。
扫描二维码可添加本人微信：
目 录
1. 简介 1
2. 三次样条插值的基本数学原理 1
2.1. 插值的问题的提出 1
2.2. 插值函数的待定系数变量和约束方程 2
2.3. 插值函数的最终表达式 3
2.4. 添加边界条件求解mi 3
2.5. 根据插值函数进行插值计算 4
3. 三次样条插值函数的C#实现 4
4. SPLine类的使用方法 6
5. 附录1-SPLine类的源代码 7
6. 附录2-Chase类的源代码 11
三次样条曲线插值的数学原理及其C#实现
作者: simxpert
简介
在工程实践中经常会用到插值，最简单的插值是在相邻的两个样本点之间线性插值，有时候希望精度高一点，一般都会用三次样条曲线插值。
本文先从理论上讲解三次样条插值的数学原理，然后讲解使用C#实现三次样条插值算法的基本过程，最后以实例对算法进行验证。
三次样条插值的基本数学原理
本文并不打算详细讲解完整的三次样条曲线的推导过程，只是梳理一下三次样条曲线差值的大概的脉络，详细的推导过程可以参考随意一本数值计算教程。
插值的问题的提出
在给定的区间[a,b]上有N个点，已知其横坐标为：a=x0<x1<……<xn-1<xn=b，对应纵坐标为：y0,y1,……yn-1,yn。
这N个点的每两个相邻点xi,xi+1构成一个子区间[xi,xi+1]，[a,b]之间共计有n个子区间：[x0,x1],…,[xn-1,xn] 。(N=n+1)
问题：根据上述已知条件，任意给定一个在xt，xt∈[a,b]，如何近似求得一个合理的yt值？
一个最简单易行的办法就是线性插值：在每个子区间[xi,xi+1]中，把两个端点用直线连接起来，可以由这段直线的方程Si(x)来求这个子区间上任意一点的纵坐标值。我们只需要确定给定的点xt落在哪个子区间就可以了。
线性插值的缺点很明显：首先是精度低。其次，整个曲线不光滑。因为在整个区间的插值函数S(x)是由多段的直线Si(x)连接而成，在每个拐点处连接是不平滑，从数学上看也就是在拐点处的一阶导数不连续，更别提二阶段导数的连续性了。
插值函数的待定系数变量和约束方程
为了得到精度比较高而且过渡比较平滑的曲线，最常见的就是使用三次样条插值。如果我们把前面讲的子区间的线性函数提升为3次多项式：
(1)
为了保证曲线的连续性和光滑性，必须要对每个区间的插值函数提出一些要求，这些要求包括：
1).连续性。除了两个边界节点外，任意一个节点xi同时属于[xi-1,xi]和[xi,xi+1]这两个相邻的子区间。我们必须要保证在用这两个区间上的插值函数计算出来的xi的纵坐标值相等，而且要等于已知的yi值。即：
Si-1(xi)=Si(xi)=yi,i=1,2,…,n-1。 (2)
式(2)包含了2n-2个方程，再加两个边界点上S0(x0)=y0,Sn-1(xn)=yn，共计得到2n个方程。
2).一阶导数和二阶导数连续。这个约束条件是为了保证插值函数的曲线足够的光滑。除了两头的节点外，任意一个中间节点的一阶导数和二阶导数必须连续。即：
,i=1,2,…,n-1 (3)
,i=1,2,…,n-1 (4)
式(3),(4)可以得到2n-2个方程。
由于上述约束条件总计可以得到4n-2个方程。从(1)式可以看出，每个插值函数有4个未知的待定系数，[a,b]区间内有n个子区间，故一共有4n个未知待定系数。未知的变量的数目比约束方程的数目多2，也就是说我们还缺少两个约束条件。
缺少的这两个约束条件可以由边界条件来指定。也就是在两个边界点上各施加一个边界条件。边界条件是是人为指定的，常见的边界条件有4种，这个我们在后面会提到。
插值函数的最终表达式
根据前面的论述，我们只要求出每个子区间的插值函数的4个系数a0, a1, a2,a3就能完全确定这n个插值函数。但是实际中，为了数学推导的方便，我们并不会直接求出这4个变量，而是采取了间接的表达方式来描述每个插值函数。
不同的推导方法，得到的表达形式也各不相同，常见的有三转角法，三弯矩法，B样条基函数法等方法。本文并不打算进行详细的推导。因为没必要，任何一本数值计算的书上都可以找到非常详细的。在这里我直接给出推导结果:
(5)
其中 。
(6)
式(5)就是插值函数的最终表达式。式(5)看起来比(1)复杂多了，其实仔细看，如果全部展开，仍然是一个三次多项式。在这个公式中，只有 是未知的，其他都可以根据已知条件直接算出来的。剩下的任务就是如何求出 。
我们把(5)所代表的插值函数代入到之前说的约束条件中(2),(3),(4)中，经过整理，可以得到式(7)：
，i=1,2,…n-1. (7)
其中 ， ， 。
添加边界条件求解mi
前面讲过，我们还需要添加两个边界条件才能求解出位置系数。常见的边界条件大概有三、四种。本文不打算涉及所有的边界条件，只选用比较常用的第二类边界条件(又叫”自然边界条件”)：两个边界点的二阶导数为给定的值，其中最常见给定值为0，即在两个边界点上有：
， (8)
根据式(7)，结合边界条件(8)，可以推导出求解系数mi的矩阵：
= (9)
从上面的公式中λ, μ, g都是可以根据已知条件算出来的，只有mi是未知量。式(9)中的矩阵是一个典型的严格对角占优的对角方程组，可以用追赶法很快得到解。
根据插值函数进行插值计算
在根据2.4章中的(9)式求得mi之后，代入到2.3 章中的式(5)，即得到所有的子区间上的插值函数Si(x)。
对于一个给定的待插值点xt，我们先要确定xt落入到哪个子区间，然后调用该子区间的插值函数Si(x)即可求出xt对应的插值yt。
三次样条插值函数的C#实现
在C#中实现了三次样条插值函数，实现过程封装在SPLine类中。SPLine类的成员变量为：
class SPLine
{
double[] Xi;//存储已知点的x值
double[] Yi;//存储已知点的y值；
double[] A;//存储追赶法矩阵中的A
double[] B; //存储追赶法矩阵中的B
double[] C; //存储追赶法矩阵中的C
double[] H;//存储前面公式中的hi
double[] Lamda;//存储前面公式中的λi
double[] Mu;//存储μi
double[] G;//存储gi
double[] M;//存储待求的未知量mi
int N;//已知点的总数
int n;//n=N-1。 n为区间数，N为点数。
}
成员函数介绍：
private void GetH()
功能是求取前文中的h。hi=xi-xi-1。
private void GetLamda_Mu_G()
功能是根据已知条件求出λ, μ, g。计算公式见第3页中间的公式。
private void GetABC()
功能是求出追赶法中的矩阵里面的A,B,C，至于追赶法里面的A,B,C是怎么回事，这里简单给一个截图：
数组A，B，C分别保存的是上图中的a,b,c。
上图中的待求变量x对应本文中的待求变量m。上图中的d对应本文中的数组G。
private double fai0(double x)
private double fai1(double x)
这两个函数对应2.3章中的式(6)中的两个三次多项式。
private int GetSection(double x)
判定待求值的x属于哪个子区间，并返回区间编号。
有一种特殊情况：如果给定的插值点x不在区间[a,b]范围内，是不能应用插值函数的。基于实际项目的需要，本文中的程序对这种情况的处理方式：取距离x最近的边界节点的值作为其插值结果。当然也可以根据需要改成抛出”超出边界”的异常。
当x比最小节点x0还要小，则GetSection返回的区间编号为-1，如果x比最大节点xn还要大，则GetSection返回的区间编号为-999。相应地，如果返回-1，则S(x0)=f(x0)=y0,如果返回-999则S(x0)=f(xn-1)=yn-1
public double Interpolate(double x)
该函数为计算x的插值的插值函数。计算公式依据的是第2.3章中的式(6)。
public bool Init(double[] Xi, double[] Yi)
该函数为初始化函数。把[a,b]区间的(xi,yi)作为输入数据，初始化SPLine类的成员变量，并调用追赶法Chase.Solve()求出待定系数mi。
在SPLine类中引用了追赶法求解代数方程组的Chase类，Chase类也是自己编写的专门求解形如式(9)的方程组的。
SPLine类和Chase类的代码在附录中给出。
SPLine类的使用方法
首先当然是把SPLine.cs和Chase.cs这两个个代码源文件加入到solution中，然后在要使用该类的地方添加对命名空间的引用using SPLINE; using MyAlgbra;
下面是一段验证代码，通过这段代码也同时说明了该如何使用SPLine类。
private void button2_Click(object sender, EventArgs e)
{
double[] X = {1,2,4,5};
double[] Y = {1,3,4,2};
double s;
SPLine sp = new SPLine();
sp.Init(X, Y);
s = sp.Interpolate(0.5);
}
附录1-SPLine类的源代码
完整的SPLine.cs的代码如下：
using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using System.Collections;
using MyAlgbra;
namespace SPLINE
{
class SPLine
{
double[] Xi;
double[] Yi;
double[] A;
double[] B;
double[] C;
double[] H;
double[] Lamda;
double[] Mu;
double[] G;
double[] M;
int N;//
int n;//
public SPLine()
{
N = 0;
n = 0;
}
public bool Init(double[] Xi, double[] Yi)
{
if (Xi.Length != Yi.Length)
return false;
if (Xi.Length == 0)
return false;
//根据给定的Xi,Yi元素的数目来确定N的大小。
this.N = Xi.Length;
n = N - 1;
//根据实际大小给各个成员变量分配内存
A = new double[N-1];
B = new double[N];
C = new double[N-1];
this.Xi = new double[N];
this.Yi = new double[N];
H = new double[N-1];
Lamda = new double[N-1];
Mu = new double[N-1];
G = new double[N];
M = new double[N];
//把输入数据点的数值赋给成员变量。
for (int i = 0; i <= n; i++)
{
this.Xi[i] = Xi[i];
this.Yi[i] = Yi[i];
}
/************ 求出hi,Lamda(i),Mu(i),gi *******************/
GetH();
GetLamda_Mu_G();
GetABC();
/***************** 调用追赶法求出系数矩阵M *********************/
Chasing chase = new Chasing();
chase.Init(A, B, C, G);
chase.Solve(out M);
return true;
}
private void GetH()
{
//Get H first;
for (int i = 0; i <=n-1; i++)
{
H[i] = Xi[i + 1] - Xi[i];
}
}
private void GetLamda_Mu_G()
{
double t1, t2;
for (int i = 1; i <= n - 1; i++)
{
Lamda[i] = H[i] / (H[i] + H[i - 1]);
Mu[i] = 1 - Lamda[i];
t1 = (Yi[i] - Yi[i - 1]) / H[i - 1];
t2 = (Yi[i + 1] - Yi[i]) / H[i];
G[i]=3*(Lamda[i]*t1+Mu[i]*t2);
}
G[0] = 3*(Yi[1] - Yi[0])/H[0];
G[n] = 3*(Yi[n] - Yi[n-1])/H[n-1];
Mu[0] = 1;
Lamda[0] = 0;
}
private void GetABC()
{
for (int i = 1; i <= n - 1; i++)
{
A[i-1] = Lamda[i];
C[i] = Mu[i];
}
A[n-1] = 1; C[0] = 1;
for (int i = 0; i <= n; i++)
{
B[i] = 2;
}
}
private double fai0(double x)
{
double t1, t2;
double s;
t1 = 2 * x + 1;
t2 = (x - 1) * (x - 1);
s = t1 * t2;
return s;
}
private double fai1(double x)
{
double s;
s=x*(x - 1) * (x - 1);
return s;
}
public double Interpolate(double x)
{
double s = 0;
double P1, P2;
double t = x;
int iNum;
iNum=GetSection(x);
if (iNum == -1) //
{
iNum = 0;
t = Xi[iNum];
return Yi[0];
}
if (iNum == -999)//
{
iNum = n - 1;
t = Xi[iNum+1];
return Yi[n];
}
P1 = (t - Xi[iNum]) / H[iNum];
P2 = (Xi[iNum + 1] - t) / H[iNum];
s = Yi[iNum] * fai0(P1) + Yi[iNum + 1] * fai0(P2) + M[iNum] * H[iNum] * fai1(P1) - M[iNum+1] * H[iNum] * fai1(P2);
return s;
}
private int GetSection(double x)
{
int iNum = -1;
// double EPS = 1.0e-6;
if (x < Xi[0])
{
return -1;
}
if (x > Xi[N - 1])
{
return -999;
}
for (int i = 0; i <= n - 1; i++)
{
if (x >= Xi[i] && x <= Xi[i + 1])
{
iNum = i;
break;
}
}
return iNum;
}
}
}
附录2-Chase类的源代码
追赶法求线性代数方程组的Chase类完整代码如下：
using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using System.Windows.Forms;
namespace MyAlgbra
{
class Chasing
{
protected int N;//Dimension of Martrix Ax=d;
protected double[] d;//Ax=d;
protected double[] Aa;//a in A;
protected double[] Ab; //b in A:
protected double[] Ac;// c in A;
protected double[] L;//LU-->L;
protected double[] U;//LU-->U;
public double[] S;//store the result;
//constructor without parameters;
public Chasing()
{
}
public bool Init(double[] a ,double[] b,double[] c,double[] d)
{
//check validation of dimentions;
int na = a.Length;
int nb = b.Length;
int nc = c.Length;
int nd = d.Length;
if(nb<3)
return false;
N = nb;
if (na != N - 1 || nc != N - 1 ||nd!=N)
return false;
S = new double[N];
L = new double[N - 1];
U = new double[N];
Aa = new double[N - 1];
Ab = new double[N];
Ac = new double[N - 1];
this.d = new double[N];
//init Aa,Ab,Ac,Ad;
for (int i = 0; i <= N - 2; i++)
{
Ab[i] = b[i];
this.d[i] = d[i];
Aa[i] = a[i];
Ac[i] = c[i];
}
Ab[N - 1] = b[N - 1];
this.d[N - 1] = d[N - 1];
return true;
}
public bool Solve(out double[] R)
{
R = new double[Ab.Length]; /*********************A=LU***********************************/
U[0] = Ab[0];
for (int i = 2; i <= N; i++)
{
// L[i] = Aa[i] / U[i - 1];
L[i - 2] = Aa[i - 2] / U[i - 2];
//U[i]=Ab[i]-Ac[i-1]*L[i];
U[i-1] = Ab[i-1] - Ac[i - 2] * L[i-2];
}
/*************************END of A=LU **********************/
/**************** Ly=d ******************************/
double[] Y = new double[d.Length];
Y[0] = d[0];
for (int i = 2; i <= N; i++)
{
//Y[k]=d[k]-L[k]*Y[k-1];
Y[i - 1] = d[i - 1]-(L[i-2]) * (Y[i - 2]);
}
/**************** End of Ly=d ****************************/
/*************** Ux=Y ********************************/
//X[n]=Y[n]/U[n];
R[N - 1] = Y[N - 1] / U[N - 1];
//X[k]=(Y[k]-C[k]*X[k+1])/U[k];(n-1,,.....1)
for (int i =N-1; i >= 1; i--)
{
R[i - 1] = (Y[i - 1] - Ac[i - 1] * R[i]) / U[i - 1];
}
/*************** End of Ux=Y *************************/
for (int i = 0; i < R.Length; i++)
{
if (double.IsInfinity(R[i]) || double.IsNaN(R[i]))
return false;
}
return true;
}
}
}
三次样条曲线插值的基本原理及其C#实现相关推荐
- 三次样条曲线插值（cubic spline)实例应用
目标 工作需要,需要达成这样得一个需求,给一系列得三维点,三维点按照顺序连接,形成一条折线.需要依照这条折线,进行曲线1m等距离插值.具体如下图 其中,红色圆圈点为原始点集OrigPoints(原始点 ...
- MATLAB---构造一个插值三次样条曲线
function InterpCubicSplineCurv() %本程序的功能是构造一个插值三次样条曲线N = 12; SamplPs = CollectSPFergusonCi2(N);%采集型值 ...
- 数值计算 --- 三次样条函数插值(Cubic spline function interpolation)
三次样条函数插值(Cubic spline function interpolation) Part I 插值 预备知识: 什么是插值? 已知部分离散的数据,但不知道满足这些数据 ...
- 机械臂规划----三次样条曲线
机械臂规划----三次样条曲线 原理讲解 源代码 三次样条曲线将稀疏点变成稠密点,是常用的一种规划方法. 原理讲解 源代码 #!/usr/bin/env python #-*-coding:utf-8 ...
- matlab导数曲线怎样画,matlab三次样条曲线的绘制（spline和csape函数详解）
matlab三次样条函数的绘制(spline和csape函数详解) 样条函数是工程中常用的插值函数.早期工程师制图时,把富有弹性的细长木条(所谓样条)用压铁固定在样点上,在其他地方让它自由弯曲,然后沿 ...
- 三次样条曲线CubicSpline
本文参考老张在上海轨迹规划 之 三次样条曲线(概念+性质) - 知乎 (zhihu.com) 什么是三次样条曲线 之 三次 样条是一种数据插值的方式,在多项式插值中,多项式是给出的单一公式来尽可能满足 ...
- 基于曲线插值的规划方法(Interpolating Curve Planners)
在二维运动规划中,曲线扮演着重要的角色. 对于给定一系列路点 (way-points),计算机辅助几何技术 (Computer Aided Geometric Design,CAGD) 在路径平滑方面 ...
- matlab三次样条曲线的绘制（spline和csape函数详解）
matlab三次样条函数的绘制(spline和csape函数详解) 前言 1.spline函数详解 1.一维非节点边界 2.第二边界条件 3.高维无约束 4.高维第二边界 5.利用第二边界条件绘制圆 ...
- matlab 里catmull rom,Unity中的曲线插值CatmullRom
ps:博客从typecho换成jekyll后,文章复制到简书来因为图片链接原因变得麻烦了.- -! 正文 之前写了个插件,有个需要曲线插值的功能.给定一些点的位置,物体成一条平滑曲线依次通过这些点. ...
最新文章
- nginx负载均衡和lvs负载均衡的比较分析
- php中系统函数的特征,php 常用的系统函数
- inline-block的兼容性问题
- Ubuntu增加Swap分区大小
- windows 和 ubuntu服务器之间用Xshell互传文件
- MVCC(Multiversion concurrency control)
- 用idea创建vue项目
- LTE下行物理层传输机制（9）-集中式和分布式资源映射
- linux配置vsftp红帽子,linux红帽子VSFTPD的配置.doc
- （第八天）记忆系统训练软件3.0
- vue组件库(Element UI)
- 什么是微信公众平台？
- 页面设计如何进行颜色搭配
- 实习面试感悟-阿里云
- 将JBoss启动做成Windows的系统服务
- 依赖性检查工具Snyk与Dependency-check对比
- 卷积神经网络CNN——使用keras识别猫咪
- 解决群晖 “由于系统可用存储空间不足,您将无法登录“ 的问题
- JavaScript 数据结构与算法（二）哈希表
- 查找mac系统下的隐藏文件以及隐藏文件夹的方法