声振论坛

 找回密码
 我要加入

QQ登录

只需一步,快速开始

查看: 2599|回复: 1

[求助]询问matlab中的spline函数

[复制链接]
发表于 2006-5-18 09:35 | 显示全部楼层 |阅读模式

马上注册,结交更多好友,享用更多功能,让你轻松玩转社区。

您需要 登录 才可以下载或查看,没有账号?我要加入

x
<BR>matlab中的spline是三次样条拟合函数吗?<BR>如果是的话,是按哪种边界条件的三次样条?<BR>怎么我用C语言编写的三次样条插值和matlab的spline得出的结果差别那么大?
回复
分享到:

使用道具 举报

发表于 2006-5-18 11:06 | 显示全部楼层

回复:(xcsyb)[求助]询问matlab中的spline函数

spline三次样条多项式拟合,数据点处光滑--左导等于右导<BR><BR>你试一下下面的c语言写的函数的结果是否一致?<BR><BR>
<P>/***************************************************************  */<BR>/*     具体步骤:                                                  */<BR>/*     ①根据给出的样本点x(i),y(i)计算边界条件系数r(i),t(i)       */<BR>/*     ②在边界条件下求c(1),c(2)......c(n)                        */<BR>/*     ③计算个子区间[x(i),x(i+1)]上各点的函数值然后用画线函数绘出*/<BR>/*                    具体过程请参考计算方法                      */<BR>/******************************************************************/<BR>#include&lt;stdio.h&gt;<BR>#include&lt;graphics.h&gt;<BR>main()<BR>{<BR>    int x[7]={0,150,260,370,450,520,600};<BR>    int y[7]={0,80,120,20,60,100,50};<BR>    int gdriver=DETECT,gmode;<BR>    initgraph(&amp;gdriver,&amp;gmode,"");<BR>    Wbox(x,y,6);<BR>    spline(x,y,7,3);<BR>    getch();<BR>    closegraph();<BR>}<BR>spline(x,y,n,color)<BR>int x[],y[],n,color;<BR>{<BR>   float a[50],b[50],r[50],t[50],c[50],h[50],u2,u3,v2,v3;<BR>   register int i,x0,y0,x1,y1,j;<BR>   for(i=0;i&lt;n-1;++i)<BR>      h=(float)(x[i+1]-x);<BR>   r[0]=1.0;<BR>   r[n-1]=0.0;<BR>   t[0]=(float)(y[1]-y[0])*3.0/h[0];<BR>   t[n-1]=(float)(y[n-1]-y[n-2])*3.0/h[n-2];<BR>   for(i=1;i&lt;n-1;++i)<BR>       {<BR>          r=h[i-1]/(h[i-1]+h);<BR>          t=3.0*((1-r)*(float)(y-y[i-1])/h[i-1]+r*(float)(y[i+1]-y)/h);<BR>       }<BR>   a[0]=-r[0]/2.0;<BR>   b[0]=t[0]/2.0;<BR>   for(i=1;i&lt;n;++i)<BR>       {<BR>          a=-r/(2.0+(1.0-r)*a[i-1]);<BR>          b=(t-(1.0-r)*b[i-1])/(2.0+(1.0-r)*a[i-1]);<BR>       }<BR>   c[n-1]=b[n-1];<BR>   for(i=1;i&lt;n;++i)<BR>       {<BR>          j=n-1-i;<BR>         c[j]=a[j]*c[j+1]+b[j];<BR>       }<BR>   for(i=0;i&lt;n-1;++i)<BR>       {<BR>          if(x+1&gt;=x[i+1])<BR>             {<BR>                  setcolor(color);<BR>                  line(x,y,x[i+1],y[i+1]);<BR>             }<BR>           else<BR>               {<BR>                   x0=x;<BR>                   y0=y;<BR>                   for(j=x+1;j&lt;=x[i+1];++j)<BR>                       {<BR>                            u2=((float)(x[i+1]-j)/h)*((float)(x[i+1]-j)/h);<BR>                           u3=u2*((float)(x[i+1]-j)/h);<BR>                           v2=((float)(j-x)/h)*((float)(j-x)/h);<BR>                           v3=v2*((float)(j-x)/h);<BR>                           y1=(int)((3.0*u2-2.0*u3)*y+(3.0*v2-2.0*v3)*y[i+1]<BR>                                  +h*c*(u2-u3)-h*c[i+1]*(v2-v3));<BR>                           x1=j;<BR>                           setcolor(color);<BR>                           line(x0,y0,x1,y1);<BR>                           x0=x1;<BR>                           y0=y1;<BR>                       }<BR>                }<BR>         }<BR>}<BR> Wbox(int x[],int y[],int n)<BR>/*int x[],y[],n;*/<BR>{<BR>   int i;<BR>   for(i=0;i&lt;n;++i)<BR>       {<BR>           putpixel((x-1),(y-1),6);<BR>           putpixel((x-1),y,6);<BR>           putpixel(x-1,y+1,6);<BR>           putpixel(x,y-1,6);<BR>           putpixel(x,y+1,6);<BR>           putpixel(x+1,y-1,6);<BR>           putpixel(x+1,y,6);<BR>           putpixel(x+1,y+1,6);<BR>       }<BR>}</P>
您需要登录后才可以回帖 登录 | 我要加入

本版积分规则

QQ|小黑屋|Archiver|手机版|联系我们|声振论坛

GMT+8, 2024-9-25 17:14 , Processed in 0.062550 second(s), 18 queries , Gzip On.

Powered by Discuz! X3.4

Copyright © 2001-2021, Tencent Cloud.

快速回复 返回顶部 返回列表