三次样条插值求拟合曲线 python 三次样条插值函数程序

1.简介

三次样条插值（Cubic Spline Interpolation）简称Spline插值，是通过一系列形值点的一条光滑曲线，数学上通过求解三弯矩方程组得出曲线函数组的过程。

源代码里阐述了所有的计算公式及其流程，在这里讲述的是整体的设计思想。
利用已知数据计算H[k],再计算λ和μ，利用追赶法求解矩阵M，结合第二边界条件，根据S(x)函数求解公式，构建函数S(x)，根据已知x值求解函数值，最后利用VS调用matlab中的函数，进行曲线和求解点的绘制。

3.流程图

4.计算公式及源代码

在此公布***所有计算公式和部分源代码***，所有源代码请见最下方的下载链接

（1）计算公式

1. 计算参数Hk的值

计算公式
	  Hk[k] = x[k+1] - x[k]

	  注:x[k] = Data[0][k]

2.计算λ和μ的值

计算公式
	   λ[k] = Hk[i - 1] / (Hk[i] + Hk[i - 1])
	   μ[k] = 1 - λ[k]

	   注:λ = Lu[0][i],μ = Lu[1][i]

3.计算d[i]值（求解M等式右侧矩阵Ck ）

计算公式
		已知
		y[0]'' = 0;
		g[0] = 2y[0]'';
		y[n]'' = 0;
		g[n] = 2y[0]'';
	    g[i] = (6 / (Hk[i] + Hk[i - 1]))*(\
			   ((y[i + 1] - y[i]) / Hk[i])\
			   - (y[i] - y[i - 1]) / Hk[i - 1]);   i = 1,2,...(n-1)

	    注:g[i] = Ck[i];
		   y[i] = Data[1][i]

4.计算M矩阵值（追赶法）

求解M过程
	   1.构造M左侧原始矩阵,包含λ,μ和2
	   2.利用追赶法求解M
	      由于M[0]和M[n]已知,只需求解M[1] - M[n-1]即可，M左侧为17*17矩阵
	      计算公式
		  (1)消元过程
		  u[0] = b[0],y[0] = d[0];
		  对于i = 1,2,3,...,n
		  l[i] = a[i]/u[i-1]
		  u[i] = b[i] - l[i]c[i-1]
		  y[i] = d[i] - l[i]y[i-1]
		  (2)回代过程
		  x[n] = y[n]/u[n]
		  对于i = n-1,...,1,0
		  x[i] = (y[i] - c[i]x[i+1])/u[i]
		  最终求出参数M
	   注：
	      l[i] = CLu[0][i]
	      u[i] = CLu[1][i]
		  a[i] = MatSource[i][i - 1]
	      b[i] = MatSource[i][i]
		  d[i] = Ck[i-1]
		  x[i] = Data[0][i]
		  y[i] = Data[1][i]

5.构造S(x)函数

计算公式
	   S(x) = -((powf(x - Data[0][i+1],3.0f)/(6*Hk[i]))*Mk[i])\
	   + (powf(x - Data[0][i], 3.0f) / (6 * Hk[i])*Mk[i + 1])\
	   - (Data[1][i] - ((Hk[i]*Hk[i]/6)*Mk[i])) * ((x - Data[0][i+1])/Hk[i])\
	   + (Data[1][i+1] - (Hk[i]*Hk[i]* Mk[i+1])/6)* ((x - Data[0][i]) / Hk[i]);

（2）部分源代码

#include <iostream>
#include <stdio.h>
#include <math.h>             //数学相关函数库
#include <iomanip>            //指定格式库
#include <cmath>
#include <cstdlib>            //matlab相关函数库
#include <cstdio>
#include <cstring>
#include "engine.h"

using namespace std;

//调用matlab绘图函数所需参数
const int BUFFER_SIZE = 1024;
char buffer[BUFFER_SIZE];

class Spline {
public:
	Spline(void);             //构造函数
	~Spline() {};             //析构函数
	void HkCal(void);         //计算参数hk
	void LuCal(void);         //计算参数L和u
	void CkCal(void);         //计算参数Ck
	void ChasingMat(void);    //追赶法求解M的值
	double SFcn(float x);     //构造S(x)函数
	void DrawSpline(void);    //画出图像线段
private:
	int n;                             //数据的组数
	float Data[2][19];                 //存储数据x = Data[0][i],y = Data[1][i]
	float Datax[4];                    //存储所求点x的值
	float y_1[19];                     //存储y的一阶导
	float Hk[19];                      //存储参数Hk
	float Lu[2][19];                   //存储λ和μ λ = Lu[0][i],μ = Lu[1][i]
	float Ck[19];                      //储存参数Ck为求解M等式右侧矩阵
	float Mk[19];                      //参数M

};

Spline::Spline(void) {                 //构造函数
	//存储参数列表
	n = 18;                            //共有18 + 1组数据
	float data[2][19] = { { 0.520,    3.1,     8.0,     17.95,   28.65,   39.62,   50.65,   78,    104.6,   156.6,\
		                    208.6,    260.7,   312.5,   364.4,   416.3,   468,     494,     507,   520 },\
	                       {5.288,    9.4,     13.84,   20.20,   24.90,   28.44,   31.10,   35,    36.9,    36.6,\
		                    34.6,     31.0,    26.34,   20.9,    14.8,    7.8,     3.7,     1.5,   0.2} };
	float datax[4] = { 2.30,130,350,515 };

	//进行数组的复制
	memcpy(Data, data, sizeof(Data)); 
	memcpy(Datax, datax, sizeof(Datax));

	//对已知参数进行赋值
	y_1[0] = 1.86548;
	y_1[n] = -0.046115;
}


void Spline::HkCal(void) {                            //计算参数hk
	/*
	  计算公式
	  Hk[k] = x[k+1] - x[k]

	  注:x[k] = Data[0][k]
	*/
		for (int i = 0; i < n; i++) {
			Hk[i] = Data[0][i + 1] - Data[0][i];
			cout <<"H["<< i << "]:"<< Hk[i] << endl;
		}
		cout << endl;
}

void Spline::LuCal(void) {             	//计算L和u
	/*
	计算公式
	λ[k] = Hk[i - 1] / (Hk[i] + Hk[i - 1])
	μ[k] = 1 - λ[k]

	注:λ = Lu[0][i],μ = Lu[1][i]
	*/
	//遍历计算
	cout.setf(std::ios::left);             //所有数据左对齐

	for (int i = 1; i < n; i++) {
		Lu[0][i] = Hk[i - 1] / (Hk[i] + Hk[i - 1]);
		Lu[1][i] = 1 - Lu[0][i];
		cout << "λ["<< setw(3) << i <<"] = " << setw(15) << Lu[0][i] << "   " << "μ[" << setw(3) << i << "]" << Lu[1][i] << endl;
	 }

	cout << endl;
}

void Spline::CkCal(void) {                              //求解M等式右侧矩阵Ck  
	/*
	    计算公式
		已知
		y[0]'' = 0;
		g[0] = 2y[0]'';
		y[n]'' = 0;
		g[n] = 2y[0]'';
	    g[i] = (6 / (Hk[i] + Hk[i - 1]))*(\
			   ((y[i + 1] - y[i]) / Hk[i])\
			   - (y[i] - y[i - 1]) / Hk[i - 1]);   i = 1,2,...(n-1)

	    注:g[i] = Ck[i];
		   y[i] = Data[1][i]
	*/
	Ck[0] = 0;
	cout << "Ck[" << 0 << "]" << ":" << Ck[0] << endl;
	for (int i = 1; i < n; i++) {
		Ck[i] = (6 / (Hk[i] + Hk[i - 1]))*(\
			((Data[1][i + 1] - Data[1][i]) / Hk[i])\
			- (Data[1][i] - Data[1][i - 1]) / Hk[i - 1]);
		cout << "Ck[" << i << "]:" << Ck[i] << endl;
	}
	Ck[n] = 0;
	cout << "Ck[" << n << "]:" << Ck[n] << endl;
	cout << endl;
}

5.运行结果

数据信息

调用matlab绘出图像