c10.cpp

数值分析最常用的四十种算法

//C10.cpp
//Creat Clamped Cubic Spline

#include <iostream.h>

const int N=6;

void CCS(double x[],double a[],double b[],double c[],double d[],double df0,double dfn,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];
	}
	alpha[0]=2*(a[1]-a[0])/h[0]-3*df0;
	alpha[n-1]=3*dfn-3*(a[n-1]-a[n-2])/h[n-2];
	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]=2*h[0];
	u[0]=0.5;
	z[0]=alpha[0]/l[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]=h[n-2]*(2-u[n-2]);
	z[n-1]=(alpha[n-1]-h[n-2]*z[n-2])/l[n-1];
	c[n-1]=z[n-1];
	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,
		df0=2.5,dfn=1.4,
		s,ds,dds,is;
	int i;
    //Compute b[N],c[N],d[N]
	CCS(x,a,b,c,d,df0,dfn,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+4.75*(x-x[i])+4.33994*(x-x[i])^2+1.18606*(x-x[i])^3
//From 1 to 3.5 :2.5+-0.439232*(x-x[i])+0.0701166*(x-x[i])^2+-0.0217695*(x-x[i])^3
//From 3.5 to 5 :1.5+-0.496827*(x-x[i])+-0.0931545*(x-x[i])^2+-0.0133812*(x-x[i])^3
//From 5 to 7.5 :0.5+-0.866614*(x-x[i])+-0.15337*(x-x[i])^2+0.136006*(x-x[i])^3
//From 7.5 to 10 :-0.5+0.916654*(x-x[i])+0.866677*(x-x[i])^2+-0.205335*(x-x[i])^3

//在x=2.5处，S: 1.92544
//dS: -0.375826
//ddS: -0.0556922
//IS: 17.912