(Spline)三次样条求解
三次样条(Spline)算法在当前的主流数据计算的方法选择时,是一种常用的方法,下面这种算法是在晚上找到的一种算法,
一、样条函数的定义
样条函数属于分段光滑插值,他的基本思想是,在由两相邻节点所构成的每一个小区间内用低次多项式来逼近,并且在各结点的连接处又保证是光滑的(即导数连续)。
设在区间[a,b]上给定一组结点X:
,和一组对应的函数值。若函数S(x)满足下列条件:
(1)在每一个子区间(k=1,2,...n)上,S(x)是一个不超过三次的多项式。
(2)在每一个结点上满足S(xi)=yi,i=0,1,...,n。
(3)S(x)在区间[a,b]上为二次连续可微函数。
则称S(x)为结点X上插值与Y的三次样条插值。
二、三次样条函数的构造
在工程上,构造三次样条插值函数通常有两种方法:一是以给定插值结点处得二阶导数值作为未知数来求解,而工程上称二阶导数为弯矩,因此,这种方法成为三弯矩插值。二是以给定插值结点处得一阶导数作为未知数来求解,而一阶导数右称为斜率,因此,这种方法称为三斜率插值。
三斜率插值法
根据定义,三次样条函数在插值结点处一阶导数应存在。因此设各结点处的一阶导数为:
。利用两点埃尔米特插值公式,就可以得到样条插值函数S(x)在子区间上的表达式为:
(1)
其中:。由此可知只要确定各结点的一阶导数值mk(k=0,1,2....n),则各子区间上的三次样条插值函数S(x)也就确定了。
由于S(x)在区间[a,b]上的二阶导数是连续的,即在各结点的左右两子区间上的S(x)虽然不同,但在连接点的二阶导数存在,即在连接点处的二阶左导数与二阶右导数相等:。分别求子区间[xk-1,xk]右端点xk上的二阶左导数,以及子区间[xk,xk+1]左端点xk上的右导数,即可得到:
,(2)与 ,(3) 整理后可得:,k=1,2,...n-1,(4)。
其中,并且ak+bk=1。由于有n+1个插值点,因此需要确定n+1个m(m0,m1,...,mn)。而现在方程组只有n-1个方程(公式(4)),因此还需要另外补充两个条件财能唯一确定n+1个m。在实际应用中,这两个条件可以根据实际的物理意义所给出的边界条件来得到。实际问题中常用的三种边界条件如下:
(1)给定区间两个端点处的一阶导数值,即
则根据方程组(4)可得到新的方程组:
方程组的系数矩阵为,可知此矩阵为三对角矩阵。因此方程组可用追赶法求解。
解出mk后,即得到各结点的一阶导数值,将mk带入各结点的二阶导数值得表达式公式(2)或(3)可求得各结点的二阶导数值,将mk带入各区间上S(x)的表达式公式(1)即可得到各子区间上的三次样条插值函数S(x)。
(2)给定区间两个端点处得二阶到数值,即
,由此可得两个补充方程,其中。与公式(4)联立可得如下方程组:
(3)插值函数为周期函数
yn=y0,且,由此可得两个补充方程为,其中
,并且an+bn=1。最后可得方程组:
在此给出第一种边界条件下三次样条插值的C#算法实现:
1、实现三次样条插值封装的类:
/// <summary>
/// 插值
/// </summary>
public static class SplineMath
{
/// <summary>
/// 三次样条插值
/// </summary>
/// <param name="points">排序好的数</param>
/// <param name="xs">需要计算的插值点</param>
/// <param name="chf">写1</param>
/// <returns>返回计算好的数值</returns>
public static double[] SplineInsertPoint(PointClass[] points, double[] xs, int chf)
{
int plength = points.Length;
double[] h = new double[plength];
double[] f = new double[plength];
double[] l = new double[plength];
double[] v = new double[plength];
double[] g = new double[plength];
for (int i = 0; i < plength - 1; i++)
{
h[i] = points[i + 1].x - points[i].x;
f[i] = (points[i + 1].y - points[i].y) / h[i];
}
for (int i = 1; i < plength - 1; i++)
{
l[i] = h[i] / (h[i - 1] + h[i]);
v[i] = h[i - 1] / (h[i - 1] + h[i]);
g[i] = 3 * (l[i] * f[i - 1] + v[i] * f[i]);
}
double[] b = new double[plength];
double[] tem = new double[plength];
double[] m = new double[plength];
double f0 = (points[0].y - points[1].y) / (points[0].x - points[1].x);
double fn = (points[plength - 1].y - points[plength - 2].y) / (points[plength - 1].x - points[plength - 2].x);
b[1] = v[1] / 2;
for (int i = 2; i < plength - 2; i++)
{
// Console.Write(" " + i);
b[i] = v[i] / (2 - b[i - 1] * l[i]);
}
tem[1] = g[1] / 2;
for (int i = 2; i < plength - 1; i++)
{
//Console.Write(" " + i);
tem[i] = (g[i] - l[i] * tem[i - 1]) / (2 - l[i] * b[i - 1]);
}
m[plength - 2] = tem[plength - 2];
for (int i = plength - 3; i > 0; i--)
{
//Console.Write(" " + i);
m[i] = tem[i] - b[i] * m[i + 1];
}
m[0] = 3 * f[0] / 2.0;
m[plength - 1] = fn;
int xlength = xs.Length;
double[] insertRes = new double[xlength];
for (int i = 0; i < xlength; i++)
{
int j = 0;
for (j = 0; j < plength; j++)
{
if (xs[i] < points[j].x)
break;
}
j = j - 1;
Console.WriteLine(j);
if (j == -1 || j == points.Length - 1)
{
if (j == -1)
throw new Exception("插值下边界超出");
if (j == points.Length - 1 && xs[i] == points[j].x)
insertRes[i] = points[j].y;
else
throw new Exception("插值下边界超出");
}
else
{
double p1;
p1 = (xs[i] - points[j + 1].x) / (points[j].x - points[j + 1].x);
p1 = p1 * p1;
double p2; p2 = (xs[i] - points[j].x) / (points[j + 1].x - points[j].x);
p2 = p2 * p2;
double p3; p3 = p1 * (1 + 2 * (xs[i] - points[j].x) / (points[j + 1].x - points[j].x)) * points[j].y + p2 * (1 + 2 * (xs[i] - points[j + 1].x) / (points[j].x - points[j + 1].x)) * points[j + 1].y;
double p4; p4 = p1 * (xs[i] - points[j].x) * m[j] + p2 * (xs[i] - points[j + 1].x) * m[j + 1];
// Console.WriteLine(m[j] + " " + m[j + 1] + " " + j);
p4 = p4 + p3;
insertRes[i] = p4;
//Console.WriteLine("f(" + xs[i] + ")= " + p4);
}
}
//Console.ReadLine();
return insertRes;
}
}
2、插值前需要对数据进行排序,需要使用PointClass类:
public class PointClass
{
public double x = 0;
public double y = 0;
public PointClass()
{
x = 0; y = 0;
}
//-------写一个排序函数,使得输入的点按顺序排列,是因为插值算法的要求是,x轴递增有序的---------
public static PointClass[] DeSortX(PointClass[] points)
{
int length = points.Length;
double temx, temy;
for (int i = 0; i < length - 1; i++)
{
for (int j = 0; j < length - i - 1; j++)
if (points[j].x > points[j + 1].x)
{
temx = points[j + 1].x;
points[j + 1].x = points[j].x;
points[j].x = temx;
temy = points[j + 1].y;
points[j + 1].y = points[j].y;
points[j].y = temy;
}
}
return points;
}
}
3、具体实现:
private void btnCalcSpline_Click(object sender, EventArgs e)
{
double[] x = {-100,-90,-80,-70,-60,-50,-40,-30,-20,-10,0,10 ,20,30,40,50,60,70,80,90,100};
double[] y = {9802,7922,6242,4762,3482,2402,1522,842,362,82,2,122,442,962,1682,2602,3722,5142,6562,8282,10202};
PointClass[] points = new PointClass[x.Length];
for (int i=0;i<x.Length;i++)
{
points[i] = new PointClass();
points[i].x = x[i];
points[i].y = y[i];
}
PointClass.DeSortX(points); //排序
//double[] xs = { 0,1,80 };
double[] xs = { Convert.ToDouble(txtX.Text) };
try
{
double[] Y = SplineMath.SplineInsertPoint(points, xs, 1);
txtY.Text = Y[0].ToString();
}
catch(Exception ex)
{
MessageBox.Show(ex.Message);
}
}
(Spline)三次样条求解相关推荐
- Spline(三次样条插值)
关于三次样条插值,计算方法比较复杂,但是静下心来仔细研究也是可以理解的. 本文借鉴文章来源:http://www.cnki.com.cn/Article/CJFDTotal-BGZD200611035 ...
- matlab spline三次样条插值x,Spline(三次样条插值)
关于三次样条插值,计算方法比较复杂,但是静下心来仔细研究也是可以理解的. 本文借鉴文章来源:http://www.cnki.com.cn/Article/CJFDTotal-BGZD200611035 ...
- aitken插值方法的c++代码_无人驾驶路径规划技术-三次样条插值曲线及Python代码实现...
自动驾驶运动规划(Motion Planning)是无人驾驶汽车的核心模块之一,它的主要任务之一就是如何生成舒适的.碰撞避免的行驶路径和舒适的运动速度.生成行驶路径最经典方法之一就是是Sampling ...
- 滑动轨迹 曲线 python_无人驾驶路径规划技术-三次样条插值曲线及Python代码实现...
自动驾驶运动规划(Motion Planning)是无人驾驶汽车的核心模块之一,它的主要任务之一就是如何生成舒适的.碰撞避免的行驶路径和舒适的运动速度.生成行驶路径最经典方法之一就是是Sampling ...
- bezier + spline
中国农业大学-赵明老师-计算机图形学-moochttps://www.bilibili.com/video/BV1Dt411f7Qj?p=11 目录 bezier bezier 曲线定义 bernst ...
- Spline算法实现
0 前言 本文总结不同工具对Spline算法的实现情况. 1 Matlab 1.1 spline() 具体案例详见<Spline导数及曲率计算>.<函数逼近方法>第1.3小节. ...
- matlab中如何画三次样条,Matlab之三次样条画图和表达式
这一题是得到数据点(0,3),(1,5),(2,4),(3,1)并得到它的三次样条表达式和画出三次样条后的图图形. 以及对数据点(-1,3),(0,5),(3,1),(4,1),(5,1)并得到它的三 ...
- Matlab参考函数
附录1 常用命令 附录1.1 管理用命令 函数名 功能描述 函数名 功能描述 addpath 增加一条搜索路径 rmpath 删除一条搜索路径 demo 运行Matlab演示程序 type 列出.M文 ...
- Matlab命令汇总
Matlab命令汇总 都是从网上转的,贴到一起方便查,使用的时候直接Ctrl+F搜索. 一.常用对象操作:除了一般windows窗口的常用功能键外. 1.!dir 可以查看当前工作目录的文件. ! ...
- matlab库函数大全
附录 MATLAB函数参考 附录1 常用命令 附录1.1 管理用命令 函数名 功能描述 函数名 功能描述 addpath 增加一条搜索路径 rmpath 删除一条搜索路径 demo 运行Matlab演 ...
最新文章
- angular6 mysql_angular6之路由
- 数据挖掘导论读书笔记9聚类分析
- 有关C++多态的一些测试
- CUDA从入门到精通(三):必备资料
- 导致集群重启_解析 Elasticsearch 棘手问题,集群的 RED 与 YELLOW
- 3520a新板做内存测试
- 线程属性 pthread_attr_t
- (idea)设置鼠标移到类、方法、变量上时,显示相关提示信息
- FD.io VPP配置文件详解
- javascript中的XML
- atitit 编程语言概念与原理
- 桌面云之深信服VMP管理
- 分析CHE矢量变频器在数控雕刻机床上应用
- 商品归类查询服务_喜报 | 东泽国际获批进出口商品归类服务单位资质
- 菱形的常见图案_菱形图案,简约而不简单
- 淘宝SDK高级模板,设计师模块开放接口详解
- python爬虫豆瓣电影按电影类型_Python爬虫入门 | 7 分类爬取豆瓣电影,解决动态加载问题...
- 计算机任务管理器无法响应,Win7系统电脑在任务管理器中关闭进程时总是未响应的解决方法...
- 基于 Web 的 Java Swing Kiosk 应用程序
- test题目:袋鼠过河
热门文章
- CSP 2014-03-1 相反数(C++)
- Centos7下安装Relion
- java 判断当前时间节气,请问js获取阴历节气后根据节气判断春夏秋冬四季?该怎么写呢...
- python正则表达式与re模块
- [转载]千古真人张三丰
- Udacity 传感器融合笔记 (一)lidar
- 音乐推荐,持续收集中
- 联想thinkpad E430C硬盘位换为固态,硬盘放于光驱位(win7+win10+ubuntu三系统安装教程)
- 王者荣耀在android目录下的名字,王者荣耀手q区有哪些 王者荣耀安卓手Q区名称...
- 广东科技学院计算机学院院长,陈强-广东科技学院-计算机学院