三次自然样条插值编程.doc

/ .cpp : Defines the entry point for the console application./#include stdafx.h#include stdio.hvoid GetValueSpline(int n,double x,double y,double dy,double ddy,int m,double t,double z,double dz,double ddz)int i,j;double h0,h1,alpha,beta,g,*s;h0=x1-x0;s0=3.0*(y1-y0)/2.0*h0-ddy0*h0/4.0;for(j=1;j=0;j-)dyj=dyj*dyj+1+sj;for(j=0;j=n-2;j+)sj=xj+1-xj;for(j=0;j=n-2;j+)h1=sj*sj;ddyj=6.0*(yj+1-yj)/h1-2.0*(2.0*dyj+dyj+1)/sj;h1=sn-2*sn-2;ddyn-1=6.0*(yn-2-yn-1)/h1+2.0*(2.0*dyn-1+dyn-2)/sn-2;g=0.0;for(i=0;i=n-2;i+)h1=0.5*si*(yi+yi+1);h1=h1-si*si*si*(ddyi+ddyi+1)/24.0;g=g+h1;for(j=0;j=xn-1)i=n-2;elsei=0;while(tjxi+1)i=i+1;h1=(xi+1-tj)/si;h0=h1*h1;zj=(3.0*h0-2.0*h0*h1)*yi;zj=zj+si*(h0-h0*h1)*dyi;dzj=6.0*(h0-h1)*yi/si;dzj=dzj+(3.0*h0-2.0*h1)*dyi;ddzj=(6.0-12.0*h1)*yi/(si*si);ddzj=ddzj+(2.0-6.0*h1)*dyi/si;h1=(tj-xi)/si;h0=h1*h1;zj=zj+(3.0*h0-2.0*h0*h1)*yi+1;zj=zj-si*(h0-h0*h1)*dyi+1;dzj=dzj-6.0*(h0-h1)*yi+1/si;dzj=dzj+(3.0*h0-2.0*h1)*dyi+1;ddzj=ddzj+(6.0-12.0*h1)*yi+1/(si*si);ddzj=ddzj-(2.0-6.0*h1)*dyi+1/si;for(i=0;i4;i+);printf( n );printf( t(i) z(i) dz(i) ddz(i) n);printf( %f %f %f %fn,ti,zi,dzi,ddzi);/delete s;/return (g);void main()int i;double x5,y5,dy5,ddy5;double t4,z4,dz4,ddz4;printf(:n);for(i=0;i5;i+)printf(x(%d)=,i);scanf(%f ,&xi);printf( y(%d)=,i);scanf(%fn,&yi);printf(:n); for(i=0;i4;i+)printf(t(%d)=,i);scanf(%fn,&ti);GetValueSpline(5,x5,y5,dy5,ddy5,4,t4,z4,dz4,ddz4);/for(i=0;i4;i+);/printf( n );/printf( t(i) z(i) dz(i) ddz(i) n);/for(i=0;i4;i+)/printf( %f %f %f %fn,ti,zi,dzi,ddzi);/*/ aa.cpp : Defines the entry point for the console application./#include stdafx.h#includestdio.hvoid main()int i,j;double h0,h1,alpha,beta,*s;double x5,y5,dy5,ddy5;double t4,z4,dz4,ddz4;int m,n;m=5;n=4; x0=-3.0;x1=-1.0;x2=2.0;x3=4.0;x4=7.0;y0=3.0;y1=6.0;y2=8.0;y3=2.0;y4=5.0;t0=-2.5;t1=0.0;t2=2.5;t3=5.0;dy0=-0.5;ddy0=0;ddy4=0;h0=x1-x0;s0=3.0*(y1-y0)/2.0*h0-ddy0*h0/4.0;for(j=1;j=0;j-) dyj=dyj*dyj+1+sj;for(j=0;j=n-2;j+) sj=xj+1-xj;for(j=1;j=n-2;j+)h1=sj*sj; ddyj=6.0*(yj+1-yj)/h1-2.0*(2.0*dyj+dyj+1)/sj;for(j=0;j=xn-1)i=n-2;elsei=0;while(tjxi+1) i=i+1;h1=(xi+1-tj)/si;h0=h1*h1;zj=(3.0*h0-2.0*h0*h1)*yi;zj=zj+si*(h0-h0*h1)*dyi;dzj=6.0*(h0-h1)*yi/si;dzj=dzj+(3.0*h0-2.0*h1)*dyi;ddzj=(6.0-12.0*h1)*yi/(si*si);ddzj=ddzj+(2.0-6.0*h1)*dyi/si;h1=(tj-xi)/si;h0=h1*h1;zj=zj+(3.0*h0-2.0*h0*h1)*yi+1;zj=zj-si*(h0-h0*h1)*dyi+1;dzj=dzj-6.0*(h0-h1)*yi+1/si;dzj=dzj+(3.0*h0-2.0*h1)*dyi+1;ddzj=ddzj+(6.0-12.0*h1)*yi+1/(si*si);ddz