二次样条插值及其C语言的实现文档格式.docx
《二次样条插值及其C语言的实现文档格式.docx》由会员分享，可在线阅读，更多相关《二次样条插值及其C语言的实现文档格式.docx（13页珍藏版）》请在冰豆网上搜索。
//v[11]用于存放λ1~λ10，g[12]用于存放g0~g11
inti,j=0;
intAllD=1;
//三角对阵的行列式
doubleMv=1,Mu=1;
//Mv表示λ1*λ2*…*λ10，Mu表示μ1*μ2*…*μ10
doubleM[12],X3[11],X2[11],X1[11],X0[11];
//M[12]用于存放M0~M11，X3[11]用于存放x^3的系数，X2[11]用于存放x^2的系数,
//X1[11]用于存放x的系数,X0[11]用于存放常数项的系数
printf（"
Pleaseinputx（i）:
\n"
）;
for（i=0;
i<
12;
i++）x[i]=0;
i++）scanf（"
%lf"
&
x[i]）;
Pleaseinputy（i）:
y[i]）;
for（i=1;
11;
i++）
{
u[i]=（x[i]-x[i-1]）/（x[i+1]-x[i-1]）;
//μ1~μ10的求解公式
v[i]=1-u[i];
//λ1~λ10的求解公式
g[i]=（6/（x[i+1]-x[i-1]））*（（y[i+1]-y[i]）/（x[i+1]-x[i]）-（y[i]-y[i-1]）/（x[i]-x[i-1]））;
//g1~g10的求解公式
}
g[0]=（6/（x[1]-x[0]））*（（y[1]-y[0]）/（x[1]-x[0]）-0.75）;
//g0的求解公式
g[11]=（6/（x[11]-x[10]））*（18-（y[11]-y[10]）/（x[11]-x[10]））;
//g11的求解公式
u
（1）tou（10）:
printf（"
%f"
u[i]）;
j++;
if（j%4==0）printf（"
j=0;
v
（1）tov（10）:
v[i]）;
g（0）tog（11）:
g[i]）;
=12;
AllD=AllD*2;
=10;
Mv=Mv*v[i];
Mu=Mu*u[i];
//采用克拉默法则进行线性方程组的求解x[i]=D[i]/D
doubleD[12];
D[0]=g[0]*AllD/2+Mv*g[11];
D[11]=g[11]*AllD/2+Mu*g[0];
i++）D[i]=g[i]*AllD/2;
=11;
i++）M[i]=D[i]/AllD;
X3[i]=（M[i]-M[i-1]）/（6*（x[i]-x[i-1]））;
//计算想x^3系数的公式
X2[i]=（x[i]*M[i-1]-x[i-1]*M[i]）/（2*（x[i]-x[i-1]））;
//计算想x^2系数的公式
X1[i]=（-x[i]*x[i]*M[i-1]+x[i-1]*x[i-1]*M[i]）/（2*（x[i]-x[i-1]））+（（y[i]-（M[i]*（x[i]-x[i-1]）*（x[i]-x[i-1]））/6）-（y[i-1]-（M[i-1]*（x[i]-x[i-1]）*（x[i]-x[i-1]））/6））/（x[i]-x[i-1]）;
//计算想x系数的公式
X0[i]=x[i]*x[i]*x[i]-x[i-1]*x[i-1]*x[i-1]-（（y[i]-（M[i]*（x[i]-x[i-1]）*（x[i]-x[i-1]））/6）*x[i-1]）/（x[i]-x[i-1]）+（（y[i-1]-（M[i-1]*（x[i]-x[i-1]）*（x[i]-x[i-1]））/6）*x[i]）/（x[i]-x[i-1]）;
//计算常数项的公式
[%.2f,%.2f]:
x[i],x[i+1]）;
S=%fx^3+%fx^2+%fx+%f\n"
X3[i+1],X2[i+1],X1[i+1],X0[i+1]）;
}
、最小二乘拟合函数：
由散点图，设想y=f（x）是双曲型的，并且具有下面的形式y=
①。
做变量替换y1=
，x1=
则式①变为y1=a+b*x1。
i
1
2
3
4
5
6
7
x1（i）=
0.05263
0.04000
0.03226
0.02632
0.02273
0.01819
0.01852
y1（i）=
0.03096
0.02041
0.01364
0.01012
0.00928
0.00835
解方程组
即
由克拉默法则
解得a=-0.024，b=1.4878
代入①式，得经验方程
y=
#include<
stdlib.h>
malloc.h>
math.h>
Smooth（double*x,double*y,double*a,intn,intm,double*dt1,double*dt2,double*dt3）;
voidmain（）
{
inti,n,m;
double*x,*y,*a,dt1,dt2,dt3,b;
n=7;
m=3;
b=0;
x=（double*）calloc（n,sizeof（double））;
if（x==NULL）
printf（"
内存分配失败\n"
exit（0）;
}
y=（double*）calloc（n,sizeof（double））;
if（y==NULL）
a=（double*）calloc（n,sizeof（double））;
if（a==NULL）
x[0]=19;
;
x[1]=25;
x[2]=31;
x[3]=38;
x[4]=40;
x[5]=50;
x[6]=54;
y[0]=19.0;
y[1]=32.3;
y[2]=49.0;
y[3]=73.3;
y[4]=98.8;
y[5]=107.8;
y[6]=119.7;
Smooth（x,y,a,n,m,&
dt1,&
dt2,&
dt3）;
for（i=1;
=m;
i++）
a[%d]=%.10f\n"
（i-1）,a[i-1]）;
拟合多项式与数据点偏差的平方和为：
%.10e\n"
dt1）;
拟合多项式与数据点偏差的绝对值之和为：
dt2）;
拟合多项式与数据点偏差的绝对值最大值为：
dt3）;
free（x）;
free（y）;
free（a）;
Smooth（double*x,double*y,double*a,intn,intm,double*dt1,double*dt2,double*dt3）
inti,j,k;
double*s,*t,*b,z,d1,p,c,d2,g,q,dt;
s=（double*）calloc（n,sizeof（double））;
if（s==NULL）
t=（double*）calloc（n,sizeof（double））;
if（t==NULL）
b=（double*）calloc（n,sizeof（double））;
if（b==NULL）
z=0;
=n;
z=z+x[i-1]/n;
b[0]=1;
d1=n;
p=0;
c=0;
p=p+x[i-1]-z;
c=c+y[i-1];
c=c/d1;
p=p/d1;
a[0]=c*b[0];
if（m>
1）
t[1]=1;
t[0]=-p;
d2=0;
g=0;
q=x[i-1]-z-p;
d2=d2+q*q;
c=y[i-1]*q+c;
g=（x[i-1]-z）*q*q+g;
c=c/d2;
p=g/d2;
q=d2/d1;
d1=d2;
a[1]=c*t[1];
a[0]=c*t[0]+a[0];
for（j=3;
j<
j++）
s[j-1]=t[j-2];
s[j-2]=-p*t[j-2]+t[j-3];
if（j>
=4）
for（k=j-2;
k>
=2;
k--）
s[k-1]=-p*t[k-1]+t[k-2]-q*b[k-1];
s[0]=-p*t[0]-q*b[0];
q=s[j-1];
for（k=j-1;
=1;
q=q*（x[i-1]-z）+s[k-1];
a[j-1]=c*s[j-1];
t[j-1]=s[j-1];
a[k-1]=c*s[k-1]+a[k-1];
b[k-1]=t[k-1];
t[k-1]=s[k-1];
*dt1=0;
*dt2=0;
*dt3=0;
q=a[m-1];
for（k=m-1;
q=q*（x[i-1]-z）+a[k-1];
dt=q-y[i-1];
if（fabs（dt）>
*dt3）
*dt3=fabs（dt）;
*dt1=*dt1+dt*dt;
*dt2=*dt2+fabs（dt）;
free（s）;
free（t）;
free（b）;
return
（1）;
4、上机调试说明：
、三次样条差值函数：
5、心得体会：
通过此次试验我加深了对三次样条插值和曲线拟合的最小二乘法的理解，明白到计算机科学对数值计算的作用，即计算机能极大地减轻人的工作量，但这必须建立在能熟练编程的基础上。
所以二者能否完美地结合取决于我们是否弄清问题本质并用计算机语言表达出来。
升级会员 