三次样条插值

水下修卫星

已于 2024-03-24 20:05:42 修改

阅读量1k

点赞数 16

分类专栏：数学 C++ 文章标签：线性代数矩阵

于 2024-03-24 15:22:36 首次发布

本文链接：https://blog.csdn.net/qq_26860485/article/details/136921118

版权

C++ 同时被 2 个专栏收录

9 篇文章 0 订阅

订阅专栏

数学

4 篇文章 2 订阅

订阅专栏

三次样条

三次样条曲线拟合主要是为解决三次样条函数不能解决的问题而提出的。三次样条函数要求自变量满足单调递增，为解决此限制条件，三次样条拟合曲线通常采用参数方程的形式表示。

1、原理

三次样条插值法是一种常用的数值分析方法，可以通过给定的一组散点数据来拟合出一条光滑的连续函数曲线。其基本思想是用低次多项式逼近一段小区间内的数据，并利用这些多项式的连接处衔接条件来保证整个曲线的光滑性。
定义：对于已知的 $n$ 个数据点 $x_{i},y_{i})$ ，其中 $x_{i}<x_{i+1}$ ，三次样条插值法可以找到一组函数 $S (x)$ ，并满足以下条件：
1、 $S (x)$ 在每个小区间上均为三次多项式，即：
$S_{i}(x_{i})=a_{i}+b_{i}(x-x_{i})+c_{i}(x-x_{i})^2+d_{i}(x-x_{i})^3$
2、 $S (x)$ 经过所有给定的数据点，即：
$S_{1}(x_{1})=y_{1},S_{n-1}(x_{n})=y_{n},S_{i}(x_{i+1})=S_{i+1}(x_{i+1})=y_{i+1},i=1,2,...,n-2$ ；
3、 $S (x)$ 在两个相邻的小区间的连接处拥有相同的一阶、二阶导数，即：
$S_{i}'(x_{i+1})=S_{i+1}'(x_{i+1}),S_{i}''(x_{i+1})=S_{i+1}''(x_{i+1}),i=1,2,...,n-2$
根据以上条件，可以得到以下等式组：
$\begin{aligned}&S_{1}(x_{1})=y_{1}\\ S_{1}(x_{2})=&S_{2}(x_{2})=y_{2}\\ S_{2}(x_{3})=&S_{3}(x_{3})=y_{3}\\ &...\\ S_{n-2}(x_{n-1})=&S_{n-1}(x_{n-1})=y_{n-1}\\ S_{n-1}(x_{n})=&y_{n}\\ S'_{1}(x_{2})&=S'_{2}(x_{2})\\ S'_{2}(x_{3})&=S'_{3}(x_{3})\\ &...\\ S'_{n-2}(x_{n-1})&=S'_{n-1}(x_{n-1})\\ S''_{1}(x_{2})&=S''_{2}(x_{2})\\ S''_{2}(x_{3})&=S''_{3}(x_{3})\\ &...\\ S''_{n-2}(x_{n-1})&=S''_{n-1}(x_{n-1}) \end{aligned}$
第一型边界条件：
$S''_{1}(x_{1})=acc_{1}\\S''_{n-1}(x_{n})=acc_{n-1}$
当 $acc_{1}=acc_{n-1}=0$ 时，称为自然边界条件
第二型边界条件：
$S'_{1}(x_{1})=vel_{1}\\S'_{n-1}(x_{n})=vel_{n-1}$
自然边界条件下的数学推导：
由 $\begin{aligned}S_{1}(x_{1})&=y_{1}\\ S_{2}(x_{2})&=y_{2}\\ S_{3}(x_{3})&=y_{3}\\ &...\\ S_{n-2}(x_{n-2})&=y_{n-2}\\ S_{n-1}(x_{n-1})&=y_{n-1}\\ \end{aligned}$
可得：
$\begin{aligned} a_{i}&=y_{i}，i=1,2,...,n-1 \end{aligned}$
由
$\begin{aligned}S''_{1}(x_{2})&=S''_{2}(x_{2})\\ S''_{2}(x_{3})&=S''_{3}(x_{3})\\ &...\\ S''_{n-2}(x_{n-1})&=S''_{n-1}(x_{n-1})\\ S''_{n-1}(x_{n})&=0 \end{aligned}$
可得：
$\begin{aligned} d_{i}&=\frac{c_{i+1}-c_{i}}{3h_{i}},i=1,2,...,n-2\\ d_{n-1}&=\frac{-c_{n-1}}{3h_{n-1}},i=n-1 \end{aligned}$
由 $\begin{aligned}S_{1}(x_{2})&=y_{2}\\ S_{2}(x_{3})&=y_{3}\\ &...\\ S_{n-2}(x_{n-1})&=y_{n-1}\\ S_{n-1}(x_{n})&=y_{n}\\ \end{aligned}$
可得：
$\begin{aligned} b_{i}&=\frac{y_{i+1}-y_{i}}{h_{i}}-\frac{c_{i+1}+2c_{i}}{3}h_{i},i=1,2,...,n-2\\ b_{n-1}&=\frac{y_{n}-y_{n-1}}{h_{n-1}}-\frac{2c_{n-1}}{3}h_{n-1},i=n-1 \end{aligned}$
由
$\begin{aligned}S''_{1}(x_{1})&=0\\ S'_{1}(x_{2})&=S'_{2}(x_{2})\\ S'_{2}(x_{3})&=S'_{3}(x_{3})\\ &...\\ S'_{n-2}(x_{n-1})&=S'_{n-1}(x_{n-1})\\ S''_{n-1}(x_{n})&=0 \end{aligned}$
可得：
$\begin{aligned}c_{1}&=0\\ h_{i}c_{i}+2(h_{i}+h_{i+1})c_{i+1}+h_{i+1}c_{i+2}&=3(\frac{y_{i+2}-y_{i+1}}{h_{i+1}}-\frac{y_{i+1}-y_{i}}{h_{i}}),i=1,2,...,n-2\\ c_{n}&=0 \end{aligned}$
将上式写成矩阵 $\boldsymbol{AX}=\boldsymbol{B}$ 形式：
$\boldsymbol{A}=\begin{bmatrix}1&0&0&0&...&0\\ h_{1}&2(h_{1}+h_{2})&h_{2}&0&...&0\\ 0&h_{2}&2(h_{2}+h_{3})&h_{3}&...&0\\ \vdots&\vdots&\vdots&\vdots&\vdots&\vdots\\ 0&...&0&h_{n-2}&2(h_{n-2}+h_{n-1})&h_{n-1}\\ 0&...&0&0&0&1 \end{bmatrix}$
$\boldsymbol{X}=\begin{bmatrix}c_{1}\\ c_{2}\\ c_{3}\\ \vdots\\ c_{n-1}\\ c_{n} \end{bmatrix} \boldsymbol{B}=\begin{bmatrix}0\\ 3(\frac{y_{3}-y_{2}}{h_{2}}-\frac{y_{2}-y_{1}}{h_{1}})\\ 3(\frac{y_{4}-y_{3}}{h_{3}}-\frac{y_{3}-y_{2}}{h_{2}})\\ \vdots\\ 3(\frac{y_{n}-y_{n-1}}{h_{n-1}}-\frac{y_{n-1}-y_{n-2}}{h_{n-2}})\\ 0 \end{bmatrix}$

2、追赶法求方程组

对于系数矩阵是三对角矩阵的特殊方程组，若用原有的一般方法来求解，势必造成存储和计算的浪费，可以采用“追赶法”来求解。
$\begin{bmatrix}b_{1}&a_{1}&0&0&...&0\\ c_{2}&b_{2}&a_{2}&0&...&0\\ 0&c_{3}&b_{3}&a_{3}&...&0\\ \vdots&\vdots&\vdots&\vdots&\vdots&\vdots\\ 0&...&0&c_{n-1}&b_{n-1}&a_{n-1}\\ 0&...&0&0&c_{n}&b_{n} \end{bmatrix} \begin{bmatrix}x_{1}\\ x_{2}\\ x_{3}\\ \vdots\\ x_{n-1}\\ x_{n} \end{bmatrix}= \begin{bmatrix}d_{1}\\ d_{2}\\ d_{3}\\ \vdots\\ d_{n-1}\\ d_{n} \end{bmatrix}$
将上述三角阵进行 $LU$ 分解可得：
$\begin{bmatrix}b_{1}&a_{1}&0&0&...&0\\ c_{2}&b_{2}&a_{2}&0&...&0\\ 0&c_{3}&b_{3}&a_{3}&...&0\\ \vdots&\vdots&\vdots&\vdots&\vdots&\vdots\\ 0&...&0&c_{n-1}&b_{n-1}&a_{n-1}\\ 0&...&0&0&c_{n}&b_{n} \end{bmatrix}= \begin{bmatrix}l_{1}&0&0&0&...&0\\ r_{2}&l_{2}&0&0&...&0\\ 0&r_{3}&l_{3}&0&...&0\\ \vdots&\vdots&\vdots&\vdots&\vdots&\vdots\\ 0&...&0&r_{n-1}&l_{n-1}&0\\ 0&...&0&0&r_{n}&l_{n} \end{bmatrix} \begin{bmatrix}1&u_{1}&0&0&...&0\\ 0&1&u_{2}&0&...&0\\ 0&0&1&u_{3}&...&0\\ \vdots&\vdots&\vdots&\vdots&\vdots&\vdots\\ 0&...&0&0&1&u_{n-1}\\ 0&...&0&0&0&1 \end{bmatrix}$
由上式可得：
当 $i = 1$ 时：
$\begin{aligned}l_{1}&=b_{1}\\u_{1}&=a_{1}/b_{1}\end{aligned}$
当 $i = 2, 3, ..., n$ 时：
$\begin{aligned}l_{i}&=b_{i}-r_{i}u_{i-1}\\r_{i}&=c_{i}\\u_{i}&=a_{i}/l_{i}\end{aligned}$
即：
$\begin{aligned}l_{i}&=b_{i}-c_{i}a_{i-1}/l_{i-1}\\r_{i}&=c_{i}\\u_{i}&=a_{i}/(b_{i}-c_{i}u_{i-1})\end{aligned}$
令：
$\begin{bmatrix}1&u_{1}&0&0&...&0\\ 0&1&u_{2}&0&...&0\\ 0&0&1&u_{3}&...&0\\ \vdots&\vdots&\vdots&\vdots&\vdots&\vdots\\ 0&...&0&0&1&u_{n-1}\\ 0&...&0&0&0&1 \end{bmatrix}\begin{bmatrix}x_{1}\\ x_{2}\\ x_{3}\\ \vdots\\ x_{n-1}\\ x_{n} \end{bmatrix}= \begin{bmatrix}y_{1}\\ y_{2}\\ y_{3}\\ \vdots\\ y_{n-1}\\ y_{n} \end{bmatrix}$
则：
$\begin{bmatrix}l_{1}&0&0&0&...&0\\ r_{2}&l_{2}&0&0&...&0\\ 0&r_{3}&l_{3}&0&...&0\\ \vdots&\vdots&\vdots&\vdots&\vdots&\vdots\\ 0&...&0&r_{n-1}&l_{n-1}&0\\ 0&...&0&0&r_{n}&l_{n} \end{bmatrix} \begin{bmatrix}y_{1}\\ y_{2}\\ y_{3}\\ \vdots\\ y_{n-1}\\ y_{n} \end{bmatrix}= \begin{bmatrix}d_{1}\\ d_{2}\\ d_{3}\\ \vdots\\ d_{n-1}\\ d_{n} \end{bmatrix}$
由 $L y = d$ 可得：
当 $i = 1$ 时：
$y_{1}=d_{1}$
当 $i = 2, 3, ..., n$ 时：
$\begin{aligned}y_{i}&=(d_{i}-c_{i}y_{i-1})/l_{i}\\ &=(d_{i}-c_{i}y_{i-1})/(b_{i}-c_{i}u_{i-1})\end{aligned}$
再由 $Ux = y$ 可得：
当 $i = n$ 时：
$x_{n}=y_{n}$
当 $i = n - 1, n - 2, ..., 1$ 时：
$\begin{aligned}x_{i}&=y_{i}-u_{i}x_{i+1}\end{aligned}$

#ifndef _CUBIC_SPLINE_H
#define _CUBIC_SPLINE_H

#include <iostream>

class CSpline
{
public:
	CSpline();
	~CSpline();

	void Splint(double *xa, double *ya, double *m, int n, double &x, double &y);

	void Spline(double *xa, double *ya, int n, double *&m);
};

#endif

#include "math.h"
#include "cubic_spline.h"

CSpline::CSpline() {};

CSpline::~CSpline() {};

//================================================================
// 函数功能： 利用求出的二阶偏导数，求给定点值处的函数值
// 输入参数： *xa 为横坐标值，ya为纵坐标值，n为点个数，m为二阶偏导数
//            x为给定点，y接收插出来的值
// 返回值：   无返回值
void CSpline::Splint(double *xa, double *ya, double *m, int n, double &x, double &y)
{
	int klo, khi, k;
	klo = 0; khi = n - 1;
	double a, b, c, d, h, deltaX;

	while (khi - klo > 1)            //  二分法查找x所在区间段
	{
		k = (khi + klo) >> 1;
		if (xa[k] > x)  khi = k;
		else klo = k;
	}
	h = fabs(xa[khi] - xa[klo]);
	a = ya[klo];
	c = m[klo];
	if (klo == n - 2) {
		b = (ya[khi] - ya[klo]) / h - (2 * m[klo])*h / 3;
		d = -m[klo] / (3 * h);
	}
	else {
		b = (ya[khi] - ya[klo]) / h - (m[khi] + 2 * m[klo]) * h / 3;
		d = (m[khi] - m[klo]) / (3 * h);
	}
	deltaX = x - xa[klo];
	y = a + b * deltaX + c * pow(deltaX, 2) + d * pow(deltaX, 3);
}

//===========================================================================
// 函数功能： 根据自然边界条件对一系列点求二阶偏导数，点横坐标单调递增
// 输入参数： *xa 为横坐标值，ya为纵坐标值，n为点个数，m为二阶偏导数(输出值）
// 返回值：   无返回值
void CSpline::Spline(double *xa, double *ya, int n, double *&m)
{
	double *a = new double[n];			//  a:稀疏矩阵最上边一串数
	double *b = new double[n];			//  b:稀疏矩阵最中间一串数
	double *c = new double[n];			//  c:稀疏矩阵最下边一串数
	double *d = new double[n];			//  d:结果数据
	double *h = new double[n - 1];			//  h:步长

	double *l = new double[n];
	double *u = new double[n];

	m = new double[n];

	//  计算步长
	for (int i = 0; i < n - 1; i++) {
		h[i] = xa[i + 1] - xa[i];
	}
	//  计算三对角矩阵最上边的一串数据
	a[0] = 0;
	for (int i = 1; i < n - 1; i++) {
		a[i] = h[i];
	}
	a[n - 1] = 0;//  无用数据
	//  计算三对角矩阵最中间的一串数据
	b[0] = 1;
	for (int i = 1; i < n - 1; i++) {
		b[i] = 2 * (h[i - 1] + h[i]);
	}
	b[n - 1] = 1;
	//  计算三对角矩阵最下边的一串数据
	c[0] = 0;//  无用数据
	for (int i = 1; i < n - 1; i++) {
		c[i] = h[i - 1];
	}
	c[n - 1] = 0;
	//  计算L对角线上数据
	l[0] = b[0];
	for (int i = 1; i < n; i++) {
		l[i] = b[i] - c[i] * a[i - 1] / l[i - 1];
	}
	//  计算U最上边的一串数据
	for (int i = 0; i < n; i++) {
		u[i] = a[i] / l[i];
	}
	//  计算结果数据
	d[0] = 0;
	for (int i = 1; i < n; i++) {
		d[i] = 3 * ((ya[i + 1] - ya[i]) / h[i] - (ya[i] - ya[i - 1]) / h[i - 1]);
	}
	d[n - 1] = 0;
	//  更新结果数据
	d[0] = d[0];
	for (int i = 1; i < n; i++) {
		d[i] = (d[i] - c[i] * d[i - 1]) / l[i];
	}
	//  计算系数m
	m[n - 1] = 0; //  无用数据
	m[n - 2] = d[n - 2];
	for (int i = n - 3; i >= 0; i--) {
		m[i] = d[i] - u[i] * m[i + 1];
	}

	//  输出系数m
	for (int i = 0; i < n; i++) {
		std::cout << m[i] << "  ";
	}
	std::cout << std::endl;

	delete a;
	delete b;
	delete c;
	delete d;
	delete h;

	delete l;
	delete u;
}

示例：

#include <fstream>
#include "cubic_spline.h"

int main()
{
	int num = 8;
	// initialize data
	double x[8] = { 9.59,60.81,105.57,161.59,120.5,100.1,50.0,10.0 };
	double y[8] = { 61.97,107.13,56.56,105.27,120.5,150.0,110.0,180.0 };

	//  generate parameterized variable
	double *s = new double[num];
	s[0] = 0;
	double tmp_s = 0;
	for (int i = 0; i < num - 1; i++) {
		tmp_s = tmp_s + sqrt((x[i + 1] - x[i])*(x[i + 1] - x[i]) + (y[i + 1] - y[i])*(y[i + 1] - y[i]));
		s[i + 1] = tmp_s;
	}

	//  smoothed spline points
	int N = tmp_s;
	double *sd = new double[N];
	double *xd = new double[N];
	double *yd = new double[N];

	for (int i = 0; i < N - 1; i++)
		sd[i] = i;
	sd[N - 1] = tmp_s;

	// 自然边界条件
	CSpline spline;
	double *mx, *my;
	spline.Spline(s, x, num, mx);
	spline.Spline(s, y, num, my);

	for (int i = 0; i < N; i++)
	{
		spline.Splint(s, x, mx, num, sd[i], xd[i]);
		spline.Splint(s, y, my, num, sd[i], yd[i]);
	}
	
	delete sd;
	delete xd;
	delete yd;
	delete mx;
	delete my;

	return 0;
}