c9.cpp

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

//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