已知函数在区间上个互异节点,处的函数值为,若构造函数,知足:ios
则称为的三次样条插值函数。算法
根据定义知道规律为:函数
已知:spa
根据定点,求出每段样条曲线方程中的系数,便可获得每段曲线的具体表达式。.net
插值和连续性:3d
其中.code
微分连续性:blog
其中.教程
样条曲线的微分式:ci
将步长 带入样条曲线的条件:
得
至此,能够获得一个关于m为未知数的线性方程组,经过解这个方程组,肯定其系数。
边界条件有三种:自由边界条件、固定边界条件、非节点边界条件。
两端没有任何限制,即。求解方程为:
数据两端节点的微分值是已知的,令分别为A 和B。则有
则有
指定样条曲线的三次微分匹配,
由 .
则有
一、matlab
在matlab中,spline函数能够实现三次样条插值。例:
clc; clear; //x = [ 0, 3, 5, 7, 9, 11, 12, 13, 14, 15 ]; //y = [0,1.2,1.7,2.0,2.1,2.0,1.8,1.2,1.0,1.6 ]; //x_ = 0:.1:15; x = -pi:pi; y = sin(x); x_ = -pi:.1:pi; p1 = spline(x,y,x_); % 分段三次样条插值 plot(x,y,'ko',x_,p1,'b.') legend('插值节点','分段三次样条插值','location','southeast')
能够看出,由点组成得曲线比较平滑。
二、C++
程序我有用到矩阵库Armadillo,没有安装的请点这个安装教程连接
#include<iostream> #include<armadillo> #include<vector> using namespace std; using namespace arma; int main() { double x[] = { 0, 3, 5, 7, 9, 11, 12, 13, 14, 15 }; double y[] = { 0,1.2,1.7,2.0,2.1,2.0,1.8,1.2,1.0,1.6 }; int size = sizeof(x) / sizeof(double); const double min = x[0]; const double max = x[size - 1]; vector<double> xx, yy; //插值点与计算的函数值 for (double i = x[0]; i <= x[size - 1]; i = i + 0.1) { xx.push_back(i); } int size_xx = xx.size(); vector<double> dx, dy; //差分,即步长 for (int i = 0; i < size - 1; i++) { double temp_dx = x[i + 1] - x[i]; dx.push_back(temp_dx); double temp_dy = y[i + 1] - y[i]; dy.push_back(temp_dy); } mat H, Y, M; // H * M = Y H.zeros(size, size); Y.zeros(size,1); //M.zeros(1, size); H(0, 0) = 1; H(size - 1, size - 1) = 1; for (int i = 1; i < size - 1; i++) { H(i, i - 1) = dx[i - 1]; H(i, i) = 2 * (dx[i - 1] + dx[i]); H(i, i + 1) = dx[i]; Y(i) = 3 * (dy[i] / dx[i] - dy[i - 1] / dx[i - 1]); } M = solve(H,Y); vector<double> ai, bi, ci, di; //系数 for (int i = 0; i < size - 1; i++) { ai.push_back(y[i]); di.push_back((M(i + 1) - M(i)) / (3 * dx[i])); bi.push_back(dy[i] / dx[i] - dx[i] * (2 * M(i) + M(i + 1)) / 3); ci.push_back(M(i)); } vector<int> x_, xx_; for (int i = 0; i < size; i++) { int temp = x[i] / 0.1; x_.push_back(temp); } for (int i = 0; i < size_xx; i++) { int temp = xx[i] / 0.1; xx_.push_back(temp); } for (int i = 0; i < size_xx; i++) { int k = 0; for (int j = 0; j < size - 1; j++) { if (xx_[i] >= x_[j] && xx_[i] < x_[j + 1]) { k = j; break; } else if (xx[i] == x[size - 1]) { k = size - 1; } } //yy(i) = y[i] + bi(k) * (xx[i] - x[k]) + 1 / 2.0 * M(i) * pow((xx[i] - x[k]) , 2) + di(k) * pow((xx[i] - x[k]),3); double temp = ai[k] + bi[k] * (xx[i] - x[k]) + M(k) * pow((xx[i] - x[k]), 2) + di[k] * pow((xx[i] - x[k]), 3); yy.push_back(temp); } std::ofstream output; output.open("Spline.txt"); for (unsigned i = 0; i < size_xx; i++) { output << xx[i] << '\t' << yy[i] << std::endl; } output.close(); cout << "Hello World !" << endl; }
结果用matlab显示为:
一、算法流程
二、难点
三次样条插值的总体推导,以及 H 矩阵的填充。