python_matlab_样条插值求未知曲线的函数解析式

一、问题引入

    对于给出如下的离散的数据点，现在想根据如下的数据点来推测x=5时的值，我们应该采用什么方法呢？

用于拟合样条函数的数据：
x          f ( x)
3.0 2.5
4.5 1.0
7.0 2.5
9.0 0.5
   我们知道在平面上两个点确定一条直线，三个点确定一条抛物线（假设曲线的类型是抛物线），那么现在有四个点，我们很自然的会想到，既然两个点确定一条直线，那么最简单的方法就是，两个点之间连一条线，两个点之间连一条线，最后得到的一种折线图如下：这样我们只要确定x=5时的直线，把自变量的值带进去，就显然会得到预测的函数值。但是，这种方法显然具有很明显的缺陷：曲线不够光滑，连接点处的斜率变化过大。可能会导致函数的一阶导数不连续。那么我们应该如何解决这个问题呢？

   

 

二次样条的原理

   我们会想到既然直线不行，那么我们就用曲线来近似的代替和描述。最简单的曲线是二次函数，如果我们用二次函数：aX^2+bx+c来描述曲线，最后的结果可能会好一点，下图中一共有4个点，可以分成3个区间。每一个区间都需要一个二次函数来描述，一共需要9个未知数。下面的任务就是找出9个方程。

   如下图所示：一共有x0,x1,x2,x3四个点，三个区间，每个区间上都有一个方程。

   1>曲线方程在节点处的值必须相等，即函数在x1,x2两个点处的值必须符合两个方程，这里一共是4个方程：

       a1*x1^2+b1*x1+c1=f(x1)

       a2*x1^2+b2*x1+c2=f(x1)

       a2*x2^2+b2*x2+c2=f(x2)

       a3*x2^2+b3*x2+c3=f(x2)

   2>第一个端点和最后一个端点必须过第一个和最后一个方程：这里一共是2个方程

   3>节点处的一阶导数的值必须相等。这里为两个方程。

         2*a1*x1+b1=2*a2*x1+b2

         2*a2*x2+b2=2*a3*x2+b3

  4>在这里假设第一个方程的二阶导数为0：这里为一个方程：

          a1=0

 上面是对应的9个方程，现在只要把九个方程联立求解，最后就可以实现预测x=5处节点的值。

下面是写成矩阵的形式,由于a1=0，所以未知数的个数少了一个。

 

 

下面是二次样条的python实现，最后将结果绘制在图上。

  1 import numpy as np
  2 import matplotlib.pyplot as plt
  3 from pylab import mpl
  4 """
  5 三次样条实现：
  6 函数x的自变量为:3, 4.5, 7, 9
  7 因变量为：2.5, 1 2.5, 0.5
  8 """
  9 x = [3, 4.5, 7, 9]
 10 y = [2.5, 1, 2.5, 0.5]
 11 
 12 
 13 """
 14 功能：完后对三次样条函数求解方程参数的输入
 15 参数：要进行三次样条曲线计算的自变量
 16 返回值：方程的参数
 17 """
 18 def calculateEquationParameters(x):
 19 #parameter为二维数组，用来存放参数，sizeOfInterval是用来存放区间的个数
 20 parameter = []
 21 sizeOfInterval=len(x)-1;
 22 i = 1
 23 #首先输入方程两边相邻节点处函数值相等的方程为2n-2个方程
 24 while i < len(x)-1:
 25 data = init(sizeOfInterval*4)
 26 data[(i-1)*4] = x[i]*x[i]*x[i]
 27 data[(i-1)*4+1] = x[i]*x[i]
 28 data[(i-1)*4+2] = x[i]
 29 data[(i-1)*4+3] = 1
 30 data1 =init(sizeOfInterval*4)
 31 data1[i*4] =x[i]*x[i]*x[i]
 32 data1[i*4+1] =x[i]*x[i]
 33 data1[i*4+2] =x[i]
 34 data1[i*4+3] = 1
 35 temp = data[2:]
 36 parameter.append(temp)
 37 temp = data1[2:]
 38 parameter.append(temp)
 39 i += 1
 40 # 输入端点处的函数值。为两个方程, 加上前面的2n - 2个方程，一共2n个方程
 41 data = init(sizeOfInterval * 4 - 2)
 42 data[0] = x[0]
 43 data[1] = 1
 44 parameter.append(data)
 45 data = init(sizeOfInterval * 4)
 46 data[(si