二次样条插值及其C语言的实现.docx
《二次样条插值及其C语言的实现.docx》由会员分享，可在线阅读，更多相关《二次样条插值及其C语言的实现.docx（13页珍藏版）》请在冰豆网上搜索。
二次样条插值及其C语言的实现
数值计算实验报告
实验名称：
三次样条差值及曲线拟合的最小二乘法
姓名：
陈明（09064836）
专业：
生物医学工程（09064811）
1、实验要求：
用C语言编程实现三次样条插值函数及曲线拟合的最小二乘法的求解。
2、实验目的：
、通过编程加深对三次样条插值及曲线拟合的最小二乘法的理解；
、学习用计算机解决工程问题，主要包括数据处理与分析。
3、程序说明：
、三次样条插值函数的求法：
分析：
在解线性方程组确定M值时，因为矩阵是三对角阵，所以我用克拉默法则
=
进行求解，其中D为三对角阵的行列式，值为2^12。
程序如下：
#include
voidmain（）
{
doublex[12],y[12],u[11],v[11],g[12];//x[12]用于存放点的横坐标，y[12]用于存放点的纵坐标，u[11]用于存放μ1~μ10
//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;
for（i=0;i<12;i++）scanf（"%lf",&x[i]）;
printf（"Pleaseinputy（i）:
\n"）;
for（i=0;i<12;i++）scanf（"%lf",&y[i]）;
for（i=1;i<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的求解公式
printf（"u
（1）tou（10）:
\n"）;
for（i=1;i<11;i++）
{
printf（"%f",u[i]）;
j++;
if（j%4==0）printf（"\n"）;
}
printf（"\n"）;
j=0;
printf（"v
（1）tov（10）:
\n"）;
for（i=1;i<11;i++）
{
printf（"%f",v[i]）;
j++;
if（j%4==0）printf（"\n"）;
}
printf（"\n"）;
j=0;
printf（"g（0）tog（11）:
\n"）;
for（i=0;i<12;i++）
{
printf（"%f",g[i]）;
j++;
if（j%4==0）printf（"\n"）;
}
for（i=1;i<=12;i++）
AllD=AllD*2;
for（i=1;i<=10;i++）
{
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];
for（i=1;i<=10;i++）D[i]=g[i]*AllD/2;
for（i=0;i<=11;i++）M[i]=D[i]/AllD;
for（i=1;i<=11;i++）
{
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]）;//计算常数项的公式
}
for（i=0;i<11;i++）
{
printf（"[%.2f,%.2f]:
\n",x[i],x[i+1]）;
printf（"S=%fx^3+%fx^2+%fx+%f\n",X3[i+1],X2[i+1],X1[i+1],X0[i+1]）;
printf（"\n"）;
}
}
、最小二乘拟合函数：
由散点图，设想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.05263
0.03096
0.02041
0.01364
0.01012
0.00928
0.00835
解方程组
=
即
=
由克拉默法则
=
解得a=-0.024，b=1.4878
代入①式，得经验方程
y=
程序如下：
#include
#include
#include
#include
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）
{
printf（"内存分配失败\n"）;
exit（0）;
}
a=（double*）calloc（n,sizeof（double））;
if（a==NULL）
{
printf（"内存分配失败\n"）;
exit（0）;
}
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;i<=m;i++）
printf（"a[%d]=%.10f\n",（i-1）,a[i-1]）;
printf（"拟合多项式与数据点偏差的平方和为：
\n"）;
printf（"%.10e\n",dt1）;
printf（"拟合多项式与数据点偏差的绝对值之和为：
\n"）;
printf（"%.10e\n",dt2）;
printf（"拟合多项式与数据点偏差的绝对值最大值为：
\n"）;
printf（"%.10e\n",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）
{
printf（"内存分配失败\n"）;
exit（0）;
}
t=（double*）calloc（n,sizeof（double））;
if（t==NULL）
{
printf（"内存分配失败\n"）;
exit（0）;
}
b=（double*）calloc（n,sizeof（double））;
if（b==NULL）
{
printf（"内存分配失败\n"）;
exit（0）;
}
z=0;
for（i=1;i<=n;i++）
z=z+x[i-1]/n;
b[0]=1;
d1=n;
p=0;
c=0;
for（i=1;i<=n;i++）
{
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;
c=0;
g=0;
for（i=1;i<=n;i++）
{
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<=m;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];
d2=0;
c=0;
g=0;
for（i=1;i<=n;i++）
{
q=s[j-1];
for（k=j-1;k>=1;k--）
q=q*（x[i-1]-z）+s[k-1];
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[j-1]=c*s[j-1];
t[j-1]=s[j-1];
for（k=j-1;k>=1;k--）
{
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;
for（i=1;i<=n;i++）
{
q=a[m-1];
for（k=m-1;k>=1;k--）
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、心得体会：
通过此次试验我加深了对三次样条插值和曲线拟合的最小二乘法的理解，明白到计算机科学对数值计算的作用，即计算机能极大地减轻人的工作量，但这必须建立在能熟练编程的基础上。
所以二者能否完美地结合取决于我们是否弄清问题本质并用计算机语言表达出来。
升级会员 