/*********************下面程序是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
*/
问题补充，因字数限制，挪到这
1.拉格朗日插值简介：
对给定的n个插值节点x1,x2,…,xn,及其对应的函数值y1=f(x1), y2=f(x2),…, yn=f(xn);使用拉格朗日插值公式，计算在x点处的对应的函数值f(x);
2.一维拉格朗日插值c语言程序：
Int lagrange(x0, y0, n, x, y)
Float xo[], yo[], x;
Int n;
Float *y
{
Int i, j;
Float p;
*y=0;
If (n1)
{
For(i=0;in;i++)
{
P=1;
For(j=1;jn;j++)
{
If(i!=J)
P=p*(x-x0[j]/x0[i]-x0[j]);
}
*y=*y+p*y0[i];
Return(0);
}
Else
Return(-1);
}
3.例题。已知函数如下表所示，求x=0.472处的函数值：
X 0.46 0.47 0.48 0.49
Y 0.484655 0.4903745 0.502750 0.511668
计算这个问题的c语言程序如下：
#minclude stdio
#includeMnath.h
Main()
{
Float x0[4]={ 0.46, 0.47,0.48,0.49};
Float y0[4]={ 0.484655 ,0.4903745 ,0.502750 ,0.511668};
Float x, y;
Int n, rtn;
N=4;
X=0.472;
Rth=lagrange(x0,y0,n,x,y);
If(rtn=0)
{
Prinf(“Y(0.472)=:%f\n”,y);
}
Else
{
Prinf(“n must be larger than 1.\n”);
}
}
计算结果：Y(0.472)=:0.495553
4.问题补充
我的问题与上面的例子类似，计算三维空间一点(x,y,z)对应的函数值（Vx,Vy,Vz）.不同的是自变量(point_coordinate.txt)为三维空间散乱点（不是正方体的顶点），因变量(point_data.txt)为矢量（向量 ）。插值算法比较多，常数法，拉格朗日插值，埃特金插值，三阶样条插值等。最简单的就是常数法，查找离目标点（x,y,z）距离最近的已知自变量(Xi,Yi,Zi)，把该点的函数值赋给目标点做函数值，求高手帮忙写写。
#includestdio.h
#includestdlib.h
#includeiostream.h
typedef struct data
{
float x;
float y;
}Data;//变量x和函数值y的结构
Data d[20];//最多二十组数据
float f(int s,int t)//牛顿插值法，用以返回插商
{
if(t==s+1)
return (d[t].y-d[s].y)/(d[t].x-d[s].x);
else
return (f(s+1,t)-f(s,t-1))/(d[t].x-d[s].x);
}
float Newton(float x,int count)
{
int n;
while(1)
{
cout"请输入n值(即n次插值):";//获得插值次数
cinn;
if(n=count-1)// 插值次数不得大于count－1次
break;
else
system("cls");
}
//初始化t，y，yt。
float t=1.0;
float y=d[0].y;
float yt=0.0;
//计算y值
for(int j=1;j=n;j++)
{
t=(x-d[j-1].x)*t;
yt=f(0,j)*t;
//coutf(0,j)endl;
y=y+yt;
}
return y;
}
float lagrange(float x,int count)
{
float y=0.0;
for(int k=0;kcount;k++)//这儿默认为count－1次插值
{
float p=1.0;//初始化p
for(int j=0;jcount;j++)
{//计算p的值
if(k==j)continue;//判断是否为同一个数
p=p*(x-d[j].x)/(d[k].x-d[j].x);
}
y=y+p*d[k].y;//求和
}
return y;//返回y的值
}
void main()
{
float x,y;
int count;
while(1)
{
cout"请输入x[i],y[i]的组数，不得超过20组:";//要求用户输入数据组数
cincount;
if(count=20)
break;//检查输入的是否合法
system("cls");
}
//获得各组数据
for(int i=0;icount;i++)
{
cout"请输入第"i+1"组x的值:";
cind[i].x;
cout"请输入第"i+1"组y的值:";
cind[i].y;
system("cls");
}
cout"请输入x的值:";//获得变量x的值
cinx;
while(1)
{
int choice=3;
cout"请您选择使用哪种插值法计算："endl;
cout" (0):退出"endl;
cout" (1):Lagrange"endl;
cout" (2):Newton"endl;
cout"输入你的选择：";
cinchoice;//取得用户的选择项
if(choice==2)
{
cout"你选择了牛顿插值计算方法，其结果为：";
y=Newton(x,count);break;//调用相应的处理函数
}
if(choice==1)
{
cout"你选择了拉格朗日插值计算方法，其结果为：";
y=lagrange(x,count);break;//调用相应的处理函数
}
if(choice==0)
break;
system("cls");
cout"输入错误!!!!"endl;
}
coutx" , "yendl;//输出最终结果
}
void
SPL(int
n,
double
*x,
double
*y,
int
ni,
double
*xi,
double
*yi)；
是你所要。
已知
n
个点
x,y;
x
必须已按顺序排好。要插值
ni
点，横坐标
xi[],
输出
yi[]。
程序里用double
型，保证计算精度。
SPL调用现成的程序。
现成的程序很多。端点处理方法不同，结果会有不同。想同matlab比较，你需
尝试
调用
spline（）函数
时，令
end1
为
1,
设
slope1
的值，令
end2
为
1
设
slope2
的值。
#include
stdio.h
#include
math.h
int
spline
(int
n,
int
end1,
int
end2,
double
slope1,
double
slope2,
double
x[],
double
y[],
double
b[],
double
c[],
double
d[],
int
*iflag)
{
int
nm1,
ib,
i,
ascend;
double
t;
nm1
=
n
-
1;
*iflag
=
0;
if
(n
2)
{
/*
no
possible
interpolation
*/
*iflag
=
1;
goto
LeaveSpline;
}
ascend
=
1;
for
(i
=
1;
i
n;
++i)
if
(x[i]
=
x[i-1])
ascend
=
0;
if
(!ascend)
{
*iflag
=
2;
goto
LeaveSpline;
}
if
(n
=
3)
{
d[0]
=
x[1]
-
x[0];
c[1]
=
(y[1]
-
y[0])
/
d[0];
for
(i
=
1;
i
nm1;
++i)
{
d[i]
=
x[i+1]
-
x[i];
b[i]
=
2.0
*
(d[i-1]
+
d[i]);
c[i+1]
=
(y[i+1]
-
y[i])
/
d[i];
c[i]
=
c[i+1]
-
c[i];
}
/*
----
Default
End
conditions
*/
b[0]
=
-d[0];
b[nm1]
=
-d[n-2];
c[0]
=
0.0;
c[nm1]
=
0.0;
if
(n
!=
3)
{
c[0]
=
c[2]
/
(x[3]
-
x[1])
-
c[1]
/
(x[2]
-
x[0]);
c[nm1]
=
c[n-2]
/
(x[nm1]
-
x[n-3])
-
c[n-3]
/
(x[n-2]
-
x[n-4]);
c[0]
=
c[0]
*
d[0]
*
d[0]
/
(x[3]
-
x[0]);
c[nm1]
=
-c[nm1]
*
d[n-2]
*
d[n-2]
/
(x[nm1]
-
x[n-4]);
}
/*
Alternative
end
conditions
--
known
slopes
*/
if
(end1
==
1)
{
b[0]
=
2.0
*
(x[1]
-
x[0]);
c[0]
=
(y[1]
-
y[0])
/
(x[1]
-
x[0])
-
slope1;
}
if
(end2
==
1)
{
b[nm1]
=
2.0
*
(x[nm1]
-
x[n-2]);
c[nm1]
=
slope2
-
(y[nm1]
-
y[n-2])
/
(x[nm1]
-
x[n-2]);
}
/*
Forward
elimination
*/
for
(i
=
1;
i
n;
++i)
{
t
=
d[i-1]
/
b[i-1];
b[i]
=
b[i]
-
t
*
d[i-1];
c[i]
=
c[i]
-
t
*
c[i-1];
}
/*
Back
substitution
*/
c[nm1]
=
c[nm1]
/
b[nm1];
for
(ib
=
0;
ib
nm1;
++ib)
{
i
=
n
-
ib
-
2;
c[i]
=
(c[i]
-
d[i]
*
c[i+1])
/
b[i];
}
b[nm1]
=
(y[nm1]
-
y[n-2])
/
d[n-2]
+
d[n-2]
*
(c[n-2]
+
2.0
*
c[nm1]);
for
(i
=
0;
i
nm1;
++i)
{
b[i]
=
(y[i+1]
-
y[i])
/
d[i]
-
d[i]
*
(c[i+1]
+
2.0
*
c[i]);
d[i]
=
(c[i+1]
-
c[i])
/
d[i];
c[i]
=
3.0
*
c[i];
}
c[nm1]
=
3.0
*
c[nm1];
d[nm1]
=
d[n-2];
}
else
{
b[0]
=
(y[1]
-
y[0])
/
(x[1]
-
x[0]);
c[0]
=
0.0;
d[0]
=
0.0;
b[1]
=
b[0];
c[1]
=
0.0;
d[1]
=
0.0;
}
LeaveSpline:
return
0;
}
double
seval
(int
n,
double
u,
double
x[],
double
y[],
double
b[],
double
c[],
double
d[],
int
*last)
{
int
i,
j,
k;
double
w;
i
=
*last;
if
(i
=
n-1)
i
=
0;
if
(i
0)
i
=
0;
if
((x[i]
u)
||
(x[i+1]
u))
{
i
=
0;
j
=
n;
do
{
k
=
(i
+
j)
/
2;
if
(u
x[k])
j
=
k;
if
(u
=
x[k])
i
=
k;
}
while
(j
i+1);
}
*last
=
i;
w
=
u
-
x[i];
w
=
y[i]
+
w
*
(b[i]
+
w
*
(c[i]
+
w
*
d[i]));
return
(w);
}
void
SPL(int
n,
double
*x,
double
*y,
int
ni,
double
*xi,
double
*yi)
{
double
*b,
*c,
*d;
int
iflag,last,i;
b
=
(double
*)
malloc(sizeof(double)
*
n);
c
=
(double
*)malloc(sizeof(double)
*
n);
d
=
(double
*)malloc(sizeof(double)
*
n);
if
(!d)
{
printf("no
enough
memory
for
b,c,d\n");}
else
{
spline
(n,0,0,0,0,x,y,b,c,d,iflag);
if
(iflag==0)
printf("I
got
coef
b,c,d
now\n");
else
printf("x
not
in
order
or
other
error\n");
for
(i=0;ini;i++)
yi[i]
=
seval(ni,xi[i],x,y,b,c,d,last);
free(b);free(c);free(d);
};
}
main(){
double
x[6]={0.,1.,2.,3.,4.,5};
double
y[6]={0.,0.5,2.0,1.6,0.5,0.0};
double
u[8]={0.5,1,1.5,2,2.5,3,3.5,4};
double
s[8];
int
i;
SPL(6,
x,y,
8,
u,
s);
for
(i=0;i8;i++)
printf("%lf
%lf
\n",u[i],s[i]);
return
0;
}
