s1.cpp

数值分析,三次样条插值的实验代码,可工大家参考

#include<iostream>
#define N 4
using namespace std;
double x[N+1],fx[N+1],h[N],d[N+1],u[N+1],l[N+1],M[N+1],alpha[N+1],beta[N+1],gama[N+1],y[N+1];
double tiaojian1,tiaojian2;
double juncha2(int i,int j)
{
	return (fx[j]-fx[i])/(x[j]-x[i]);
}
double juncha3(int i,int j,int k)
{
	return (juncha2(j,k)-juncha2(i,j))/(x[k]-x[i]);
}
void Init1()
{
	int i;
	for(i=0;i<=N-1;i++)
		h[i]=x[i+1]-x[i];
	u[N]=1;
	l[0]=1;
	d[0]=(6./h[0])*(juncha2(0,1)-tiaojian1);
	d[N]=(6./h[N-1])*(tiaojian2-juncha2(N-1,N));
}
void Init2()
{
	int i;
	for(i=0;i<=N-1;i++)
		h[i]=x[i+1]-x[i];
	u[N]=l[0]=d[0]=d[N]=0;
}
void Init(int s)
{
	int i;
	if(s==1)
		Init1();
	else
		Init2();
	for(i=1;i<=N-1;i++)
	{
		u[i]=h[i-1]/(h[i-1]+h[i]);
		l[i]=h[i]/(h[i-1]+h[i]);
		d[i]=6*juncha3(i-1,i,i+1);
	}
}
void Zhuigan()
{
	int i;
	beta[0]=l[0]/2.;
	y[0]=d[0]/2.;
	for(i=1;i<=N-1;i++)
		beta[i]=l[i]/(2.-u[i]*beta[i-1]);
	for(i=1;i<=N;i++)
		y[i]=(d[i]-u[i]*y[i-1])/(2.-u[i]*beta[i-1]);
	M[N]=y[N];
	for(i=N-1;i>=0;i--)
		M[i]=y[i]-beta[i]*M[i+1];
}
void Prints()
{
	int i;
	for(i=0;i<N;i++)
	{
		cout<<M[i]/6./h[i]<<"("<<x[i+1]<<"-x)^3";
		cout.setf(ios::showpos);
		cout<<M[i+1]/6./h[i]<<"(x-";
		cout.unsetf(ios::showpos);
		cout<<x[i]<<")^3";
		cout.setf(ios::showpos);
		cout<<(fx[i]-(M[i]*h[i]*h[i])/6.)/h[i]<<"(";
		cout.unsetf(ios::showpos);
		cout<<x[i+1]<<"-x)";
		cout.setf(ios::showpos);
		cout<<(fx[i+1]-(M[i+1]*h[i]*h[i])/6.)/h[i]<<"(x-";
		cout.unsetf(ios::showpos);
		cout<<x[i]<<"),   ["<<x[i]<<","<<x[i+1]<<"]"<<endl;
	}
}
int main()
{
	int i;
	x[0]=0.2;fx[0]=0.9798652;
	x[1]=0.4;fx[1]=0.9177710;
	x[2]=0.6;fx[2]=0.8080348;
	x[3]=0.8;fx[3]=0.6386093;
	x[4]=1.0;fx[4]=0.3843735;
	tiaojian1=0.20271;
	tiaojian2=1.55741;
	Init(1);
	Zhuigan();
	for(i=0;i<=N;i++)
		cout<<M[i]<<endl;
	Prints();
	return 0;
}