189 8069 5689
大部分油气藏的数据是散乱分布的，因此称为散乱数据。散乱数据指的是在二维平面上或三维空间中，无规则的、随机分布的数据。利用散乱数据建模就要对散乱数据进行插值或拟合。
网站建设哪家好，找成都创新互联公司！专注于网页设计、网站建设、微信开发、微信小程序开发、集团企业网站建设等服务项目。为回馈新老客户创新互联还提供了五华免费建站欢迎大家使用！
设在二维平面上有n个点 （xi，yi） （i ＝1，…，n），并有Zi ＝f（xi，yi）。插值问题就是要构造一个连续的函数F （x，y），使其在 （xk，yk） （后＝1，…，n） 点的函数值为Zk，即Zk ＝f（xk，yk） （k＝1 ，…，n）。
早在20世纪60年代，散乱数据的插值问题就已引起人们的注意。近50年来，已经有多种算法被提了出来。但是，由于应用问题千差万别，数据量大小不同，对连续性的要求也不同等等，没有一种算法适用于所有的场合。而且大多数算法只能适用于具有中、小规模数据量的散乱点插值问题。大规模散乱数据 （例如，10000个点以上） 的插值问题还正在研究之中。
据散乱数据的复杂程度，其可分为单自变量、双自变量及多自变量3种类型。下面将主要讨论双自变量散乱数据的插值问题。
（一） 插值的一般概念
插值的概念最早可追溯到 “控制论之父” 诺伯特·维纳的不朽著作 《平稳时间序列的外推、插值和光滑及其工程应用》 （Wiener，1949）。随着计算机技术的发展，插值的概念已广泛地应用于数据处理、自动控制、数值分析、地球物理及数学地质等领域。
1. 插值方法
计算机插值方法一般可分为两大类：拟合函数方法和加权平均方法。它们的原理都是来自手工方法。Crain （1970） 把这两种方法得到的曲面分别称为数学曲面和数值曲面。Alfeld ＆ Barnhill （1984） 分别称它们为分片方法和点方法，而Cuyt （1987） 则把这两种方法分别称为系数问题和数值问题。
利用拟合函数方法进行插值，就是利用二元多项式来表示一个插值曲面，插值问题化为确定这个二元多项式的系数的问题。一般来说，往往可通过求解一个线性代数方程组来获取这些系数。这个方程组的系数可由观测数据来确定，它们代表了观测数据的影响，而其方程的次数则表示了控制多项式拟合程度。方程次数越高，拟合的程度就越高。当这个多项式的系数确定以后，将空间某一点处的坐标代入该多项式，即可求得该点处的值。
这个方法的特点是可以制服畸变的原始数据或带有噪声的原始数据。所以，用一个函数进行拟合是一个光滑的过程，一些局部的细节可能消失。所得的插值曲面的复杂程度取决于多项式的次数，即所求解的线性方程的数目。
加权平均插值方法把求插值的问题化为求取观测数据的加权平均。每个观测数据点对应的加权系数恰恰反映了该数据点对插值点的影响大小。为了求取一个插值，必须要计算出一组加权系数。
加权平均方法的一个主要优点是，可以获取在观测数据点附近变量的小尺度趋势。而利用一个适当次数的多项式是无法获取曲面的这种局部细节的。
从原则上讲，这两种方法的差别就在于：加权平均方法强调了曲面的局部细节，而拟合函数方法则概括了曲面的整体性质。从计算时间上看，前者花费的时间比后者要多得多。
2. 插值效果评判
从理论上讲，一个插值方法的效果如何应通过插值结果和客观存在的原始曲面的比较，按以下3条标准来进行判断：
（1） 原始曲面和插值结果之间差异的最大值为最小。
（2） 原始曲面和插值结果之差的平方和为最小。
（3） 在每个观测数据点处，插值结果本身的数值及其1阶到k阶导数和原始曲面的相等。
由于原始曲面本身是未知的，所以在以上3个标准中，第一个和第二个标准是无法检验的，仅有第三个标准是在一定的模型假设之下可以进行检验。在一些实际应用中，当原始曲面可用解析函数来表达时，仅利用第三个标准来检验插值的效果也是可行的。然而，在地质建模中，观测数据点往往不够多，且还有一定的观测误差，不可能断定原始曲面是否可用解析函数表达。这时，插值技术的合理性必须从直观的几何和人们的经验等方面进行评价。
利用计算机进行插值所遇到的困难，主要来自观测数据点数目不足和观测误差。如果观测数据充分多且精确，那么几乎所有的插值方法都会给出良好的效果。另一方面，对于圆形或狭长的隆起，凹陷和鞍点等变量的空间变化几何特征，在数据点分布较稀的情况下，用任何插值方法都是难以推断出它们的存在的。所以，插值方法需要考虑曲面的局部斜率的影响。
下面主要介绍加权平均方法和拟合函数方法中最常用的插值算法。这些算法能解决大部分油气藏建模问题。
（二） 与距离成反比的加权法
距离成反比加权插值方法是基于如下的模型：每个数据点都有局部影响，这个影响随着数据点和插值点距离的增加而减弱，且在一定的范围以外，可以忽略不计；这个影响是以该数据点为中心；而在任一点处的插值恰是各数据点影响之和。
1. 与距离成反比加权插值公式
这一方法首先是由气象学及地质学工作者提出来的，后来由于D. Shepard的工作被称为Shepard方法。其基本思想是将插值函数F （x，y）定义为各数据点函数值f i的加权平均，即：
油气田开发地质学
在 （xk，yk） 点处函数值可写成：
油气田开发地质学
式中： 表示由第i个（xi，yi）点到插值点（xk，yk）的距离；Wi（xk，yk）——权函数；μ——功率因素，通过改变值来调整权函数与距离的关系，与距离成反比加权和与距离成平方反比加权分别是μ＝1和μ＝2时的特殊情形。
与距离成反比加权插值方法是最早使用的计算机插值方法，至今仍被广泛地应用着。在大多数商业性的等值线图绘制软件包中被用来形成网格化数据。这种方法较为直观：一个数据点对于插值点的影响模型化为与这两点之间的距离成反比。
2. 与距离成反比加权插值改进
距离成反比加权算法中的功率因素μ应该取为μ≥0，否则表明距离越远的点作用越突出，这违背了普通常识。考虑以下最极端的情况是：
油气田开发地质学
假设有n个数据点，一个插值点 （xk，yk） 位于第i个数据点附近，相应的权函数可写成：
油气田开发地质学
式中：dj （xk，yk），di （xk，yk）——插值点 （xk，yk）到各数据点的距离。
当插值点 （xk，yk） 和第i个数据点很靠近时，可以认为其他数据点对WD的影响是一个常数，即C为常数。当该插值点和第i个数据点的距离趋于零时，对于不同的μ，WD会有不同的性质。首先看WD对D的导数，有：
油气田开发地质学
然后再有：
油气田开发地质学
可见，当μ＝1时，随着D趋于0，W′D近似为不随D变化的常数，这可用图6-8中左端的图形来表示。
当μ1时，随着D趋于零，W′D趋于零。这说明在该数据点附近，加权系数的变化为零，即可用图6-8中间的图形来表示。
当μ1时，随着D趋于零，W′D趋于无穷大。这说明加权系数在该数据点附近还有一个尖点，如图6-8右端的图形所示。
图6-8 与距离反比加权的权数随参数μ的变化
（1）功率因素μ越小时，近距离点和远距离点的作用越接近，生成的平面网格数据越平滑。随着μ增大，平面网格数据的光滑性越差。同时，μ直接影响网格数据的极值和均值，μ越小网格数据的均值越接近原始数据的均值，但极值相差越大。因此，当要求插值结果尽可能接近原始数据的均值时，功率因素不能选择过大。例如对于开发早期的油气藏，因为仅有少数探井控制，此时网格化得到的各类物性参数应该在总体上符合井点的统计结果，均值是比极值更有价值的参数，因而通常将功率因素取为1。相反，μ取值越大，网格数据越能恢复原始数据的极值，但也容易使均值误差增大。原因是当μ取较大的值时，近距离点的作用越突出，原始数据点分布的不均匀性使得部分点在网格节点上发挥了更大的作用，而另外一些点的作用则受到屏蔽。因此，当要求突出数据的局部特征，体现储层的非均质性特征时，功率因素应选择得大些，一般取为2。
（2） 利用与距离成反比加权法进行插值时，当增加、删除或改变一个点时，权函数Wi （xk，yk） 均需重新计算，因而该方法是一个全局插值算法。
为了克服Shepard方法的上述缺陷，Franke及Nielson提出了MQS （Modified QuadraticShepard） 方法，它仍然是一个与距离成反比的加权方法。对它的改进如下：
插值点 （xk，yk） 到已知数据点的距离di作适当修改，使其只能在局部范围内起作用，以改变Shepard方法的全局插值性质。这时重新定义距离函数：
油气田开发地质学
式中：rw为一个常数。而
油气田开发地质学
因此，当 （xk，yk） 点与某一点的距离大于rw时，权值就为零。
（3） 与距离成反比加权法的权函数Wi（xk，yk）始终满足Wi（xk，yk）≤1/n，因此插值结果不会大于或小于原始数据的最大值、最小值。当已知数据点过少时 （这种情况在早期地质研究中是最为常见），使用具有外推能力的曲面样条或趋势面分析，得到的结果往往背离实际。其原因是这两种方法在远点不具有控制能力，它将沿趋势无限发展下去。特别是在仅有少数井资料可用的情况下，与距离成反比加权法应是优先选择的方法。
但是，隐含在原始数据中的尖峰会被淹没而无法显示出来。因为距离成反比加权进行插值时，每个数据点所发挥的作用基本上是中心对称的，因此对山脊和山谷等非各向同性的几何形状的显示不利。为了克服这种状况，需要考虑变量的局部变化趋势，为此需要对梯度进行估计。
当已知数据点用与距离成反比加权方法形成数据插值曲面时，该曲面被称数据曲面。与距离成反比加权方法也可用于各个数据点处的切平面，把任一点处的插值值取成各切平面在该点处取值的一个加权平均，所形成的曲面称为与距离成反比加权梯度插值曲面，简称梯度曲面。如果在数据点以外的一个点是变量的局部高点，那么梯度曲面在该点的值容易大于变量的真实值，即呈现 “过估计” 的状态。如果数据曲面在该点的值小于变量的真实值，则呈现 “欠估计” 的状态。因此，数据曲面可以通过和梯度曲面的相互结合来克服本身的缺陷。
可以用数据曲面和梯度曲面之差乘以一个系数作为一个修正量，对数据曲面进行修正，这样，可以用下述的曲面来代替单纯的数据曲面：
油气田开发地质学
式中：L （x，y） ——距离反比加权插值；Wi——和第i个数据点的距离反比加权系数；τi——曲面在该点处的粗糙度指数；Si（x，y）——第i个数据点 （xi，yi） 处的切平面在点（x，y） 处的值；H（Wi，τi）——混合函数。
如此得到的曲面称为混合曲面。它通过所有的数据点，具有连续的坡度，其变化在空间的分布更均匀。数据曲面和梯度曲面是混合曲面的两种极端情况。由于数据曲面和梯度曲面之差在各数据点处为零，还因为混合函数的变化范围为0～1，且当混合函数等于0或1时，其一阶导数为零，故混合曲面和梯度曲面相切于各数据点处，且其高阶导数在各数据点处亦为零。
（4） 与距离成反比加权法仅考虑了插值点与数据点之间距离的影响，没有考虑到各数据点之间的关系，物性参数分布的趋势性没有得到充分的体现。为此，Franke及Nielson进行了改进。
用节点函数Qi （x，y） 代替fi，Qi （x，y） 是一个插值于 （xi，yi） 点的二次多项式，即有Qi （xi，yi） ＝f，i＝1，…，n。Qi可由下式表示：
Qi（x，y）＝fi＋a1（x-xi）＋a2（y-yi）＋a3（x-xi）2＋a4（x-xi）（y-yi）＋a5（y-yi）2
式中：a1，a2，…，a5是按下式最小二乘法得出的优化解：
油气田开发地质学
式中：fi，fj分别为 （xi，yi）和 （xj，yj）点的函数值，而ρj可按下式选取：
油气田开发地质学
其中rq为一常数，而
油气田开发地质学
求出Qi （x，y） 后，插值函数可表示为：
油气田开发地质学
上述方法消除了Shepard方法中的一些缺陷，因而在散乱点插值中得到广泛的应用。但是，为了求得Qi （x，y） （i＝1，…，n），需要多次求解线性方程组，计算量大，因此，一般只用于中、小规模散乱点的插值运算。
（三） 多项式趋势面法
由计算机产生的曲面一般不会总是和原始的观测数据一致。如果两者的差别在给定的尺度之下不是很明显，那么产生的曲面可被认为是插值曲面，否则就被认为是近似曲面。如果观测数据含有明显的观测误差，近似曲面就显得更合理。这时，和插值曲面相比，近似曲面由于数据的各种误差所产生的扰动不太容易看得清，但是近似曲面空间变化的一些主要性质还是能清晰地被体现出来的。
近似曲面和每个数据点之间的差称为残差，可视为每一个数据点上的一种误差表示。然而，计算出来的这种残差是意味着对未知的观测误差的一种度量，还是意味着一种允许的插值误差，或者意味两者都是，这要依变量的空间性质和观测数据的获取方法而定。
确定近似曲面的方法可分为3种。第一种方法是以残余的平方和最小为条件，确定多项式的系数，以获取曲面。第二种方法是利用观测数据误差的附加信息，并满足最小曲率的原则以确定曲面。最后一种方法是利用观测误差和插值误差的附加信息，以满足最小平方差或最小曲率为条件确定曲面。以下主要讨论多项式构造趋势面法。
多项式构造趋势面是目前最常用的方法，一次多项式表示的趋势面是空间的一个平面，二次趋势面是抛物面，椭球面或双曲面，三次及三次以上的趋势面是形态复杂的空间曲面，随着趋势面的次数增高，曲面的形态就越复杂。
如果有一组总共n个观测数据，其观测点的平面坐标为 （xi，yi），地质变量的观测值为fi（xi，yi）。对于这组观测数据的多项式趋势面方程表示成如下形式：
油气田开发地质学
式中： ——第i个观测点的趋势值；a1，a2，…，a5——待定系数，它们的个数m与所选用趋势面方程的多项式次数n存在下列关系式：
m＝［n（n＋3）＋2］/2
为使趋势面最大限度地逼近原始观测数据，可采用最小二乘法使每个观测点的观测值与趋势值之差 （残差） 的平方和最小，即：
油气田开发地质学
得到需求解的m阶正规方程组。当系数矩阵满秩时，趋势面方程也就被唯一确定了。
由于趋势面分析不具有过点性，使得局部井点上误差可能很大。同时参数场在三维空间中的分布过于复杂，无论从理论上还是实验中都无法确证某类参数场能较好地符合某确定次数的曲面，多项式次数过高，会导致趋势面发生频繁振动，多项式次数过低，得到的趋势面又过于光滑，丧失许多细节。再者，趋势面方程在外推过程中容易使参数场发生畸变，产生无意义的结果。因而现代地质建模研究中已经很少将趋势面分析单独作为插值方法使用。但对于下列两种情况趋势面分析仍然能达到较为理想的效果，一是对小范围内具有显著趋势性分布的数据点；二是数据点分布过于密集，而且可能存在若干异常数据点时，趋势面分析会自动削弱异常数据点的影响。为提高趋势面分析的精度，可采用残差来校正趋势面分析的结果。具体过程如下：
（1） 由趋势面分析得到任一网格节点 （xk，yk）处的趋势值
油气田开发地质学
（2）计算n个已知点 （xi，yi）处的残差△fi（xi，yi）＝fi（xi，yi）
油气田开发地质学
（3） 调用某种插值方法将残差分配到每个网格节点上，对网格节点 （xk，yk） 有△fk（xk，yk）；
（4） 网格节点 （xk，yk）经校正后的最终结果为
油气田开发地质学
yk）。
特别地，如果残差分配的插值算法也是趋势面分析，就形成所谓的多级趋势面分析。此时，网格节点上的值是多次趋势面分析的结果，多级趋势面分析在不断减少观测点插值误差的同时，整张参数场曲面仍然保持其连续性和光滑性，原因是该算法同样符合线性迭加原理。
（四） 径向基函数插值法
径向基函数的名字来源于这样一种情况，即基函数是由单个变量的函数构成的。一个点 （x，y） 的这种基函数的形式往往是hk（x，y）＝h（dk），这里的dk表示由点 （x，y） 至第k个数据点的距离。一般说来，这种方法不具有多项式精度，但只要稍加改进，即可获得具有多项式精度的插值公式：
油气田开发地质学
式中：qk（x，y）是一个多项式基，其阶次小于m。
上式中的系数ak和bk应满足下面的联立方程组：
油气田开发地质学
油气田开发地质学
第一式中的n个方程式满足了插值要求，而第二式中的m个方程式则保证了多项式精度。两式中共有m＋n个未知数，同时存在m＋n个方程式，联立求解，即可得出待定系数。
下面，介绍两种主要的径向基函数插值法。
1. Multiquadric方法
Multiquadric方法是由R. L. Hardy在1971年提出来的。它是最早被提出并且应用得最为成功的一种径向基函数插值法。它采用的插值函数，即 （x，y）处的值F （x，y）：
油气田开发地质学
式中： 为基函数；ei——非负常数；ai——加权系数，满足如下方程组：
MVa＝Vz
式中：Va=(a1,a2,…,an)T,Vz=[f(x1,y1),f(x2,y2),…,f(xn,yn)]T,f(xi,yi)是(xi,yi)处的数据点的值。
油气田开发地质学
由于M和Vz都不依赖于插值点的坐标 （x，y），所以ai也不依赖于 （x，y）。然而，基函数C（X-Xi）则是以Xi为参数的 （x，y）的函数。所以说，插值曲面F（x，y）是n个基函数C（X-Xi）所构成的n个空间曲面配置而成的。
Arther提出了如下形式的基函数：C（d）＝1-d2/e2，其中d是数据点到插值点之间距离，而e则是一个常数。显然，基函数C（d）是d的一个衰减函数。当d＝0时，C（d）取得最大值1，而当d≤e时，有0≤C（d）≤1。这时，基函数呈现为椭圆抛物面，而加权系数ai（i＝1，2，…，n）所满足的线性代数方程组的矩阵M应作相应的改动，其对角线元素应改成1。Hardy（1971）引入如下的基函数：C（d）＝（d2＋e2）1/2，该基函数呈现为椭圆双曲面。Hardy还建议将e2取成0.815乘以数据点间距离的平均值。
2. 薄板样条法
样条 （Spline）本来是绘图员用来绘制光滑曲线的工具，是一种用木材或金属等弹性材料做成的细条。在绘图时，沿着通过图纸各已知点的样条，便可绘出一条光滑曲线。数学上所说的样条 （多项式样条） 实质上是分段多项式曲线的光滑连接。当函数为分段的m次多项式，在分段点上有直至m-1阶连续导数，那么该函数则称为m次样条函数，简称为样条。一般来说，研究和应用得比较多的是三次样条。零次和一次样条函数分别是台阶状函数和折线状函数。以上所述的是关于一维样条函数，对于二维样条函数也可作为类似的考虑。
样条函数的主要功能是进行插值，其主要优点在于，能在插值多项式的次数尽可能低的条件下，使插值曲线或插值曲面取得较高的光滑度，且只需要利用函数本身的值，而不需要提供函数的各阶导数的值。
三次样条曲面包含有三种不同的类型：双三次样条、伪三次样条及薄片样条。这些样条曲面以m和s为其两个参数，使得希氏空间Hs中元素的m阶导数的范数所构成的一个泛函达到最小。此外，这个泛函具有旋转不变性。
对于薄片样条曲面，m＝2，s＝0。这一方法是由R.L. Harder及R. N. Desmarais在1972年提出来的，后来由J. Duchon及J. Meinguet等人予以发展。薄板样条法得名于如下事实，即用此方法求出的散乱点的插值函数使下面这一泛函表达式具有最小值：
油气田开发地质学
在这里，I（F）表示受限于插值点的无限弹性薄板的弯曲能量。因此，这一方法的实质从力学观点看是使插值函数所代表的弹性薄板受限于插值点，并且具有最小的弯曲能量。这是一个泛函求极值的问题。这一变分问题的解即为我们所需要的插值函数，具有径向基函数插值法的一般形式。
R.L. Harder及R. N. Desmarais提出解析形式如下：
油气田开发地质学
且有 其中t（x，y）和ti（xi，yi）是二维空间中的点，而fi是ti处的观测值。还有，K（ti，t） 这里，ri代表点t和ti之间的距离，
由解析表达式及其约束条件，可给出用以确定系数的线性代数方程组：
油气田开发地质学
其中，K=[kij]n×n,kij=K(ti,tj),kii=0,FT[f1,f2,…,fn],αT=[b,a1,a2],AT［λ1，λ2，…，λn］，且有：
油气田开发地质学
求解上述n＋3阶方程组则得到待定系数b，a1，a2，λi，然后可插值出平面任一位置的函数值F（x，y）。
上述方程与Enriguez等人给出的方程相类似，其微小的差异就是基函数中的自然对数（In） 变成了常用对数 （lg）。
在构造薄片样条曲面的过程中，Franke （1982）提出了以r2lgr作为基函数，Sandwell（1987）则提出了以双调和格林函数r2（lgr-1）作为基函数。另外，Ayeni （1979）也对不同的基函数进行了讨论。
使用平面上n个已知点进行曲面样条插值时，实质上是求解一个n＋3阶线性方程组以确定n＋3个系数。为了保证解的存在性和唯一性，系数矩阵应该是满秩的。对下列3种情况必须避免：
（1） 在给定的n个已知点中存在着距离过近的两点 （xi，yi） 与 （xj，yj） 极端的情形是同一个数据点的重复输入，此时系数矩阵的第i行与第j行对应各元素非常接近，导致线性方程组的系数矩阵是奇异的。因此在数据预处理过程中必须消除沉余数据。
（2） 给定的已知数据点过少，此时系数矩阵的后3行线性相关，矩阵是不满秩的。
（3） n个已知点数据分布在一条直线上，显然由这样的n个点不能唯一决定一张曲面，系数矩阵表现为不满秩。
针对情况 （1），通常在做曲面样条插值时，首先对数据进行预处理，通过给定一个适当的距离下限rmin来滤掉那些相距过近的点，研究发现rmin＝（△x＋△y）/8是一个较合理选择 （△x和△y分别表示x和y方向的步长）。情况 （2） 和 （3） 实际上意味着不能进行曲面样条插值，除非通过数据均整来改变数据分布状态。
曲面样条插值方法是一种严格的过点插值法，即由生成的样条曲面必定通过给定的n个已知数据点，这样井点数据的控制作用自然得到体现。同时，曲面样条方法充分考虑了数据间的相对位置，其插值精度很高，在外推过程中，总是沿数据点的分布趋势外推，因此曲面样条法是具有一定外推能力的插值方法。同时曲面样条方程得到的是一张连续光滑的曲面。
曲面样条插值方法特别适合于地层层面的生成和地层厚度的插值。如果从曲面样条法严格的过点性、良好的光滑性及外推性看，适用于那些光滑、趋势性明显、变化连续的储层物性参数诸如油气饱和度的插值。
曲面样条法插值的精度很大程度上取决于数据点分布的均匀程度，稀疏区域主要由邻近区域的外推得到，其插值结果可能偏差较大。同时，曲面样条法的严格过点性使得它不能分别对待不同精度点的数据，因此它不具备数据的校正能力。
#includestdio.h
double fun(double x)
{//多项式
double a[]={2,-4,0,3};
double sum;
int n=4;
n--;
sum=a[n];
while(--n = 0)
sum=sum*x+a[n];
return sum;
}
typedef double (*Function)(double);
double cal(Function f,double a,double b,double e)
{
double a1,a2,a3,t_a2,t;
a1=a;a3=b;a2=(a+b)/2;
while(1)
{
t=((f(a3)-f(a2))/(a3-a2)-(f(a2)-f(a1))/(a2-a1))/(a3-a1);
t_a2=((a1+a2)*t-(f(a2)-f(a1))/(a2-a1))/(2*t);
if(t_a2-a2 e t_a2-a2 -e)break;
if(t_a2 a2)a1=a2;
else a3=a2;
a2=t_a2;
}
return f(t_a2);
}
void main()
{
printf("%lf\n",cal(fun,0,2,0.01));
}希望能帮到你 望采纳 谢谢
/*********************下面程序是C语言程序（标准C）******************/
/* 计算给定M0,Mn值的三次样条插值多项式 */
/*给定离散点(1.1,0.4),(1.2,0.8),(1.4,1.65),(1.5,1.8),M0=Mn=0,*/
/*用M关系式构造三次样条插值多项式S(x),计算S(1.25)。 */
/*************************************************************/
#include stdio.h
#define Max_N 20
main()
{int i,k,n;
double h[Max_N+1],b[Max_N+1],c[Max_N+1],d[Max_N+1],M[Max_N+1];
double u[Max_N+1],v[Max_N+1],yy[Max_N+1],x[Max_N+1],y[Max_N+1];
double xx,p,q,S;
printf("\nPlease input n value:"); /*输入插值点数n*/
do
{ scanf("%d",n);
if(nMax_N)
printf("\nplease re-input n value:");
}
while(nMax_N||n=0);
printf("Input x[i],i=0,...%d:\n",n-1);
for(i=0;in;i++) scanf("%lf",x[i]);
printf("Input y[i],i=0,...%d:\n",n-1);
for(i=0;in;i++) scanf("%lf",y[i]);
printf("\nInput the M0,Mn value:");
scanf("%lf%lf",M[0],M[n]);
printf("\nInput the x value:"); /*输入计算三次样条插值函数的x值*/
scanf("%lf",xx);
if((xxx[n-1])||(xxx[0]))
{printf("Please input a number between %f and %f.\n",x[0],x[n-1]);
return;
}
/*计算M关系式中各参数的值*/
h[0]=x[1]-x[0];
for(i=1;in;i++)
{h[i]=x[i+1]-x[i];
b[i]=h[i]/(h[i]+h[i-1]);
c[i]=1-b[i];
d[i]=6*((y[i+1]-y[i])/h[i]-(y[i]-y[i-1])/h[i-1])/(h[i]+h[i-1]);
}
/*用追赶法计算Mi,i=1,...,n-1*/
d[1]-=c[1]*M[0];
d[n-1]-=b[n-1]*M[n];
b[n-1]=0; c[1]=0; v[0]=0;
for(i=1;in;i++)
{u[i]=2-c[i]*v[i-1];
v[i]=b[i]/u[i];
yy[i]=(d[i]-c[i]*y[i-1])/u[i];
}
for(i=1;in;i++)
{M[n-i]=yy[n-i]-v[n-i]*M[n-i+1];
}
/*计算三次样条插值函数在x处的值*/
k=0;
while(xx=x[k]) k++;
k=k-1;
p=x[k+1]-xx;
q=xx-x[k];
S=(p*p*p*M[k]+q*q*q*M[k+1])/(6*h[k])+(p*y[k]+q*y[k+1])/h[k]-h[k]*(p*M[k]+q*M[k+1])/6;
printf("S(%f)=%f\n",xx,S); /*输出*/
getch();
}
/*----------------------------------- End of file ------------------------------------*/
/*程序输入输出：
Please input n value:4
Input x[i],i=0,...3:
1.1 1.2 1.4 1.5
Input y[i],i=0,...3:
0.4 0.8 1.65 1.8
Input the M0,Mn value: 0 0
Input the x value:1.25
S(1.250000)=1.033171
*/