Cubic Spline三次样条曲线

//C9.cpp
//Creat Natuarl Cubic Spline

#include <iostream.h>

const int N=6;

void NCS(double x[],double a[],double b[],double c[],double d[],int n)
{
double h[N-1],alpha[N-1],l[N],z[N],u[N];
int i;
for (i=0;i<n-1;i++)
{
h[i]=x[i+1]-x[i];
}
for(i=1;i<n-1;i++)
{
alpha[i]=3*(a[i+1]-a[i])/h[i]-3*(a[i]-a[i-1])/h[i-1];
}
//Now solve the linear system,using Court Factorazation
l[0]=u[0]=z[0]=0;
for(i=1;i<n-1;i++)
{
l[i]=2*(x[i+1]-x[i-1])-h[i-1]*u[i-1];
u[i]=h[i]/l[i];
z[i]=(alpha[i]-h[i-1]*z[i-1])/l[i];
}
l[n-1]=1;
z[n-1]=c[n-1]=0;
for(i=n-2;i>-1;i--)
{
c[i]=z[i]-u[i]*c[i+1];
b[i]=(a[i+1]-a[i])/h[i]-h[i]*(c[i+1]+2*c[i])/3;
d[i]=(c[i+1]-c[i])/(3*h[i]);
}
}

void main()
{
double a[N]={4.0,2.5,1.5,0.5,-0.5,4.0},
x[N]={0.0,1.0,3.5,5.0,7.5,10.0},
c[N],d[N],b[N],
xp=2.5,
s,ds,dds,is;
int i;
    //Compute b[N],c[N],d[N]
NCS(x,a,b,c,d,N);
    //locate xp
for(i=0;i<N-1;i++)
{
if (xp>x[i]&&xp<x[i+1])
break;
}
s=a[i]+b[i]*(xp-x[i])+c[i]*(xp-x[i])*(xp-x[i])+d[i]*(xp-x[i])*(xp-x[i])*(xp-x[i]);
ds=b[i]+2*c[i]*(xp-x[i])+3*d[i]*(xp-x[i])*(xp-x[i]);
dds=2*c[i]+6*d[i]*(xp-x[i]);
is=0;
for (i=0;i<N-1;i++)
{
is+=a[i]*(x[i+1]-x[i])+b[i]/2*(x[i+1]-x[i])*(x[i+1]-x[i])+c[i]/3*(x[i+1]-x[i])*(x[i+1]-x[i])*(x[i+1]-x[i])+d[i]/4*(x[i+1]-x[i])*(x[i+1]-x[i])*(x[i+1]-x[i])*(x[i+1]-x[i]);
}
cout<<"The cubic spline is:"<<endl;
for(i=0;i<N-1;i++)
{
cout<<"From "<<x[i]<<" to "<<x[i+1]<<" :"<<a[i]<<"+"<<b[i]<<"*(x-x[i])+"<<c[i]<<"*(x-x[i])^2+"<<d[i]<<"*(x-x[i])^3 "<<endl;
}
cout<<"S: "<<s<<endl;
cout<<"dS: "<<ds<<endl;
cout<<"ddS: "<<dds<<endl;
cout<<"IS: "<<is<<endl;
cin>>i;
}
//The cubic spline is:
//From 0 to 1 :4+-1.68875*(x-x[i])+0*(x-x[i])^2+0.188746*(x-x[i])^3
//From 1 to 3.5 :2.5+-1.12251*(x-x[i])+0.566237*(x-x[i])^2+-0.110893*(x-x[i])^3
//From 3.5 to 5 :1.5+-0.370573*(x-x[i])+-0.265462*(x-x[i])^2+0.0453779*(x-x[i])^3
//From 5 to 7.5 :0.5+-0.86066*(x-x[i])+-0.0612619*(x-x[i])^2+0.0982103*(x-x[i])^3
//From 7.5 to 10 :-0.5+0.674474*(x-x[i])+0.675315*(x-x[i])^2+-0.0900421*(x-x[i])^3

//在x=2.5处，S: 1.716
//dS: -0.172328
//ddS: 0.134434
//IS: 12.0992