sanmu.cpp

三次样条函数算法C++源代码,用于工程上计算平滑曲线的插值计算.

// sanmu.cpp : Defines the entry point for the console application.


#include "stdafx.h"
#include <math.h>
#include <iostream>
using namespace std;

#define J 11
 
int main(int argc, char* argv[])
{
	//start my program
	double x[11]={0,20,40,60,80,100,120,140,160,180,280};
	double f[11]={4.0,4.001,4.047,4.245,4.584,4.874,5.045,5.184,5.315,5.442,5.442};
	double h[J],lanbd[J],u[J],ff[J][J],fff[J],d[J];
	int i,j; //lanbd = h[j+1]/(h[j]+h[j+1])
	// declare finish
	for(j=1;j<J;j++)
		h[j]=x[j]-x[j-1];
	for(j=1;j<J-1;j++)
	{
		lanbd[j]=h[j+1]/(h[j]+h[j+1]);
		u[j]=1-lanbd[j];
	}
	for(j=1;j<J-1;j++)
	{	cout<<h[j]<<endl;
		cout<<lanbd[j]<<endl;
		cout<<u[j]<<endl;
	}
//input data fjnjhsh,first f[xj x+1]
	for(j=0;j<J-2;j++)
		for(i=1;i<J;i++)
	{
		ff[j][i]=(f[j]-f[i])/(x[j]-x[i]);	
		cout<<ff[j][i]<<" ";
	}
		cout<<endl;
	for(j=1;j<J-1;j++)
	{
		fff[j]=(ff[j-1][j+1]-ff[j-1][j])/(x[j+1]-x[j]);
		cout<<fff[j]<<endl;
	}
	for(j=1;j<J-1;j++)
	{
		d[j]=6*fff[j];
	cout<<d[j]<<endl;
	}
//finish right value of equotion d[j],start the computer equotion,zhui gan method
//prepare for next
// a[i]=u[i] (i=2,n-1),b[i]=2;c[i]=lanbd[i](i=1,n-2)
//(1) beta i  for next step  
	double beta[J];
	beta[1]=lanbd[1]/2;
	for(i=2;i<J-2;i++)
		beta[i]=lanbd[i]/(2-u[i]*beta[i-1]);
//(2) y[i]   for next step
	double y[J];
	y[1]=d[1]/2;
	for(i=2;i<J-1;i++)
	y[i]=(d[i]-u[i]*y[i-1])/(2-u[i]*beta[i-1]);
//(3) solution for next
	double M[J];
	M[J-1]=0;
	M[0]=0;
	M[J-2]=y[J-2];
	for(i=J-3;i>=1;i--)
	M[i]=y[i]-beta[i]*M[i+1];
	for(i=0;i<J;i++)
	cout<<M[i]<<endl;
//here  is  s(x) with diffrent range
cout<<"s(x)=(x[j+1]-x)^3/(6*h[j+1]*M[j]+(x-x[j])^3/(6*h[j+1])*M[j+1]+(x-x[j])/h[j+1]*\
(f[j+1]-f[j]-h[j+1]^2/6*(M[J+1]-M[j]))+f[j]-h[j+1]^2/6*M[j]"<<endl;
////////////////////////////////////////////
	double xx=0,sx=0;
	cout<<"please input value between "<<x[0]<<"and"<<x[J-1]<<endl;
		cin>>xx;
/*	if(x[0]>xx||xx>x[J-1]) 
		cout<<"error,pls check your data input!"<<endl;
	else if(xx==x[0])
		cout<<"sx="<<f[0]<<endl;
	else if(xx==x[J-1])
		cout<<"sx="<<f[J-1]<<endl;
	else
*/
	for(i=0;i<J;i++)
		if(xx==x[i])
		{	sx=f[i];
			cout<<"sx="<<f[i]<<endl;
		}
		else 
			continue;
	for(i=0;i<J-1;i++)
	   if(x[i]<xx&&xx<x[i+1])
	   {
	    	j=i;
   sx=pow((x[j+1]-xx),3)/(6*h[j+1])*M[j]+pow((xx-x[j]),3)/(6*h[j+1])*M[j+1]+(xx-x[j])/h[j+1]*\
(f[j+1]-f[j]-pow(h[j+1],2)/6*(M[J+1]-M[j]))+f[j]-pow(h[j+1],2)/6*M[j];
   cout<<sx<<endl;
	}
	   
//xx=0.5,sx=0.7069032821, pass! god ,take one day! it's 21:17 now
//the end

	return 0;
}