三次样条插值函数

[tianfeng]

有谁知道三次样条插值函数的计算代码？

[yonk]

参见 三次样条插值函数 的讲解
《计算方法论》里面有。
或者你看看matalab中有没有源程序（当初学它的时候尽是偷懒了，所以不清楚有没有）
我的《计算方法论》带回家了，现在不能帮你。
如果你实在不能解决，等我回来找找书咱们再搓商。
现在要去bathing了。

[tianfeng]

急！！！

[tianfeng]

请问哪里有求三次样条插值函数的源代码下载

[X_Ben2000]

随便找本计算方法的书里面都有！！

[feixiaer]

Matlab中查函数spline或csapi、csape等函数的源代码，应该有的。一般的计算数学书上都会讲这个东东的。

[yonk]

刚bathing 回来。拿同学的《计算方法引论》第二版，徐萃薇 孙绳武编著，抄来如下：

2.7 样条插值函数
。。。。
设在区间[a, b]上取n+1个节点
a = X0 < X1 < X2 <.....<Xn-1 < Xn =b
给定这些点上的函数值 f(Xi) = Yi (i = 0, 1,...n)。构造一个三次样条插值函数
S(x),使得满足下列条件：
(1) S(Xi) = Yi (i = 0, 1, 2 , ..., n)
(2) 在每个小区间[Xi, Xi+1]上是一个三次多项式
(3) S(x) ∈ C^2 [a, b]
假设在区间[a, b] 上三次样条插值函数S(x）存在，并利用Mi来表示S(x)在点Xi处的
微商值。由于曲线通过点 (Xi, Yi) （i = 0, 1, ..., n)并且在每一个小区间上满足
条件
S(Xi) = Yi, S(Xi+1) = Yi+1 ， S'(Xi) = Mi, S'(Xi+1) = Mi+1
使用Hermite插值公式写出小区间[Xi, Xi+1]上的三次样条插值函数S(x)的计算公式：
S(x) = (1 + 2 * (X - Xi)/(Xi+1 - Xi)) * ((X - Xi+1)/(Xi+1 - Xi))^2 * Yi
+ (1 + 2 * (X - Xi+1)/(Xi-Xi+1)) * ((X - Xi)/(Xi+1 - Xi))^2 * Yi+1
+ (X - Xi) * (((X - Xi+1)/(Xi+1 - Xi))^2) * Mi
+ (X - Xi+1) * (((X - Xi)/(Xi+1 - Xi))^2) * Mi+1
记作 式 （2.55）
注意，在样点Xi上的微商值Mi并不知道。因此如果要利用上述公式，必须先求出Mi
利用 二阶微商的性质，对上式作二阶微商，并令 Hi = Xi+1 - Xi,便得到
S''(x) = (6/Hi^2 - 12 * (Xi+1 - X)/Hi^3) * Yi +
(6/Hi^2 - 12 * (X - Xi)/Hi^3) * Yi+1 +
(2/Hi - 6 * (Xi+1 - X)/Hi^2) * Mi -
(2/Hi - 6 * (X - Xi)/Hi^2) * Mi+1
从而我们可以得到Xi的右微商， 记作 S''(Xi+):
S''(Xi+) = - 6 * Yi / Hi^2 + 6 * Yi+1 / Hi^2 - 4 * Mi / Hi - 2 * Mi+1 / Hi
记作 式 （2.56）
同样得到Xi的左微商：
S''(Xi-) = 6 * Yi-1 / Hi-1 ^ 2 - 6 * Yi / Hi-1 ^2 +
2 * Mi-1 / Hi-1 + 4 * Mi / Hi-1
记作 式 （2.57）

由于 在点Xi处微商连续，所以 S'' (Xi+) == S'' (Xi-)
将上式整理，并令
Ai = Hi-1 / (Hi-1 + Hi)
Bi = 3 * ((1 - Ai) * (Yi - Yi-1) / Hi-1 + Ai * (Yi+1 - Yi)/Hi)
记作 （2.58）
对每一个内点建立方程，得到方和组
(1 - Ai) * Mi-1 + 2 * Mi + Ai * Mi+1 = Bi, i = 1, 2, .. , n-1
记作 （2.59）
这是关于n+1 个未知量M0, M1, .. Mn的n-1个线性方程组。有无穷多个解。但实际问题
只能取特定的一个解。因而补充两个附加条件，称作边界条件。
常见的边界条件有：
（！） 曲线在两端点X0, Xn处的切线斜率已知，即S'(X0) = M0以及
S'(Xn) = Mn已知，方程组（2.59）就成为具有n-1个未知数M1, M2, ... Mn-1的n-1个方程
的方程组，有唯一解。
（2） 函数 Y = f(x)在两端点X0, Xn处二阶微商为0，从（2.56）和（2.57）得：
2 * M0 + M1 = 3 * (Y1 - Y0)/H0
Mn-1 + 2 * M0 = 3 * (Yn - Yn-1)
记作 式 （2.60）把2.60和2.59联立，可唯一解出Mi
（3）函数 Y = f(x)是周期函数，基本周期是 b - a = Xn - X0,这时，
Y0 = Yn,相应地要求样条插值函数 S(X）也是周期函数。在端点满足S'(X0) = S'(Xn)
S''(X0) = S''(Xn)，于是有
M0 = Mn,

哎呀，累死我了。说明一下，由于无法按数学方式表示出来，Mi+1表示的下标为i+1的M
其它也一样。
我建议你还是找相关的书籍看看吧。这些公式虽说可以帮你做出来，但是，我想也不是很容易就能实现的。
3 * (Y1 - Y0)/H062 - （2*M0 + M1)/H0 =
3 * (Yn-1 - Yn)/Hn-1^2 + (Mn-1 + 2 * Mn)/Hn-1 记作2。61
联立2.59, 2.61，也能解出Mi

[yonk]

不知怎么回事，本来应该在下面的话（哎呀。。。）却不在最下面。呵呵

[tianfeng]

你虽然写了很多关于这个算法的思想，但我更需要的是要能运行的源代码呀！

[tianfeng]

我写了一个函数，但运行起来有误，谁能帮我检查一下

/* 三次样条插值计算算法 */
#include "math.h"
#include "stdio.h"
#include "stdlib.h"


/*
N：已知节点数N+1
R：欲求插值点数R+1
x，y为给定函数f(x)的节点值{x(i)} (x(i)<x(i+1)) ,以及相应的函数值{f(i)} 0<=i<=N
P0=f(x0)的二阶导数；Pn=f(xn)的二阶导数
u:存插值点{u(i)} 0<=i<=R
求得的结果s(ui)放入s[R+1] 0<=i<=R
返回0表示成功，1表示失败
*/
int SPL(int N,int R,double x[],double y[],double P0,double Pn,double u[],double s[])
{
/*声明局部变量*/
double *h; /*存放步长:{hi} 0<=i<=N-1 */
double *a; /*存放系数矩阵{ai} 1<=i<=N ； 分量0没有利用 */
double *c; /*先存放系数矩阵{ci} 后存放{Bi} 0<=i<=N-1 */
double *g; /*先存放方程组右端项{gi} 后存放求解中间结果{yi} 0<=i<=N */
double *af; /*存放系数矩阵{a(f)i} 1<=i<=N ； */
double *ba; /*存放中间结果 0<=i<=N-1*/
double *m; /*存放方程组的解{m(i)} 0<=i<=N ； */

int i,k;
double p1,p2,p3,p4;

/*分配空间*/
if(!(h=(double*)malloc(N*sizeof(double)))) exit(1);
if(!(a=(double*)malloc((N+1)*sizeof(double)))) exit(1);
if(!(c=(double*)malloc(N*sizeof(double)))) exit(1);
if(!(g=(double*)malloc((N+1)*sizeof(double)))) exit(1);
if(!(af=(double*)malloc((N+1)*sizeof(double)))) exit(1);
if(!(ba=(double*)malloc((N)*sizeof(double)))) exit(1);
if(!(m=(double*)malloc((N+1)*sizeof(double)))) exit(1);

/*第一步：计算方程组的系数*/
for(k=0;k<N;k++)
h[k]=x[k+1]-x[k];
for(k=1;k<N;k++)
a[k]=h[k]/(h[k]+h[k-1]);
for(k=1;k<N;k++)
c[k]=1-a[k];
for(k=1;k<N;k++)
g[k]=3*(c[k]*(y[k+1]-y[k])/h[k]+a[k]*(y[k]-y[k-1])/h[k-1]);
c[0]=a[N]=1;
g[0]=3*(y[1]-y[0])/h[0]-P0*h[0]/2;
g[N]=3*(y[N]-y[N-1])/h[N-1]+Pn*h[N-1]/2;

/*第二步：用追赶法解方程组求{m(i)} */
ba[0]=c[0]/2;
g[0]=g[0]/2;
for(i=1;i<N;i++)
{
af[i]=2-a[i]*ba[i-1];
g[i]=(g[i]-a[i]*g[i-1])/af[i];
ba[i]=c[i]/af[i];
}
af[N]=2-a[N]*ba[N-1];
g[N]=(g[N]-a[N]*g[N-1])/af[N];

m[N]=g[N]; /*P110 公式：6.32*/
for(i=N-1;i>=0;i--)
m[i]=g[i]-ba[i]*m[i+1];

/*第三步：求值*/
for(i=0;i<=R;i++)
{
/*判断u(i)属于哪一个子区间，即确定k */
if(u[i]<x[0] || u[i]>x[N])
{
/*释放空间*/
free(h);
free(a);
free(c);
free(g);
free(af);
free(ba);
free(m);
return 1;
}
k=0;
while(u[i]>x[k+1])
k++;


p1=(h[k]+2*(u[i]-x[k])*pow((u[i]-x[k+1]),2)*y[k])/pow(h[k],3);
p2=(h[k]-2*(u[i]-x[k+1])*pow((u[i]-x[k]),2)*y[k+1])/pow(h[k],3);
p3=(u[i]-x[k])*pow((u[i]-x[k+1]),2)*m[k]/pow(h[k],2);
p4=(u[i]-x[k+1])*pow((u[i]-x[k]),2)*m[k+1]/pow(h[k],2);
s[i]=p1+p2+p3+p4;
}

/*释放空间*/
free(h);
free(a);
free(c);
free(g);
free(af);
free(ba);
free(m);

return 0;
}


void main()
{
int N,R;
double *x,*y,*u,*s;
double P0,Pn;
int i;

/*验证算法：*/
N=7;
R=6;

/*分配空间*/
if(!(x=(double*)malloc((N+1)*sizeof(double))))
{
printf("malloc error!\n");
exit(1);
}
if(!(y=(double*)malloc((N+1)*sizeof(double))))
{
printf("malloc error!\n");
exit(1);
}
if(!(u=(double*)malloc((R+1)*sizeof(double))))
{
printf("malloc error!\n");
exit(1);
}
if(!(s=(double*)malloc((R+1)*sizeof(double))))
{
printf("malloc error!\n");
exit(1);
}

x[0]=0.5;x[1]=0.7;x[2]=0.9;x[3]=1.1;x[4]=1.3;x[5]=1.5;x[6]=1.7;x[7]=1.9;
y[0]=0.4794;y[1]=0.6442;y[2]=0.7833;y[3]=0.8912;y[4]=0.9636;y[5]=0.9975;y[6]=0.9917;y[7]=0.9463;
u[0]=0.6;u[1]=0.8;u[2]=1.0;u[3]=1.2;u[4]=1.4;u[5]=1.6;u[6]=1.8;

P0=-0.4794;
Pn=-0.9463;

if(!SPL( N, R, x, y, P0, Pn, u, s))
{
/*打印结果*/
printf("\nx= ");
for(i=0;i<=N;i++)
printf("%8.1f",x[i]);

printf("\ny= ");
for(i=0;i<=N;i++)
printf("%8.4f",y[i]);

printf("\n\nu= ");
for(i=0;i<=R;i++)
printf("%9.2f",u[i]);
printf("\ns= ");
for(i=0;i<=R;i++)
printf("%9.5f",s[i]);
printf("\nsin= ");
for(i=0;i<=R;i++)
printf("%9.5f",sin(u[i]));
}

/*释放空间*/
free(x);
free(y);
free(u);
free(s);
}

[tianfeng]

我已经检查出来了。
p1=(h[k]+2*(u[i]-x[k]))*pow((u[i]-x[k+1]),2)*y[k]/pow(h[k],3);
p2=(h[k]-2*(u[i]-x[k+1]))*pow((u[i]-x[k]),2)*y[k+1]/pow(h[k],3);

[yonk]

机器中毒了，整了一天。呵呵，你已经搞定了。祝贺！

[chisoft]

Mark