612 三次样条函数插值(第三种边界条件)与微分
一、 功能
给出 n 个节点Xi 及节点上函数值Yi = f (Xi),i=1,2,…,n, 关于第三种边界条件构造三次样条插值函数S (X),对 m 个插值点Xj ( j=1,2,…,m),作成组插值,并可计算插值点上的一、二阶导数;还可计算节点上的一、二阶导数;也可用求出 S (X) 在插值区间 [a,b] 上的积分,近似代替 f (X) 在[a,b]上的积分。
二、算法简介
设在 [a,b] 上取 n 个节点
a = X1<X2<…<Xn-1<Xn = b
给定这些节点上的函数值Yi = f (Xi) (i=1,2,…,n),要求构造一个三次样条插值函数S (X),使得满足下列条件:
(1)S (Xi) = Yi, i=1,2,…,n ;
(2)在每个小区间 [Xi, Xi+1]上是一个三次多项式;
(3)S (X )∈C2[a,b]。
第三种边界条件有几种给法,我们这里给的第三种边界条件,是指给定的函数f(X)是周期函数,因而三次样条函数S(X)亦为周期函数,其基本周期为Xn-X1,此时就是Y1=f(X1)=Yn=f(Xn),在插值区间的端点上满足:
S’(X1)=S’(Xn), S”(X1)=S”(Xn)
由此边界条件,可由递推法求出 n 个节点上S (X)的一阶导数mj (j=2,3,…,n-1),计算公式为:
a1=-0, b1=0, c1=0
h1=hn-1, Y2-Y1=Yn-Yn-1
hi=Xi+1-Xi
αi=hi-1 / (hi-1+hi)
βi=3[(1-αi)(Yi-Yi-1)/ hi-1+αi(Yi+1-Yi)/ hi] (i=1,2,…,n-1)
ai=-αi-1/[2+(1-αi)αi-1]
bi=-(1-αi) bi-1/[2+(1-αi-1)ai-1]
ci=[βi-1-(1-αi-1)ci-1]/[2+(1-αi-1)ai-1] (i=2,3,…,n-1)
tn=1, vn=0
ti=aiti+1+bi
vi=aivi+1+cI (i=n-1,n-2,…,2)
mn-1=[βn-1-αn-1v2-(1-αn-1)vn-1]/[2+αn-1t2+(1-αn-1)tn-1]
mi=ti+1mn-1+ vI+1 (i=1,2,…,n-2))
将计算出的mi (i=1,2,…,n) 和节点及节点上的函数值,在小区间[Xi,Xi+1] 上构造三次样条插值函数S(x)及S’(x)与S”(x):
S (X)=[3/ h2i(Xi+1-X)2-2/h3i(Xi+1-X)3]Yi-[3/h2i(X-Xi)2-2/h3i(X-Xi)3]Yi+1+hi [1/h2i(Xi+1-X)2-1/h3i(Xi+1-X)3]mi - hi[3/h2i(X-Xi)2-2/h3i(X-Xi)3]mi+1
S’(X)=6/hi[1/h2i(Xi+1-X)2-1/hi(Xi+1-X)]Yi-6/hi[1/h2i(X-Xi)2-1/hi(X-Xi)]Yi+1+ [3/h2i(Xi+1-X)2-2/hi(Xi+1-X)]mi - [3/h2i(X-Xi)2-2/hi(X-Xi)]mi+1
S”(X)=1/h2i[6-12/hi(Xi+1-X)]Yi-1/h2i[6-12/hi(X-Xi)]Yi+1+1/hi[2-6/hi(Xi+1-X)]mi- 1/hi[2-6/hi(X-Xi)]mi+1
在节点Xi与Xi+1上(I=1,2,…n-1),二阶导数
Y”i=S”(Xi)=6/h2i(Yi+1-Yi)-2/hi(
Y”i=S”(Xi)=6/h2i(Yi -Yi+1)-2/hi(mi+
利用节点上的二阶导数,可得到改进的梯形求积公式,计算出S (X)在插值区间[a,b]上积分S,近似代替在上的积分:
S=∫baS(X)dx=1/2∑
由上述知,三次样条函数插值,是一种改进的分段插值,在分段处具有连续的二阶导数,且连续点保持光滑。但因为三次多项式计算简单,且满足一般实际问题的要求,故用的最多。由于样条函数明显的优点,目前被广泛应用。
三、程序使用说明
1、输入参数
n,m分别为节点个数与插值点个数
w 控制计算内容的变量,当w=1时,计算节点上的一,二阶导数; 当w=2时,除计算节点上的一,二阶导数外,还计算插值点上的函数值及一、二阶导数;当w<>1且w<>2时,如w=3,除计算w=2时的内容外,还计算S (X) 在插值区间[a,b]上的积分,近似代替f(X)在[a,b]上的积分。
X(N),Y(N)是一维实数组,分别存放节点及节点上的函数值
C(M)是一维实数组,存放插值点的数值
X(N),Y(N),C(M),分别存放在程序的70-150语句中
2、输出参数
n,m分别为节点个数及插值点个数
w 控制计算内容的变量[见1]
部分节点及节点上的函数值
部分插值点及插值点上函数值的一、二阶导数
部分节点及节点的一、二阶导数
S S (X)在插值区间[a,b]上的积分
请大家帮帮忙啊!!!